脑网络分析NBS方法：从功能连接到组间差异统计检验的MATLAB实践-育师

简介：本资源是面向神经科学与认知研究者的MATLAB脑网络分析实践工具包，聚焦脑功能网络建模、连接属性计算与疾病相关网络差异检测，特别适用于fMRI/DTI数据驱动的脑连接组学研究初学者与进阶用户。压缩包含68个文件（1.59MB），涵盖13个核心MATLAB函数（如NBS.m、NBSfdr.m、get_components.m）、8个预处理与结果数据文件（.mat格式，含SchizophreniaExample设计矩阵、AAL脑区模板及COG坐标）、44个说明与帮助文档（.txt），以及.nii结构模板和.asv备份脚本，构成从数据加载、NBS统计检验、FDR校正到可视化分析的完整闭环流程。已有2712人学习下载，资源内置精神分裂症案例实操路径、模块化帮助文档体系（help_*.txt）及图标化GUI支持（NBSicon.m、readUI.m），可直接复用代码框架开展临床组间比较、小世界属性分析或社区结构识别，显著降低BCT工具箱上手门槛与调试成本。

1. 项目概述：从数据到洞察的神经科学桥梁

如果你正在处理功能磁共振成像（fMRI）数据，试图从一堆时间序列中找出大脑不同区域之间“谁和谁在聊天”，并且正在为如何从统计上验证这些连接的显著性而头疼，那么你很可能已经听说过或正在寻找NBS（Network-Based Statistic）。这个项目标题“NBS_脑功能网络_脑网络分析_脑连接工具箱_脑网络_matlab_”精准地指向了神经科学，特别是脑影像数据分析中的一个核心痛点：群体水平的脑网络差异统计检验。简单来说，NBS不是一个单一的软件，而是一套在MATLAB环境中实现的、用于分析脑功能或结构连接网络组间差异的统计方法框架。它解决了一个非常具体但至关重要的问题——当我们比较两组人（例如患者组 vs. 健康对照组）的大脑连接网络时，如何判断哪些连接上的差异不是随机噪声，而是具有统计学意义的。

传统的脑网络分析，比如基于图论的度、聚类系数等全局或局部指标，虽然能给出网络整体的属性差异，但往往丢失了“差异具体发生在哪些连接上”的空间信息。而如果对成千上万个连接（一个包含90个脑区的网络就有4005个可能的连接）逐一进行统计检验，则会面临严峻的多重比较校正问题，可能导致大量假阴性（发现不了真实差异）。NBS巧妙地提供了一种折中方案：它不单独检验每条边，而是寻找在组间存在显著差异的“连接子网络”。其核心思想是，疾病或某种条件对大脑的影响可能不是孤立的某一条连接，而是一个相互关联的“电路”或“模块”。NBS通过基于网络的置换检验来识别这些连续的、空间上聚集的差异连接集合，并控制家族错误率。

对于研究者而言，拥有NBS工具箱意味着你获得了一把利器。它通常与脑网络构建的上下游流程紧密集成：从原始的fMRI时间序列预处理、去噪、头动校正，到定义脑区（如使用AAL、Desikan-Killiany等脑图谱），计算区域间的相关性（如皮尔逊相关、偏相关）以构建个体水平的连接矩阵，再到使用NBS进行组间统计比较，最后可视化结果。整个过程高度依赖MATLAB及其强大的矩阵运算和统计工具箱。因此，这个标题背后隐含的是一条完整的分析流水线，而NBS是这条流水线上最关键、也最具方法论挑战性的一环。它适合神经科学、心理学、生物医学工程等领域的研究生、博士后以及临床研究人员，帮助他们从复杂的大脑数据中提取出可靠且可解释的生物学发现。

2. NBS方法的核心原理与统计逻辑拆解

要正确使用NBS，绝不能把它当作一个黑箱。理解其背后的统计逻辑，对于正确设置参数、合理解读结果至关重要。NBS的本质是一种基于置换检验的非参数统计方法，专门用于处理高维、稀疏且相互依赖的脑网络数据。

2.1 为什么是“基于网络”的统计？

想象一下，我们比较两组大脑网络，得到了一个差异矩阵，其中每个元素代表对应连接在两组间的t值或F值。如果直接对矩阵中所有元素（边）进行独立检验，并采用Bonferroni等传统校正方法，阈值会变得极其严格（例如，4005次比较，校正后p<0.05意味着单边p值需小于0.0000125）。这会导致统计效力极低，许多真实的微弱但一致的连接差异会被淹没。NBS提出了一个不同的假设：有意义的差异往往不是随机散布的，而是形成连通的子图。比如，与特定认知功能相关的几个脑区，它们之间的连接可能协同变化。

