简介面向脑影像与神经科学研究者的NBS脑网络分析工具包内含MATLAB环境下用于脑连接组统计比较的完整工具箱源码与示例数据重点解决功能/结构脑网络在组间差异分析中的多重比较校正问题。资源共68个文件以txt说明文档、m脚本和mat数据矩阵为主另有少量asv自动保存文件与nii脑图谱文件压缩包仅1.59MB便于下载后直接加入MATLAB路径使用。包中提供精神分裂症研究示例、AAL脑区模板及标签坐标以及NBS流程、统计建模等源码使用者可对照帮助文档快速理解基于网络统计的算法原理并在自己的fMRI、DTI或EEG数据上迁移应用。已有2722人学习下载适合掌握一定MATLAB基础、希望深入脑网络组分析方法的研究生或临床科研人员也可为撰写数据分析方法学部分提供参考。1. 一句话说清 NBS 在脑网络分析里的位置它解决的是多重比较不是连接矩阵本身NBS 全称 Network-Based Statistic是一个跑在 MATLAB 里的脑网络分析工具箱专门处理「上千条连接同时做统计检验」时的多重比较问题。做功能连接或结构连接的人都会撞上这个场景90 个脑区得到 4005 条边逐条边做 t 检验假阳性会多到没法看。NBS 的思路是不逐条边下结论而是先把边聚类成连通分量再对整个分量做置换检验最后给你一组校正过 p 值的显著连接子网络。它的适用场景很明确你有两组或多组被试的脑网络矩阵想知道组间在哪些连接上出现差异。GRETNA、DPABI、SPM 的 DCM 流程做完数据落在你手里之后NBS 就是那个把「网络差异」变成「可发表结论」的统计出口。2. 先看懂 NBS 的统计逻辑再做数据准备分量、置换检验与输入矩阵2.1 每条边都检验一遍会发生什么从单边 t 检验到连通分量脑网络分析的原始数据是连接矩阵以 AAL90 模板为例去掉对角线后共有 4005 条边。如果你对每条边单独做 t 检验显著性水平取 0.05那么即使两组数据完全是同一批被试复制出来的也会有大约 200 条边被误报为显著。这是多重比较问题最直白的形态。常见的应对方式有三种。第一种是 Bonferroni 校正把阈值除以检验次数4005 条边时要求 p 小于 1.2e-5功效低到离谱真实差异也会被压掉。第二种是 FDR控制错误发现率比 Bonferroni 温和但它把每条边当作独立假设处理完全不管脑网络的空间结构。第三种才是 NBS它假设真正的差异不会孤立出现在某一条边上而是集中在某个子网络里。于是先对每条边算一个统计量超过阈值的边被保留下来保留的边在脑区拓扑上会聚成若干个连通分量每个分量的大小或强度作为新的统计量再拿到置换检验里去校正。这样校正的是「分量」而不是「边」功效高很多代价是只能告诉你哪些脑区组成的子网络有差异说不清具体哪一条边铁定不同。这个逻辑决定了 NBS 的两个性格特点。第一个它对输入矩阵的对称性和节点顺序极其敏感因为连通分量的计算依赖真实的邻接关系。第二个它对阈值很敏感阈值太高可能一个分量都没有阈值太低整个网络变成一个巨大的分量置换检验又回到无差异。这两个点后面都会展开。2.2 输入矩阵的硬性要求对称矩阵、节点顺序、上三角向量化NBS 的输入不是任意形态的表格而是一批被试的矩阵每个被试一个节点数乘节点数的方阵。假设你有两批被试一批健康对照一批患者每个被试都算好了一张 90×90 的功能连接矩阵那么最终的数据形态就是 90×90×NN 是总被试数。这里最容易翻车的是节点顺序。NBS 不知道你的第 7 个节点是左侧额中回还是右侧楔前叶它只认矩阵下标。如果对照组用 AAL 模板顺序导出患者组用自定义 ROI 顺序导出两组同一个位置放的根本不是同一个脑区那么后续所有统计都是对着错误的边在做。连接矩阵里没有「名字」这个维度保存成 mat 文件之前一定要用同样的模板、同样的顺序重排。对称矩阵还有一个默认约定。90×90 的矩阵里第 i 行第 j 列和第 j 行第 i 列是同一个连接NBS 内部通常只取上三角或者下三角来做边级检验。你在准备数据时最好也遵循这个约定把矩阵显式对称化之后再取上三角。用代码摊平这一步很直接但容易被漏掉的是对角线。AAL 模板的对角线本身不是连接摊平之前要把对角线置零这部分我在第 3 章会给出完整的转换脚本。2.3 在动手之前把假设写成 GLM组别、协变量与对比向量NBS 的统计内核是通用线性模型所以在跑之前你脑子里得有一个明确的模型方程。最常见的是两组比较设计矩阵里有一列组别编码比如健康对照记为 0患者记为 1如果你还想控制年龄和性别就再放两列协变量。对比向量则告诉 NBS 你要检验的是哪一项比如只关心组别效应对比向量就是 [1 0 0]。这里有一个新手容易误解的地方NBS 的「对比」不是给你一个交互作用的便利菜单它就是你学线性模型时写的那个 C 向量。协变量建议先做去均值或者 zscore否则截距项和协变量之间会出现共线性置换检验跑出来要么报设计矩阵奇异要么结果不稳定。先把所有连续协变量中心化这比事后排查玄学问题省钱得多。另一个值得提前想清楚的问题是你的假设到底适不适合用 NBS。NBS 擅长的是无向网络上的子网络差异检验。如果你的数据是有向连接比如有效连接或 Granger 因果连通分量的定义会变得含糊常见做法是先对称化但这就丢掉了方向信息。如果研究问题本身要求保留方向NBS 不是最合适的工具倒不如回到逐边 FDR。3. 在 MATLAB 里跑通 NBS从连接矩阵到结果文件3.1 准备数据与 Fisher Z两行代码把相关矩阵变成正态可检验的连接值功能连接矩阵通常是皮尔逊相关系数取值范围在 -1 到 1 之间直接拿去检验并不理想。Fisher Z 变换能把相关系数映射到近似正态的尺度上是脑网络分析里的常规预处理。下面这段代码处理一批被试的原始相关矩阵输出一个被试数乘连接数的矩阵这是 NBS 最常吃的输入格式。% all_r: 90 x 90 x N 的原始相关矩阵N 为总被试数 % 初始化每条连接占一行总连接数为上三角个数 90*89/2 N size(all_r, 3); nnode size(all_r, 1); nedge nnode * (nnode - 1) / 2; allZ zeros(N, nedge); % 取上三角下标用于把所有被试摊平成 N x nedge [I, J] ind2sub([nnode nnode], find(triu(ones(nnode), 1))); for sub 1:N mat all_r(:, :, sub); % 对角线不是连接强制置零后再取上三角 mat(1:(nnode1):end) 0; vec mat(I (J - 1) * nnode); % Fisher Zatanh 即反双曲正切 allZ(sub, :) atanh(vec); end这段代码的核心是triu(ones(nnode), 1)产生一个只保留严格上三角的掩膜ind2sub把掩膜里每个 1 的位置换成行列下标再用线性索引I (J - 1) * nnode从矩阵里取出对应的边值。这样每个被试得到的是相同顺序的 4005 个值顺序完全由掩膜决定不会出现两组被试边顺序不一致的问题。atanh在 MATLAB 里就是 Fisher Z 变换。如果连接矩阵里出现绝对值等于 1 的相关值atanh 会得到无穷大通常是因为某个脑区时间序列完全相同时长太短。这种情况下先把该被试的矩阵检查一遍必要的时候删除异常脑区或做正则化不要带着 inf 进统计。3.2 最小批处理脚本NBS.run 的参数怎么填NBS 工具箱自带图形界面在 MATLAB 命令行输入NBS会打开主窗口。但批处理做多组对比、多阈值扫描时脚本化调用更靠谱。以下是 NBS 批处理接口最常见的调用方式参数名在不同小版本之间略有出入跑之前先help NBS.run核对一下本机字段名。% allZ: N x nedge每个被试一条连接向量 % group: N x 10 表示健康对照1 表示患者 % covariates: N x kk 个连续协变量使用前已做 zscore % 组装设计矩阵组别列 协变量列 design [group, covariates]; % 对比向量只检验组别效应长度为设计矩阵列数 con [1, zeros(1, size(covariates, 2))]; % 调用 NBS置换 10000 次边级阈值取 t3.2 % sizeextent 表示用分量包含的边数作为分量统计量 [nbs, res] NBS.run(con, design, allZ, [], ... test, ttest, ... thresh, 3.2, ... perm, 10000, ... size, extent);第四个参数我在上面传了空数组它在有些版本里用来放额外的被试级协变量矩阵。如果你只用设计矩阵里的协变量这个位置保持空即可。test指定检验类型两组比较用ttest多组或交互设计用ftest。thresh是边级 t 阈值决定哪些边先进入候选集合3.2 是我在 60 到 100 人样本里常用的起点。perm是置换次数正式分析至少 5000我一般直接上 10000。size决定分量怎么比大小extent数边数intensity把分量内所有边的统计量加权求和后者对弱效应更敏感但计算也更重。跑完后的res结构里主要看两样东西校正后的 p 值以及显著分量对应的连接矩阵。很多版本里res.nbs.pval存着 p 值res.nbs.con_mat里每个单元存一个分量矩阵。字段名在不同版本里确实有差异拿到结果先fieldnames(res)和fieldnames(res.nbs)扫一眼别凭记忆硬取字段。3.3 先在小样本上验证流程1000 次置换的快速跑通正式置换检验跑 10000 次在 90 节点、几十个被试的数据上通常要几分钟到十几分钟这取决于机器。第一次跑通流程时不建议直接上 10000先用 1000 次置换跑一遍确认数据维度、设计矩阵、对比向量都没有问题再改成 10000 挂机。我一般会在正式跑之前做一个最小自检把两组的连接矩阵随机打乱分组理论上不应该得到显著的子网络。如果打乱后 p 值小于 0.05说明流程有 bug常见于设计矩阵顺序和数据行顺序没对齐。NBS 的置换检验会反复重排分组标签如果你的allZ行顺序和group顺序不是同一个被试排列结果就是一团乱麻。这一点用随机分组的自检能很快暴露。另一个快速验证是 GUI 对照。脚本能跑通之后把同样一组变量加载到 NBS 界面里手动点一遍对照输出。脚本和 GUI 填的字段对不上时通常不是 NBS 坏了而是你脚本里某个参数的语义理解成另一个意思。花十分钟做这个对照能省下后面整晚的排查时间。4. NBS 最常见的 5 个坑现象、原因与解决办法4.1 跑出来的 p 值全是一节点顺序没对齐现象统计跑完所有分量的 p 值都是 1或者结果里根本找不到显著分量。原因两组数据虽然都是 90×90 的矩阵但节点顺序不一致。最常见的来源是一个组用 AAL 模板按文件名字母排序导出另一个组按 ROI 编号排序导出两边矩阵第 17 行对应的脑区不一样。NBS 完全不知道你在说什么它只会觉得两组在所有连接上都没差异。解决在拼接数据之前用同一份 ROI 编号表重排所有被试的矩阵。写一段脚本读入每个 ROI 的编号和名称确认两组第一个被试矩阵的行列顺序一致后再进入第 3 章的摊平流程。如果项目里有人换过模板这种坑特别隐蔽因为矩阵尺寸一模一样肉眼根本看不出来。4.2 阈值设太低整片网络全连通分量检验的玄学现象阈值取 1.5 时跑出一个覆盖几乎所有脑区的巨大分量p 值反而不显著阈值取 3.5 时又什么都没有。原因NBS 的分量统计量是「超过阈值的边连通成多大一片」。阈值太低时绝大多数边都超过阈值整个大脑连成一个分量置换检验里随机打乱分组也经常出现类似大小的分量校正后自然不显著。阈值本身就是检验的一部分这个跟逐边检验很不一样。解决不要只跑一个阈值就下结论。常见做法是取 2.5 到 4.0 之间几个阈值做扫描看显著子网络的核心连接是否稳定存在。如果只在某个很窄的阈值区间显著结论要写得更谨慎。分量显著但边数很少时也要回头检查是不是阈值卡在临界点附近。4.3 协变量没去均值设计矩阵共线性报警现象NBS 报设计矩阵秩亏或者加入年龄、性别协变量之后结果和完全不加协变量时差异巨大。原因设计矩阵里同时放了未中心化的连续协变量和组别编码连续变量均值和截距项高度相关导致设计矩阵近似奇异。置换检验里每次打乱分组都在和这个病态矩阵打交道统计量不稳定。解决所有连续协变量在进入设计矩阵前做 zscore 或至少减去均值。离散协变量比如性别用 0/1 编码即可。处理完之后检查一下设计矩阵的秩rank(design)应该等于列数。4.4 置换结果不可复现随机种子没固定现象同一个脚本跑两次p 值有时是 0.042有时是 0.038甚至显著与否都会变。原因NBS 的置换检验需要生成随机置换序列不同版本的工具箱对随机数种子的处理不一致。第二天重跑、或者换台机器跑结果自然有波动。置换次数越少波动越大。解决脚本开头固定随机种子例如rng(2024)。置换次数上到 10000 后p 值的随机波动会小很多但定种子仍然是个好习惯尤其当你需要复算审稿人要求的变体分析时。发表时把随机种子写进方法部分很多人会忽略这一点。4.5 显著分量只有一两条边结果的生物学解释要谨慎现象p 值小于 0.05但显著分量只包含 2 到 3 条边脑区也很零散。原因NBS 的功效来自聚合分量越小它对「网络层面差异」的支撑越弱。一个只有两条边的分量在统计上被校正为显著往往说明这两条边的效应确实很强但它并不代表一个稳定的子网络。解决报告时不要把「分量显著」和「整个子网络差异显著」划等号。建议同时报告分量的边数、涉及的脑区、效应方向并在讨论里承认这是局部连接差异。如果审稿人问起来一份变阈值扫描结果比任何话术都管用。5. 解读 NBS 结果并可视化p 值、显著分量与脑网络图5.1 结果对象里有哪些字段pval、con_mat 与分量编号拿到res之后第一件事不是画图而是把结果结构摸清楚。多数版本里结果会嵌在res.nbs下面几个核心字段几乎都有字段含义常见用途pval校正后的分量 p 值判断每个分量是否显著con_mat每个显著分量对应的连接矩阵后续可视化、定量描述test_stat每条边的原始检验统计量导出加权连接图sizes每个分量的边数或强度报告分量规模先跑一句fieldnames(res.nbs)确认你的版本里这些字段叫什么。这里有个容易踩的细节con_mat可能是一个 cell第几个 cell 对应第几个分量每个 cell 是一个节点数乘节点数的稀疏矩阵非零位置就是该分量经过的边。取出来之后用full()转成普通矩阵再写文件免得后面处理稀疏结构出错。p 值的解读也要小心。NBS 给的是 FWE 校正后的分量 p 值不是每条边的 p 值。一个分量显著说明「这个连通分量在置换分布下不太可能随机出现」不代表分量里每条边单独检验都显著。写论文时建议明确写「经过 NBS 校正发现一个包含 X 条边的显著子网络p0.03」不要写成「两条边显著」。5.2 把显著分量导出成 .edge/.node给 BrainNet Viewer 用的脚本BrainNet Viewer 是脑网络可视化的常用工具它需要两个文件.node文件描述脑区节点坐标.edge文件描述节点之间的连接。下面这段代码把显著分量矩阵写成 edge 文件权重可以填原始 t 统计量也可以只填 1 表示存在连接。% conMat: 90x90显著分量矩阵非零位置为显著连接 % tstat: 90x90每条边的 t 统计量用于显示权重 edgeMat zeros(size(conMat)); edgeMat(conMat 0) tstat(conMat 0); % 写出 .edge 文件BrainNet 要求空格分隔的矩阵 writematrix(edgeMat, sig_network.edge, Delimiter, space); % 写出 .node 文件需要先准备 90x3 的坐标矩阵 coords % 列顺序为 x y z模板不同坐标也不同通常来自模板自带模板文件 nodeData [(1:90), coords]; writematrix(nodeData, sig_network.node, Delimiter, space);writematrix从 R2019a 开始可用老版本用dlmwrite替代。.node文件的实际格式比这复杂一点BrainNet 还支持第五、第六列控制节点大小和颜色如果模板导出时已经带好这些列直接复用模板的原生 node 文件最省事。.edge文件只包含数字矩阵对角线必须为 0BrainNet 在读取时会忽略自连接。画图时有一个显示层面的建议把显著分量里的脑区按模块或半球着色比所有节点一个颜色更容易看出差异网络的拓扑。BrainNet 在 node 文件的颜色列里填 RGB 或者调色板索引都可以先出一张全默认图再慢慢调视觉细节。5.3 写进论文或报告前必须交代的信息阈值、置换次数、节点模板NBS 结果的可复现性严重依赖参数报告方法部分至少要写清楚五件事节点模板和节点数、边级阈值、置换次数、分量统计量是 extent 还是 intensity、是否对协变量做了中心化。这五条缺一条别人就没法复算你的结果。我在写方法部分时通常这样组织先说明连接矩阵是怎么算出来的皮尔逊相关还是偏相关有没有做 Fisher Z然后写 NBS 的配置比如「在边级 t 阈值 3.2、10000 次置换下采用基于分量大小的 NBS 进行 FWE 校正」最后补一句随机种子。阈值扫描的结果作为补充材料也很常见尤其当主阈值结果只有边缘显著时。结果部分建议放一张网络图加一张表格。表格列出显著分量的编号、p 值、边数、涉及的主要脑区比在正文里数着脑区名字挨个念清楚得多。如果同时做了多个阈值把阈值作为表格的列变量一眼能看出结果稳定性。6. 把 NBS 跑得更稳的最后一招变阈值扫描与置换次数验证6.1 用循环做变阈值扫描从 2.5 到 4.0 试一遍单独跑一个阈值的结果总让人心里没底。把阈值从 2.5 到 4.0 每 0.5 扫一次记录每个阈值下显著性分量数量和覆盖的边能在十几分钟内判断结论稳不稳。扫描结果如果只有 3.2 这一个点显著那这个结论大概率敏感写讨论时得多加限制如果 2.5 到 3.5 都指向同一组核心连接那就比较踏实。thrList 2.5:0.5:4.0; for k 1:numel(thrList) [~, res] NBS.run(con, design, allZ, [], ... test, ttest, ... thresh, thrList(k), ... perm, 5000, ... size, extent); fprintf(thresh%.1f pval, thrList(k)); % 打印第一个显著分量的 p 值没有显著分量就显示 NaN if ~isempty(res.nbs.pval) fprintf(%.3f\n, res.nbs.pval(1)); else fprintf(NaN\n); end end这个循环里每次置换次数用 5000 就够扫描阶段不需要追求精确的 p 值看趋势为主。扫描结果能直观回答一个常见审稿问题结果对阈值的选择是否敏感。6.2 我每次都会做的一件事固定种子并重跑一次 10000 置换正式结果出来之后我不会直接拿去写材料。先把脚本开头的rng固定住再把置换次数改成 10000重跑一遍看 p 值和分量组成有没有变化。这一步虽然耗时但能拦住绝大部分「换台机器结果就变」的翻车现场。置换检验本身随机固定种子保证可复现10000 次保证 p 值稳定到小数点后两位。我个人的习惯是最后把随机种子、版本号、阈值扫描结果一起整理进分析日志这个习惯已经帮我避免过不止一次「审稿人要复现数据但我自己也说不清当初怎么跑的」的尴尬。希望帮到你。本文还有配套的精品资源点击获取