Nilearn从大脑图谱中提取时间序列及制作功能连接

from nilearn import datasets
dataset = datasets.fetch_atlas_harvard_oxford('cort-maxprob-thr25-2mm')
atlas_filename = dataset.maps
labels = dataset.labels
print('Atlas ROIs are located in nifti image (4D) at: %s' %
      atlas_filename)  # 4D data
# One subject of brain development fmri data
data = datasets.fetch_development_fmri(n_subjects=1)
fmri_filenames = data.func[0]

解释一下：首先是获取了哈佛的地图所有信息（dataset），atlas_filename就是获取哈佛地图，其次获取label，label就是地图的标签。len(lables)=49，哈佛地图一共有49个标签。但是，需要注意的是，第0个标签是background，背景呈现用的，所以第0个标签是不作为ROI去提取时间序列的，所以masker默认是48。

from nilearn.input_data import NiftiLabelsMasker
masker = NiftiLabelsMasker(labels_img=atlas_filename, standardize=True,
                           memory='nilearn_cache', verbose=5)
# Here we go from nifti files to the signal time series in a numpy
# array. Note how we give confounds to be regressed out during signal
# extraction
time_series = masker.fit_transform(fmri_filenames, confounds=data.confounds)

NiftiLabelsMasker 是个重要的函数，该函数可以根据地图标签，制作mask。可惜的是，这个mask无法被预览，是个特殊的格式，NiftiLabelsMasker数据格式。
fmri_filenames是一个4D的时间序列。

time_series导师很普通的，是一个（168，48）的矩阵。一共有有48个标签，每个标签有168个时间点组成的时间序列，所以就是（168，48）。

计算并显示相关矩阵

from nilearn.connectome import ConnectivityMeasure
correlation_measure = ConnectivityMeasure(kind='correlation')
correlation_matrix = correlation_measure.fit_transform([time_series])[0]

# Plot the correlation matrix
import numpy as np
from nilearn import plotting
# Make a large figure
# Mask the main diagonal for visualization:
np.fill_diagonal(correlation_matrix, 0)
# The labels we have start with the background (0), hence we skip the
# first label
# matrices are ordered for block-like representation
plotting.plot_matrix(correlation_matrix, figure=(10, 8), labels=labels[1:],
                     vmax=0.8, vmin=-0.8, reorder=True)

这段的意思是：计算的方法是相关（即dparsf的full correlation），然后生成correlation_matrix，相关性矩阵，矩阵维度是（48，48），最后呈现出的就是矩阵图。

同样的事情没有混淆，强调混淆的重要性

time_series = masker.fit_transform(fmri_filenames)
# Note how we did not specify confounds above. This is bad!

correlation_matrix = correlation_measure.fit_transform([time_series])[0]

# Mask the main diagonal for visualization:
np.fill_diagonal(correlation_matrix, 0)

plotting.plot_matrix(correlation_matrix, figure=(10, 8), labels=labels[1:],
                     vmax=0.8, vmin=-0.8, title='No confounds', reorder=True)

plotting.show()

这个就不展示了，图像基本都是强相关，没啥意义。
这一步的主要意思是，计算时间序列相关性时，需要把白质、脑脊液等信号给回归掉，这个在dparsf里就是nuisance这一步，是自动处理的，就是不用太关注这个，了解就行。

总结：总结一下提取时间序列的过程：
1、导入4D的图像（因为3D的图像没有时间维度，所以必须是4D的volumes）
2、选择带有标签的地图集，map and labels
3、制作labels的masks，用NiftiLabelsMasker函数
4、使用masker.fit_transform功能在4D的volumes上提取出masks的时间序列
5、计算并呈现相关性矩阵
重要补充：对于举一反三来说，假如我想提取masks中特定的几个labels，其实很简单，就是把labels和time_series根据索引值去掉那几项，因为这两个都是numpy格式，所以去掉很简单，这样就能得到期望中的几个labels的时间序列。很完美。

补充代码：

# 导入labels和altas
from nilearn import datasets
dataset = datasets.fetch_atlas_harvard_oxford('cort-maxprob-thr25-2mm')
atlas_filename = dataset.maps
labels = dataset.labels
# 制作ROI的masker
from nilearn.input_data import NiftiLabelsMasker
masker = NiftiLabelsMasker(labels_img=atlas_filename, standardize=True)
# 根据masker提取时间序列，混淆矩阵就是csv_file
# 当然，对于已经去除掉混淆矩阵的nii_files，就不用啦
time_series = masker.fit_transform(frmi_files, confounds=csv_file)

如果是概率图谱提取时间序列，只需要把NiftiLabelsMasker换成NiftiMapsMasker就行了，很简单记住，确定的就是labels，因为规定好了，概率的就是maps，因为只是地图，没有规定好。

from nilearn.input_data import NiftiMapsMasker
masker = NiftiMapsMasker(maps_img=atlas_filename, standardize=True)
time_series = masker.fit_transform(frmi_files, confounds=csv_file)

功能连接

