传统的多变量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_dirassert 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 = 1TR = 1.5scan_duration = 540
#将数据平移一定量shift_size = 2from utils import latatt_dirassert 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 = 1TR = 1.5scan_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()
"""绘制几个体素的时间序列"""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")
# 绘制ROIsplotting.plot_roi(atlas_filename);print('Harvard-Oxford cortical atlas')
# 打印标签
# 创建一个图谱数据的Pandas数据框架，以便于检查atlas_pd = pd.DataFrame(atlas)print(atlas_pd['labels'])
# 创建一个掩膜对象，可以用它来选择ROIsmasker_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%')
往期推荐：
小伙伴们点个“在看”，加
编辑
(星标)关注，将第一时间收到精彩内容推送哦～
