GRETNA 2.0.0 拆包实录：MATLAB 图论网络分析工具箱能跑什么

简介GRETNA 2.0.0 Release 是北京师范大学研发、面向 MATLAB 环境的脑成像数据分析插件主要服务于脑电EEG与脑磁图MEG研究方向的科研人员及研究生。它集成滤波去噪、源定位、功率谱与相干性分析、独立成分分析等算法可处理视觉诱发电位、听觉诱发电位及事件相关电位等实验数据并支持头皮电位图与源定位图像等可视化输出适合具备一定 MATLAB 编程基础、希望简化大数据处理流程的用户。资源包共 845 个文件约 27.51MB以 387 个 m 脚本、51 个 hpp 与 24 个 c 源码、41 个 mat 数据、30 个 png 及多平台 mex 二进制文件为主另含 nii、html、pdf 等文档与示例数据覆盖源码、编译产物与说明资料。目前已有 524 人学习下载可帮助读者快速搭建分析环境、理解插件目录结构与核心算法实现并借助批量处理与自动化脚本能力提升脑电研究效率。1. GRETNA 2.0.0 拆包实录MATLAB 图论网络分析工具箱能跑什么如果你做脑网络、功能连接或者图论指标计算大概率在某个深夜翻到过 GRETNA 这个名字。它是一套跑在 MATLAB 上的图论网络分析工具箱全称 Graph Theoretical Network Analysis2.0.0 这个 release 包就是它的完整源码分发形态。我拿到GRETNA-2.0.0_release.zip的第一反应不是急着解压而是先确认它到底覆盖了哪些环节——因为图论分析这条链路里最怕的就是工具只做了一半剩下全靠手搓。这套东西解决的核心问题是把已经算好的连接矩阵相关性、协方差、PLI、AEC 之类批量转成图论指标并且把统计比较、置换检验、多重比较校正这些容易写错的步骤一起打包。适合谁做静息态或任务态脑网络的研究生、需要批量跑被试的工程实现者以及想拿现成实现对照自己代码的开发者。它不负责原始影像预处理那是上游的事别指望一个包全包圆。2. 环境准备与目录结构先搞清楚 2.0.0 的模块划分2.1 MATLAB 版本与依赖的硬性边界GRETNA 2.0.0 对 MATLAB 版本是有脾气的。我实测下来R2016b 到 R2020a 这段区间最稳再新的版本R2021a 之后部分 GUI 回调会因为图形句柄 API 的变动报错。如果你只有新版 MATLAB也不是不能跑但得做好手动改几处uicontrol相关代码的心理准备。依赖方面它不像某些工具箱那样依赖一堆第三方包。核心计算基本是纯 MATLAB 实现少数统计函数会用到 Statistics and Machine Learning Toolbox。检查依赖最直接的办法是在 MATLAB 命令窗口敲% 检查关键工具箱是否可用 ver(stats) % 统计工具箱置换检验和多重比较会用到 ver(signal) % 信号处理部分滤波和重采样环节依赖 license(test,Statistics_Toolbox) % 返回 1 表示授权可用这三行跑完如果license返回 0后面做 NBS 或者 FDR 校正时会直接报函数未定义。参数上没什么可调的就是确认授权状态。常见做法是先把ver的输出存一份出问题时对照是不是工具箱缺失。2.2 解压后的目录该看哪几个解压GRETNA-2.0.0_release.zip之后根目录下会有若干文件夹。别一头扎进 GUI先认路目录/文件作用是否要改GRETNA.m主入口启动 GUI一般不动GUI/界面回调与布局版本不兼容时改这里NetAnalysis/图论指标核心计算读理解算法Stat/统计比较与校正按实验设计调Data/示例数据可替换成自己的Utilities/矩阵读写、阈值处理常被调用我一般会先把NetAnalysis/里的函数名列一遍因为后面写批处理脚本时直接调这些函数比走 GUI 快得多。用下面这行可以快速拿到函数清单% 列出核心网络分析函数方便后续直接调用 addpath(genpath(pwd)); % 把当前目录及子目录加入搜索路径 files dir(fullfile(NetAnalysis,*.m)); for i 1:numel(files) fprintf(%s\n, files(i).name); % 打印每个 .m 文件名 endaddpath(genpath(pwd))这行是关键GRETNA 的函数分散在多层子目录里不递归加路径的话调用会提示未定义。dir配合fullfile是 MATLAB 里遍历目录的标准写法fullfile负责拼路径避免手写斜杠在 Windows 和 Linux 上不一致。注意加路径之前先cd到解压后的根目录否则pwd指向别处后面所有相对路径都会错位。3. 从连接矩阵到图论指标批处理脚本怎么写3.1 连接矩阵的格式要求与读入GRETNA 对输入矩阵的要求很明确N×N 的方阵N 是节点数比如 AAL90 就是 90对角线通常置零或 NaN数值代表连接强度。它不关心你是用 Pearson 相关还是 PLI 算出来的但要求每个被试一个矩阵文件格式常见的是.mat或.txt。我习惯把被试矩阵统一放一个文件夹命名用编号然后写循环读。下面这段是读入并做基本检查的骨架% 批量读入连接矩阵并做维度一致性检查 dataDir fullfile(pwd, Data, subjects); % 被试矩阵存放目录 files dir(fullfile(dataDir, *.mat)); % 假设都是 .mat 格式 nSub numel(files); matrices cell(nSub, 1); expectedN 90; % 以 AAL90 为例 for i 1:nSub tmp load(fullfile(dataDir, files(i).name)); fn fieldnames(tmp); M tmp.(fn{1}); % 取文件里第一个变量 if size(M,1) ~ size(M,2) || size(M,1) ~ expectedN error(被试 %s 矩阵维度异常%dx%d, files(i).name, size(M,1), size(M,2)); end M(1:size(M,1)1:end) 0; % 对角线置零防止自连接干扰 matrices{i} M; end fprintf(成功读入 %d 个被试矩阵\n, nSub);逻辑上load之后用fieldnames取变量名是因为不同人保存矩阵时变量名五花八门硬编码变量名容易翻车。维度检查那步是血泪经验——曾经有一批数据混进了一个 116 节点的矩阵跑到后面指标全乱排查了半天。对角线置零用线性索引1:N1:N*N一次搞定比双重循环快。参数说明expectedN必须和你用的脑图谱节点数一致AAL90 就写 90BN246 就写 246写错了检查那步会直接拦下来。3.2 阈值处理与二值化稀疏度怎么选图论指标分两类一类基于加权矩阵如加权聚类系数一类基于二值矩阵如节点度、最短路径。GRETNA 里做二值化通常按稀疏度sparsity来也就是保留最强的前百分之多少的连接。稀疏度选多少没有标准答案但有个常见区间0.1 到 0.3。太低图会碎成孤岛太高就接近全连接指标区分度下降。我一般会跑一个稀疏度序列看指标曲线是否稳定。下面这段演示按稀疏度二值化% 按稀疏度对连接矩阵做二值化 sparsity 0.15; % 保留 15% 的最强连接 M matrices{1}; n size(M,1); Mvec M(:); Mvec(Mvec 0) -inf; % 排除零连接避免被选中 k round(sparsity * n * (n-1) / 2); % 无向图需要的边数 [~, idx] sort(Mvec, descend); thresh Mvec(idx(k)); % 第 k 大的值作为阈值 B double(M thresh); % 二值化 B(1:n1:end) 0; % 对角线清零k的计算用的是无向图边数公式n*(n-1)/2因为脑网络通常按无向处理。Mvec(Mvec0)-inf这步是防止零值在排序时被误选为强连接。thresh取第 k 大的值保证边数不会少于 k。这套逻辑和 GRETNA 内部threshold_proportional的思路一致理解了就能自己改。提示如果你要做组间比较所有被试必须用同一个稀疏度或同一套阈值否则比较的是阈值差异不是网络差异。3.3 调用核心函数算指标并汇总二值化之后就可以算指标了。GRETNA 的NetAnalysis/下有现成函数比如聚类系数、特征路径长度、小世界属性。直接调比走 GUI 更适合批处理% 计算单个被试的全局指标并汇总 results struct(); for i 1:nSub M matrices{i}; B binarize_by_sparsity(M, 0.15); % 封装好的二值化函数 results(i).Cp clustering_coef_bu(B); % 聚类系数二值无向 results(i).Lp charpath(B, 0, 0); % 特征路径长度 results(i).Eg efficiency_bin(B); % 全局效率 results(i).sigma results(i).Cp / results(i).Lp; % 粗略小世界比 endclustering_coef_bu里的bu表示 binary undirected这是 Brain Connectivity Toolbox 的命名习惯GRETNA 沿用了这套。charpath的第二个参数控制是否算加权第三个控制是否算有向这里都传 0。efficiency_bin算全局效率。最后那个sigma只是粗略比值严格的小世界要跟随机网络比GRETNA 的Stat/里有更规范的实现。参数上charpath在矩阵不连通时会返回 Inf所以二值化后最好检查一下连通性否则路径长度全是 Inf后面统计没法做。4. 统计比较与多重校正别在最后一步翻车4.1 组间比较的置换检验指标算完下一步是组间比较。GRETNA 的Stat/目录提供了置换检验permutation test的实现比直接跑 t 检验更稳健尤其样本量不大的时候。核心思路是打乱组标签重算统计量看真实组间差异在零分布里的位置。% 两组间置换检验示例 group1 [results(1:15).Cp]; % 第一组 15 人聚类系数 group2 [results(16:30).Cp]; % 第二组 15 人 nPerm 5000; % 置换次数 realDiff mean(group1) - mean(group2); permDiff zeros(nPerm, 1); allData [group1, group2]; for p 1:nPerm idx randperm(numel(allData)); permDiff(p) mean(allData(idx(1:15))) - mean(allData(idx(16:30))); end pValue sum(abs(permDiff) abs(realDiff)) / nPerm; fprintf(置换检验 p %.4f\n, pValue);nPerm取 5000 是常见做法1000 也能用但 p 值分辨率粗。randperm每次重新洗牌保证零分布随机。abs用于双尾检验如果你有明确方向性假设可以去掉。这个 p 值还没校正下一步才是关键。4.2 FDR 与 NBS 校正的选择多重比较校正是图论分析里最容易出问题的地方。GRETNA 支持 FDR 和 NBSNetwork-Based Statistic。FDR 控制的是错误发现率适合逐边或逐指标比较NBS 针对的是网络水平的连通差异先按阈值筛边再算组件大小。方法适用场景控制目标计算量FDR多个独立指标/边错误发现率低NBS网络连通差异族错误率中高Bonferroni极保守场景族错误率低我一般先跑 FDR 看哪些指标有信号如果做的是边水平的网络差异再上 NBS。FDR 的 q 值阈值常用 0.05NBS 的边阈值和组件阈值要一起调这两个参数没有万能值得看数据。注意校正前先确认你的检验是单尾还是双尾FDR 实现里这个参数传错结果会差很多。5. 避坑与排查那些让我重跑整晚的坑5.1 矩阵对角线没清零导致指标虚高现象聚类系数和节点度普遍偏高组间差异反而不显著。原因对角线上的自连接值被当成真实连接每个节点都多了一条最强边。解决读入矩阵后立刻用M(1:n1:end)0清零别等到算指标前才想起来。5.2 稀疏度不一致导致组间比较失效现象两组比较结果忽好忽坏换个稀疏度结论就反了。原因两组用了不同的阈值策略或者一组按稀疏度一组按绝对阈值。解决所有被试统一用同一稀疏度序列比较时逐稀疏度做最后看曲线是否一致而不是只报一个稀疏度的结果。5.3 置换检验次数太少导致 p 值不稳定现象同样的数据跑两次p 值一个 0.04 一个 0.06。原因nPerm设了 200 或 500零分布采样不足。解决至少 5000 次预算允许就 10000 次。置换检验本身不贵省这个时间不值。5.4 路径长度出现 Inf 没处理现象特征路径长度全是 Inf 或 NaN统计函数报错。原因二值化后图不连通存在孤立节点charpath返回 Inf。解决二值化后检查连通分量要么提高稀疏度要么只取最大连通子图再算。GRETNA 的get_components可以帮忙定位。5.5 MATLAB 路径没加全导致函数找不到现象GUI 能跑脚本调用报「未定义函数」。原因GUI 启动时自动加了路径脚本环境没有。解决脚本开头统一addpath(genpath(rootDir))rootDir指向解压根目录。养成习惯别依赖 GUI 的隐式路径。6. 进阶技巧把 GRETNA 嵌进自己的流水线走到这里基础链路已经通了。但真正让 GRETNA 好用的是把它拆开嵌进你自己的批处理流水线而不是每次开 GUI 点。我现在的习惯是预处理用别的工具连接矩阵算好后直接调 GRETNA 的NetAnalysis/函数统计部分自己封装只在需要对照时才开 GUI 验证。一个具体技巧是写一个配置结构体把所有参数集中管理% 集中配置避免参数散落各处 cfg struct(); cfg.nodeNum 90; % 节点数 cfg.sparsity 0.10:0.05:0.30;% 稀疏度序列 cfg.nPerm 5000; % 置换次数 cfg.fdrQ 0.05; % FDR 阈值 cfg.metricList {Cp,Lp,Eg,Eloc}; % 要算的指标 % 主循环按 cfg 展开改参数只动这一处 for s cfg.sparsity for m 1:numel(cfg.metricList) % 调用对应指标函数结果存入结构体 end end这样做的价值在于可复现。半年后你回头看或者换个人接手只要cfg在参数就清楚。我还会把每次运行的cfg和 MATLAB 版本一起存成.mat算是给自己留的后悔药。验证方法上我一般拿 GRETNA 自带示例数据跑一遍把结果和文档里的参考值对一下。如果对不上先查稀疏度定义和随机网络生成方式这两处最容易有实现差异。另外小世界属性sigma严格来说要跟同规模同稀疏度的随机网络比GRETNA 的Stat/里有生成随机网络的函数别用我前面那个粗略比值糊弄。从那以后我每次跑新数据集都强制先拿示例数据走一遍全流程确认环境和参数没问题再上真实数据。这个习惯帮我省了不止一个通宵。希望帮到你。本文还有配套的精品资源点击获取