教程 | fMRI连接性分析

传统的多变量fMRI技术侧重于局部区域（ROI或探照灯）活动模式中的信息。有时，相关信息可能会跨脑区网络表示，因此无法通过ROI分析或探照灯识别。功能连接测量有助于在全局范围内检查相距较远的脑区的信息，重点关注网络交互而不是空间定位。在执行连接性分析时，将跨区域比较BOLD时间序列（通常使用相关性指标），并且关系的大小决定了它们的功能连接强度。通过包含或排除刺激/任务变量，我们可以研究不同认知状态对连接性的调节。本文接下来将描述如何运行基于种子的连接分析，以及使用图谱来进行分割和定义种子点(基于Python)。

import warnings
import sys 
if not sys.warnoptions:
    warnings.simplefilter("ignore")

import numpy as np
import os 
import nibabel as nib
from nilearn.input_data import NiftiMasker, NiftiLabelsMasker
from nilearn import plotting
from nilearn import datasets
from nilearn.connectome import ConnectivityMeasure
from scipy import stats
from scipy.ndimage.measurements import center_of_mass
import matplotlib.pyplot as plt
import seaborn as sns 
import pandas as pd
import brainiak.utils.fmrisim as sim
from brainiak.fcma.util import compute_correlation
from nilearn import input_data
import time
from utils import shift_timing

%autosave 5
%matplotlib inline
sns.set(style = 'white', context='poster', rc={"lines.linewidth": 2.5})
sns.set(palette="colorblind")

## 每5秒自动保存一次。

使用的数据集来自Hutchinson等人(2016)的研究，要求被试(在每个block中)注意左侧或右侧的一个场景。下面是描述数据集的README文件。该数据集经过运动校正和线性去趋势预处理。

README：注意和连接

在神经层面，注意由顶叶和额叶皮层控制，它们调节感觉系统的加工，增强关注信息并抑制未关注的信息。因此，为了研究注意对感知加工的影响，不仅需要检查局部大脑区域，还需要检查这些区域如何相互作用。此时，检查局部大脑区域活动模式的传统MVPA技术是不够的，因而需要更全面的分析。在这种情况下，两个或多个大脑区域的共变反应变得至关重要，功能连接测量对于研究这些现象是很关键的。

加载数据

加载被试1的预处理数据。

from utils import latatt_dir
assert os.path.exists(latatt_dir)

dir_time = os.path.join(latatt_dir, 'onsets/fsl')
dir_motion = os.path.join(latatt_dir, 'processed_data/motionnuisance/')
dir_motion_background = os.path.join(latatt_dir, 'processed_data/background/motionnuisance/')

sub = 'sub01'
num_runs = 1
TR = 1.5
scan_duration = 540

#将数据平移一定量
shift_size = 2

a、创建刺激标签和时移

从刺激时序文件中，我们需要为BOLD数据中的每个TR创建标签。因为有该数据集的FSL起始文件，所以可以在BrainIAK中使用fmrisim为每个TR创建标签。这涉及调用在utils中定义的generate_stimfunctionthen和shift_timingutils。

## 时移与否？在传统的MVPA分析中，对BOLD信号进行时移以考虑血流动力学的时滞情况。在功能连接分析中，有时不进行时移以确保被比较的体素在刺激呈现之前、期间和之后是相似的。如果不对数据进行时移，可以设置“shift_size = 0”。注意：这里设置了‘shift_size = 2’。

# Use the utilities from the simulator to create an event time course based on an FSL onset file

right_stimfunction = sim.generate_stimfunction(
    onsets='', 
    event_durations='', 
    total_time=scan_duration,
    temporal_resolution=1/TR, 
    timing_file=(dir_time + '/%s/right.txt' % (sub))
)

left_stimfunction  = sim.generate_stimfunction(
    onsets='', 
    event_durations='', 
    total_time=scan_duration,
    temporal_resolution=1/TR, 
    timing_file=(dir_time + '/%s/left.txt' % (sub))
)

# 移动时间进程以考虑血流动力学的时滞情况
right_stim_lag = shift_timing(right_stimfunction, shift_size)
left_stim_lag = shift_timing(left_stimfunction, shift_size)

b、检查头文件

# 读取nifti对象
nii = nib.load(dir_motion + '%s.nii.gz' % (sub))
hdr=nii.get_header()
print(hdr)
print('Voxel size in %s, Time in %s' %hdr.get_xyzt_units())

c、绘制刺激标签

该实验有两个条件，每个条件都有一个单独的计时文件。

