{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/oesteban/workspace/niworkflows/niworkflows/__init__.py:24: UserWarning: \n", "This call to matplotlib.use() has no effect because the backend has already\n", "been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n", "or matplotlib.backends is imported for the first time.\n", "\n", "The backend was *originally* set to 'pgf' by the following code:\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/runpy.py\", line 193, in _run_module_as_main\n", " \"__main__\", mod_spec)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/runpy.py\", line 85, in _run_code\n", " exec(code, run_globals)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py\", line 16, in \n", " app.launch_new_instance()\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/traitlets/config/application.py\", line 658, in launch_instance\n", " app.start()\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/kernelapp.py\", line 477, in start\n", " ioloop.IOLoop.instance().start()\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/zmq/eventloop/ioloop.py\", line 177, in start\n", " super(ZMQIOLoop, self).start()\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/tornado/ioloop.py\", line 888, in start\n", " handler_func(fd_obj, events)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/tornado/stack_context.py\", line 277, in null_wrapper\n", " return fn(*args, **kwargs)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py\", line 440, in _handle_events\n", " self._handle_recv()\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py\", line 472, in _handle_recv\n", " self._run_callback(callback, msg)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py\", line 414, in _run_callback\n", " callback(*args, **kwargs)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/tornado/stack_context.py\", line 277, in null_wrapper\n", " return fn(*args, **kwargs)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 283, in dispatcher\n", " return self.dispatch_shell(stream, msg)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 235, in dispatch_shell\n", " handler(stream, idents, msg)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 399, in execute_request\n", " user_expressions, allow_stdin)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/ipkernel.py\", line 196, in do_execute\n", " res = shell.run_cell(code, store_history=store_history, silent=silent)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/zmqshell.py\", line 533, in run_cell\n", " return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2717, in run_cell\n", " interactivity=interactivity, compiler=compiler, result=result)\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2821, in run_ast_nodes\n", " if self.run_code(code, result):\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2881, in run_code\n", " exec(code_obj, self.user_global_ns, self.user_ns)\n", " File \"\", line 14, in \n", " from matplotlib import pyplot as plt\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/matplotlib/pyplot.py\", line 71, in \n", " from matplotlib.backends import pylab_setup\n", " File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/matplotlib/backends/__init__.py\", line 16, in \n", " line for line in traceback.format_stack()\n", "\n", "\n", " matplotlib.use('Agg')\n", "/home/oesteban/.anaconda3/lib/python3.6/importlib/_bootstrap.py:219: ImportWarning: can't resolve package from __spec__ or __package__, falling back on __name__ and __path__\n", " return f(*args, **kwds)\n" ] } ], "source": [ "# %load_ext autoreload\n", "# %autoreload 2\n", "# %matplotlib inline\n", "import os\n", "from pathlib import Path\n", "import warnings\n", "\n", "import numpy as np\n", "import nibabel as nb\n", "import pandas as pd\n", "\n", "import matplotlib as mpl\n", "mpl.use('pgf')\n", "from matplotlib import pyplot as plt\n", "from matplotlib import gridspec\n", "import seaborn as sn\n", "import palettable\n", "\n", "from niworkflows.data import get_template\n", "\n", "from nilearn.image import concat_imgs, mean_img\n", "from nilearn import plotting\n", "\n", "warnings.simplefilter('ignore')\n", "\n", "DATA_HOME = Path(os.getenv('FMRIPREP_DATA_HOME', os.getcwd())).resolve()\n", "DS030_HOME = DATA_HOME / 'ds000030' / '1.0.3'\n", "DERIVS_HOME = DS030_HOME / 'derivatives'\n", "ATLAS_HOME = get_template('MNI152NLin2009cAsym')\n", "ANALYSIS_HOME = DERIVS_HOME / 'fmriprep_vs_feat_2.0-oe'\n", "\n", "fprep_home = DERIVS_HOME / 'fmriprep_1.0.8' / 'fmriprep'\n", "feat_home = DERIVS_HOME / 'fslfeat_5.0.10' / 'featbids'\n", "\n", "out_folder = Path(os.getenv('FMRIPREP_OUTPUTS') or '').resolve()\n", "\n", "# Load MNI152 nonlinear, asymmetric 2009c atlas\n", "atlas = nb.load(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_T1w.nii.gz'))\n", "mask1mm = nb.load(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_brainmask.nii.gz')).get_data() > 0\n", "mask2mm = nb.load(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-02_brainmask.nii.gz')).get_data() > 0\n", "data = atlas.get_data()\n", "data[~mask1mm] = data[~mask1mm].max()\n", "atlas = nb.Nifti1Image(data, atlas.affine, atlas.header)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# sn.set_style(\"whitegrid\", {\n", "# 'ytick.major.size': 5,\n", "# 'xtick.major.size': 5,\n", "# })\n", "# sn.set_context(\"notebook\", font_scale=1.5)\n", "\n", "# pgf_with_custom_preamble = {\n", "# 'ytick.major.size': 0,\n", "# 'xtick.major.size': 0,\n", "# 'font.size': 30,\n", "# 'font.sans-serif': ['HelveticaLTStd-Light'],\n", "# 'font.family': 'sans-serif', # use serif/main font for text elements\n", "# 'text.usetex': False, # use inline math for ticks\n", "# }\n", "# mpl.rcParams.update(pgf_with_custom_preamble)\n", "\n", "\n", "pgf_with_custom_preamble = {\n", " 'text.usetex': True, # use inline math for ticks\n", " 'pgf.rcfonts': False, # don't setup fonts from rc parameters\n", " 'pgf.texsystem': 'xelatex',\n", " 'verbose.level': 'debug-annoying',\n", " \"pgf.preamble\": [\n", " r\"\"\"\\usepackage{fontspec}\n", "\\setsansfont{HelveticaLTStd-Light}[\n", "Extension=.otf,\n", "BoldFont=HelveticaLTStd-Bold,\n", "ItalicFont=HelveticaLTStd-LightObl,\n", "BoldItalicFont=HelveticaLTStd-BoldObl,\n", "]\n", "\\setmainfont{HelveticaLTStd-Light}[\n", "Extension=.otf,\n", "BoldFont=HelveticaLTStd-Bold,\n", "ItalicFont=HelveticaLTStd-LightObl,\n", "BoldItalicFont=HelveticaLTStd-BoldObl,\n", "]\n", "\\setmonofont{Inconsolata-dz}\n", "\"\"\",\n", " r'\\renewcommand\\familydefault{\\sfdefault}',\n", "# r'\\setsansfont[Extension=.otf]{Helvetica-LightOblique}',\n", "# r'\\setmainfont[Extension=.ttf]{DejaVuSansCondensed}',\n", "# r'\\setmainfont[Extension=.otf]{FiraSans-Light}',\n", "# r'\\setsansfont[Extension=.otf]{FiraSans-Light}',\n", " ]\n", "}\n", "mpl.rcParams.update(pgf_with_custom_preamble)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def mean_std_map(pipe_home, meanmask, force=False, lazy=False, maskval=1000):\n", " pipe_std = pipe_home / 'summary_stdev.nii.gz'\n", " pipe_mean = pipe_home / 'summary_means.nii.gz'\n", "\n", " if force or not pipe_mean.is_file():\n", " print('Forced or %s not found' % pipe_mean)\n", " all_mus = []\n", " if lazy:\n", " all_mus = [nb.load(str(f)) for f in pipe_home.glob(\n", " 'sub-*/func/sub-*_task-stopsignal_bold_space-MNI152NLin2009cAsym_avgpreproc.nii.gz')]\n", " \n", " if not all_mus:\n", " print('Generating means file')\n", " pipe_files = list(pipe_home.glob(\n", " 'sub-*/func/sub-*_task-stopsignal_bold_space-MNI152NLin2009cAsym_preproc.nii.gz'))\n", " all_mus = []\n", " for f in pipe_files:\n", " mean = mean_img(str(f))\n", " data = mean.get_data()\n", " sigma = np.percentile(data[meanmask], 50) / maskval\n", " data /= sigma\n", " all_mus.append(nb.Nifti1Image(data, mean.affine, mean.header))\n", " \n", " meannii = concat_imgs(all_mus, auto_resample=False)\n", " meannii.to_filename(str(pipe_mean))\n", " force = True\n", "\n", " if force or not pipe_std.is_file():\n", " print('Generating standard deviation map')\n", " meannii = nb.load(str(pipe_mean))\n", " nb.Nifti1Image(meannii.get_data().std(3), meannii.affine, meannii.header).to_filename(str(pipe_std))\n", " \n", " return pipe_mean, pipe_std" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Use the WM mask to normalize intensities of EPI means\n", "meanmask = nb.load(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-02_class-WM_probtissue.nii.gz')).get_data() > 0.9\n", "\n", "# Calculate average and std\n", "fprep_mean, fprep_std = mean_std_map(fprep_home, meanmask)\n", "feat_mean, feat_std = mean_std_map(feat_home, meanmask)\n", "\n", "# Trick to avoid nilearn zooming in\n", "feat_nii = nb.load(str(feat_std))\n", "\n", "fd = feat_nii.get_data()\n", "fd[0, 0, :] = 50\n", "fd[0, -1, :] = 50\n", "fd[-1, 0, :] = 50\n", "fd[-1, -1, :] = 50\n", "\n", "nb.Nifti1Image(fd, feat_nii.affine, feat_nii.header).to_filename('newfeat.nii.gz')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "df = pd.read_csv(str(ANALYSIS_HOME / 'smoothness.csv'))\n", "\n", "plt.clf()\n", "fig = plt.gcf()\n", "_ = fig.set_size_inches(15, 2 * 3.1)\n", "# gs = gridspec.GridSpec(2, 4, width_ratios=[38, 7, 60, 10], height_ratios=[1, 1], hspace=0.0, wspace=0.03)\n", "gs = gridspec.GridSpec(2, 3, width_ratios=[42, 9, 64], height_ratios=[1, 1], hspace=0.0, wspace=0.03)\n", "\n", "a_ax1 = plt.subplot(gs[0, 0])\n", "a_ax2 = plt.subplot(gs[1, 0])\n", "\n", "\n", "fmriprep_smooth = df[df.pipeline.str.contains('fmriprep')][['fwhm_pre', 'fwhm_post']]\n", "feat_smooth = df[df.pipeline.str.contains('feat')][['fwhm_pre', 'fwhm_post']]\n", "\n", "cols = palettable.tableau.ColorBlind_10.hex_colors\n", "sn.distplot(fmriprep_smooth.fwhm_post, color=cols[0], ax=a_ax2,\n", " axlabel='Smoothing', label='fMRIPrep')\n", "sn.distplot(feat_smooth.fwhm_post, color=cols[1], ax=a_ax2,\n", " axlabel='Smoothing', label=r'\\texttt{feat}')\n", "sn.distplot(fmriprep_smooth.fwhm_pre, color=cols[0], ax=a_ax1,\n", " axlabel='Smoothing', label='fMRIPrep')\n", "sn.distplot(feat_smooth.fwhm_pre, color=cols[1], ax=a_ax1,\n", " axlabel='Smoothing', label=r'\\texttt{feat}')\n", "\n", "a_ax2.set_xlim([3, 8.8])\n", "a_ax2.set_xticks([])\n", "a_ax2.set_xticklabels([])\n", "a_ax2.xaxis.tick_bottom()\n", "a_ax2.grid(False)\n", "a_ax2.set_xlabel('')\n", "\n", "a_ax2.set_ylim([-1.1, 9.9])\n", "a_ax2.set_yticks([])\n", "a_ax2.spines['left'].set_visible(False)\n", "\n", "a_ax1.set_ylabel(r'\\noindent\\parbox{4.8cm}{\\centering\\textbf{Before smoothing} fraction of images}',\n", " size=13)\n", "a_ax1.yaxis.set_label_coords(-0.1, 0.4)\n", "a_ax2.set_ylabel(r'\\noindent\\parbox{4.8cm}{\\centering\\textbf{After smoothing} fraction of images}',\n", " size=13)\n", "a_ax2.yaxis.set_label_coords(-0.1, 0.6)\n", "\n", "# ax4.spines['bottom'].set_position(('outward', 20))\n", "a_ax2.invert_yaxis()\n", "a_ax2.spines['top'].set_visible(False)\n", "a_ax2.spines['bottom'].set_visible(False)\n", "a_ax2.spines['left'].set_visible(False)\n", "a_ax2.spines['right'].set_visible(False)\n", "\n", "a_ax1.set_xlim([3, 8.8])\n", "a_ax1.set_ylim([-0.6, 10.4])\n", "a_ax1.grid(False)\n", "\n", "a_ax1.set_xlabel('(mm)')\n", "a_ax1.xaxis.set_label_coords(0.95, 0.1)\n", "a_ax1.set_yticks([])\n", "a_ax1.set_yticklabels([])\n", "a_ax1.set_xticks([3, 4, 5, 6 , 7, 8])\n", "a_ax1.set_xticklabels([3, 4, 5, 6 , 7, 8])\n", "a_ax1.tick_params(axis='x', zorder=100, direction='inout')\n", "a_ax1.spines['left'].set_visible(False)\n", "a_ax1.spines['right'].set_visible(False)\n", "a_ax1.spines['top'].set_visible(False)\n", "a_ax1.spines['bottom'].set_visible(True)\n", "a_ax1.zorder = 100\n", "\n", "# a_ax2.xaxis.set_label_position('top')\n", "# a_ax2.xaxis.set_label_coords(0.45, 0.95)\n", "a_ax1.annotate(\n", " r'\\noindent\\parbox{6.8cm}{\\centering\\textbf{Estimated smoothness} full~width~half~maximum~(mm)}',\n", " xy=(0.15, 0.8), xycoords='axes fraction', xytext=(.0, .0),\n", " textcoords='offset points', va='center', color='k', size=13,\n", ")\n", "\n", "legend = a_ax2.legend(ncol=2, loc='upper center', bbox_to_anchor=(0.5, 0.45), prop={'size': 15})\n", "legend.get_frame().set_facecolor('w')\n", "legend.get_frame().set_edgecolor('none')\n", "\n", "# a_ax2.annotate(\n", "# r'\\noindent\\parbox{15cm}{Panels A, B present statistics derived from N=257 biologically independent participants}',\n", "# xy=(-0.2, 0.03), xycoords='axes fraction', xytext=(.0, .0),\n", "# textcoords='offset points', va='center', color='k', size=11,\n", "# )\n", "\n", "###### PANEL B\n", "\n", "b_ax1 = plt.subplot(gs[0, 2])\n", "b_ax2 = plt.subplot(gs[1, 2])\n", "\n", "\n", "thres = 20\n", "vmin = 50\n", "vmax = 200\n", "\n", "disp = plotting.plot_anat(str(fprep_std), display_mode='z', annotate=False,\n", " cut_coords=[-5, 10, 20], cmap='cividis', threshold=thres, vmin=vmin, vmax=vmax,\n", " axes=b_ax1)\n", "disp.annotate(size=12, left_right=True, positions=True)\n", "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_class-CSF_probtissue.nii.gz'), colors=['k'], levels=[0.8])\n", "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_class-WM_probtissue.nii.gz'), colors=['w'], levels=[0.8], linewidths=[1], alpha=0.7)\n", "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_brainmask.nii.gz'), colors=['k'], levels=[0.8], linewidths=[3], alpha=.7)\n", "\n", "\n", "disp = plotting.plot_anat('newfeat.nii.gz', display_mode='z', annotate=False,\n", " cut_coords=[-5, 10, 20], cmap='cividis', threshold=thres, vmin=vmin, vmax=vmax,\n", " axes=b_ax2)\n", "disp.annotate(size=12, left_right=False, positions=False)\n", "disp.annotate(size=12, left_right=False, positions=False, scalebar=True,\n", " loc=3, size_vertical=2, label_top=False, frameon=True, borderpad=0.1)\n", "\n", "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_class-CSF_probtissue.nii.gz'), colors=['k'], levels=[0.8])\n", "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_class-WM_probtissue.nii.gz'), colors=['w'], levels=[0.8], linewidths=[1], alpha=0.7)\n", "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_brainmask.nii.gz'), colors=['k'], levels=[0.8], linewidths=[3], alpha=.7)\n", "\n", "b_ax1.annotate(\n", " 'fMRIPrep',\n", " xy=(0., .5), xycoords='axes fraction', xytext=(-20, .0),\n", " textcoords='offset points', va='center', color='k', size=15,\n", " rotation=90)\n", "\n", "b_ax2.annotate(\n", " r'\\texttt{feat}',\n", " xy=(0., .5), xycoords='axes fraction', xytext=(-20, .0),\n", " textcoords='offset points', va='center', color='k', size=12,\n", " rotation=90)\n", "\n", "\n", "# inner_grid = gridspec.GridSpecFromSubplotSpec(1, 2, width_ratios=[1, 15],\n", "# subplot_spec=gs[:, -1], wspace=0.01)\n", "\n", "# b_ax3 = fig.add_subplot(inner_grid[0])\n", "# gradient = np.hstack((np.zeros((50,)), np.linspace(0, 1, 120), np.ones((130,))))[::-1]\n", "# gradient = np.vstack((gradient, gradient))\n", "# b_ax3.imshow(gradient.T, aspect='auto', cmap=plt.get_cmap('cividis'))\n", "# b_ax3.xaxis.set_ticklabels([])\n", "# b_ax3.xaxis.set_ticks([])\n", "# b_ax3.yaxis.set_ticklabels([])\n", "# b_ax3.yaxis.set_ticks([])\n", "\n", "# b_ax4 = fig.add_subplot(inner_grid[1])\n", "# sn.distplot(nb.load(str(fprep_std)).get_data()[mask2mm], label='fMRIPrep',\n", "# vertical=True, ax=b_ax4, kde=False, norm_hist=True)\n", "# sn.distplot(nb.load(str(feat_std)).get_data()[mask2mm], label=r'\\texttt{feat}', vertical=True,\n", "# color='darkorange', ax=b_ax4, kde=False, norm_hist=True)\n", "\n", "# # plt.gca().set_ylim((0, 300))\n", "# plt.legend(prop={'size': 15}, edgecolor='none')\n", "\n", "# b_ax4.xaxis.set_ticklabels([])\n", "# b_ax4.xaxis.set_ticks([])\n", "# b_ax4.yaxis.set_ticklabels([])\n", "# b_ax4.yaxis.set_ticks([])\n", "\n", "# plt.axis('off')\n", "# b_ax3.axis('off')\n", "# b_ax4.axis('off')\n", "\n", "\n", "a_ax1.set_title('A', fontdict={'fontsize': 24}, loc='left', x=-0.2);\n", "b_ax1.set_title('B', fontdict={'fontsize': 24}, loc='left');\n", "\n", "plt.savefig(str(out_folder / 'figure03.pdf'),\n", " format='pdf', bbox_inches='tight', pad_inches=0.2, dpi=300)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "coords = [-27, 0, 7]\n", "\n", "thres = 20\n", "vmin = 50\n", "vmax = 200\n", "\n", "\n", "# Plot\n", "plt.clf()\n", "fig = plt.figure(figsize=(20,10))\n", "\n", "plotting.plot_anat('newfeat.nii.gz', cut_coords=coords, colorbar=True, cmap='cividis',\n", " threshold=thres, vmin=vmin, vmax=vmax, title='feat',\n", " axes=plt.subplot(2,2,1)\n", ");\n", "plotting.plot_anat(str(fprep_std), cut_coords=coords, colorbar=True, cmap='cividis',\n", " threshold=thres, vmin=vmin, vmax=vmax, title='fmriprep',\n", " axes=plt.subplot(2,2,3)\n", ");\n", "plotting.plot_glass_brain(str(feat_std), threshold=200, colorbar=True, title='feat',\n", " axes=plt.subplot(2,2,2));\n", "plotting.plot_glass_brain(str(fprep_std), threshold=200, colorbar=True, title='fmriprep',\n", " axes=plt.subplot(2,2,4));" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }