9.4.5 基于ICA和字典学习的群体fMRI空间图提取

本人能力有限，如有错误，欢迎指出

---

9.4.5 基于ICA和字典学习的群体fMRI空间图提取

原文见：
https://nilearn.github.io/auto_examples/03_connectivity/plot_compare_decomposition.html#sphx-glr-auto-examples-03-connectivity-plot-compare-decomposition-py

存在各种方法从组fmr数据中导出空间图或网络。这些方法提取在一段时间内有相似血氧水平波动的区域。分解方法（Decomposition methods ）允许同时生成许多独立的映射，而不需要提供先验信息（例如：生成种子或先验信息）。
这个例子将使用两种常用的分解方法，ICA（独立成分分析）和字典学习（Dictionary Learning），来测量儿童和年轻人看电影时的fMRI数据。生成的图将使用atlas绘图工具进行可视化。
CanICA是一种ICA方法，用于fMRI数据的组水平分析。与其他策略相比，该算法不仅具有良好的控制群体模型，而且还具有阈值算法，通过信号的显式模型来控制特异性和敏感性。参考文献为：
G. Varoquaux et al. “A group model for stable multi-subject ICA on fMRI datasets”, NeuroImage Vol 51 (2010), p. 288-299 preprint

9.4.5.1. 加载大脑发育功能磁共振数据集

from nilearn import datasets

rest_dataset = datasets.fetch_development_fmri(n_subjects=30)
func_filenames = rest_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' %
      rest_dataset.func[0])  # 4D data

9.4.5.2. 在数据上使用CanICA

我们使用“template”作为计算掩模（mask）的策略，因为这会导致更快和更可重复的结果。但是，图像需要在MNI模板空间中。

from nilearn.decomposition import CanICA

canica = CanICA(n_components=20,
                memory="nilearn_cache", memory_level=2,
                verbose=10,
                mask_strategy='template',
                random_state=0)
canica.fit(func_filenames)

#得到大脑空间中的独立成分。直接使用属性‘components_img_’
# Retrieve the independent components in brain space. Directly accesible through attribute `components_img_`.
canica_components_img = canica.components_img_
#得到的类型是Nifti Image object，可以使用以下方法存储：
# components_img is a Nifti Image object, and can be saved to a file with the following line:
canica_components_img.to_filename('canica_resting_state.nii.gz')

为了可视化，我们在一个图形上绘制所有组件的轮廓

from nilearn.plotting import plot_prob_atlas

# Plot all ICA components together
plot_prob_atlas(canica_components_img, title='All ICA components')

最后，我们分别绘制每个ICA组件

from nilearn.image import iter_img
from nilearn.plotting import plot_stat_map, show

for i, cur_img in enumerate(iter_img(canica_components_img)):
    plot_stat_map(cur_img, display_mode="z", title="IC %d" % i,
                  cut_coords=1, colorbar=False)

9.4.5.3. 比较CanICA和字典学习（dictionary learning）

字典学习是一种基于稀疏性的空间图分解方法。它提取的图自然稀疏，通常比ICA更干净。在这里，我们将比较用CanICA构建的网络和用字典学习构建的网络。
Arthur Mensch et al. Compressed online dictionary learning for fast resting-state fMRI decomposition, ISBI 2016, Lecture Notes in Computer Science

创建字典学习估计器（estimator）

from nilearn.decomposition import DictLearning

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

print('[Example] Fitting dicitonary learning model')
dict_learning.fit(func_filenames)
print('[Example] Saving results')
#抓取提取的组件，并将其还原到Nifti图像。
# Grab extracted components umasked back to Nifti image.
# Note: For older versions, less than 0.4.1. components_img_
# is not implemented. See Note section above for details.
dictlearning_components_img = dict_learning.components_img_
dictlearning_components_img.to_filename('dictionary_learning_resting_state.nii.gz')

将结果可视化
首先，将所有字典学习的组件（DictLearning components）绘制在一起

plot_prob_atlas(dictlearning_components_img,
                title='All DictLearning components')

每个组件的一个图

for i, cur_img in enumerate(iter_img(dictlearning_components_img)):
    plot_stat_map(cur_img, display_mode="z", title="Comp %d" % i,
                  cut_coords=1, colorbar=False)

估计每个组成部分的解释方差（explained variance），使用matplotlib绘制
dict_learning可用于计算每个组件的得分

scores = dict_learning.score(func_filenames, per_component=True)

# Plot the scores
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.ticker import FormatStrFormatter

plt.figure(figsize=(4, 4))
positions = np.arange(len(scores))
plt.barh(positions, scores)
plt.ylabel('Component #', size=12)
plt.xlabel('Explained variance', size=12)
plt.yticks(np.arange(20))
plt.gca().xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
plt.tight_layout()

show()

Note

To see how to extract subject-level timeseries’ from regions created using Dictionary Learning, see example Regions extraction using Dictionary Learning and functional connectomes.

为了从使用字典学习创建的区域提取时间序列，看这个例子 Regions extraction using Dictionary Learning and functional connectomes.