因此，NBS的检验单位不是“单一边”，而是“连通的边集合”。它首先需要一个用户定义的初级阈值（通常是一个t值或F值阈值，如t> 3.0），用于筛选出那些初步看来组间差异较大的边。然后，在这些超过初级阈值的边中，寻找所有连通的成分（即子网络）。每一个连通成分的大小（包含的边数或节点数）被作为一个检验统计量。最后，通过置换检验（随机打乱组别标签成千上万次）来评估观察到的连通成分大小的极端程度，从而得到一个校正后的p值。

2.2 置换检验：构建零分布的关键

置换检验是NBS的引擎。它的基本步骤如下：

计算原始统计量：基于真实组别标签，计算每个连接的组间统计量（如双样本t检验的t值），应用初级阈值，找出最大的连通成分大小S_obs。
构建零假设分布：将两组被试的标签随机打乱（例如，100次观测中随机分配50个为“组A”，50个为“组B”），但保持每个被试的数据不变。然后，基于这个随机标签，重复步骤1，计算出一个“在零假设（组间无差异）下”可能出现的最大连通成分大小S_perm。
重复置换：将步骤2重复数千次（如5000次或10000次），得到一组S_perm的集合，这就是零分布。
计算校正p值：比较S_obs与这个零分布。校正后的p值等于零分布中大于或等于S_obs的置换次数所占的比例。例如，在5000次置换中，有50次得到的S_perm>=S_obs，那么校正p值就是 50/5000 = 0.01。

这个过程控制了由于在连通成分水平进行检验而产生的家族错误率。它回答的问题是：“在组间根本没有真实差异的情况下，随机出现一个规模如我们观察到的连通差异子网络的概率有多大？”

2.3 初级阈值的选择：艺术与科学的平衡

初级阈值（t或F阈值）是NBS分析中最需要经验和谨慎对待的参数之一，它没有绝对的金标准。

阈值过高：可能过滤掉所有真实的差异边，导致找不到任何连通成分，统计效力为零。
阈值过低：会导致大量噪声边超过阈值，这些边可能连接起来形成一个巨大的、但生物学意义模糊的“伪网络”，使得置换检验的零分布向右偏移，最终难以获得显著结果（因为大网络在零假设下也容易出现）。

实操心得：初级阈值的选择往往需要结合先验知识和探索性分析。一个常见的策略是参考文献中类似研究使用的阈值（例如，t=2.5到3.5）。更稳健的做法是进行一个初步的、宽松的阈值分析，观察差异边的空间分布模式，然后选择一个能使差异模式看起来“合理”（例如，集中在特定功能网络内，如默认模式网络）且连通成分大小适中的阈值。也可以尝试一个阈值范围，观察结果的稳定性。记住，这个阈值本质上是一个“聚类形成”的过滤器，而不是最终的显著性判断标准，最终的显著性由置换检验的p值决定。

3. 脑网络构建的前置流程详解

NBS是分析的最后一步，但它的输入——个体水平的脑连接矩阵——的质量直接决定了分析的成败。一个标准的脑功能网络构建流程包括以下核心环节，每个环节都有其“坑点”。

3.1 数据预处理与质量控制在MATLAB中的实现

原始的fMRI数据必须经过严格的预处理，才能用于计算有意义的连接。这一系列操作通常在SPM、DPARSF、CONN或fMRIPrep等工具中完成，但理解和在MATLAB中检查关键步骤至关重要。

