neuroimaging_in_python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with neuroimaging data in Python\n",
    "\n",
    "We've already learned about general Python tools, but there is also a growing ecosystem of Python-based tools for working with neuroimaging data.\n",
    "\n",
    "In this training, we will focus on `pybids`, `nibabel`, and `nilearn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os.path as op  # for navigating the system\n",
    "from pprint import pprint  # for pretty-printing dictionaries\n",
    "import numpy as np\n",
    "\n",
    "dataset_folder = '/scratch/madlab/condassoc_bids/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PyBIDS\n",
    "[PyBIDS](https://bids-standard.github.io/pybids/) is a package for **working with BIDS datasets**.\n",
    "\n",
    "With PyBIDS, you can:\n",
    "1. Validate datasets to make sure they're BIDS compliant.\n",
    "2. Search datasets for specific files.\n",
    "3. Read in metadata for files within datasets.\n",
    "4. Generate citable reports summarizing datasets.\n",
    "5. Apply transformations to variables in datasets for analysis.\n",
    "\n",
    "In this training we will cover 1-3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bids import BIDSLayout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in BIDS dataset as BIDSLayout object\n",
    "layout = BIDSLayout(dataset_folder, validate=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First let's look at the dataset description file to learn about the dataset\n",
    "pprint(layout.get_dataset_description())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can use the BIDSLayout object to search the data\n",
    "# for different kinds of information\n",
    "subjects = layout.get_subjects()\n",
    "print('{} subjects in dataset'.format(len(subjects)))\n",
    "print('Subject IDs: {}'.format(', '.join(subjects)))\n",
    "print()\n",
    "\n",
    "tasks = layout.get_tasks()\n",
    "print('{} tasks in dataset'.format(len(tasks)))\n",
    "print('Tasks: {}'.format(', '.join(tasks)))\n",
    "print()\n",
    "\n",
    "data_types = layout.get_datatypes()\n",
    "print('{} data types in dataset'.format(len(data_types)))\n",
    "print('Data types: {}'.format(', '.join(data_types)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you're confused about what you can search for,\n",
    "# check the 'entities' attribute\n",
    "# Each entity can be search with a \"get_<entity>s\" function\n",
    "layout.entities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Each search function can take other entities as arguments\n",
    "# For example, what repetition times exist in the dataset?\n",
    "print(layout.get_RepetitionTimes())\n",
    "\n",
    "# Ah, but what about for functional scans only?\n",
    "print(layout.get_RepetitionTimes(datatype='func'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# There is also the general search method `get` that returns BIDS-type file object\n",
    "sub_runs = layout.get(subject='001', task='condassoc', extension='nii.gz')\n",
    "pprint(sub_runs)\n",
    "\n",
    "# Notice that these files are stored as BIDSImageFiles, a custom class created for pybids.\n",
    "# We generally just want one thing directly from these objects- the file path.\n",
    "path = sub_runs[0].path\n",
    "print()\n",
    "print('The file path for the first file:')\n",
    "print(path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's get the metadata for one of the functional scans\n",
    "func_file = sub_runs[0].path\n",
    "pprint(layout.get_metadata(func_file))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "Load the dataset located in `/scratch/cis-training/ds001491` and print the name of the file with the **extension** `'nii.gz'` for **echo** `1` of the `'images'` **task** of **subject** `'01'`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "Use `BIDSLayout.get_metadata` to determine the EchoTime for that file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NiBabel\n",
    "[NiBabel](https://nipy.org/nibabel/) allows you to **read and write neuroimaging files**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nibabel as nib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab a random nifti file from the dataset\n",
    "f = layout.get(extension='nii.gz', task='condassoc')[0].path\n",
    "print(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the file as an image object\n",
    "img = nib.load(f)\n",
    "print(type(img))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Each image has an affine matrix that contains information\n",
    "# about the orientation and dimensionality of the data\n",
    "print(img.affine)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The image header contains a lot of its own information\n",
    "hed = img.header\n",
    "print(type(hed))\n",
    "\n",
    "# get_zooms() prints the voxel sizes and\n",
    "# repetition time (fourth dimension) for the image\n",
    "hed.get_zooms()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The data are a 3D-to-4D array\n",
    "data = img.get_data()\n",
    "print(type(data))\n",
    "print(data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can also create images from an affine matrix and a numpy array\n",
    "random_data = np.random.random((96, 96, 42))\n",
    "random_img = nib.Nifti1Image(random_data, img.affine)\n",
    "print(type(random_img))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can also save image objects\n",
    "random_img.to_filename('example_random_image.nii.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3\n",
    "Compute the mean of the data from the loaded functional scan and generate a new 3D image with the mean data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nilearn\n",
    "[Nilearn](http://nilearn.github.io) is primarily for doing machine learning on fMRI data, but it also provides a wide range of functionality, including:\n",
    "- Applying masks to data\n",
    "- Visualizing fMRI data\n",
    "- Simple denoising and analysis\n",
    "- Performing simple math on images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nilearn as nl\n",
    "from nilearn import plotting, image, masking\n",
    "\n",
    "first_volume = image.index_img(f, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### nilearn.masking\n",
    "This submodule can be used to create masks, apply masks to data, and unmask masked data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import masking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab a random functional file from the dataset\n",
    "f = layout.get(suffix='bold', extension='nii.gz')[0].path\n",
    "print(f)\n",
    "print('The shape of the data in the file is {}'.format(nib.load(f).shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate a general brain mask from the functional data\n",
    "mask = masking.compute_epi_mask(f)\n",
    "\n",
    "# Let's plot it against the first volume in the functional scan\n",
    "fig, ax = plt.subplots(figsize=(16, 6))\n",
    "plotting.plot_roi(mask, bg_img=first_volume, axes=ax)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the mask to the file\n",
    "data = masking.apply_mask(f, mask)\n",
    "\n",
    "print('Shape of the functional data: {}'.format(nib.load(f).shape))\n",
    "print('Shape of the mask: {}'.format(mask.shape))\n",
    "print('Number of voxels within mask: {}'.format(np.sum(mask.get_data())))\n",
    "print('Shape of the masked data: {}'.format(data.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots(figsize=(16, 5))\n",
    "ax.imshow(data.T, aspect='auto')\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_ylabel('Voxel')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's plot one voxel's time series\n",
    "fig, ax = plt.subplots(figsize=(16, 5))\n",
    "ax.plot(data[:, 0])\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_ylabel('BOLD signal (arbitrary unit)')\n",
    "ax.set_xlim(0, data.shape[0])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### nilearn.image\n",
    "This submodule can perform operations (e.g., math, thresholding, resampling) on image files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can use index_img to grab a specific\n",
    "# volume within a 4D dataset\n",
    "first_volume = image.index_img(f, 0)\n",
    "print(first_volume.shape)\n",
    "\n",
    "# Plot the first volume of the functional data\n",
    "fig, ax = plt.subplots(figsize=(16, 6))\n",
    "plotting.plot_epi(first_volume, axes=ax, cut_coords=[0, 0, 20])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also compute the mean image\n",
    "mean_img = image.mean_img(f)\n",
    "print(mean_img.shape)\n",
    "\n",
    "# Plot the mean image\n",
    "fig, ax = plt.subplots(figsize=(16, 6))\n",
    "plotting.plot_epi(mean_img, axes=ax, cut_coords=[0, 0, 20])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Or apply a threshold\n",
    "mean_img_thresh = image.threshold_img(mean_img, np.mean(masking.apply_mask(mean_img, mask)))\n",
    "\n",
    "# Plot the first volume of the functional data\n",
    "fig, ax = plt.subplots(figsize=(16, 6))\n",
    "plotting.plot_epi(mean_img_thresh, axes=ax, cut_coords=[0, 0, 20])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Or smooth an image\n",
    "smooth_mean_img_thresh = image.smooth_img(mean_img_thresh, 8)\n",
    "\n",
    "# Plot the first volume of the functional data\n",
    "fig, ax = plt.subplots(figsize=(16, 6))\n",
    "plotting.plot_epi(smooth_mean_img_thresh, axes=ax, cut_coords=[0, 0, 20])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4\n",
    "Use [nilearn.image.math_img](http://nilearn.github.io/modules/generated/nilearn.image.math_img.html#nilearn.image.math_img) to square the last volume in the functional dataset, then plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### nilearn.datasets\n",
    "Nilearn has the ability to download a range of datasets and useful images from the internet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's show what functions there are in this module\n",
    "from inspect import getmembers, isfunction\n",
    "\n",
    "functions_list = [o[0] for o in getmembers(datasets) if isfunction(o[1])]\n",
    "print('\\n'.join(functions_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can fetch atlases\n",
    "aal = datasets.fetch_atlas_aal()\n",
    "\n",
    "# Let's plot this atlas\n",
    "fig, ax = plt.subplots(figsize=(16, 6))\n",
    "plotting.plot_roi(aal['maps'], axes=ax, cut_coords=[0, 0, 20])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also fetch and plot *probabilistic* atlases, like MSDL\n",
    "atlas = datasets.fetch_atlas_msdl()\n",
    "\n",
    "# Let's plot this atlas\n",
    "fig, ax = plt.subplots(figsize=(16, 6))\n",
    "plotting.plot_prob_atlas(atlas['maps'], axes=ax, cut_coords=[0, 0, 20])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5\n",
    "Fetch and plot the `'sub-maxprob-thr25-2mm'` version of the Harvard-Oxford atlas. Determine if it's a probabilistic atlas (atlas is a 4D image) or a deterministic one (atlas is a 3D image).\n",
    "Hint: Use the [documentation](http://nilearn.github.io/modules/generated/nilearn.datasets.fetch_atlas_harvard_oxford.html#nilearn.datasets.fetch_atlas_harvard_oxford) to figure out what arguments are required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other packages\n",
    "Some other neuroimaging-related Python packages to know include:\n",
    "- nistats: Modeling and GLMs for fMRI\n",
    "- Nipype: Python wrappers for non-Python neuroimaging tools\n",
    "- heudiconv: Automated conversion of scanner data to BIDS format\n",
    "- DiPy: Diffusion-weighted imaging processing and analysis in Python"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}