{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Clustering methods to learn a brain parcellation from fMRI\n\nWe use spatially-constrained Ward-clustering, KMeans, and Recursive Neighbor\nAgglomeration (ReNA) to create a set of parcels.\n\nIn a high dimensional regime, these methods can be interesting\nto create a 'compressed' representation of the data, replacing the data\nin the fMRI images by mean signals on the parcellation, which can\nsubsequently be used for statistical analysis or machine learning.\n\nAlso, these methods can be used to learn functional connectomes\nand subsequently for classification tasks.\n\n## References\n\nWhich clustering method to use, an empirical comparison can be found in this\npaper\n\n * Bertrand Thirion, Gael Varoquaux, Elvis Dohmatob, Jean-Baptiste Poline.\n `Which fMRI clustering gives good brain parcellations ?\n `_ Frontiers in Neuroscience,\n 2014.\n\nThis parcellation may be useful in a supervised learning, see for\ninstance\n\n * Vincent Michel, Alexandre Gramfort, Gael Varoquaux, Evelyn Eger,\n Christine Keribin, Bertrand Thirion. `A supervised clustering approach\n for fMRI-based inference of brain states.\n `_.\n Pattern Recognition, Elsevier, 2011.\n\nThe big picture discussion corresponding to this example can be found\nin the documentation section `parcellating_brain`.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Download a brain development fmri dataset and turn it to a data matrix\n\nWe download one subject of the movie watching dataset from Internet\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from nilearn import datasets\ndataset = datasets.fetch_development_fmri(n_subjects=1)\n\n# print basic information on the dataset\nprint('First subject functional nifti image (4D) is at: %s' %\n dataset.func[0]) # 4D data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Brain parcellations with Ward Clustering\n\nTransforming list of images to data matrix and build brain parcellations,\nall can be done at once using `Parcellations` object.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from nilearn.regions import Parcellations\n\n# Computing ward for the first time, will be long... This can be seen by\n# measuring using time\nimport time\nstart = time.time()\n\n# Agglomerative Clustering: ward\n\n# We build parameters of our own for this object. Parameters related to\n# masking, caching and defining number of clusters and specific parcellations\n# method.\nward = Parcellations(method='ward', n_parcels=1000,\n standardize=False, smoothing_fwhm=2.,\n memory='nilearn_cache', memory_level=1,\n verbose=1)\n# Call fit on functional dataset: single subject (less samples).\nward.fit(dataset.func)\nprint(\"Ward agglomeration 1000 clusters: %.2fs\" % (time.time() - start))\n\n# We compute now ward clustering with 2000 clusters and compare\n# time with 1000 clusters. To see the benefits of caching for second time.\n\n# We initialize class again with n_parcels=2000 this time.\nstart = time.time()\nward = Parcellations(method='ward', n_parcels=2000,\n standardize=False, smoothing_fwhm=2.,\n memory='nilearn_cache', memory_level=1,\n verbose=1)\nward.fit(dataset.func)\nprint(\"Ward agglomeration 2000 clusters: %.2fs\" % (time.time() - start))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualize: Brain parcellations (Ward)\n\nFirst, we display the parcellations of the brain image stored in attribute\n`labels_img_`\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ward_labels_img = ward.labels_img_\n\n# Now, ward_labels_img are Nifti1Image object, it can be saved to file\n# with the following code:\nward_labels_img.to_filename('ward_parcellation.nii.gz')\n\nfrom nilearn import plotting\nfrom nilearn.image import mean_img, index_img\n\nfirst_plot = plotting.plot_roi(ward_labels_img, title=\"Ward parcellation\",\n display_mode='xz')\n\n# Grab cut coordinates from this plot to use as a common for all plots\ncut_coords = first_plot.cut_coords" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compressed representation of Ward clustering\n\nSecond, we illustrate the effect that the clustering has on the signal.\nWe show the original data, and the approximation provided by the\nclustering by averaging the signal on each parcel.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Grab number of voxels from attribute mask image (mask_img_).\nimport numpy as np\nfrom nilearn.image import get_data\noriginal_voxels = np.sum(get_data(ward.mask_img_))\n\n# Compute mean over time on the functional image to use the mean\n# image for compressed representation comparisons\nmean_func_img = mean_img(dataset.func[0])\n\n# Compute common vmin and vmax\nvmin = np.min(get_data(mean_func_img))\nvmax = np.max(get_data(mean_func_img))\n\nplotting.plot_epi(mean_func_img, cut_coords=cut_coords,\n title='Original (%i voxels)' % original_voxels,\n vmax=vmax, vmin=vmin, display_mode='xz')\n\n# A reduced dataset can be created by taking the parcel-level average:\n# Note that Parcellation objects with any method have the opportunity to\n# use a `transform` call that modifies input features. Here it reduces their\n# dimension. Note that we `fit` before calling a `transform` so that average\n# signals can be created on the brain parcellations with fit call.\nfmri_reduced = ward.transform(dataset.func)\n\n# Display the corresponding data compressed using the parcellation using\n# parcels=2000.\nfmri_compressed = ward.inverse_transform(fmri_reduced)\n\nplotting.plot_epi(index_img(fmri_compressed, 0),\n cut_coords=cut_coords,\n title='Ward compressed representation (2000 parcels)',\n vmin=vmin, vmax=vmax, display_mode='xz')\n# As you can see below, this approximation is almost good, although there\n# are only 2000 parcels, instead of the original 60000 voxels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Brain parcellations with KMeans Clustering\n\nWe use the same approach as with building parcellations using Ward\nclustering. But, in the range of a small number of clusters,\nit is most likely that we want to use standardization. Indeed with\nstandardization and smoothing, the clusters will form as regions.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# class/functions can be used here as they are already imported above.\n\n# This object uses method='kmeans' for KMeans clustering with 10mm smoothing\n# and standardization ON\nstart = time.time()\nkmeans = Parcellations(method='kmeans', n_parcels=50,\n standardize=True, smoothing_fwhm=10.,\n memory='nilearn_cache', memory_level=1,\n verbose=1)\n# Call fit on functional dataset: single subject (less samples)\nkmeans.fit(dataset.func)\nprint(\"KMeans clusters: %.2fs\" % (time.time() - start))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualize: Brain parcellations (KMeans)\n\nGrab parcellations of brain image stored in attribute `labels_img_`\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "kmeans_labels_img = kmeans.labels_img_\n\nplotting.plot_roi(kmeans_labels_img, mean_func_img,\n title=\"KMeans parcellation\",\n display_mode='xz')\n\n# kmeans_labels_img is a Nifti1Image object, it can be saved to file with\n# the following code:\nkmeans_labels_img.to_filename('kmeans_parcellation.nii.gz')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Brain parcellations with ReNA Clustering\n\nOne interesting algorithmic property of ReNA (see References) is that\nit is very fast for a large number of parcels (notably faster than Ward).\nAs before, the parcellation is done with a Parcellations object.\nThe spatial constraints are implemented inside the Parcellations object.\n\n### References\n\nMore about ReNA clustering algorithm in the original paper\n\n * A. Hoyos-Idrobo, G. Varoquaux, J. Kahn and B. Thirion, \"Recursive\n Nearest Agglomeration (ReNA): Fast Clustering for Approximation of\n Structured Signals,\" in IEEE Transactions on Pattern Analysis and\n Machine Intelligence, vol. 41, no. 3, pp. 669-681, 1 March 2019.\n https://hal.archives-ouvertes.fr/hal-01366651/\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "start = time.time()\nrena = Parcellations(method='rena', n_parcels=5000, standardize=False,\n smoothing_fwhm=2., scaling=True)\n\nrena.fit_transform(dataset.func)\nprint(\"ReNA 5000 clusters: %.2fs\" % (time.time() - start))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualize: Brain parcellations (ReNA)\n\nFirst, we display the parcellations of the brain image stored in attribute\n`labels_img_`\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rena_labels_img = rena.labels_img_\n\n# Now, rena_labels_img are Nifti1Image object, it can be saved to file\n# with the following code:\nrena_labels_img.to_filename('rena_parcellation.nii.gz')\n\nplotting.plot_roi(ward_labels_img, title=\"ReNA parcellation\",\n display_mode='xz', cut_coords=cut_coords)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compressed representation of ReNA clustering\n\nWe illustrate the effect that the clustering has on the signal.\nWe show the original data, and the approximation provided by\nthe clustering by averaging the signal on each parcel.\n\nWe can then compare the results with the compressed representation\nobtained with Ward.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Display the original data\nplotting.plot_epi(mean_func_img, cut_coords=cut_coords,\n title='Original (%i voxels)' % original_voxels,\n vmax=vmax, vmin=vmin, display_mode='xz')\n\n# A reduced data can be created by taking the parcel-level average:\n# Note that, as many scikit-learn objects, the ReNA object exposes\n# a transform method that modifies input features. Here it reduces their\n# dimension.\n# However, the data are in one single large 4D image, we need to use\n# index_img to do the split easily:\nfmri_reduced_rena = rena.transform(dataset.func)\n\n# Display the corresponding data compression using the parcellation\ncompressed_img_rena = rena.inverse_transform(fmri_reduced_rena)\n\nplotting.plot_epi(index_img(compressed_img_rena, 0), cut_coords=cut_coords,\n title='ReNA compressed representation (5000 parcels)',\n vmin=vmin, vmax=vmax, display_mode='xz')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even if the compressed signal is relatively close\nto the original signal, we can notice that Ward Clustering\ngives a slightly more accurate compressed representation.\nHowever, as said in the previous section, the computation time is\nreduced which could still make ReNA more relevant than Ward in\nsome cases.\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.8.5" } }, "nbformat": 4, "nbformat_minor": 0 }