Nilearn_tutorial_2018_2_Extracting_Features.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Section 2: Extracting signals from a brain parcellation\n",
    "============================================\n",
    "\n",
    "Here we show how to extract signals from a brain parcellation and compute\n",
    "a correlation matrix.\n",
    "\n",
    "We also show the importance of defining good confounds signals: the\n",
    "first correlation matrix is computed after regressing out simple\n",
    "confounds signals: movement regressors, white matter and CSF signals, ...\n",
    "The second one is without any confounds: all regions are connected to\n",
    "each other.\n",
    "\n",
    "\n",
    "One reference that discusses the importance of confounds is `Varoquaux and\n",
    "Craddock, Learning and comparing functional connectomes across subjects,\n",
    "NeuroImage 2013\n",
    "<http://www.sciencedirect.com/science/article/pii/S1053811913003340>`_.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the cell below in order to activate the tutorial download script (otherwise, ignore it!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import os.path\n",
    "\n",
    "from sklearn.datasets.base import Bunch\n",
    "from nilearn.datasets.utils import _get_dataset_dir, _fetch_files\n",
    "\n",
    "\n",
    "def fetch_data(n_subjects=30, data_dir=None, url=None, resume=True,\n",
    "               verbose=1):\n",
    "    \"\"\"Download and load the dataset.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    n_subjects: int, optional\n",
    "        The number of subjects to load from maximum of 40 subjects.\n",
    "        By default, 30 subjects will be loaded. If None is given,\n",
    "        all 40 subjects will be loaded.\n",
    "\n",
    "    data_dir: string, optional\n",
    "        Path of the data directory. Used to force data storage in a specified\n",
    "        location. Default: None\n",
    "\n",
    "    url: string, optional\n",
    "        Override download URL. Used for test only (or if you setup a mirror of\n",
    "        the data). Default: None\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    data: sklearn.datasets.base.Bunch\n",
    "        Dictionary-like object, the interest attributes are :\n",
    "         - 'func': Paths to functional resting-state images\n",
    "         - 'phenotypic': Explanations of preprocessing steps\n",
    "         - 'confounds': CSV files containing the nuisance variables\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    :Download:\n",
    "        https://openneuro.org/datasets/ds000228/versions/00001\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    if url is None:\n",
    "        url = 'https://openneuro.org/crn/datasets/ds000228/snapshots/00001/files/'\n",
    "\n",
    "    # Preliminary checks and declarations\n",
    "    dataset_name = 'ds000228'\n",
    "    data_dir = _get_dataset_dir(dataset_name, data_dir=data_dir,\n",
    "                                verbose=verbose)\n",
    "    max_subjects = 155\n",
    "    if n_subjects is None:\n",
    "        n_subjects = max_subjects\n",
    "    if n_subjects > max_subjects:\n",
    "        warnings.warn('Warning: there are only %d subjects' % max_subjects)\n",
    "        n_subjects = max_subjects\n",
    "    ids = range(1, n_subjects + 1)\n",
    "\n",
    "    # First, get the metadata\n",
    "    phenotypic = (\n",
    "            'participants.tsv',\n",
    "            url + 'participants.tsv', dict())\n",
    "\n",
    "    phenotypic_file = _fetch_files(data_dir, [phenotypic], resume=resume,\n",
    "                              verbose=verbose)[0]\n",
    "\n",
    "    # Load the csv file\n",
    "    phenotypic = np.genfromtxt(phenotypic_file, names=True, delimiter='\\t',\n",
    "                               dtype=None)\n",
    "\n",
    "    # Keep phenotypic information for selected subjects\n",
    "    int_ids = np.asarray(ids, dtype=int)\n",
    "    phenotypic = phenotypic[[i - 1 for i in int_ids]]\n",
    "\n",
    "    # Download dataset files\n",
    "\n",
    "    functionals = [\n",
    "        'derivatives:fmriprep:sub-pixar%03i:sub-pixar%03i_task-pixar_run-001_swrf_bold.nii.gz' % (i, i)\n",
    "        for i in ids]\n",
    "    urls = [url + name for name in functionals]\n",
    "    functionals = _fetch_files(\n",
    "        data_dir, zip(functionals, urls, (dict(),) * n_subjects),\n",
    "        resume=resume, verbose=verbose)\n",
    "\n",
    "    confounds = [\n",
    "        'derivatives:fmriprep:sub-pixar%03i:sub-pixar%03i_task-pixar_run-001_ART_and_CompCor_nuisance_regressors.mat'\n",
    "        % (i, i)\n",
    "        for i in ids]\n",
    "    confound_urls = [url + name for name in confounds]\n",
    "\n",
    "    confounds = _fetch_files(\n",
    "        data_dir, zip(confounds, confound_urls, (dict(),) * n_subjects),\n",
    "        resume=resume, verbose=verbose)\n",
    "\n",
    "    return Bunch(func=functionals, confounds=confounds,\n",
    "                 phenotypic=os.path.join(data_dir, phenotypic_file),\n",
    "                 description='ds000228')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retrieve the atlas and the data\n",
    "--------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we're using rs-fmri data, it makes sense to use an atlas defined using rs-fmri data\n",
    "\n",
    "Let's use the MIST atlas, created here in Montreal using the BASC method. This atlas has multiple resolutions, for larger networks or finer-grained ROIs. Let's use a 64-ROI atlas to allow some detail, but to ultimately keep our connectivity matrices manageable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import datasets\n",
    "\n",
    "parcellations = datasets.fetch_atlas_basc_multiscale_2015(version='sym')\n",
    "atlas_filename = parcellations.scale064\n",
    "\n",
    "\n",
    "print('Atlas ROIs are located in nifti image (4D) at: %s' %\n",
    "       atlas_filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at that atlas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import plotting\n",
    "\n",
    "plotting.plot_roi(atlas_filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, let's load an example 4D fmri time-series for one subject\n",
    "\n",
    "We have prepared some data especially for this tutorial. It is based on an open dataset of children and young adults. More details can be found here:\n",
    "\n",
    "https://openneuro.org/datasets/ds000228/versions/00001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have already downloaded the data, this next line should happen instantly. Otherwise, the data download will begin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# One subject of resting-state data\n",
    "data = fetch_data(n_subjects=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmri_filenames = data.func[0]\n",
    "print('fmri timeseries are located in nifti image (4D) at: %s' %\n",
    "       fmri_filenames)  # 4D data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's have a look at that 4D resting-state image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_stat_map(fmri_filenames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uh oh! An error! What's the problem here?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's just a look at a single frame\n",
    "# First we load the image into memory\n",
    "from nilearn import image\n",
    "\n",
    "myImg = image.load_img(fmri_filenames)\n",
    "print(myImg.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now, let's isolate the first frame\n",
    "first_frame = image.index_img(myImg, 0)\n",
    "print(first_frame.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's try to plot it\n",
    "plotting.plot_stat_map(first_frame)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first frame can be a bit wonky. What about an average image?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "averaged_Img = image.mean_img(image.mean_img(myImg))\n",
    "averaged_Img.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_stat_map(averaged_Img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract signals on a parcellation defined by labels\n",
    "Using the NiftiLabelsMasker\n",
    "\n",
    "So we've loaded our atlas and 4D data for a single subject. Let's practice extracting features!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# here is a short script to deal with the confounds in this particular dataset\n",
    "# don't worry about the details\n",
    "# but here is an example of how you can create your own function!\n",
    "\n",
    "import h5py\n",
    "import numpy as np\n",
    "\n",
    "def prepare_confounds(conf, key = 'R', transpose=True):\n",
    "    arrays = {}\n",
    "    f = h5py.File(conf)\n",
    "    for k, v in f.items():\n",
    "        arrays[k] = np.array(v)\n",
    "    \n",
    "    if transpose:\n",
    "        output = arrays[key].T\n",
    "    else:\n",
    "        output = arrays[key]\n",
    "    \n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiLabelsMasker\n",
    "\n",
    "masker = NiftiLabelsMasker(labels_img=atlas_filename, standardize=True, \n",
    "                           memory='nilearn_cache', verbose=5)\n",
    "\n",
    "# Here we go from nifti files to the signal time series in a numpy\n",
    "# array. Note how we give confounds to be regressed out during signal\n",
    "# extraction\n",
    "conf = prepare_confounds(data.confounds[0])\n",
    "time_series = masker.fit_transform(myImg, confounds=conf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what did we just create here?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(time_series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Its a 168 (timeframes) x 64 (region) array!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are these \"confounds\" and how are they used?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.confounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conf.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import DataFrame\n",
    "DataFrame(conf).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute and display a correlation matrix\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.connectome import ConnectivityMeasure\n",
    "\n",
    "correlation_measure = ConnectivityMeasure(kind='correlation')\n",
    "correlation_matrix = correlation_measure.fit_transform([time_series])[0]\n",
    "correlation_matrix.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the correlation matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# Mask the main diagonal for visualization:\n",
    "np.fill_diagonal(correlation_matrix, 0)\n",
    "\n",
    "# The labels we have start with the background (0), hence we skip the\n",
    "# first label\n",
    "plotting.plot_matrix(correlation_matrix, figure=(10, 8), \n",
    "                     labels=range(time_series.shape[-1]),\n",
    "                     vmax=0.8, vmin=-0.8, reorder=False)\n",
    "\n",
    "# matrices are ordered for block-like representation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Same thing without confounds, to stress the importance of confounds\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series = masker.fit_transform(myImg)\n",
    "# Note how we did not specify confounds above. This is bad!\n",
    "\n",
    "correlation_matrix = correlation_measure.fit_transform([time_series])[0]\n",
    "\n",
    "# Mask the main diagonal for visualization:\n",
    "np.fill_diagonal(correlation_matrix, 0)\n",
    "\n",
    "plotting.plot_matrix(correlation_matrix, figure=(10, 8), \n",
    "                     labels=range(time_series.shape[-1]),\n",
    "                     vmax=0.8, vmin=-0.8,\n",
    "                    title='No confounds', reorder=False)\n",
    "\n",
    "plotting.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise** We generated a connectivity matrix with one subject. Try doing the same thing with a second subject. Are the connectomes similar? What about the influence of confounds?\n",
    "\n",
    "Hint, to get a second subject, change to n=2 on the download function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## WORKSPACE\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}