Giving credit

Previous topic

8.4.6. Extracting signals from a brain parcellation

Next topic

8.4.8. Extract signals on spheres from an atlas and plot a connectome

8.4.7. Dictionary Learning and ICA for doing group analysis of resting-state fMRIΒΆ

This example applies dictionary learning and ICA to resting-state data, visualizing resulting components using atlas plotting tools.

Dictionary learning is a sparsity based decomposition method for extracting spatial maps. It extracts maps that are naturally sparse and usually cleaner than ICA

  • Gael Varoquaux et al. Multi-subject dictionary learning to segment an atlas of brain spontaneous activity Information Processing in Medical Imaging, 2011, pp. 562-573, Lecture Notes in Computer Science

Available on https://hal.inria.fr/inria-00588898/en/

Load ADHD rest dataset

from nilearn import datasets

adhd_dataset = datasets.fetch_adhd(n_subjects=30)
func_filenames = adhd_dataset.func  # list of 4D nifti files for each subject

# print basic information on the dataset
print('First functional nifti image (4D) is at: %s' %
      adhd_dataset.func[0])  # 4D data

Out:

First functional nifti image (4D) is at: /home/parietal/gvaroqua/nilearn_data/adhd/data/0010042/0010042_rest_tshift_RPI_voreg_mni.nii.gz

Create two decomposition estimators

from nilearn.decomposition import DictLearning, CanICA

n_components = 40

Dictionary learning

dict_learning = DictLearning(n_components=n_components,
                             memory="nilearn_cache", memory_level=2,
                             verbose=1,
                             random_state=0,
                             n_epochs=1)

CanICA

canica = CanICA(n_components=n_components,
                memory="nilearn_cache", memory_level=2,
                threshold=3.,
                n_init=1,
                verbose=1)

Fit both estimators

estimators = [dict_learning, canica]
names = {dict_learning: 'DictionaryLearning', canica: 'CanICA'}
components_imgs = []

for estimator in estimators:
    print('[Example] Learning maps using %s model' % names[estimator])
    estimator.fit(func_filenames)
    print('[Example] Saving results')
    # Decomposition estimator embeds their own masker
    masker = estimator.masker_
    # Drop output maps to a Nifti   file
    components_img = masker.inverse_transform(estimator.components_)
    components_img.to_filename('%s_resting_state.nii.gz' %
                               names[estimator])
    components_imgs.append(components_img)

Out:

[Example] Learning maps using DictionaryLearning model
[MultiNiftiMasker.fit] Loading data from [/home/parietal/gvaroqua/nilearn_data/adhd/data/0010042/0010042_rest_tshift_RPI_voreg_mni.nii.gz, /home/parietal/gvaroqua/nilearn_data/adhd/data/0010064/0010064_rest_tshift_RPI_voreg_mni.nii.gz, /home
[MultiNiftiMasker.fit] Computing mask
[MultiNiftiMasker.transform] Resampling mask
[DictLearning] Loading data
[DictLearning] Learning initial components
[DictLearning] Computing initial loadings
[DictLearning] Learning dictionary
[Example] Saving results
[Example] Learning maps using CanICA model
[MultiNiftiMasker.fit] Loading data from [/home/parietal/gvaroqua/nilearn_data/adhd/data/0010042/0010042_rest_tshift_RPI_voreg_mni.nii.gz, /home/parietal/gvaroqua/nilearn_data/adhd/data/0010064/0010064_rest_tshift_RPI_voreg_mni.nii.gz, /home
[MultiNiftiMasker.fit] Computing mask
[MultiNiftiMasker.transform] Resampling mask
________________________________________________________________________________
[Memory] Calling sklearn.utils.extmath.randomized_svd...
randomized_svd(array([[-0.002998, ..., -0.001381],
       ...,
       [-0.002863, ...,  0.007679]]), n_iter=3, random_state=None, transpose=True, n_components=40)
__________________________________________________randomized_svd - 14.9s, 0.2min
________________________________________________________________________________
[Memory] Calling sklearn.decomposition.fastica_.fastica...
fastica(array([[ 0.002945, ...,  0.002152],
       ...,
       [ 0.003551, ..., -0.007815]]), fun='cube', random_state=1754728016, whiten=True)
_________________________________________________________fastica - 17.0s, 0.3min
[Example] Saving results
________________________________________________________________________________
[Memory] Calling nilearn.masking.unmask...
unmask(array([[ 0.006561, ...,  0.      ],
       ...,
       [-0.      , ..., -0.      ]]),
<nibabel.nifti1.Nifti1Image object at 0x2aed8d4df750>)
___________________________________________________________unmask - 0.9s, 0.0min

Visualize the results

from nilearn.plotting import (plot_prob_atlas, find_xyz_cut_coords, show,
                              plot_stat_map)
from nilearn.image import index_img

# Selecting specific maps to display: maps were manually chosen to be similar
indices = {dict_learning: 1, canica: 31}
# We select relevant cut coordinates for displaying
cut_component = index_img(components_imgs[0], indices[dict_learning])
cut_coords = find_xyz_cut_coords(cut_component)
for estimator, components in zip(estimators, components_imgs):
    # 4D plotting
    plot_prob_atlas(components, view_type="filled_contours",
                    title="%s" % names[estimator],
                    cut_coords=cut_coords, colorbar=False)
    # 3D plotting
    plot_stat_map(index_img(components, indices[estimator]),
                  title="%s" % names[estimator],
                  cut_coords=cut_coords, colorbar=False)
show()
  • ../../_images/sphx_glr_plot_compare_resting_state_decomposition_001.png
  • ../../_images/sphx_glr_plot_compare_resting_state_decomposition_002.png
  • ../../_images/sphx_glr_plot_compare_resting_state_decomposition_003.png
  • ../../_images/sphx_glr_plot_compare_resting_state_decomposition_004.png

Total running time of the script: ( 4 minutes 12.694 seconds)

Generated by Sphinx-Gallery