Nilearn：多个被试连接组的组稀疏逆协方差

总代码

import numpy as np
from nilearn import plotting
n_subjects = 4  # subjects to consider for group-sparse covariance (max: 40)
def plot_matrices(cov, prec, title, labels):
    """Plot covariance and precision matrices, for a given processing. """
    prec = prec.copy()  # avoid side effects
    # Put zeros on the diagonal, for graph clarity.
    size = prec.shape[0]
    prec[list(range(size)), list(range(size))] = 0
    span = max(abs(prec.min()), abs(prec.max()))
    # Display covariance matrix
    plotting.plot_matrix(cov, cmap=plotting.cm.bwr,
                         vmin=-1, vmax=1, title="%s / covariance" % title,
                         labels=labels)
    # Display precision matrix
    plotting.plot_matrix(prec, cmap=plotting.cm.bwr,
                         vmin=-span, vmax=span, title="%s / precision" % title,
                         labels=labels)

第一步 获取数据集

from nilearn import datasets
n_subjects = 4  
msdl_atlas_dataset = datasets.fetch_atlas_msdl()
rest_dataset = datasets.fetch_development_fmri(n_subjects=n_subjects)
# print basic information on the dataset
print('First subject functional nifti image (4D) is at: %s' %
      rest_dataset.func[0])  # 4D data

n_subjects就是多少个被试
msdl_atlas_dataset就是地图的标签，格式是bunch格式。

rest_dataset功能像图像，另外还附带了混淆矩阵。

提取区域信号

from nilearn import input_data
# A "memory" to avoid recomputation
from joblib import Memory
mem = Memory('nilearn_cache')
masker = input_data.NiftiMapsMasker(
    msdl_atlas_dataset.maps, resampling_target="maps", detrend=True,
    high_variance_confounds=True, low_pass=None, high_pass=0.01,
    t_r=2, standardize=True, memory='nilearn_cache', memory_level=1,
    verbose=2)
masker.fit()
subject_time_series = []
func_filenames = rest_dataset.func
confound_filenames = rest_dataset.confounds
for func_filename, confound_filename in zip(func_filenames,
                                            confound_filenames):
    print("Processing file %s" % func_filename)
    region_ts = masker.transform(func_filename,
                                 confounds=confound_filename)
    subject_time_series.append(region_ts)

mem设置内存运行文件

nilearn.input_data.NiftiLabelsMasker和我们平常使用的mask（掩码）差不多，这里的设置的意思是：使用msdl_atlas_dataset.maps地图集。
resampling_target{“data”, “mask”, “maps”, None}, optional.
表示对输入的图像重采样，这里是变成和maps地图集一样的尺寸形状。它有多个选择，可以是data或者mask，默认是none。
detrend去除线性趋势。
high_variance_confounds，在指定的maps图像上计算高方差混淆并回归。
standardize标准化。
subject_time_series 空列表

func_filenames就是功能像，4D的

confound矩阵文件
然后通过for循环，提取了每一个功能像的ROI的时间序列，形成被试的时间序列。


这就是subject_time_series文件，可以看出原始文件是python原始数据格式list，然后它的形状是（4，168，39），4个被试，每个ROI有168个时间序列，一共有39个ROI。
默默吐槽一点，python自带的list还是不如numpy好用，好像nilearn以后会自全部变成使用numpy格式，这样更方便点。

计算群稀疏精度矩阵

from nilearn.connectome import GroupSparseCovarianceCV
gsc = GroupSparseCovarianceCV(verbose=2)
gsc.fit(subject_time_series)
try:
    from sklearn.covariance import GraphicalLassoCV
except ImportError:
    # for Scitkit-Learn < v0.20.0
    from sklearn.covariance import GraphLassoCV as GraphicalLassoCV
gl = GraphicalLassoCV(verbose=2)
gl.fit(np.concatenate(subject_time_series))

gsc就是 具有交叉验证的参数选择稀疏的逆协方差矩阵估算器。
group 组 sparse 稀疏的 covariance 协方差，cv就是交叉验证。如果不指定cv，默认是3. verbose就是打印出处理内容的方式，没啥用，可以设置可以不设置。

np.concatenate就是连接、拼接函数，就是把3维数据拼接成二维数据，所以被试1有168个时间序列，被试2有168个时间序列…4个被试就是672个时间序列，行是时间序列，列是39个地图的ROI。
下面这个是计算协方差矩阵估算器fit的用法：
estimator.fit([time_series_1, time_series_2, …])
estimator.covariances_[0] 协方差
estimator.precisions_[0] 逆协方差

注意，这里有两个估算器，lasso那个其实是参考，真正要用的就是gsc估算器，这个估算器是可以输入三维数据的，不需要去降维。

显示结果