# 绘制刺激的时间进程。
plt.figure(figsize=(14, 4))
plt.plot(right_stim_lag)
plt.plot(left_stim_lag)
plt.yticks([0,1])
plt.xlabel('Timepoints')

plt.legend(('Attend Right', 'Attend Left'), loc='upper right')
plt.ylim(0, 1.5)
sns.despine()

d、掩膜和提取全脑数据

使用nilearn从数据集中创建掩膜，然后从掩膜中提取数据。此函数也可以获得数据的z分数。接下来，先绘制一些体素的时间过程。

# 初始化一个同样标准化数据的掩膜对象
masker_wb = input_data.NiftiMasker(
    standardize=True,  # Are you going to zscore the data across time?
    t_r=1.5, 
    memory='nilearn_cache',  # 将掩码缓存在此处作为字符串给出的目录中，以便加载和检索
    memory_level=1,
    verbose=0
)

# 提取体素的时间进程
bold_wb = masker_wb.fit_transform(nii)
bold_wb_r = bold_wb[(right_stim_lag==1),:]
bold_wb_l = bold_wb[(left_stim_lag==1),:]

print('shape - whole brain bold time series: ', np.shape(bold_wb))
print('shape - whole brain bold time series, attend left : ', np.shape(bold_wb_l))
print('shape - whole brain bold time series, attend right : ', np.shape(bold_wb_r))

"""
绘制几个体素的时间序列
"""
voxel_ids = [0,10,100]

plt.figure(figsize=(14, 4))
plt.title('Voxel activity for rightward trials, voxel ids = ' + str(voxel_ids));
plt.plot(bold_wb_r[:, voxel_ids]);
plt.ylabel('Evoked activity');
plt.xlabel('Timepoints');
sns.despine()

创建种子

此时，已经加载了两个实验条件的全脑数据。为了检验注意力的影响，我们将创建种子ROI，并将它们的活动与大脑中的其他体素相关联。对于与种子ROI相关的任何体素，可以推断它们在功能上是相互联系的。

a、创建球形ROI

创建一个ROI来定义海马旁区域(PPA)，这是一个场景选择区域，因为该实验数据集中呈现了场景刺激，并使用MNI空间中的坐标来识别个体被试的区域。

Nilearn有一些用于绘制ROI的强大工具。这些功能使我们能够灵活地识别任何形状的ROI，并具有多个参数，如平滑、去趋势、滤波和标准化。但是，使用这些函数时很容易出错，所以要谨慎使用这些参数。接下来将使用nilearn中最基本的Sphere ROI函数。

b、绘制掩膜的Bold信号

计算并绘制掩膜中所有体素的平均Bold信号。

# 初始化掩膜对象
masker_lPPA = input_data.NiftiSpheresMasker(
    coords_lPPA, 
    radius=8, standardize=True, t_r=2.,
    memory='nilearn_cache', memory_level=1, verbose=0
)

# 对Epi数据进行掩膜，并获取ROI的时间序列
bold_lPPA = masker_lPPA.fit_transform(nii)

# 绘制两种注意条件下种子区域的数据
bold_lPPA_r = bold_lPPA[right_stim_lag == 1, :]
bold_lPPA_l = bold_lPPA[left_stim_lag == 1, :]

plt.figure(figsize=(14, 4))
plt.title('Left PPA ROI activity for right and left attention conditions')
plt.plot(bold_lPPA_r);
plt.plot(bold_lPPA_l);
plt.legend(('Attend Right', 'Attend Left'));
plt.ylabel('Evoked activity');
plt.xlabel('Timepoints');
sns.despine()

计算相关矩阵

使用相关矩阵评估功能连接性。这些矩阵中的每个单元格反映了一对体素或区域之间BOLD时间序列的相关性，并且可以为每种条件甚至每个试次单独计算矩阵。这里，将通过一种基于循环的方法来计算大脑中每个体素与PPA种子区域的相关性。

# 将种子与每个大脑体素相关联。循环并提取每个体素的数据。
start_time = time.time()
num_voxels = bold_wb_r.shape[1]
all_corr = np.zeros((num_voxels, 1))
for v in range(num_voxels): 
    all_corr[v, 0] = np.corrcoef(bold_lPPA_r.flatten(), bold_wb_r[:, v])[0, 1]

end_time = time.time()
print('Analysis duration for %d voxels: %0.2fs' % (num_voxels, (end_time - start_time)))

因为Pearson相关性的有界性质违反了某些统计假设，所以常用的方法是将相关性转换为Fisher Z分数。这通常不会对结果产生很大影响。