时间层校正与头动校正：fMRI扫描是逐层进行的，不同层面的采集时间有微小差异，需进行插值校正。头动校正则通过刚体变换将每个时间点的图像对齐到一个参考图像（通常是第一幅或平均图像）。在MATLAB中，你可以通过检查SPM生成的rp_*.txt头动参数文件来评估数据质量。通常，我们会排除平动超过3毫米或转动超过3度的被试。
空间标准化：将个体大脑图像配准到一个标准空间（如MNI空间），使得不同被试的脑区位置具有可比性。这一步的配准精度直接影响后续脑区提取的准确性。
空间平滑：使用高斯核（如6-8mm FWHM）对图像进行平滑，可以提高信噪比，但过度平滑会降低空间分辨率。在MATLAB中，可以使用SPM的spm_smooth函数。
去噪：这是功能连接分析中最关键也最复杂的步骤之一。需要移除的信号源包括：
- 生理噪声：通过记录的心跳和呼吸信号进行回归，或使用CompCor等方法从脑脊液和白质信号中提取噪声成分。
- 头动效应：不仅回归头动参数（6个），有时还需回归其逐时间点的一阶导数，甚至采用“擦除”策略（如scrubbing）移除头动过大的时间点。
- 全局信号：是否回归全局平均信号存在巨大争议。回归它可能引入负相关，但能有效移除一些全脑范围的噪声。需根据研究问题和领域惯例谨慎选择。

注意事项：预处理流程的参数选择（如平滑核大小、去噪策略）应在整个研究的所有被试中保持一致。强烈建议在组水平分析前，对每个被试的预处理后时间序列进行质量检查，例如绘制各脑区时间序列的图，查看是否有异常的尖峰或漂移。

3.2 脑区定义与时间序列提取

构建网络需要节点。节点通常由脑图谱定义。

选择脑图谱：常用的包括AAL（90或116区）、Desikan-Killiany（84区）、Harvard-Oxford（96区）等。选择时需考虑其分区是否与你的研究假设相关（例如，是否精细区分了某个特定皮层下核团）。

提取平均时间序列：将每个被试预处理后的fMRI数据，根据脑图谱的掩模，提取每个脑区内所有体素时间序列的平均值。在MATLAB中，这通常涉及读取NIFTI图像和图谱掩模文件，然后进行矩阵运算。例如，使用spm_read_vols读取数据，再对每个脑区标签内的体素求平均。

% 伪代码示例：提取单个被试所有脑区时间序列 nii_data = spm_read_vols(‘func.nii’); % 4D数据: [x, y, z, t] atlas = spm_read_vols(‘atlas.nii’); % 3D图谱，值为脑区编号 num_regions = max(atlas(:)); num_timepoints = size(nii_data, 4); time_series = zeros(num_timepoints, num_regions); for r = 1:num_regions mask = (atlas == r); for t = 1:num_timepoints vol = nii_data(:,:,:,t); time_series(t, r) = mean(vol(mask)); end end

可能遇到的问题：部分脑区（尤其是小脑、边缘系统的小核团）在某些被试的标准化图像中可能缺失或只有极少量体素，导致时间序列信噪比极低或为NaN。需要制定规则处理，如体素数量少于10的脑区，该被试该脑区数据标记为缺失，在后续连接计算中需相应处理。

3.3 功能连接矩阵的计算

有了每个脑区的时间序列，下一步是计算它们两两之间的“连接强度”，即功能连接。最常用的度量是皮尔逊相关系数。

皮尔逊相关：计算简单，解释直观。在MATLAB中，使用corrcoef函数即可。conn_mat = corrcoef(time_series);会得到一个对称的N x N矩阵，对角线为1。
考虑其他度量：
- 偏相关：在控制其他所有脑区影响的前提下，衡量两个脑区之间的直接关联。更能反映“直接连接”，但计算更复杂且对数据长度和信噪比要求更高。可以使用partialcorr函数或基于逆协方差矩阵的方法（如Graphical Lasso）。
- 相位同步性：对于研究脑振荡同步的研究，可能需要在特定频带（如Alpha波）计算连接。
矩阵后处理：得到的相关矩阵通常需要进一步处理。
- Fisher z变换：由于相关系数的分布不是正态的，尤其在高相关时，通常将其转换为近似正态分布的Fisher‘s z值：z = 0.5 * log((1+r)/(1-r))。这在后续的组水平统计（如t检验）中更合适。
- 阈值化（二值化或加权）：对于图论分析，有时需要将连续的相关矩阵转换为二值邻接矩阵（连接存在或不存在）。这需要设定一个相关性阈值。阈值的选择同样敏感，可采用绝对阈值（如r > 0.3）、比例阈值（保留前10%最强的连接）或基于网络属性的阈值（如图的密度）。注意：NBS分析通常直接使用连续的连接强度值（如z值）作为输入，其初级阈值是基于统计量（t值），而非直接的相关性阈值。

4. 在MATLAB中实施NBS分析：一步步实操指南

