{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nTutorial 0: Preparing your data for gradient analysis\n=====================================================\nIn this example, we will introduce how to preprocess raw MRI data and how\nto prepare it for subsequent gradient analysis in the next tutorials.\n\nRequirements\n------------\nFor this tutorial, you will need to install the Python package\n`load_confounds `_. You can\ndo it using ``pip``::\n\n pip install load_confounds\n \n\nPreprocessing\n-------------\nBegin with an MRI dataset that is organized in `BIDS\n`_ format. We recommend preprocessing your data\nusing `fmriprep `_, as described below, but\nany preprocessing pipeline will work.\n\nFollowing is example code to run `fmriprep `_\nusing docker from the command line::\n\n docker run -ti --rm \\\n -v :/data:ro \\\n -v :/out poldracklab/fmriprep:latest \\\n --output-spaces fsaverage5 \\\n --fs-license-file license.txt \\\n /data /out participant\n\n

Note

For this tutorial, it is crucial to output the data onto a cortical surface\n template space.

\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import the dataset as timeseries\n++++++++++++++++++++++++++++++++\nThe timeseries should be a numpy array with the dimensions: nodes x timepoints \n\nFollowing is an example for reading in data:: \n\n import nibabel as nib\n import numpy as np\n\n filename = 'filename.{}.mgz' # where {} will be replaced with 'lh' and 'rh'\n timeseries = [None] * 2\n for i, h in enumerate(['lh', 'rh']):\n timeseries[i] = nib.load(filename.format(h)).get_fdata().squeeze()\n timeseries = np.vstack(timeseries)\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a **working example**, simply fetch timeseries:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from brainspace.datasets import fetch_timeseries_preprocessing\ntimeseries = fetch_timeseries_preprocessing()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Confound regression\n++++++++++++++++++++++++\nTo remove confound regressors from the output of the fmriprep pipeline, first\nextract the confound columns. For example::\n\n import load_confounds\n confounds_out = load_confounds(\"path to confound file\",\n strategy='minimal',\n n_components=0.95,\n motion_model='6params')\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a **working example**, simply read in confounds\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from brainspace.datasets import load_confounds_preprocessing\nconfounds_out = load_confounds_preprocessing()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do the confound regression\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from nilearn import signal\nclean_ts = signal.clean(timeseries.T, confounds=confounds_out).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And extract the cleaned timeseries onto a set of labels\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nfrom nilearn import datasets\nfrom brainspace.utils.parcellation import reduce_by_labels\n\n# Fetch surface atlas\natlas = datasets.fetch_atlas_surf_destrieux()\n\n# Remove non-cortex regions\nregions = atlas['labels'].copy()\nmasked_regions = [b'Medial_wall', b'Unknown']\nmasked_labels = [regions.index(r) for r in masked_regions]\nfor r in masked_regions:\n regions.remove(r)\n\n# Build Destrieux parcellation and mask\nlabeling = np.concatenate([atlas['map_left'], atlas['map_right']])\nmask = ~np.isin(labeling, masked_labels)\n\n# Distinct labels for left and right hemispheres\nlab_lh = atlas['map_left']\nlabeling[lab_lh.size:] += lab_lh.max() + 1\n\n# extract mean timeseries for each label\nseed_ts = reduce_by_labels(clean_ts[mask], labeling[mask], axis=1, red_op='mean')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate functional connectivity matrix\n++++++++++++++++++++++++++++++++++++++++\nThe following example uses\n`nilearn `_:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from nilearn.connectome import ConnectivityMeasure\n\ncorrelation_measure = ConnectivityMeasure(kind='correlation')\ncorrelation_matrix = correlation_measure.fit_transform([seed_ts.T])[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the correlation matrix:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from nilearn import plotting\n\n# Reduce matrix size, only for visualization purposes\nmat_mask = np.where(np.std(correlation_matrix, axis=1) > 0.2)[0]\nc = correlation_matrix[mat_mask][:, mat_mask]\n\n# Create corresponding region names\nregions_list = ['%s_%s' % (h, r.decode()) for h in ['L', 'R'] for r in regions]\nmasked_regions = [regions_list[i] for i in mat_mask]\n\n\ncorr_plot = plotting.plot_matrix(c, figure=(15, 15), labels=masked_regions,\n vmax=0.8, vmin=-0.8, reorder=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run gradient analysis and visualize\n+++++++++++++++++++++++++++++++++++\n\nRun gradient analysis\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from brainspace.gradient import GradientMaps\n\ngm = GradientMaps(n_components=2, random_state=0)\ngm.fit(correlation_matrix)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visualize results\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from brainspace.datasets import load_fsa5\nfrom brainspace.plotting import plot_hemispheres\nfrom brainspace.utils.parcellation import map_to_labels\n\n# Map gradients to original parcels\ngrad = [None] * 2\nfor i, g in enumerate(gm.gradients_.T):\n grad[i] = map_to_labels(g, labeling, mask=mask, fill=np.nan)\n\n\n# Load fsaverage5 surfaces\nsurf_lh, surf_rh = load_fsa5()\n\n# sphinx_gallery_thumbnail_number = 2\nplot_hemispheres(surf_lh, surf_rh, array_name=grad, size=(1200, 400), cmap='viridis_r',\n color_bar=True, label_text=['Grad1', 'Grad2'], zoom=1.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This concludes the setup tutorial. The following tutorials can be run using\neither the output generated here or the example data.\n\n" ] } ], "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.7.6" } }, "nbformat": 4, "nbformat_minor": 0 }