def seed_correlation(wbBold, seedBold):
    """计算种子体素与其他体素之间的相关性 
    Parameters
    ----------
    wbBold [2d array], n_stimuli x n_voxels 
    seedBold, 2d array, n_stimuli x 1

    Return
    ----------    
    seed_corr [2d array], n_stimuli x 1
    seed_corr_fishZ [2d array], n_stimuli x 1
    """
    num_voxels = wbBold.shape[1]
    seed_corr = np.zeros((num_voxels, 1))
    for v in range(num_voxels):    
        seed_corr[v, 0] = np.corrcoef(seedBold.flatten(), wbBold[:, v])[0, 1]
    # 将相关值转换为Fisher z分数    
    seed_corr_fishZ = np.arctanh(seed_corr)
    return seed_corr, seed_corr_fishZ

# 使用函数并输出结果的范围
corr_lPPA_r, corr_fz_lPPA_r = seed_correlation(bold_wb_r, bold_lPPA_r)

print("Seed-based correlation Fisher-z transformed: min = %.3f; max = %.3f" % (
    corr_fz_lPPA_r.min(), corr_fz_lPPA_r.max()
))

# A histogram is always a useful first way of looking at your data.
plt.hist(corr_fz_lPPA_r)
plt.ylabel('Frequency');
plt.xlabel('Fisher-z score');
sns.despine()

## 基于种子相关的Fisher-z转换：min=-0.561；max=1.875

# 将相关数组转换回一个可以保存的Nifti图像对象
img_corr_lPPA_r= masker_wb.inverse_transform(corr_fz_lPPA_r.T)
# img_corr_lPPA_r.to_filename('seed_rtstim.nii.gz')

a、绘制种子相关性

绘制与其他所有体素的种子相关性。为了更好地可视化，这里将设置一个阈值，只显示高于阈值的体素。通常，阈值是根据统计显著性来选择的。这里随意选择了以下阈值。

# 将种子与每个体素的相关性进行可视化
threshold = .8

# Nilearn提供了可将结果绘制成地图的工具
r_map_ar = plotting.plot_stat_map(
    img_corr_lPPA_r, 
    threshold=threshold,
    cut_coords=coords_lPPA[0],
)
# 添加种子
r_map_ar.add_markers(
    marker_coords=coords_lPPA, 
    marker_color='g',
    marker_size=50
)

以下是绘制相同信息的另一种方法。

# 绘制一个玻璃大脑
plotting.plot_glass_brain(
    img_corr_lPPA_r, 
    threshold=threshold,
    colorbar=True, 
    plot_abs=False,
    display_mode='lyrz', 
);

从图谱创建种子

除了创建我们自己的种子ROI，还可以使用可用的图谱（哈佛-牛津皮质图谱）来提取ROI。Nilearn提供了一种简单的方法来实现这一点。

atlas = datasets.fetch_atlas_harvard_oxford('cort-maxprob-thr25-2mm')
atlas_filename = atlas.maps

# 图谱就保存在这里。
print("Atlas path: " + atlas_filename + "\n\n")

# 绘制ROIs
plotting.plot_roi(atlas_filename);
print('Harvard-Oxford cortical atlas')

# 打印标签
# 创建一个图谱数据的Pandas数据框架，以便于检查
atlas_pd = pd.DataFrame(atlas)
print(atlas_pd['labels'])

# 创建一个掩膜对象，可以用它来选择ROIs
masker_ho = NiftiLabelsMasker(labels_img=atlas_filename)
print(masker_ho.get_params())

# 将图谱应用到Nifti对象，这样就可以从单个parcels/ROIs中提取数据
bold_ho = masker_ho.fit_transform(nii)
print('shape: parcellated bold time courses: ', np.shape(bold_ho))

前面计算了单个种子区域随时间变化的平均活动。但是，Nilearn有一些工具可以轻松计算提供给掩膜对象的所有ROI的活动时间进程。

# 仅获取用于右侧注意的数据
bold_ho_r = bold_ho[(right_stim_lag==1),:]

# 我们的数据结构是什么样的？
print("Parcellated data shape (time points x num ROIs)")
print("All time points  ", bold_ho.shape)
print("Rightward attention trials: ", bold_ho_r.shape)

# 提取一个对应于后海马旁皮层的ROI
# 意标签#35是海马旁回，后部。 
roi_id = 34
bold_ho_pPHG_r = np.array(bold_ho_r[:, roi_id])
bold_ho_pPHG_r = bold_ho_pPHG_r.reshape(bold_ho_pPHG_r.shape[0],-1)
print("Posterior PPC (region 35) rightward attention trials: ", bold_ho_pPHG_r.shape)