假设你已经准备好了两组被试（如HC组和PAT组）的所有个体连接矩阵（均为N x N的对称矩阵，已进行Fisher z变换），并存储在MATLAB工作区中。我们将使用NBS工具箱（可从其官网获取）进行实操。

4.1 环境准备与数据组织

首先，确保NBS工具箱路径已添加到MATLAB。addpath(genpath(‘/your/path/to/NBS’))。数据组织是关键，NBS通常期望一种简单的结构。

组1数据：一个单元格数组HC_mats，其中HC_mats{i}是第i个健康对照被试的N x N连接矩阵。
组2数据：类似地，PAT_mats。
设计矩阵与对比：你需要构建一个描述所有被试组别的设计矩阵design。例如，有20个HC和20个PAT，则design是一个40 x 2的矩阵。第一列通常是全1的截距项，第二列是组别指示变量（如HC为0，PAT为1）。对比向量contrast则指定你要检验的效应，对于组间比较，contrast = [0, 1]表示检验设计矩阵第二列的系数（即组别效应）。

4.2 配置NBS并运行分析

NBS提供了一个图形用户界面（GUI），但对于可重复研究和批量处理，更推荐使用脚本调用。核心函数是nbs_bct（如果使用Brain Connectivity Toolbox的格式）或直接使用NBS的主函数。

下面是一个典型的脚本示例：

% 1. 准备数据 % 假设我们已经将40个矩阵加载到两个元胞数组中 group{1} = HC_mats; % 组1: 健康对照 group{2} = PAT_mats; % 组2: 患者 % 2. 设置参数 NBS.thresh = 3.1; % 初级t值阈值，这是一个需要调整的关键参数 NBS.k = 10000; % 置换检验的次数，建议至少5000次 NBS.tail = ‘both’; % 检验方向：‘both’（双尾）， ‘left’， ‘right’ NBS.alpha = 0.05; % 显著性水平 NBS.exchange = []; % 置换块定义，对于独立样本t检验，留空即可 NBS.contrast = [0, 1]; % 对比向量，检验组别差异 NBS.design = design; % 设计矩阵 NBS.node_coor = coor; % (可选) 节点的三维坐标（Nx3矩阵），用于可视化 NBS.node_label = labels; % (可选) 节点标签 % 3. 运行NBS % 注意：NBS的输入数据格式要求可能因版本而异。常见的是将所有矩阵堆叠成一个3D矩阵 (N x N x Subject) % 假设我们已将group{1}和group{2}合并成一个3D矩阵‘all_mats’ [nbs_stats, nbs_net, nbs_mat] = nbs_stats(all_mats, design, contrast, NBS.thresh, NBS.k, NBS.tail, NBS.exchange); % 4. 查看结果 % nbs_stats 包含每个显著子网络的信息：大小、p值等 % nbs_net 是一个元胞数组，每个元胞包含一个显著子网络的边信息 % nbs_mat 是经过NBS分析后得到的显著性矩阵（0/1，表示边是否属于某个显著子网络） disp(nbs_stats);

运行后，NBS会输出在给定的初级阈值下，通过置换检验发现的任何显著连通子网络。nbs_stats会列出每个子网络包含的边数、节点数以及经过置换检验校正的族系错误率p值。这个p值才是最终报告的依据。

4.3 结果可视化与解读

发现显著子网络后，可视化至关重要。

连接矩阵可视化：使用imagesc或heatmap函数显示nbs_mat，可以直观看到差异连接集中在哪些脑区对之间。
脑网络图可视化：这是更直观的方式。你需要节点的3D坐标（MNI坐标）和标签。
- 使用BrainNet Viewer：这是一个非常流行的MATLAB脑网络可视化工具。你可以将NBS输出的显著边列表（nbs_net）和节点信息保存为.node和.edge文件，然后用BrainNet Viewer加载并渲染出漂亮的3D大脑图形。
- 自定义绘图：你也可以使用scatter3绘制节点，用plot3绘制连接边，并通过边的颜色或粗细来编码差异的强度（如t值大小）。
结果解读：
- 报告内容：必须报告初级阈值（t=3.1）、置换次数（10000）、检验方向（双尾）以及每个显著子网络的校正p值、包含的边和节点。
- 生物学解释：结合子网络涉及的脑区，从神经科学角度进行解释。例如，“发现一个主要涉及前额叶和顶叶脑区的子网络在患者组中连接减弱，这可能与执行功能缺损有关”。
- 谨慎因果推断：功能连接差异不代表结构损伤或直接的因果影响。它反映的是脑区活动模式的协同性变化。

