ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

基于Pearson相关和LOFC的静息态fMRI功能连接分析全流程

基于Pearson相关和LOFC的静息态fMRI功能连接分析全流程 预处理跑完头动曲线、白质信号、脑脊液信号都看过了配准结果也逐层检查过没问题那接下来干什么如果你做的是静息态fMRI最经典、最不容易出错、审稿人最买账的下一步就是功能连接分析。而功能连接里最基础、最稳定、也最适合作为论文第一张图的就是低阶功能连接LOFC计算方式就是Pearson相关。这篇文章我打算把整个流程串一遍从DPABI预处理完的数据开始到用Nilearn提取时间序列、算相关矩阵、做Fisher Z变换再到最后能拿去做统计分析的特征矩阵。中途会讲清楚每一步为什么要这么做以及我实际跑数据时踩过的坑。适合刚入门脑影像、或者已经跑完预处理但对后续特征提取还不太有把握的同学参考。1. 整体思路拆解为什么偏偏是Pearson相关和LOFC1.1 功能连接到底在算什么功能连接Functional Connectivity这个概念听起来玄其实本质很简单大脑不同区域的血氧信号BOLD信号在时间维度上是不是同步变化。如果两个脑区在扫描过程中信号涨落高度同步就认为它们之间存在功能上的连接关系。这里要特别强调它不是解剖上的白质纤维连接纯粹是统计意义上的时间同步性。那怎么衡量这种同步性最直接的办法就是算相关系数。你把每个脑区当成一个时间序列两两配对算Pearson相关得到一个N×N的相关矩阵N是脑区数量。这个矩阵就是低阶功能连接LOFC的核心产物。之所以叫“低阶”是因为还有高阶功能连接——比如把每个脑区当成一个节点计算它和其他所有脑区连接模式的相似性再对这个相似性矩阵做相关这叫做功能连接强度相关FCS或者功能连接模式相似性FCPS。低阶就是直接算原始信号的相关不额外做嵌套计算。为什么低阶连接这么常用两个原因。第一它稳定。Pearson相关是一个非常成熟的统计量计算简单可解释性强不同研究之间容易横向对比。第二它是高阶分析的基础。你后面要做图论分析、做分类特征、做动态功能连接起点都是这个低阶连接矩阵。1.2 用ROI做分析为什么不是体素级功能连接有两种粒度体素级Voxel-wise和ROI级Region of Interest。体素级就是把全脑每个体素的时间序列都提出来选一个种子点算种子点和所有体素的相关最后得到一张相关图。ROI级则是先定义一个脑区模板把模板里每个脑区内部所有体素的时间序列取平均得到一个代表该脑区整体活动的平均时间序列然后脑区两两之间算相关。我建议新手直接做ROI级。原因很实际体素级计算结果是一个统计图需要做多重比较校正、做团块阈值流程繁琐而且结果解读依赖你对脑解剖的熟悉程度。ROI级则是得到一个人数×人×脑区×脑区的矩阵套一个分类器或者做一个组间比较都特别方便对于批量处理、机器学习特征提取来说比体素级友好得多。当然ROI级也有代价——空间分辨率低。你定义一个脑区假设这个区域内部的功能是均一的但实际情况不一定如此。不过这是目前主流做法SCN、Dosenbach、Power这些常用图谱都是这么用的做出来结果也相当稳健不用担心这个。1.3 工具选型DPABI预处理nilearn算连接这里的选型逻辑是DPABI和nilearn各有自己最擅长的事情。DPABI是一个基于SPM的图形界面工具包写论文引用率高、教程丰富、处理流程规范。它最强大的是可视化——预处理每一步的质控图、头动报告、配准效果全都做得很完善特别适合需要严格质控的场景也适合新手找到标准化的处理模板。它内部跑的还是SPM的算法意味着生成的文件格式、命名规则和SPM完全兼容后面接任何SPM生态的工具都没问题。Nilearn则是Python生态里做神经影像的利器玩的就是灵活性和可编程性。它的NiftiLabelsMasker可以一句话完成“读取图谱提取每个ROI的时间序列可选的平滑/去噪”整套操作配合Pandas和Scikit-learn后面做机器学习特征输入非常顺滑。而且Nilearn的代码可复现性远胜图形界面操作跑批处理的时候尤其明显。所以我个人习惯是DPABI管预处理出干净数据Nilearn管后处理做特征提取、统计分析和可视化。两个工具各管一段发挥各自长处。2. 数据准备与预处理关键点这一步没做好后面全白搭2.1 从原始数据到DPABI能识别的格式假设你已经有了DICOM原始数据。第一步是格式转换DPABI推荐用DICOM to NIfTI工具实际上调用的是SPM的dicom2nifti接口界面里直接选“DICOM to NIfTI”指定输入和输出文件夹就能跑。转换之后注意检查生成的nii文件头信息尤其是TR值必须正确后面时间层校正要用。还有一种情况是拿到手的数据已经是NIfTI格式比如从OpenfMRI、ABCD、UK Biobank这类公开数据库下载的数据那就省了转换环节。不过即使已经是NIfTI也建议先跑一遍DPABI的检查功能Check Data把每个被试的所有文件列出来确认没有缺时序、缺文件的情况。我遇到过下载数据时某个被试的run-2只下了一半的情况DPABI的检查功能一下就标红了。2.2 DPABI预处理流程十个步骤按顺序跑DPABI的预处理界面打开后默认有十个步骤按顺序是删除前几帧DICOM to NIfTI之后、时间层校正Slice Timing、头动校正Realign、T1配准到MNI空间Coregister Normalize、分割灰白质脑脊液Segmentation、去除协变量Nuisance Covariate Regression、滤波Filter、去线性漂移Detrend、共配准到标准空间Normalize to MNI、平滑Smooth。我强烈建议不要跳过任何一步更不要打乱顺序。有几个地方特别注意第一删除前几帧。因为扫描刚开始时磁场还没稳定被试还没适应扫描噪音前几帧信号通常不可靠。默认删5帧是通常做法但如果数据质量一般也可以删10帧。DPABI界面里有个选项是“Remove first N time points”填上数字就好。第二时间层校正。这一步要求你知道扫描的采集顺序——是从下往上还是隔层采集以及参考层是哪一层。这些信息在扫描序列参数里能看到或者问扫描技师。填错会导致时间校正无效甚至产生虚假的信号延迟。第三头动校正。DPABI会输出每个被试的头动参数平移和旋转共6个以及一张头动曲线图。一般要求平移不超过3mm、旋转不超过3度严格一点的期刊要求不超过2mm和2度。超标的被试数据要么剔除要么至少要在后续分析里把头动参数作为协变量回归掉。第四配准。SPM的配准分两步先T1像对齐到功能像的平均像Coregister再把这个对齐后的T1配准到MNI模板Normalize最后用这个变换参数把功能像也转到MNI空间。DPABI界面里这步对应的是“T1 New Segment DARTEL”或者“T1 Unified Segmentation”两种方式。我自己用DARTEL比较多配准精度比旧版Unified Segmentation好一些尤其是皮层区域。2.3 质控必须自己看一遍预处理跑完不是直接进入下一个环节而是先检查质量。DPABI最常见的质控方式是在结果目录下生成一系列质控图包括每个被试的头动曲线、配准前后对比图用checkerboard形式显示T1和功能像的重叠情况、标准化后功能像的叠加图。这里我想额外强调一下一定要亲自看不能只靠软件报错。软件没报错不代表结果就好。比如配准有些被试的结构像本身有伪影或者和模板差异过大哪怕软件显示“completed”实际效果可能就是歪的。我习惯是把每个被试的配准结果拼成大图一张一张翻过去看发现问题立刻回预处理重新调参。这一步花的时间远比后面发现问题再返工少得多。另一个容易忽略的质控点是检查预处理后的图像覆盖范围。有些被试扫描时头部位置偏了图像边缘被截断预处理后MNI空间里可能有一侧脑区数据缺失这时候该脑区的时间序列就不用算数了。判断方法就是看看MNI空间后的全脑覆盖率如果有一个ROI覆盖严重不足后续分析里需要单独标注或剔除。2.4 是否需要做全局信号回归全局信号回归Global Signal Regression是个争议话题。支持者认为它能去除全脑性的生理噪声比如心跳、呼吸引起的伪影让局部连接更准确反对者认为它引入了负相关偏倚导致组间比较出现假性差异。我的建议是如果你的研究目的是一般的静息态功能连接探索先不要做全局信号回归只回归白质、脑脊液信号和头动参数就够了。如果审稿人要求或者你确实需要就在敏感性分析里补充一组做全局信号回归的结果做对照。DPABI界面里有一个“Nuisance Covariate Regression”的模块你可以选择回归白质、脑脊液、头动参数也可以勾选全局信号灵活处理。3. Pearson相关计算与LOFC特征提取实操3.1 ROI模板选择Dosenbach还是Power还是其他计算LOFC前需要先确定用哪个脑区模板。常用的模板有几种AALAutomated Anatomical Labeling最经典90个脑区基于解剖边界优点是每个脑区有明确的解剖学意义缺点是脑区划分比较粗可能把功能不同的区域混在一起。Power图谱264个节点基于功能连接定义了覆盖全脑的网络节点在功能社区里使用极广。Dosenbach图谱160个ROI也是基于功能连接定义涵盖了六个主要的认知控制网络区域。我个人做静息态连接分析最常用的是Dosenbach 160因为这个模板和很多已发表文献的接轨度高而且里面每个ROI都自带网络标签比如“default mode network”、“fronto-parietal network”后面对比分析组间差异落在哪个网络时非常方便。但如果你觉得160个ROI太多、算出来的图太密Power 264其实也差不多AAL 90则相对稀疏。选模板的核心原则只有一个和你参考的文献保持一致这样结果才有可比性。还有一点模板选好之后不要中途换全流程统一用同一个。3.2 Nilearn提取时间序列预处理完成之后现在假设我们有一个经过标准化、平滑的功能像文件例如func_preprocessed.nii还有一个图谱文件Dosenbach_160.nii.gz。接下来用Nilearn提取时间序列核心API是NiftiLabelsMaskerfrom nilearn.input_data import NiftiLabelsMasker import pandas as pd import nibabel as nib # 初始化掩模器 masker NiftiLabelsMasker( labels_imgDosenbach_160.nii.gz, standardizeFalse, detrendTrue, low_pass0.1, high_pass0.01, t_r2.0, memorynilearn_cache, memory_level1, verbose1 ) # 读取预处理后的功能像 fmri_img nib.load(func_preprocessed.nii) # 提取时间序列形状为 (n_timepoints, n_rois) time_series masker.fit_transform(fmri_img) # 保存为 Dataframe df pd.DataFrame(time_series, columns[fROI_{i} for i in range(time_series.shape[1])]) df.to_csv(subject001_timeseries.csv, indexFalse)这里面有几个参数要解释一下detrendTrue对每个ROI的时间序列做线性去趋势去除扫描过程中可能的信号漂移。虽然DPABI预处理时已经做过一次去漂移了但Nilearn里再做一次也不会有什么坏处而且如果在Nilearn里直接处理未去趋势的数据这一步就是必须的。low_pass和high_pass带通滤波。静息态功能连接通常关注低频波动0.01~0.1 Hz因为这个频段被认为是自发神经活动的窗口。如果你在DPABI里已经做了同样的滤波这里就可以不设或者设为None避免重复滤波。t_r2.0这是重复时间要和采集参数一致。带通滤波需要这个参数来计算滤波器的系数。standardizeFalse如果要输出原始信号设为False。但如果你后面要算相关矩阵建议设置成True让每个脑区的时间序列做z-score标准化这样算相关之前信号已经归一化了。一个小技巧是memory参数。它指定了一个缓存目录Nilearn会把掩模计算、信号提取的中间结果缓存下来。当你要对几百个被试批量处理时这个缓存能省下大量重复计算的时间。3.3 用Pearson相关构建功能连接矩阵时间序列提取出来后Pearson相关就是一行代码的事。Nilearn里有两个选择一个是直接用numpy.corrcoef一个是用手动计算的方式目的是把自相关对角线置零方便后面做图论分析时忽略自环。import numpy as np # 手动计算 Pearson 相关系数矩阵 n_rois time_series.shape[1] conn_matrix np.corrcoef(time_series, rowvarFalse) # 对角线置零因为节点和自身不构成连接 np.fill_diagonal(conn_matrix, 0) # 保存 np.savetxt(subject001_connectivity_matrix.txt, conn_matrix, delimiter,)这里np.corrcoef输出的实际上就是标准的Pearson相关系数取值范围在-1到1之间。为什么要把对角线置零从数学上讲任何变量和自身的相关系数都是1但功能连接网络分析中我们关心的是节点之间的连接自己和自己不算连接。而且很多图论指标比如聚类系数、效率计算时假设邻接矩阵对角线为0不做这个处理会报错或者结果异常。3.4 Fisher Z变换从相关系数到适合统计的正态分布量Pearson相关系数有个不太方便的性质它的抽样分布不是正态的当实际相关系数接近±1时分布会被压缩到边界附近。这意味着如果你直接对相关系数做t检验、ANOVA这些参数检验容易违反正态性假设。解决办法就是Fisher Z变换公式是z 0.5 * ln((1 r) / (1 - r))这个变换把[-1, 1]区间映射到(-∞, ∞)而且变换后的z近似服从正态分布。统计功效更高组间比较也更有把握。Nilearn的NiftiLabelsMasker本身不做这个变换需要自己算z_matrix np.arctanh(conn_matrix)其实就是np.arctanh。一句话的事但这是整个流程里最容易被新手忽略的步骤。如果直接拿相关矩阵去做统计分析结果不一定错但统计严谨性肯定要打折扣。很多审稿人看到你没有做Fisher Z变换会直接提问。所以建议处理连接矩阵后统一做一个z变换再保存。3.5 可视化功能连接矩阵和连接图计算完矩阵可视化是下一步。最常用的可视化方式有热力图和连接组图。import matplotlib.pyplot as plt from nilearn import plotting # 热力图 plt.figure(figsize(8, 6)) plt.imshow(z_matrix, cmapRdBu_r, vmin-1, vmax1) plt.colorbar(labelFisher Z-transformed correlation) plt.xlabel(ROI index) plt.ylabel(ROI index) plt.title(LOFC Matrix - Subject 001) plt.tight_layout() plt.savefig(subject001_LOFC_heatmap.png, dpi300) plt.show() # 连接组图 plotting.plot_connectome( z_matrix, node_coordscoords, # 从图谱文件获取ROI坐标 edge_threshold50%, node_size20, edge_cmapRdBu_r, ) plotting.show()连接组图是从大脑的三维视角展示连接模式的图非常直观。做的时候注意edge_threshold要根据你的连接密度来调整默认阈值保留了一半边。如果边太多图会成一团乱麻太少又看不出整体架构我一般会试着调到30%~60%之间选一个好看又能看出特征的阈值。3.6 批量处理一个人跑完几十个被试上面讲的是单个被试的流程。实际做研究肯定是一批被试肯定要写成循环。我的建议是把单个被试的处理封装成一个函数然后用列表遍历import os from nilearn.input_data import NiftiLabelsMasker masker NiftiLabelsMasker( labels_imgDosenbach_160.nii.gz, standardizeTrue, detrendTrue, low_pass0.1, high_pass0.01, t_r2.0, memorynilearn_cache, memory_level1, verbose0 ) subjects [sub-001, sub-002, sub-003] # 实际用 os.listdir 读取 for sub in subjects: func_file fpath/to/{sub}/func_preprocessed.nii if not os.path.exists(func_file): print(f{sub} missing functional image) continue ts masker.fit_transform(func_file) conn np.corrcoef(ts, rowvarFalse) np.fill_diagonal(conn, 0) z np.arctanh(conn) np.save(foutput/{sub}_LOFC_z.npy, z)强烈建议把处理后的数据存成NumPy的.npy格式而不是文本.txt因为后面要把所有人拼成一个大的数组喂给模型或者做统计分析.npy加载快、省空间。这里的missing functional image就是一个重要的容错逻辑——不能假设每个被试的文件都一定存在。4. 常见问题与排查技巧实录4.1 头动伪影对功能连接的污染头动是静息态fMRI最大的敌人。哪怕被试只动了很小的幅度脑边界的体素信号也会被污染而且这种污染会通过配准和标准化扩散到更大的范围。最典型的“假连接”出现在眼眶和脑干附近——因为这些地方信号本来就弱容易被运动伪影盖过。所以处理功能连接时必须检查头动参数是否作为协变量回归掉。如果你的DPABI预处理没做头动回归或者不放心可以在Nilearn里手动加# 假设运动参数文件是6列的文本 motion_params np.loadtxt(subject001_motion.txt) from nilearn.image import clean_img clean_img clean_img( fmri_img, confoundsmotion_params, detrendTrue, low_pass0.1, high_pass0.01, t_r2.0 ) ts masker.fit_transform(clean_img)这里的clean_img就是回归掉头动参数之后的图像。注意confounds可以传多个协变量比如白质信号、脑脊液信号Nilearn会自动进行多元回归。4.2 Pearson相关的伪相关现象噪声和信号混淆Pearson相关对异常值极度敏感。如果一个脑区的时间序列里有一个巨大幅度的尖峰比如因为头动或者设备信号跳动这个尖峰可以单独拉高它和所有其他脑区的相关系数。这种现象在文献里叫“伪相关”spurious correlation。排查方法在计算相关矩阵之前先画几个关键ROI的时间序列图看看有没有明显的异常尖峰。如果发现异常可以用winsorize截尾的方式处理把超过某个百分位数的值压到那个百分位数上。Nilearn没有内置winsorize但可以用scipy.stats.mstats.winsorize实现。不过更推荐的做法是回到预处理阶段重新做头动校正因为尖峰往往意味着某个体素对齐不准而不是纯粹的统计问题。4.3 多重比较校正LOFC矩阵56万次检验的问题160个脑区算出的相关矩阵有160×159/2 12720个连接对。如果两组被试比较每个连接对都要做一个假设检验12720次检验如果不做校正假阳性的数量会非常可观。比如p0.05时纯随机数据里也会有约600个“显著”连接。所以Nilearn里提供了一些多重比较校正方法from scipy.stats import ttest_ind from nilearn.connectome import ConnectivityMeasure import numpy as np # group1 和 group2 分别是两个组的连接矩阵列表 # 假设 shape 都是 (n_subjects, n_rois, n_rois) t_stats [] for i in range(n_rois): for j in range(i1, n_rois): # 提取两组对应连接的值 vals1 group1[:, i, j] vals2 group2[:, i, j] t, p ttest_ind(vals1, vals2) t_stats.append((i, j, t, p)) # 用 FDR 校正 from statsmodels.stats.multitest import multipletests pvals [item[3] for item in t_stats] reject, pvals_corrected, _, _ multipletests(pvals, methodfdr_bh)这里用的是statsmodels的FDR校正方法Benjamini-Hochberg。更省事的方式是直接用Nilearn的ConnectivityMeasure加上permutation_testing功能做置换检验不过那个计算量大一些适合最终分析阶段用。4.4 时间序列长度不够短扫描导致的相关矩阵不稳定静息态扫描时间长短直接影响相关矩阵的可靠性。理论上为了稳定估计低频段0.01~0.1 Hz的连接至少需要包含几个完整周期也就是最少5分钟左右的静息态数据。如果你用的是某公开数据库里只有3分钟或更短的数据相关矩阵的方差会变大复现性变差。一种补救办法是使用“分半信度”来检验你的连接矩阵是否稳定把时间序列分成前一半和后一半分别计算相关矩阵然后看两个矩阵之间的相关有多高。如果这个分半相关的值低于某个阈值经验上0.3~0.5说明这个被试的数据质量不足以支撑稳定的连接估计可以考虑剔除或者标记为低质量。4.5 协变量选择哪些该回归掉哪些不该预处理时常用协变量包括头动6参数、白质信号、脑脊液信号、全局信号。前三个是几乎必做的头动参数去除运动干扰白质和脑脊液信号去除生理噪声。全局信号做不做看文献习惯前面已经讨论过。还有一个容易被忽略的协变量是“运动导数”把头动参数做差分捕捉瞬时运动幅度有些严格的预处理会把6参数和它们的一阶差分共12个变量一起回归能更好地去除个体运动的瞬时效应。Nilearn的confounds参数只需把12列都传进去即可非常方便。4.6 处理不同扫描站点/不同批次数据时Harmonization问题如果你用的是多中心数据比如多医院的合研数据不同站点之间的信号差异必须处理否则组间差异可能是站点效应混入。常用的方法是ComBat调和Python里有neuroCombat库可以直接使用。具体的做法是在所有被试算完连接矩阵之后把矩阵的每个连接值作为特征做ComBat校正把扫描站点作为批处理变量。这个方法能有效去除中心效应同时保留真实的组间差异。5. 踩坑记录与实操心得最后聊几个我在实际运行中踩过的比较有价值的坑以及一些经验性总结希望帮你少走弯路。5.1 DPABI跑完以后没有生成预期文件的排查DPABI偶尔会出现预处理中途报错的情况尤其是在Windows环境下路径中文名或者空格都可能引发问题。我建议把DPABI、SPM和数据路径全都放在不含中文、不含空格的目录下。比如D:\fMRI_project\data这种结构。如果报错信息看起来像Subscript indices must either be real positive integers or logicals十有八九是指定了不存在的文件路径。另一个我遇到的情况是DPABI里选择了“Normalize by DARTEL”但T1分割结果没生成导致配准步骤失败。解决方法是回到SPM batch里检查T1分割的输出文件是否真的存在文件名为c1T1_*.nii、c2T1_*.nii这些。如果不存在说明分割步骤没跑完需要看SPM窗口有没有红色报错信息。5.2 NiftiLabelsMasker报错IndexError的排查Nilearn在处理图谱文件时偶尔会报IndexError: index N is out of bounds这个报错信息通常是因为图谱文件里的标签编号不连续。比如Dosenbach图谱里的ROI编号是1、2、3、5、6中间有个4被跳过可能是某个ROI被废弃了。Nilearn需要知道确切的标签列表如果在NiftiLabelsMasker里没有传入labels参数它会自动读取图谱中出现的所有非零整数。这种情况下报错的话手动传入标签列表解决import nibabel as nib import numpy as np labels_img nib.load(Dosenbach_160.nii.gz) labels_data labels_img.get_fdata() unique_labels np.unique(labels_data) unique_labels unique_labels[unique_labels 0] # 去掉0 print(unique_labels) # 看看标签范围是否正确 masker NiftiLabelsMasker( labels_imgDosenbach_160.nii.gz, labelsunique_labels, # 强制指定标签 )5.3 关于ROI内信号平均Nilearn默认对每个ROI内的体素信号取平均。但如果你做的是大图谱比如AAL的90个区每个区的体素数差别很大平均值的噪声特性也不同而如果用的是解剖图谱那么ROI内信号异质性问题更明显。这时候可以考虑用PCA主成分分析或ICA的方式提取每个ROI内的主要信号成分而不是简单平均。Nilearn的NiftiLabelsMasker没有直接做PCA的选项需要手动写from sklearn.decomposition import PCA from nilearn.input_data import NiftiMasker # 先对全脑数据做mask masker NiftiMasker(mask_imgDosenbach_160_roi1_mask.nii.gz) roi1_data masker.fit_transform(fmri_img) # shape (n_timepoints, n_voxels) pca PCA(n_components1) roi1_pca_signal pca.fit_transform(roi1_data).ravel()不过这种方法毕竟复杂而且不同的方法对后续结果是否有显著影响文献结论不一作为起步做均值就够了不用强求PCA提取。5.4 滤波和去趋势的顺序DPABI里的常规顺序是先回归协变量再滤波再去线性漂移。如果你用Nilearn的clean_img它的默认顺序是先去趋势再回归协变量再带通滤波。虽然看起来顺序不同但实际结果差异很小。如果你担心可以在clean_img里设置detrendTrue, standardizeFalse它会自动做去趋势后再回归、再滤波。重点是不要在预处理里做两次滤波比如DPABI已经过了0.01-0.1HzNilearn里又设置一次低频滤波——这会让信号的有效自由度进一步下降结果反而变差。5.5 标准化设置对结果的影响Nilearn的standardize参数控制是否对时间序列做z-score标准化。我在3.2节里建议如果是要算相关矩阵就设为True如果后面要做频谱分析或者保持原始幅度就设为False。因为Pearson相关本身对尺度变化不敏感z-score标准化在计算相关之前做了相当于把每个脑区的时间序列缩放到方差为1理论上对相关系数没有影响但是对数值稳定性有帮助特别是在数据存在较大尺度差异的时候。所以做一个z-score标准化没有坏处做了也不影响相关系数。5.6 关于低阶连接和高阶连接的关系这一节是对标题里“低阶功能连接LOFC”的补充。低阶连接是两两脑区的直接相关性而高阶连接是连接模式的连接。举个例子低阶关注的是“脑区A和脑区B的时间序列是否同步”高阶关注的是“脑区A的连接分布与所有其他脑区的相关向量和脑区B的连接分布是否相似”。高阶方法在2013年左右的文献里很热门因为它能捕捉到低阶连接看不到的“模式相似性”。但你要明确一点做低阶连接并不意味着你的研究层次低。恰恰相反很多疾病研究比如阿尔茨海默、精分、抑郁症在LOFC上都有非常稳定的组间差异发现。低阶连接是基础高阶连接是延伸。如果你在LOFC上没发现明显差异再到高阶上去找这种由低到高的分析路径更容易讲清楚故事。我自己做过的研究里通常流程是先算LOFC矩阵看整体网络拓扑小世界属性、聚类系数、特征路径长度再做网络内/网络间连接强度比较比如默认网络内部、默认网络与额顶网络之间最后如果有必要再做高阶FCS。这样由浅入深逻辑通顺审稿人也容易跟。5.7 关于显著性阈值的选择如果你最后要做组间比较就要考虑连接水平上的多重比较校正。常见方案是FDR错误发现率频率学派里最常用计算量小。还有一种做法是网络水平上的置换检验Permutation testing通过随机打乱组标签看真实的组间差异在随机分布里处在什么位置这个方法在样本量小的时候比较稳。Nilearn提供了permutation_testing接口但计算量大在大样本上跑一次可能需要几分钟到几十分钟要合理安排计算资源。我个人倾向的做法是先用FDR做初筛选出p0.05经过FDR校正的连接然后再对这些连接做一个聚类看看它们是否集中在某些网络上。如果显著的连接都落在同一个网络里那么即使没有单独的网络水平校正这个结果的可信度也很高。结尾写到这里整个基于Pearson相关的LOFC特征提取流程已经完整走了一遍。对我来说这套方法最让人放心的地方在于它的“透明性”——每一步都有明确的数学定义每个结果都能在标准空间里定位回大脑结构每一个环节出问题都能顺着时间序列、相关矩阵、统计结果一层层回溯排查。相比之下很多深度学习的脑影像特征提取方法虽然在分类准确率上可能更高但要解释清楚特征到底在“哪里”、为什么重要难度要大得多。最后再分享两个小建议。第一个是存档习惯跑完预处理和特征提取之后把所有中间产物——包括去趋势失配的文件、滤波后的时间序列、标准化的功能像、连接矩阵、z变换后的矩阵——都分门别类存放在独立文件夹里并且用命名规则标记清楚处理状态。你后来说不定会需要回去看某个特定处理阶段的输出。第二个是计算资源的规划160个ROI的连接矩阵在单机上跑完全没问题但如果你要扩展到全脑体素级功能连接或者上千被试的批量处理强烈建议早点学习用集群任务提交比如SLURM别等数据堆在那里才想起来要学。
返回列表