MAIN_tutorial_intro_to_nilearn.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import nilearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's keep our notebook clean, so it's a little more readable!\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Understanding neuroimaging data](http://nilearn.github.io/manipulating_images/input_output.html)\n",
    "\n",
    "### Text files: phenotype or behavior\n",
    "\n",
    "Phenotypic or behavioral data are often provided as a text or CSV (Comma Separated Values) file. They can be loaded with the [pandas package](https://pandas.pydata.org/) but you may have to specify some options. Here, we'll specify the `sep` field, since our data is tab-delimited rather than comma-delimited.\n",
    "\n",
    "For our dataset, let's load participant level information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "\n",
    "data_dir = '/home/emdupre/Desktop/MAIN_tutorial/'  # Update this path to match your local machine !\n",
    "participants = 'participants.tsv'\n",
    "phenotypic_data = pd.read_csv(os.path.join(data_dir, participants), sep='\\t')\n",
    "phenotypic_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nifti data\n",
    "\n",
    "For volumetric data, nilearn works with data stored in the [Nifti structure](http://nipy.org/nibabel/nifti_images.html) (via the [nibabel package](http://nipy.org/nibabel/)).\n",
    "\n",
    "The NifTi data structure (also used in Analyze files) is the standard way of sharing data in neuroimaging research. Three main components are:\n",
    "\n",
    "  * data:\traw scans in form of a numpy array:  \n",
    "    `data = img.get_data()`\n",
    "  * affine:\treturns the transformation matrix that maps from voxel indices of the `numpy` array to actual real-world     locations of the brain:  \n",
    "    `affine = img.affine`\n",
    "  * header:\tlow-level informations about the data (slice duration, etc.):  \n",
    "  `header = img.header`\n",
    "\n",
    "It is important to appreciate that the representation of MRI data we'll be using is a big 4D matrix representing (3D MRI + 1D for time), stored in a single Nifti file.\n",
    "\n",
    "### Niimg-like objects\n",
    "\n",
    "Nilearn functions take as input argument what we call \"Niimg-like objects\":\n",
    "\n",
    "Niimg: A Niimg-like object can be one of the following:\n",
    "\n",
    "  * A string with a file path to a Nifti image\n",
    "  * A SpatialImage from `nibabel`, i.e., an object exposing the get_data() method and affine attribute, typically a Nifti1Image from `nibabel`.\n",
    "\n",
    "Niimg-4D: Similarly, some functions require 4D Nifti-like data, which we call Niimgs or Niimg-4D. Accepted input arguments are:\n",
    "\n",
    "  * A path to a 4D Nifti image\n",
    "  * List of paths to 3D Nifti images\n",
    "  * 4D Nifti-like object\n",
    "  * List of 3D Nifti-like objects\n",
    "\n",
    "**Note:** If you provide a sequence of Nifti images, all of them must have the same affine !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Manipulating and looking at data](http://nilearn.github.io/auto_examples/plot_nilearn_101.html#sphx-glr-auto-examples-plot-nilearn-101-py)\n",
    "\n",
    "There is a whole section of the [Nilearn documentation](http://nilearn.github.io/plotting/index.html#plotting) on making pretty plots for neuroimaging data ! But let's start with a simple one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's use a Nifti file that is shipped with nilearn\n",
    "from nilearn import datasets\n",
    "\n",
    "# Note that the variable MNI152_FILE_PATH is just a path to a Nifti file\n",
    "print('Path to MNI152 template: {}'.format(datasets.MNI152_FILE_PATH))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above, MNI152_FILE_PATH is nothing more than a string with a path pointing to a nifti image. You can replace it with a string pointing to a file on your disk. Note that it should be a 3D volume, and not a 4D volume."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import plotting\n",
    "plotting.plot_img(datasets.MNI152_FILE_PATH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also directly manipulate these images using Nilearn ! As an example, let's try smoothing this image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import image\n",
    "smooth_anat_img = image.smooth_img(datasets.MNI152_FILE_PATH, fwhm=6)\n",
    "\n",
    "# While we are giving a file name as input, the function returns\n",
    "# an in-memory object:\n",
    "print(smooth_anat_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_img(smooth_anat_img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then save this manipulated image from in-memory to disk as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smooth_anat_img.to_filename('smooth_anat_img.nii.gz')\n",
    "os.getcwd()  # We'll' check our \"current working directory\" (cwd) to see where the file was saved"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Visualizing neuroimaging volumes](https://nilearn.github.io/auto_examples/01_plotting/plot_visualization.html#visualization)\n",
    "\n",
    "What if we want to view not a structural MRI image, but a functional one ?\n",
    "No problem ! Let's try loading one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmri_filename = 'downsampled_derivatives:fmriprep:sub-pixar109:sub-pixar109_task-pixar_run-001_swrf_bold.nii.gz'\n",
    "plotting.plot_epi(os.path.join(data_dir, fmri_filename))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uh-oh, what happened ?! Let's look back at the error message:\n",
    "\n",
    "> DimensionError: Input data has incompatible dimensionality: Expected dimension is 3D and you provided a 4D image. See http://nilearn.github.io/manipulating_images/input_output.html.\n",
    "\n",
    "We can fix that ! Let's take an average of the EPI image and plot that instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.image import mean_img\n",
    "\n",
    "plotting.view_img(mean_img(os.path.join(data_dir, fmri_filename)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Convert the fMRI volumes to a data matrix](http://nilearn.github.io/auto_examples/plot_decoding_tutorial.html#convert-the-fmri-volume-s-to-a-data-matrix)\n",
    "\n",
    "These are some really lovely images, but for machine learning we want matrices so that we can use all of the techniques we learned this morning !\n",
    "\n",
    "To transform our Nifti images into matrices, we'll use the `nilearn.input_data.NiftiMasker` to extract the fMRI data from a mask and convert it to data series.\n",
    "\n",
    "First, let's do the simplest possible mask&mdash;a mask of the whole brain. We'll use a mask that ships with Nilearn and matches the MNI152 template we plotted earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brain_mask = datasets.load_mni152_brain_mask()\n",
    "plotting.plot_roi(brain_mask, cmap='Paired')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiMasker\n",
    "masker = NiftiMasker(mask_img=brain_mask, standardize=True)\n",
    "masker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We give the masker a filename and retrieve a 2D array ready\n",
    "# for machine learning with scikit-learn !\n",
    "fmri_masked = masker.fit_transform(os.path.join(data_dir, fmri_filename))\n",
    "print(fmri_masked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fmri_masked.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to think about what just happened is to look at it visually:\n",
    "\n",
    "![](http://nilearn.github.io/_images/masking.jpg)\n",
    "\n",
    "Essentially, we can think about overlaying a 3D grid on an image. Then, our mask tells us which cubes or \"voxels\" (like 3D pixels) to sample from. Since our Nifti images are 4D files, we can't overlay a single grid -- instead, we use a series of 3D grids (one for each volume in the 4D file), so we can get a measurement for each voxel at each timepoint. These are reflected in the shape of the matrix ! You can check this by checking the number of positive voxels in our brain mask.\n",
    "\n",
    "There are many other strategies in Nilearn [for masking data and for generating masks](http://nilearn.github.io/manipulating_images/manipulating_images.html#computing-and-applying-spatial-masks). I'd encourage you to spend some time exploring the documentation for these !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also [display this time series](http://nilearn.github.io/auto_examples/03_connectivity/plot_adhd_spheres.html#display-time-series) to get an intuition of how the whole brain signal is changing over time.\n",
    "\n",
    "We'll display the first three voxels by sub-selecting values from the matrix. You can also find more information on [how to slice arrays here](https://docs.scipy.org/doc/numpy-1.13.0/reference/arrays.indexing.html#basic-slicing-and-indexing)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(fmri_masked[5:150, :3])\n",
    "\n",
    "plt.title('Voxel Time Series')\n",
    "plt.xlabel('Scan number')\n",
    "plt.ylabel('Normalized signal')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Extracting signals from a brain parcellation](http://nilearn.github.io/auto_examples/03_connectivity/plot_signal_extraction.html#extracting-signals-from-a-brain-parcellation)\n",
    "\n",
    "Now that we've seen how to create a data series from a single region-of-interest (ROI), we can start to scale up ! What if, instead of wanting to extract signal from one ROI, we want to define several ROIs and extract signal from all of them  ? Nilearn can help us with that, too ! 🎉\n",
    "\n",
    "For this, we'll use `nilearn.input_data.NiftiLabelsMasker`. `NiftiLabelsMasker` which works like `NiftiMasker` except that it's for labelled data rather than binary. That is, since we have more than one ROI, we need more than one value ! Now that each ROI gets its own value, these values are treated as labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, let's load a parcellation that we'd like to use\n",
    "multiscale = datasets.fetch_atlas_basc_multiscale_2015()\n",
    "print('Atlas ROIs are located at: %s' % multiscale.scale064)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_roi(multiscale.scale064)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiLabelsMasker\n",
    "label_masker = NiftiLabelsMasker(labels_img=multiscale.scale064, standardize=True)\n",
    "label_masker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmri_matrix = label_masker.fit_transform(os.path.join(data_dir, fmri_filename))\n",
    "print(fmri_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fmri_matrix.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [Compute and display a correlation matrix](http://nilearn.github.io/auto_examples/03_connectivity/plot_signal_extraction.html#compute-and-display-a-correlation-matrix)\n",
    "\n",
    "Now that we have a matrix, we'd like to create a _connectome_. A connectome is a map of the connections in the brain. Since we're working with functional data, however, we don't have access to actual connections. Instead, we'll use a measure of statistical dependency to infer the (possible) presence of a connection.\n",
    "\n",
    "Here, we'll use Pearson's correlation as our measure of statistical dependency and compare how all of our ROIs from our chosen parcellation relate to one another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import connectome\n",
    "correlation_measure = connectome.ConnectivityMeasure(kind='correlation')\n",
    "correlation_measure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation_matrix = correlation_measure.fit_transform([fmri_matrix])\n",
    "correlation_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "correlation_matrix = correlation_matrix[0]\n",
    "# Mask the main diagonal for visualization:\n",
    "# np.fill_diagonal(correlation_matrix, 0)\n",
    "plotting.plot_matrix(correlation_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [The importance of specifying confounds](http://nilearn.github.io/auto_examples/03_connectivity/plot_signal_extraction.html#same-thing-without-confounds-to-stress-the-importance-of-confounds)\n",
    "\n",
    "In fMRI, we're collecting a noisy signal. We have artifacts like physiological noise (from heartbeats, respiration) and head motion which can impact our estimates. Therefore, it's strongly recommended that you control for these and related measures when deriving your connectome measures. Here, we'll repeat the correlation matrix example, but this time we'll control for confounds. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conf_filename = 'derivatives:fmriprep:sub-pixar109:sub-pixar109_task-pixar_run-001_ART_and_CompCor_nuisance_regressors.tsv'\n",
    "clean_fmri_matrix = label_masker.fit_transform(os.path.join(data_dir, fmri_filename),\n",
    "                                               confounds=os.path.join(data_dir, conf_filename))\n",
    "clean_correlation_matrix = correlation_measure.fit_transform([clean_fmri_matrix])[0]\n",
    "np.fill_diagonal(clean_correlation_matrix, 0)\n",
    "plotting.plot_matrix(clean_correlation_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks a little different !\n",
    "\n",
    "Looking more closely, we can see that our correlation matrix is symmetrical; that is, that both sides of the diagonal contain the same information. We don't want to feed duplicate information into our machine learning classifier, and Nilearn has a really easy way to remove this redundancy ! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorized_correlation = connectome.ConnectivityMeasure(kind='correlation',\n",
    "                                                        vectorize=True, discard_diagonal=True)\n",
    "clean_vectorized_correlation = vectorized_correlation.fit_transform([clean_fmri_matrix])[0]\n",
    "clean_vectorized_correlation.shape  # Why is this value not 64 * 64 ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Interactive connectome plotting](http://nilearn.github.io/plotting/index.html#d-plots-of-connectomes)\n",
    "\n",
    "It can also be helpful to project these connection weightings back on to the brain, to visualize these connectomes ! Here, we'll use the interactive connectome plotting in Nilearn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = plotting.find_parcellation_cut_coords(multiscale.scale064)\n",
    "plotting.view_connectome(clean_correlation_matrix, coords=coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "?plotting.view_connectome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}