# 导入地图
from nilearn import datasets
yeo = datasets.fetch_atlas_yeo_2011() #yeo 是3D图像，有17个labels
# 导入数据
data = datasets.fetch_development_fmri(n_subjects=10)
# 提取yeo地图块的坐标
from nilearn.input_data import NiftiLabelsMasker
from nilearn.connectome import ConnectivityMeasure
# ConenctivityMeasure from Nilearn uses simple 'correlation' to compute
# connectivity matrices for all subjects in a list
connectome_measure = ConnectivityMeasure(kind='correlation')
# useful for plotting connectivity interactions on glass brain
from nilearn import plotting
# create masker to extract functional data within atlas parcels
masker = NiftiLabelsMasker(labels_img=yeo['thick_17'], standardize=True,
                           memory='nilearn_cache')
# extract time series from all subjects and concatenate them
time_series = []
for func, confounds in zip(data.func, data.confounds):
    time_series.append(masker.fit_transform(func, confounds=confounds))
# calculate correlation matrices across subjects and display
correlation_matrices = connectome_measure.fit_transform(time_series)
# Mean correlation matrix across 10 subjects can be grabbed like this,
# using connectome measure object
mean_correlation_matrix = connectome_measure.mean_
# grab center coordinates for atlas labels
coordinates = plotting.find_parcellation_cut_coords(labels_img=yeo['thick_17'])
# plot connectome with 80% edge strength in the connectivity
plotting.plot_connectome(mean_correlation_matrix, coordinates,
                         edge_threshold="80%",
                         title='Yeo Atlas 17 thick (func)')

导入的yeo是一个bunch数据，有很多信息。导入的data同理，也是bunch数据。

time_series是一个（10，168，17）的三维矩阵。10代表的是10个被试，168代表的是168个时间点，17代表的是yeo[‘thick_17’]里面有17个labels，这是一个皮层的数据。我记得yeo模板有好多，有皮层模板，有7大脑功能网络模板，这里的是皮层模板。如下展示，包含17个ROI。

记住：对于多个被试，就用for循环提取labels的时间序列，并制作成list。

有意思的来了，用correlation_measure计算相关矩阵时，对于list，就是list的每一条进行相关性分析，所以就是（10，17，17），很有意思。

mean_correlation_matrix就是计算10个被试的平均时间序列，记得要在mean后面加一个_，特殊记住就行了，按照这个格式来。

前面都没有啥，这一个非常重要，plotting.find_parcellation_cut_coords函数，只需要输入labels的图像，就可以自动提取labels图像里面的位置坐标，XYZ格式，非常实用！！！对于所有的有labels的图像，都可以使用这个方法。

这个就没啥说的，看代码就行。

import numpy as np
# Define a custom function to compute lag correlation on the time series
def lag_correlation(time_series, lag):
    n_subjects = len(time_series)
    n_samples, n_features = time_series[0].shape
    lag_cor = np.zeros((n_subjects, n_features, n_features))
    for subject, serie in enumerate(time_series):
        for i in range(n_features):
            for j in range(n_features):
                if lag == 0:
                    lag_cor[subject, i, j] = np.corrcoef(serie[:, i],
                                                         serie[:, j])[0, 1]
                else:
                    lag_cor[subject, i, j] = np.corrcoef(serie[lag:, i],
                                                         serie[:-lag, j])[0, 1]
    return np.mean(lag_cor, axis=0)
# Compute lag-0 and lag-1 correlations and plot associated connectomes
for lag in [0, 1]:
    lag_correlation_matrix = lag_correlation(time_series, lag)
    plotting.plot_connectome(lag_correlation_matrix, coordinates,
                             edge_threshold="90%",
                             title='Lag-{} correlation'.format(
                                 lag))

如果是对称的矩阵，就会画一个无向的图，如果是用lag-1作为连接度量，就会画一个有向的图。
尝试学习lag-1这个指标。可以先做lag0，就是普通的（同上面的功能连接），再作lag1，就是滞后1个时间点的。两个图作为补充提交。

对于概率图谱的坐标提取，比较难，仅作了解就行，代码如下：

from nilearn.input_data import NiftiMapsMasker
# create masker to extract functional data within atlas parcels
masker = NiftiMapsMasker(maps_img=difumo.maps, standardize=True,
                         memory='nilearn_cache')
# extract time series from all subjects and concatenate them
time_series = []
for func, confounds in zip(data.func, data.confounds):
    time_series.append(masker.fit_transform(func, confounds=confounds))
# calculate correlation matrices across subjects and display
correlation_matrices = connectome_measure.fit_transform(time_series)
# Mean correlation matrix across 10 subjects can be grabbed like this,
# using connectome measure object
mean_correlation_matrix = connectome_measure.mean_
# grab center coordinates for probabilistic atlas
coordinates = plotting.find_probabilistic_atlas_cut_coords(maps_img=difumo.maps)
# plot connectome with 85% edge strength in the connectivity
plotting.plot_connectome(mean_correlation_matrix, coordinates,
                         edge_threshold="85%",
                         title='DiFuMo with {0} dimensions (probabilistic)'.format(dim))
plotting.show()