def plot_matrices(cov, prec, title, labels):
    """Plot covariance and precision matrices, for a given processing. """
    prec = prec.copy()  # avoid side effects
    # Put zeros on the diagonal, for graph clarity.
    size = prec.shape[0]
    prec[list(range(size)), list(range(size))] = 0
    span = max(abs(prec.min()), abs(prec.max()))
    # Display covariance matrix
    plotting.plot_matrix(cov, cmap=plotting.cm.bwr,
                         vmin=-1, vmax=1, title="%s / covariance" % title,
                         labels=labels)
    # Display precision matrix
    plotting.plot_matrix(prec, cmap=plotting.cm.bwr,
                         vmin=-span, vmax=span, title="%s / precision" % title,
                         labels=labels)

atlas_img = msdl_atlas_dataset.maps
atlas_region_coords = plotting.find_probabilistic_atlas_cut_coords(atlas_img)
labels = msdl_atlas_dataset.labels
plotting.plot_connectome(gl.covariance_,
                         atlas_region_coords, edge_threshold='90%',
                         title="Covariance",
                         display_mode="lzr")
plotting.plot_connectome(-gl.precision_, atlas_region_coords,
                         edge_threshold='90%',
                         title="Sparse inverse covariance (GraphicalLasso)",
                         display_mode="lzr",
                         edge_vmax=.5, edge_vmin=-.5)
plot_matrices(gl.covariance_, gl.precision_, "GraphicalLasso", labels)

title = "GroupSparseCovariance"
plotting.plot_connectome(-gsc.precisions_[..., 0],
                         atlas_region_coords, edge_threshold='90%',
                         title=title,
                         display_mode="lzr",
                         edge_vmax=.5, edge_vmin=-.5)
plot_matrices(gsc.covariances_[..., 0],
              gsc.precisions_[..., 0], title, labels)
plotting.show()

首先是定义了一个plot_matrices的函数，这个函数的功能就是画相关性矩阵的图。先看下面的代码。

atlas_img就是一个4D的地图集，前3维是XYZ坐标，第四维是ROI的label。


altas_region_coords就是ROI的坐标，需要跟地图空间保持一致。

plotting.plot_connectome这个函数就是用来画功能连接图的，如下

nilearn.plotting.plot_connectome(adjacency_matrix, node_coords, node_color=‘auto’, node_size=50, edge_cmap=<matplotlib.colors.LinearSegmentedColormap object>, edge_vmin=None, edge_vmax=None, edge_threshold=None, output_file=None, display_mode=‘ortho’, figure=None, axes=None, title=None, annotate=True, black_bg=False, alpha=0.7, edge_kwargs=None, node_kwargs=None, colorbar=False)
第一个参数adjacency_matrix是表示图的链接强度。矩阵可以是对称的，这将导致无向图，也可以是不对称的，这将导致有向图。在这里就是用的-gsc.precisions_[…, 0]。
补充说明：sklearn.covariance.GraphicalLassoCV是稀疏的逆协方差估计器，来自sklearn。需要注意的是，一般使用它的负值，就是-gl.precision_，这样就得到了偏相关，也就是直接相关。这比full correlation更加的精确。而对于这里使用的组被试，就需要用GroupSparseCovarianceCV来计算逆协方差，计算出来的值和Lasso是一样的。所以它也需要负值，-gsc.precisions_[…, 0]。covariance_是普通的协方差矩阵值。
当然，最好的还是建议使用tangent space correlation，就是切线空间的相关。

用法就是：
1.先定义估算器，measure = ConnectivityMeasure(kind=‘tangent’)
2.开始估算，connectivities = measure.fit([time_series_1, time_series_2, …]) 。注意这里是time_1, time_2，所以需要用concatenate函数降维。
（注意，GraphicalLassoCV—gl是专门用于单个主题的，对于多个主题，就需要先数据降维，变成time_1, time_2…的格式。而GroupSparseCovarianceCV—gsc是专门用于组计算的，所以不需要降维，直接3维就可以计算。这里的切线空间相关，得先降维）
3.组内ROI的连接强度，使用measure.mean_ 平均值。

edge_threshold=‘90%’：如果是数字，则仅显示值大于 edge_threshold 的边。如果是字符串，则必须以百分号结尾，例如“25.3%”，并且只会显示 abs(value) 高于给定百分位数的边。

plot_matrices就是画下面这种相关性图：

plot_matrices(cov, prec, title, labels)：
输入协方差矩阵，精确矩阵，标题，脑区标签
先来看看gsc.precisions_和gl.covariance_

和我猜的一样，就是39个ROI两两做相关，但是因为gsc是3维去做相关的，所以出来的值也是三维的，所以就是（39，39，4）.
而同理，gl是二维的，它把3维数据展平成2维，所以出来的矩阵是二维的。但是有一个问题是，它出来的（39，39），对于同一个label，它是进行平均了嘛？这个函数manual上没有说，一般是这么做。

1-复制精确矩阵-prec；2-获取精确矩阵第一个维度的长度-size；3-先生成0-38的值合并成列表，随后XX; 4-计算精确矩阵的最小值和最大值，把它们取绝对值，然后取两者的最大值-span； 5-