plt.figure(figsize=(14,4))
plt.plot(bold_ho_pPHG_r)
plt.ylabel('Evoked activity');
plt.xlabel('Timepoints');
sns.despine()

# 将整个大脑时间进程与所提取的种子相关联
corr_pPHG_r, corr_fz_pPHG_r = seed_correlation(
    bold_wb_r, bold_ho_pPHG_r
) 

# 打印相关性范围
print("PHG correlation Fisher-z transformed: min = %.3f; max = %.3f" % (
    corr_fz_pPHG_r.min(), corr_fz_pPHG_r.max())
)

# 绘制直方图
plt.hist(corr_fz_pPHG_r)
plt.ylabel('Frequency');
plt.xlabel('Fisher-z score');

## PHG相关Fisher-z转换：min=-0.656；最大值=1.088

# 映射到整个大脑图像
img_corr_pPHG_r = masker_wb.inverse_transform(
    corr_fz_pPHG_r.T
)

threshold = .8 

# 使用分割法找到该ROI的切割坐标
roi_coords = plotting.find_parcellation_cut_coords(atlas_filename,label_hemisphere='left')

# 提取这个ROI的坐标
roi_coord = roi_coords[roi_id,:]

# 在一个标准的大脑中绘制相关图 
# 为了比较，还可以绘制之前创建的球体的位置
h2 = plotting.plot_stat_map(
    img_corr_pPHG_r, 
    threshold=threshold,
    cut_coords=roi_coord,
)

# 绘制玻璃脑
plotting.plot_glass_brain(
    img_corr_pPHG_r, 
    threshold=threshold,
    colorbar=True, 
    display_mode='lyrz', 
    plot_abs=False
)

a、计算parcels之间的连接

可以计算跨多个大脑区域的相关性。当我们想研究不同大脑区域的注意力情况时就很有用。这里将使用不同的方式绘制连接矩阵。

# 设置连接对象
correlation_measure = ConnectivityMeasure(kind='correlation')

# 计算每个parcel与所有其他parcel的相关性
corr_mat_ho_r = correlation_measure.fit_transform([bold_ho_r])[0]

# 为了更好地进行可视化，移除对角线(guaranteed to be 1.0)
np.fill_diagonal(corr_mat_ho_r, np.nan)

# 绘制相关矩阵
# 这里使用的是哈佛-牛津皮质图谱的标签 
# 从背景(0)开始，因此跳过了第一个标签
fig = plt.figure(figsize=(11,10))
plt.imshow(corr_mat_ho_r, interpolation='None', cmap='RdYlBu_r')
plt.yticks(range(len(atlas.labels)), atlas.labels[1:]);
plt.xticks(range(len(atlas.labels)), atlas.labels[1:], rotation=90);
plt.title('Parcellation correlation matrix')
plt.colorbar();

# 或者，我们也可以使用Nilearn的绘图函数进行绘制
plotting.plot_matrix(
    corr_mat_ho_r, 
    cmap='RdYlBu_r', 
    figure=(11, 10), 
    labels=atlas.labels[1:], 
)

绘制连接组

Nilearn中有一些绘制连接体的工具，如plotting.plot_connectome取逐个节点的相关矩阵和逐个节点的坐标矩阵，然后创建一个连接组。阈值可用于仅显示强连接。

# 加载图谱
atlas_nii = nib.load(atlas_filename)
atlas_data = atlas_nii.get_data()
labels = np.unique(atlas_data)

# 遍历所有ROIs
coords = []
for label_id in labels:

    # 跳过背景
    if label_id == 0:
        continue

    # 提取掩膜内的ROI    
    roi_mask = (atlas_data == label_id)

    # 创建为nifti对象，这样就可以由cut coords算法读取
    nii = nib.Nifti1Image(roi_mask.astype('int16'), atlas_nii.affine)

    # 找到连接组的质心
    coords.append(plotting.find_xyz_cut_coords(nii))

# 绘制连接组
plotting.plot_connectome(correlation_matrix_mcg, coords, edge_threshold='95%')

往期推荐：

小伙伴们点个“在看”，加(星标)关注茗创科技，将第一时间收到精彩内容推送哦～

更多精彩课程推荐

MC_Brain

茗创科技工作室

觉得有帮助，欢迎转发收藏或者点个在看哦~

听说点在看的人SCI接收率都提升了18%呢！