5. 常见问题、排查技巧与高级考量

在实际操作中，你几乎一定会遇到各种问题。以下是一些典型场景及解决思路。

5.1 NBS分析没有发现任何显著子网络

这是最常见的问题。可能的原因和排查步骤：

数据质量或预处理问题：这是根源。回头检查个体连接矩阵的质量。计算每个被试连接矩阵的平均连接强度或某个已知网络（如默认模式网络）的内部连接强度，看组间是否有肉眼可见的趋势？检查时间序列的信噪比。
初级阈值设置不当：阈值可能太高了。尝试逐步降低阈值（如从3.5降到2.5，步长0.2），观察是否在某个阈值下开始出现连通成分。注意，过低的阈值会产生巨大但无意义的网络，其校正p值可能仍然不显著。
效应量本身很小：也许真实的组间差异非常微弱，需要更大的样本量才能检测到。可以进行一个事后效力分析，估算在当前样本量和数据变异下，能检测到多大效应量的差异。
置换检验次数不足：理论上，次数越多，p值估计越精确。但通常5000-10000次对于alpha=0.05已经足够。增加次数（如到20000次）主要影响p值的精度（例如，是0.048还是0.052），而不会使一个完全不显著的结果变得显著。
连接度量或去噪策略不合适：尝试不同的功能连接度量（如偏相关），或调整预处理中的去噪策略（如是否回归全局信号）。

5.2 结果不稳定：改变初级阈值，结果变化巨大

这提示结果可能不够稳健。

敏感性分析：系统地报告在一个合理的阈值范围内（如t从2.8到3.5），显著子网络的出现和基本拓扑结构是否保持相对稳定。如果只在非常狭窄的阈值下出现，结果的可靠性存疑。
聚焦先验假设：如果你有很强的先验假设（例如，只关注默认模式网络），可以定义感兴趣的子网络（ROI），仅在这个子网络内部进行NBS分析，这能提高统计效力并减少多重比较负担。

5.3 如何处理协变量（如年龄、性别）？

NBS的基本框架是双样本t检验，但研究中常常需要控制协变量。这时，你需要使用广义线性模型（GLM）框架。

构建设计矩阵：在设计矩阵design中，除了组别列，加入协变量列（如年龄、性别、头动平均FD等）。确保连续变量已标准化（z-score），便于解释。
设置对比向量：对比向量contrast需要精心定义以检验你感兴趣的效应。例如，在控制了年龄和性别后检验组别差异，如果设计矩阵是[截距，组别，年龄，性别]，那么对比向量应为[0, 1, 0, 0]。
使用支持GLM的NBS版本或函数：确保你使用的NBS工具支持GLM。在调用函数时，正确传入包含协变量的设计矩阵和对应的对比向量。

5.4 NBS与FDR、TFCE等其他校正方法的比较

vs. 基于FDR的边水平校正：FDR（错误发现率）直接对每条边进行校正。它在差异连接分散且独立时更有效。NBS则在差异连接聚集形成子网络时更有力。两者互补，可以同时进行。如果NBS发现了子网络，可以再查看该子网络内各条边的FDR校正p值。
vs. TFCE：TFCE（阈值无关的簇增强）是另一种流行的体素水平聚类校正方法，也有用于网络分析的变体。与NBS需要预设初级阈值不同，TFCE整合了所有可能的阈值信息，被认为对阈值选择更稳健。但在脑网络分析中的应用和工具支持不如NBS成熟。

5.5 从功能连接迈向有效连接

NBS分析的是静态功能连接，反映的是脑区活动的“相关性”，无法指明信息流向。如果你想探究组间在“因果”或“有效连接”上的差异，需要考虑其他模型，如动态因果模型（DCM）、格兰杰因果分析（GCA）或结构方程模型（SEM）。这些方法更为复杂，通常需要更强烈的假设和更精细的模型设定。它们的结果（如连接方向上的组间差异）可以与NBS发现的差异网络结合起来，提供更丰富的解释。

最后，一个至关重要的习惯是：代码和数据的可重复性。保存完整的MATLAB脚本，记录所有参数（初级阈值、置换次数、预处理步骤版本），并将中间数据和最终结果妥善归档。神经科学领域正日益强调研究的可重复性，清晰透明的分析流程是高质量研究的基石。

本文还有配套的精品资源，点击获取