简介这份资源是面向神经科学研究者的结构协变网络构建工具集专注于处理MRI、DTI等脑影像数据覆盖从图像预处理、脑区划分、特征提取到协方差计算、网络构建与拓扑分析的全流程。包内共239个文件以MATLAB脚本.m为主另含.mat数据文件、.nii/.img/.hdr等脑影像模板、.bmp示意图及少量辅助文档和表格压缩包约19.19MB结构清晰便于直接部署到科研工作流中。目前已有1109人学习。借助该工具集使用者能快速复现结构协变网络的标准构建步骤同时获得现成的配准模板、协方差矩阵计算与网络分析模块还可结合包内示例数据进行结果验证和后续拓展为探索大脑解剖连接模式、疾病生物标记等研究提供实用支持。1. 结构协变网络构建先从共变矩阵而不是单脑区看脑形态数据一组被试的灰质体积数据摆在我面前时我通常不会先去看单个脑区萎缩了多少而是先问另一个问题哪些脑区在分别变化的同时呈现出跨脑区的一致性这种脑区与脑区之间的形态共变就是结构协变网络——它用形态学指标之间的相关关系刻画脑区之间共享发育、共享老化或共享病理过程的程度。当你手里只有T1结构像FreeSurfer或CAT12跑完后摆在面前的是几十列的脑区体积表却不知道从哪一步开始建网络NeuroimageTools-main这类工具就能把“脑区指标 → 相关矩阵 → 邻接矩阵 → 组间比较”这段路串起来解决“指标有了却不会往下走”的尴尬。它适合刚完成数据采集、准备做形态学网络分析的研究者也适合想把单脑区结论升级成网络视角的老手。2. 先把T1变成脑区指标两条提取路线与协变量回归的细节2.1 结构协变网络的核心假设为什么“形态变化相似”可以当连接结构协变网络与功能连接网络有一个根本区别功能连接用同一被试的时间序列算相关而结构协变网络是在“群体层面”上用多个被试的脑区形态指标算相关。具体来说对每一对被试你有一个特征向量对每一个脑区你有一列跨被试的测量值。如果脑区A和脑区B在被试之间同升同降就说它们之间存在结构协变。这种协变的生物学解释常见说法是共享的神经发育轨迹、共享的营养与神经递质输入或者共同暴露于同一病理过程。它并不等于直接的白质纤维连接但多项研究观察到结构协变与功能连接、白质连接之间存在部分重叠因此常被当作“结构上的功能伙伴”来使用。做网络构建前理解一个关键边界结构协变网络是“一组被试”的共变模式被试数是样本维度脑区数是特征维度。你不可能用单个被试的T1像构建出一个脑区间的相关网络至少要有几十个被试才有意义。这个前提决定了后面所有步骤的参数选择——被试量越小后面算相关矩阵就越需要谨慎。2.2 形态学指标提取的两条主力路线FreeSurfer与CAT12怎么选结构协变网络能用的形态学指标至少包括灰质体积、皮层厚度、皮层表面积、皮层曲率、脑回化指数。它们不是随便挑一个都行不同指标反映的生物学过程不完全相同。比如皮层厚度与树突密度、神经毡体积关系更直接灰质体积则还受表面积与皮层褶皱程度影响。我的习惯是一次只构建一种指标的协变网络不要混用。从T1结构像到脑区指标主流方案就两条。第一条是FreeSurfer路线。它做完整的皮层重建输出每个被试的皮层厚度、表面积、局部体积然后用脑图谱把顶点数据映射到脑区得到“被试×脑区”的指标表。优点是可以做皮层厚度和表面积相关的精细分析缺点是处理耗时很长一个被试数小时而且对内存敏感大批量跑的时候容易出现个别被试中途崩掉。第二条是CAT12路线。它基于SPM的VBM框架做分割与空间标准化输出灰质体积图再用AAL这类图谱提取每个脑区的灰质体积。优点是单被试十几分钟搞定开源且批量处理稳定很适合样本量大的项目缺点是它对皮层厚度、表面积的支持不如FreeSurfer精细主要产出是灰质体积。给一个选型参考表对比项FreeSurfer 路线CAT12 路线核心产出皮层厚度、表面积、局部体积灰质体积调制/未调制单被试耗时数小时量级十几分钟量级批量稳定性内存敏感偶发中断相对稳定常用脑图谱映射多种图谱可选AAL、Schaefer 等适合的SCN指标厚度、表面积、体积灰质体积如果你的研究核心是灰质萎缩、灰质体积差异CAT12是省心选择如果关注皮质厚度与表面积那就只能走上半条路线。很多做老年化与神经退行性疾病的研究者直接用CAT12加AAL图谱构建灰质体积协变网络批量项目里这是最常见也最可复现的路径。2.3 提取指标后先做协变量回归顺序放在相关矩阵之前形态学指标提取完成后摆在面前的是一个形状为“被试数×脑区数”的矩阵。这个矩阵不能直接拿去做相关计算因为里面混着几个与脑区形态强相关的协变量其中最典型的是年龄、性别、总颅内体积。如果不去除总颅内体积的影响所有脑区的灰质体积都会随头骨大小同向变化结果就是相关矩阵里几乎全脑强相关网络退化成一片没有区分度的连通块。这个坑非常常见我见过不少人拿着未回归的GMV直接建网出来的网络密度动辄超过80%还以为是生物学发现。协变量回归的常规做法是先把每个脑区的指标对协变量做线性回归取残差作为该脑区的去协变量指标再用残差去算相关矩阵。注意顺序必须放在相关矩阵之前不能先算相关再回归那样协变量的影响已经进入相关结构事后清不干净。下面是一个用于回归协变量的Python片段import numpy as np from numpy.linalg import lstsq def regress_confounds(metrics, age, sex, tiv): metrics: 形状为 (n_subject, n_roi) 的脑区指标矩阵 age: (n_subject,) 年龄向量 sex: (n_subject,) 性别向量已做0/1编码 tiv: (n_subject,) 总颅内体积向量 返回: 残差矩阵形状同 metrics可直接用于构建相关矩阵 n_subj metrics.shape[0] # 设计矩阵截距 年龄 性别 总颅内体积 X np.column_stack([np.ones(n_subj), age, sex, tiv]) resid np.empty_like(metrics, dtypenp.float64) for i in range(metrics.shape[1]): # 对每个脑区独立做最小二乘回归 beta, _, _, _ lstsq(X, metrics[:, i], rcondNone) fitted X beta resid[:, i] metrics[:, i] - fitted return resid这段代码的逻辑是把年龄、性别、总颅内体积放到设计矩阵里对每个脑区单独回归把预测值扣掉留下的残差就是剔除协变量影响后的脑区形态指标。参数上需要注意两点。第一设计矩阵里保留截距项不要去掉截距的作用是吸收每个脑区自身的均值水平去掉后残差会残留均值偏移。第二性别在进入设计矩阵前必须先编码成0/1数值不能放文本列总颅内体积的单位只要全样本一致量纲不影响回归结果。回归完成后再对残差做相关计算。3. 从指标矩阵到邻接矩阵皮尔逊相关、偏相关与稀疏度阈值的完整流程3.1 数据组装与质控异常被试别进相关矩阵把回归后的残差矩阵当作网络的输入之前强烈建议先做一轮快速质控。结构协变网络最怕的就是个别被试的极端值把整个脑区的相关方向带偏。我一般做两层检查。第一层是被试层面看每个脑区跨被试的z分数如果某个被试在大量脑区上出现绝对值大于3的z值优先怀疑数据质量从T1分割结果里查该被试是否存在配准失败。第二层是脑区层面如果某个脑区有超过10%的被试缺失或为NaN直接删掉这个脑区比用插补更稳妥。结构协变网络讲的是脑区之间的关系插补出来的脑区数据容易在相关矩阵里产生伪共变。质控的小工具可以用下面这个片段import numpy as np def qc_check(metrics, z_thresh3.0, roi_na_ratio0.1): metrics: (n_subj, n_roi) 残差矩阵或原始指标矩阵 返回: 可疑被试名单、建议删除的脑区名单 n_subj, n_roi metrics.shape z (metrics - np.nanmean(metrics, axis0)) / (np.nanstd(metrics, axis0) 1e-8) bad_subj_ratio np.nanmean(np.abs(z) z_thresh, axis1) bad_subj np.where(bad_subj_ratio 0.2)[0] # 超过20%脑区出现极端z值 roi_na np.mean(np.isnan(metrics), axis0) bad_roi np.where(roi_na roi_na_ratio)[0] return bad_subj, bad_roi这里逻辑很直接计算每个脑区跨被试的z分数如果一个被试在超过20%的脑区上都是极端值就标记为可疑被试如果一个脑区的缺失比例超过10%就建议删除。注意z分数要用每一列的均值和标准差计算不要用全体脑区的整体均值和标准差因为不同脑区的形态指标数值范围差异很大。3.2 相关性计算皮尔逊相关与偏相关怎么选质控完成后进入核心步骤计算脑区之间的相关矩阵。最朴素的方法是皮尔逊相关对残差矩阵直接求两两相关系数。它简单、稳定、对样本量要求相对低是结构协变网络最常用的基线做法。但皮尔逊相关有一个特性它会把两个脑区都跟其他脑区共享的变异也算进相关里。如果全脑存在一个共同的形态变化因子皮尔逊相关矩阵里所有脑区都会被这个全局因子抬高网络看起来过度连篇。想要剔除这种全局影响可以用偏相关。偏相关的计算路径是先求脑区指标矩阵的协方差矩阵再对协方差矩阵求逆得到精度矩阵最后把精度矩阵标准化成偏相关矩阵。Python里可以用下面这段代码import numpy as np from sklearn.covariance import LedoitWolf def partial_corr(metrics, shrinkTrue): metrics: (n_subj, n_roi) 残差矩阵 shrink: True 时用 LedoitWolf 收缩估计协方差矩阵 适合被试数接近或小于脑区数的场景 返回: 偏相关矩阵形状 (n_roi, n_roi) if shrink: # LedoitWolf 收缩可以稳定协方差矩阵估计防止求逆时数值爆炸 cov LedoitWolf().fit(metrics).covariance_ else: cov np.cov(metrics, rowvarFalse) prec np.linalg.inv(cov) # 精度矩阵 diag np.sqrt(np.diag(prec)) pcorr prec / np.outer(diag, diag) np.fill_diagonal(pcorr, 0.0) return pcorr代码逻辑分三步先估计协方差矩阵再求逆得到精度矩阵再标准化成偏相关。参数上shrinkTrue 是关键选项。当被试数n小于脑区数目时协方差矩阵求逆极不稳定LedoitWolf收缩能通过把协方差矩阵往单位矩阵方向压缩让求逆结果可用。我的经验是只要被试数低于脑区数的两倍就无条件打开shrink别硬算原始协方差求逆。皮尔逊与偏相关的取舍我一般这么判断如果希望保留全局共变信号用皮尔逊如果希望突出脑区之间的特异共变、剔除全脑共同因子用偏相关。两种都做、结果交叉验证是最稳妥的做法。无论用哪种相关算出来的矩阵都是带符号的。结构协变网络通常只保留正相关因为脑区形态的负相关在解释上不明确而且多数研究把SCN定义为正协变网络。如果你保留负相关别忘了在阈值化步骤里对绝对值做阈值。3.3 稀疏化与二值化让不同组别的网络具有可比性算完相关矩阵后网络还不是网络只是一个完全连通的加权矩阵。要得到稀疏的邻接矩阵必须做阈值化。两个常见做法绝对阈值与稀疏度阈值。绝对阈值比如保留|r|0.3或|r|0.4的连边简单直观但它在两组被试间会产生不同的连边数和平均度导致后续图论指标的比较被网络密度差异污染。稀疏度阈值也叫sparsity是在两组网络里分别保留相同比例的连边让组间拓扑指标可比这是更被推荐的做法。稀疏度阈值化的实现代码如下def sparsity_threshold(mat, sparsity, weightedTrue): mat: (n_roi, n_roi) 相关矩阵对角线已被置为0 sparsity: 保留连边比例例如 0.2 表示保留强度最大的前20%连边 weighted: True 保留原始权重False 输出二值矩阵 返回: 邻接矩阵 n mat.shape[0] # 只取上三角避免对称重复与对角线 triu_idx np.triu_indices(n, k1) triu_vals mat[triu_idx] n_edge_total len(triu_vals) # 按连接强度降序确定阈值位置 k int(round(n_edge_total * sparsity)) valid_vals triu_vals[np.argsort(-triu_vals)[:k]] thresh valid_vals[-1] if valid_vals.size else 0.0 adj np.zeros_like(mat, dtypenp.float64) mask (mat thresh) (mat ! 0) if weighted: adj[mask] mat[mask] # 保留相关值作为权重 else: adj[mask] 1.0 # 二值网络 adj np.maximum(adj, adj.T) # 保证对称 return adj这段代码的关键是k的确定先提取上三角所有连边强度降序排序取前sparsity比例的连边阈值作为截断点。参数需要特别留意两点。一是weighted的选择做聚类系数、特征路径长度这类指标时保留权重往往更能体现变异信息做有人提出二值化更适应疾病网络研究的说法取决于分析目标。二是np.maximum(adj, adj.T)这步不能省因为从mask恢复邻接矩阵时上三角和下三角可能不对称图论分析里绝大部分指标都要求输入对称邻接矩阵。4. 组间比较的统计验证置换检验与NBS的正确用法4.1 为什么必须在网络层面做置换检验有了患者组和对照组各自的结构协变网络接下来常见的错误是对每个脑区或每条边做两样本t检验然后挑选显著的边声称网络存在差异。这在逻辑上站不住脚——t检验只考虑了单条边各自的效应量没有考虑整张网络的联合结构。你找到了10条显著的边它们是不是碰巧组成一个高度连通的子网络单边检验回答不了这个问题。正确的做法是把整个构建流程装进置换检验的循环里每一次置换都重新打乱组别标签重新计算协变量回归、相关矩阵、稀疏化、网络统计量然后比较真实标签下的统计量与置乱标签下的零分布。这样得到的p值含了所有前序步骤的不确定性而不是只看最后一步的输出。4.2 度中心性置换检验的实现下面这段代码演示了如何对两组被试的度中心性差异做置换检验它把前面所有步骤串起来import numpy as np from sklearn.covariance import LedoitWolf def perm_test_degree(resid, group, n_perm5000, sparsity0.2): resid: (n_subj, n_roi) 回归协变量后的残差矩阵 group: (n_subj,) 0/1 组别标签 n_perm: 置换次数 sparsity: 网络稀疏度 返回: 经验p值与真实标签下的度差异 def build_and_deg(r): cov LedoitWolf().fit(r).covariance_ prec np.linalg.inv(cov) dg np.sqrt(np.diag(prec)) pcorr prec / np.outer(dg, dg) np.fill_diagonal(pcorr, 0) triu np.triu_indices(pcorr.shape[0], k1) vals pcorr[triu] k int(round(len(vals) * sparsity)) threshold np.partition(vals, -k)[-k] adj (pcorr threshold).astype(float) adj np.maximum(adj, adj.T) return adj.sum(axis1) # 度向量 real_diff build_and_deg(resid[group 1]).mean(axis0) - \ build_and_deg(resid[group 0]).mean(axis0) real_stat np.abs(real_diff).mean() count 0 for _ in range(n_perm): perm_labels np.random.permutation(group) perm_diff build_and_deg(resid[perm_labels 1]).mean(axis0) - \ build_and_deg(resid[perm_labels 0]).mean(axis0) perm_stat np.abs(perm_diff).mean() if perm_stat real_stat: count 1 return count / n_perm, real_stat这段代码里build_and_deg把偏相关、稀疏化、度计算封装成一步置换时直接复用。统计量取的是全脑所有脑区度差异绝对值的均值这样你不用逐脑区做多重比较用一个标量概括网络层面的组间差异。参数上要注意n_perm至少要5000次少于1000次的置换在p值0.05附近误差很大把随机种子固定可以复现实际脚本里建议用np.random.default_rng(seed)显式控制。4.3 连边级多重比较NBS的基本思路如果不想只看整网统计量还想定位到差异富集的子网络主流做法是NBS。它的思路是先对每条边做组间t检验设定一个连边级阈值比如t2.0把所有超过阈值的边保留下来找它们构成的连通分量对连通分量的大小做置换检验看真实标签下最大的连通分量是否显著大于置乱标签下的连通分量。NBS比逐边FDR更贴合网络数据的特性它检测的是“一组边形成的连通结构”而不是孤立地看单边显著性。它天然需要一个完整的网络构建流程和前面偏相关、稀疏化的步骤完全兼容。对结构协变网络来说NBS输出的通常是一个异常子网络比逐脑区结果更容易解释。5. 避坑实录结构协变网络构建中五个高频翻车现场5.1 被试量太小相关矩阵整体不稳现象同一样本随机拆成两半分别构建结构协变网络两半连边的重合度不到四成结果完全没法看。 原因结构协变网络是在被试维度上算相关被试数n直接决定相关估计的置信区间宽度。n30时相关标准误会非常大算出来的相关矩阵很大程度是噪声矩阵脑区越多这个现象越明显。 解决样本量至少做到50到80最好过百。如果样本量确实受限强烈建议用LedoitWolf收缩估计协方差或在构建网络时固定sparsity而不是固定绝对相关阈值。收缩并不能创造信息但能让平均值附近的相关结构更稳定总比硬算原始协方差得到一个伪网络强。5.2 协变量回归放在相关矩阵之后现象灰质体积网络密度畸高几乎全脑连通聚类系数接近1一看就是所有脑区都在同步变化。 原因总颅内体积和年龄的效应没有剔除干净就做了相关全脑共享变异把整个相关矩阵抬高或者虽然做了回归但回归步骤放在了相关计算之后残差清理发生得太晚。 解决严格按“提取指标 → 协变量回归 → 残差标准化 → 相关矩阵 → 稀疏化”的顺序执行。先在每个脑区上对年龄、性别、总颅内体积回归只保留残差再把残差作为后续所有步骤的输入。顺序错了结果就是废的。5.3 阈值一换组间差异方向翻转现象用绝对阈值0.3、0.4、0.5分别构建患者组与对照组网络比较某个脑区的度三个阈值下差异方向居然不一致。 原因绝对阈值让两组网络的连边数和平均度不同度指标同时被网络密度和真实组间差异拉扯阈值移动时谁占上风会切换。 解决改用稀疏度阈值让两组网络保持相同的连边比例更好的是做5%到50%的稀疏度扫描看结果是否在某个连续区间内保持一致再选平台期内的阈值报告。绝对阈值最多只能作为辅助验证。5.4 不同形态学指标混进同一张网络现象把灰质体积和皮层厚度拼进同一个“被试×脑区”矩阵再算相关得到的小世界属性时好时坏网络在组间差异里翻来覆去。 原因灰质体积和皮层厚度虽然都来自T1像但生物学意义与个体差异模式不同混在一个矩阵里相关同时被量纲差异和表型差异主导根本说不清在测什么。 解决一次只构建一种形态学指标的协变网络。想报告体积网络与厚度网络就分别构建、分别比较、分别展示不要放在同一矩阵里。5.5 把单脑区t检验当成网络检验现象逐脑区t检验发现五六个脑区显著补网络分析后却找不到任何显著的连通分量两组的小世界属性也没有明显差异。 原因逐脑区检验是在每个脑区独立做假设检验完全没考虑脑区之间的相关结构网络检验关注的是连边和子网络级别的信号两者的统计效率不同结论不同是正常的。 解决设计分析方案时先把问题定清楚想知道整网拓扑是否组间不同用网络统计量加置换检验想知道差异富集在哪个子网络用NBS。不要拿逐脑区t检验结果来发布网络结论容易被审稿人直接质疑统计逻辑。6. 稀疏度扫描把阈值选择从玄学变成可复现步骤6.1 写一个5%到50%的稀疏度扫描函数阈值选择是结构协变网络里最玄的一步。为了不让自己在某个密度上挑到一个好看的数字我习惯直接写一个扫描函数把5%到50%的稀疏度全部跑一遍输出一张“稀疏度-网络属性”对照表import numpy as np def sparsity_scan(mat, sparsitiesNone): mat: (n_roi, n_roi) 已置零对角线的相关矩阵 返回: 每个稀疏度对应的网络属性汇总 if sparsities is None: sparsities np.arange(0.05, 0.51, 0.05) n mat.shape[0] triu_idx np.triu_indices(n, k1) triu_vals np.abs(mat[triu_idx]) # 用绝对值排名保留强正负相关 results [] for s in sparsities: k int(round(len(triu_vals) * s)) thresh np.partition(triu_vals, -k)[-k] adj (np.abs(mat) thresh).astype(float) adj np.maximum(adj, adj.T) np.fill_diagonal(adj, 0) degree adj.sum(axis1) results.append({ sparsity: round(s * 100, 1), density: round((adj 0).sum() / (n * (n - 1)), 4), mean_degree: float(degree.mean()), degree_std: float(degree.std()), }) return results这段代码与前面的稀疏化函数不同之处在于它用绝对值对连边强度做排名正负相关都会参与阈值判断适合你还没想好是否只保留正相关的阶段排查。输出的三个指标分别对应网络密度、平均度和度分布变异性。跑完后把它存成表格和相关矩阵一起放进补充材料审稿人再问阈值问题直接甩这份表就行。6.2 怎么读扫描结果看平台期而不是看单一点拿到扫描表后不要盯着某一个稀疏度下的数字看要看整段曲线的走向。如果某个网络指标在10%到30%的稀疏度区间内保持单调平稳那么这段区间内的结论是可信的如果指标在某个稀疏度附近突然跳变比如平均度从5跳到15说明这个密度下网络结构发生了剧烈重构任何基于单一点的结论都不可信。我一般会在曲线平稳的区间中间选一个稀疏度作为主要分析报告值同时把扫描区间里所有点都跑一遍置换检验所有点都跟着做一遍。这样做的好处是最后报告的“组间差异显著”不再被绑定在一个具体的阈值上哪怕有人质疑阈值选择你也可以用整个区间的结果回应。从那以后我每次构建结构协变网络都强制自己先跑一遍5%到50%的稀疏度扫描把这个表当成网络质量的一部分来检查提前替自己踩掉阈值玄学这个坑。希望这个习惯对你也同样有用。本文还有配套的精品资源点击获取