脑电ROI连接分析实战:wPLI、置换检验与交互HTML报告 简介这份资源是面向脑电EEG研究者的ROI连接分析工具箱以EEGLAB插件形式提供解决在感兴趣脑区之间进行连通性统计的核心需求。工具箱内含Stefan Haufe开发、Arnaud Delorme完成EEGLAB接口封装与Fieldtrip集成的代码支持TRGC与CS两类连接分析方法可加载表面或体积图集并与Fieldtrip、Brainstorm及NFT头部模型兼容适合已掌握EEGLAB基本预处理流程、希望进一步做脑区级功能连接的中高级用户。压缩包共371个文件、约46.31MB主体为210个.m函数文件与117个.html帮助文档另有少量.png示意图、.mat数据、.txt说明及面向macOS、Linux、Windows的mex编译文件构成可离线查阅的完整工具链与示例文档。目前已有655人学习下载读者可从中获得ROI连接分析的完整函数库、方法演示页面与图集配置样例便于快速搭建分析流程、理解TRGC与CS的实现细节并排查接口调用问题。1. 从一批 64 导联静息态数据说起手上这批 64 导联静息态数据合作者只关心一件事左右额叶的 alpha 频段连接到底有没有差别。真做起来路径比想象中长——电极要按脑区聚合成 ROI每个频带要单独算连接矩阵还要过统计和多重比较最后得交出一份别人点开就能看懂的东西。脑电 ROI 连接分析工具包 roiconnect 干的就是把这条链路固定下来最终产出一个自带交互图的 HTML 报告而不是散落一地的 npy 和截图。它适合两类人一是做 EEG 连接分析、需要可复现结果的研究生和工程实现者二是要给合作方交付分析结论、又不想每次重写绘图脚本的人。前提是你得知道 ROI 怎么划、连接指标怎么选否则工具跑得再顺矩阵也是错的。2. ROI 划分与连接度量8 个脑区、5 个频带的参数怎么定这一步决定后面所有数字的含义。ROI 划错连接矩阵看起来再漂亮也是噪音的排列组合指标选错容积传导会把相邻电极连成一片结论直接翻车。2.1 电极到 ROI 的映射表ROI 聚合的本质是空间降维把几十个电极压到个位数脑区换取更高的信噪比和更少的多重比较次数。代价是空间分辨率变糊两个相邻但功能不同的区域可能被并进同一个 ROI。所以映射表不能随手按字母分组常见做法是按 10-20/10-10 系统的解剖邻近性来划并且保证每个 ROI 至少有 3 个电极——单电极 ROI 等于没降维还容易受坏道影响。下面是一份 64 导联静息态的示例映射按左右半球各四个区域加中线来组织。你自己的导联表不同照着改电极名即可结构不用动。ROI 名称电极示例半球说明L_frontalFp1, AF3, F3, F7, FC3, FC5左额叶眼电污染重灾区R_frontalFp2, AF4, F4, F8, FC4, FC6右同上需重点查坏道L_centralFC1, C3, C1, CP3左运动相关注意肌电R_centralFC2, C4, C2, CP4右同上L_parietalCP5, P3, P7, PO3左alpha 源区附近R_parietalCP6, P4, P8, PO4右同上L_occipitalPO5, PO7, O1左枕区alpha 最明显R_occipitalPO6, PO8, O2右同上midlineFz, Cz, Pz中单独成组避免污染左右对比映射写成 Python 字典最省事后面所有步骤都复用它。写完之后先跑一个校验脚本别等算完连接才发现某个电极名拼错了。# roi_map.py ROI_MAP { L_frontal: [Fp1, AF3, F3, F7, FC3, FC5], R_frontal: [Fp2, AF4, F4, F8, FC4, FC6], L_central: [FC1, C3, C1, CP3], R_central: [FC2, C4, C2, CP4], L_parietal: [CP5, P3, P7, PO3], R_parietal: [CP6, P4, P8, PO4], L_occipital: [PO5, PO7, O1], R_occipital: [PO6, PO8, O2], midline: [Fz, Cz, Pz], } def check_roi_map(roi_map, ch_names): assigned [ch for chs in roi_map.values() for ch in chs] missing sorted({ch for ch in assigned if ch not in ch_names}) dup sorted({ch for ch in assigned if assigned.count(ch) 1}) small {roi: len(chs) for roi, chs in roi_map.items() if len(chs) 3} print(缺失电极:, missing) print(重复分配:, dup) print(电极数不足3的ROI:, small) return not missing and not dup and not smallassigned把所有 ROI 的电极摊平missing抓的是映射表里有、数据里没有的名字dup抓的是一个电极被分进两个 ROI——这种情况不报错但会让该电极的方差被重复计入。small是经验阈值少于 3 个电极的 ROI 建议合并或者直接剔除。校验通过再往下走能省掉后面大量返工。2.2 连接度量选型PLV、PLI、wPLI、相干差在哪这四个指标都基于相位或相位加幅值但对容积传导和参考电极的敏感度完全不同。头皮 EEG 的相位同步里混着大量零相位滞后成分主要来源是同一个源经容积传导到多个电极或者平均参考把公共信号注入每个通道。PLV 和相干对这种零滞后同步毫无抵抗力会把一个源误判成一片高度连接的网络。指标对零相位滞后的鲁棒性偏差特征计算成本适用场景相干 coherence差随电极距离下降慢低探索性分析、源空间PLV差样本量小时向上偏低短窗任务态、已知无容积传导PLI较好对噪声敏感中静息态长数据wPLI较好加权后方差更小中高默认推荐正交化 AEC好需要足够样本中源空间或 CSD 之后我一般的默认选择是 wPLI 配合 CSD电流源密度预处理两条腿走路CSD 在电极层面压低容积传导wPLI 在指标层面忽略零相位成分。如果只能选一个优先 CSD——预处理层面的修正对下游所有指标都生效而换指标只是换个角度偏。选型顺序可以按这个清单走先确认数据是静息态还是任务态静息态且时长足够用 wPLI 或 PLI任务态试次少、窗短PLV 的偏差会更明显宁可先做 CSD 再用 wPLI如果一定要用相干至少把结果和 wPLI 对比一遍两种指标结论差异大的边直接标为不可信。别迷信某一个指标的绝对数值连接值本身没有单位跨研究比较意义有限同一批数据内部的相对差异才是你能用的信息。2.3 频带与窗长的耦合为什么 1 秒的 epoch 算不好 delta频率分辨率约等于 1/TT 是分析窗长。一个 1 秒的窗只能分辨 1 Hz 以上的频率间隔而 delta 频带1–4 Hz的下限本身就是 1 Hz窗内连三个完整周期都凑不齐相位估计的方差会大到让连接矩阵在不同试次之间几乎不相关。这不是算法问题是物理限制。频带频率范围建议最短窗≥3 周期静息态常用窗delta1–4 Hz3 s2 s接受偏差theta4–8 Hz0.75 s1–2 salpha8–13 Hz0.4 s1–2 sbeta13–30 Hz0.25 s1 sgamma30–45 Hz0.1 s0.5–1 s麻烦在于多个频带共用同一个 epoch。常见做法是统一用 2 秒窗、50% 重叠delta 的结果标注为「窗长受限」或者干脆把 delta 单独用 4 秒窗重算一遍。gamma 则是另一个极端头皮 EEG 的 30 Hz 以上成分很容易被额肌和颞肌的肌电污染算出来的高连接值往往是肌肉而不是大脑除非你有 EMG 通道或者独立成分分析把肌电去掉否则 gamma 的连接结果我一般只作参考、不进结论。还有一个容易忽略的参数是多窗法的时间带宽积。用 DPSS 多窗时窗数越多方差越小但频率平滑也越强alpha 和 theta 的边界会糊在一起。静息态我通常取 3–5 个窗任务态短数据取 2–3 个具体在连接函数的n_tapers或带宽参数里调。3. 从原始文件到 ROI 时间序列预处理流水线连接分析对预处理比对功率谱更敏感。功率谱里一个坏道只会抬高一个通道的功率连接矩阵里一个坏道会污染它所在 ROI 和其他所有 ROI 的边所以这一步的处理日志要留着出问题好回溯。3.1 输入格式与目录约定roiconnect 这类工具通常吃几种中间格式但原始数据来源五花八门。比较稳的做法是先用 MNE-Python 统一读进来、存成 FIF 或 npy再交给连接计算环节。不要在每个脚本里重复写读取逻辑格式差异会把调试时间耗光。格式读取方式注意点.edf / .bdfmne.io.read_raw_edf/read_raw_bdf通道名可能带前缀需手动对齐.set / .fdtmne.io.read_raw_eeglab事件码和采样率以 .set 为准.vhdrmne.io.read_raw_brainvision参考电极信息可能缺失.npynp.load只存数据采样率和通道名要另存目录按被试分每个被试一个文件夹里面放原始文件、预处理后的 FIF、ROI 时间序列 npy 和一份 JSON 处理日志。日志至少记录滤波参数、参考方式、坏道列表、剔除的试次数。后面组水平统计如果发现某个被试异常靠这份日志能快速定位是哪一步动过手。data/ sub-01/ raw.edf clean_raw.fif roi_ts.npy roi_names.json log.json sub-02/ ...3.2 滤波、重参考、坏道处理滤波下限不要设得太高。0.5–1 Hz 的高通能压掉漂移但也会削弱 delta做 delta 连接时我一般用 0.5 Hz 高通加 45 Hz 低通再对 50 Hz 工频做陷波。高通滤波器用 FIR相位失真小IIR 虽然陡但会引入相位偏移对相位类连接指标不友好。重参考是最大的分歧点。平均参考简单但会把公共信号注入每个通道制造零滞后同步双侧乳突参考保留了更多空间信息但乳突本身可能带肌电CSD 用球面样条做空间高通能有效压低容积传导代价是边缘电极的估计不稳。我的习惯是先做 CSD再看一眼结果是否合理如果边缘 ROI 的连接值异常低就退回乳突参考做对比。import numpy as np import mne from mne.preprocessing import compute_current_source_density def load_and_prepare(path, montage_namestandard_1020, l_freq0.5, h_freq45.0): raw mne.io.read_raw_edf(path, preloadTrue, verboseFalse) raw.pick_types(eegTrue, eogFalse, miscFalse, verboseFalse) # 工频陷波50 Hz 及其谐波 raw.notch_filter(freqs[50, 100], pickseeg, verboseFalse) raw.filter(l_freql_freq, h_freqh_freq, fir_designfirwin, verboseFalse) montage mne.channels.make_standard_montage(montage_name) raw.set_montage(montage, on_missingwarn, verboseFalse) return raw def detect_bad_channels(raw, z_thresh3.0): data raw.get_data(pickseeg) logvar np.log(np.var(data, axis1) 1e-12) z (logvar - logvar.mean()) / (logvar.std() 1e-12) bads [ch for ch, zz in zip(raw.ch_names, z) if abs(zz) z_thresh] raw.info[bads] bads return bads raw load_and_prepare(sub-01/raw.edf) bads detect_bad_channels(raw) raw.interpolate_bads(reset_badsFalse, verboseFalse) csd compute_current_source_density(raw, verboseFalse)l_freq和h_freq是滤波边界做 delta 时别把l_freq抬到 1 Hz 以上否则低频连接会被自己滤掉。z_thresh控制坏道判定的松紧3 是个起点数据干净可以放到 4噪声大就降到 2.5但降太低会把正常通道误杀。interpolate_bads用的是球面样条插值插值后的通道参与 ROI 平均是可以的但要在日志里记下插了哪几个后面如果某个 ROI 结果异常先查它是不是靠插值撑起来的。compute_current_source_density要求 montage 完整缺电极位置会报错。如果原始数据的通道名和标准蒙太奇对不上先做名称映射别硬跑。3.3 ROI 聚合与几个边界ROI 时间序列有两种取法先 CSD 再对电极平均或者先平均再 CSD。前者保留空间信息后者等于把 ROI 当成一个虚拟电极。我倾向前者因为 CSD 是线性空间滤波先做再平均在数学上更干净。如果工具包默认是后一种至少要保证同一个数据集里所有被试的处理顺序一致。def roi_timeseries(raw, roi_map): data raw.get_data(pickseeg) ch_idx {ch: i for i, ch in enumerate(raw.ch_names)} out, names [], [] for roi, chs in roi_map.items(): idx [ch_idx[c] for c in chs if c in ch_idx] if len(idx) 2: print(f跳过 {roi}有效电极不足) continue out.append(data[idx].mean(axis0)) names.append(roi) return np.array(out), names arr, names roi_timeseries(csd, ROI_MAP) np.save(sub-01/roi_ts.npy, arr)这里arr的形状是(n_roi, n_times)后面要转成 epoch 才能算连接。电极不足 2 个的 ROI 直接跳过而不是报错是为了让批量脚本能跑完所有被试、最后统一看日志不至于一个被试缺电极就中断整批。聚合完之后建议画一眼各 ROI 的时域波形和功率谱。如果某个 ROI 的波形明显比其他 ROI 平多半是它包含的电极里有坏道没处理干净如果所有 ROI 在 50 Hz 都有一个尖峰说明陷波没生效或者数据里还有未标记的工频段。4. 连接矩阵、置换检验与 HTML 报告一次跑完到这一步输入是(n_roi, n_times)的数组输出应该是一份能直接发出去的 HTML 文件中间的计算和统计都在脚本里完成不依赖手工导出。4.1 逐频带算 wPLI 矩阵先把 ROI 时间序列切成固定长度的 epoch再做连接估计。切窗时长按第 2 章的结论来静息态统一 2 秒、50% 重叠比较稳。连接函数在新版 MNE 里已经拆到 mne-connectivity 包导入路径变了这是第一个容易卡住的地方。import numpy as np import mne from mne_connectivity import spectral_connectivity_epochs BANDS {delta: (1, 4), theta: (4, 8), alpha: (8, 13), beta: (13, 30)} info mne.create_info(names, sfreqcsd.info[sfreq], ch_typeseeg) roi_raw mne.io.RawArray(arr, info, verboseFalse) epochs mne.make_fixed_length_epochs( roi_raw, duration2.0, overlap1.0, preloadTrue, verboseFalse) conn spectral_connectivity_epochs( epochs, methodwpli, modemultitaper, sfreqepochs.info[sfreq], fmin[b[0] for b in BANDS.values()], fmax[b[1] for b in BANDS.values()], faverageTrue, tmin0.2, tmax1.8, n_jobs1, verboseFalse) mat conn.get_data(outputdense) # (n_roi, n_roi, n_bands) mat (mat mat.transpose(1, 0, 2)) / 2 # 强制对称压掉数值误差 for b in range(mat.shape[2]): np.fill_diagonal(mat[:, :, b], np.nan) # 对角线是无意义的自连接tmin0.2和tmax1.8是去掉窗两端的边缘效应所有连接函数都会在窗边界产生估计偏差留出 10% 左右比较稳妥。faverageTrue表示对频带内所有频点取平均得到每个频带一个标量矩阵如果想保留频率细节设成 False 会多出一维。modemultitaper比fourier方差小但频率平滑更强频带很窄时比如 alpha 只有 5 Hz 宽要注意别把相邻频带糊进来。如果你的环境里from mne_connectivity import ...报错说明装的是老版本 MNE连接函数还在mne.connectivity下面函数名是spectral_connectivity_epochs参数基本一致。两个版本的结果不保证完全一致同一批数据别中途换环境。4.2 置换检验与 FDR 校正组水平比较的对象是边不是节点。每个频带有n_roi * (n_roi - 1) / 2条独立边9 个 ROI 就是 36 条4 个频带 144 次比较不做校正必然假阳性泛滥。常见做法是置换检验拿经验零分布再用 Benjamini-Hochberg 控 FDR。import numpy as np from scipy.stats import false_discovery_control def paired_perm_test(a, b, n_perm5000, seed42): a, b: (n_subjects, n_edges)配对设计 rng np.random.default_rng(seed) d a - b obs d.mean(axis0) count np.zeros_like(obs) for _ in range(n_perm): sign rng.choice([-1.0, 1.0], size(d.shape[0], 1)) count (np.abs((d * sign).mean(axis0)) np.abs(obs)) p (count 1) / (n_perm 1) return obs, p obs, p paired_perm_test(cond_a, cond_b, n_perm5000) p_adj false_discovery_control(p, methodbh) sig (p_adj 0.05) (np.abs(obs) 0.05)sign的形状是(n_subjects, 1)广播到所有边上保证同一个被试的所有边在同一轮里被同时翻转——这是配对设计的核心打乱被试标签而不是打乱单条边的值。n_perm取 5000 时最小可分辨 p 值约 2e-4够用如果要做更严格的校正加到 10000。obs是均值差也可以换成 Cohens d 之类标准化效应量报告里更好解释。false_discovery_control在 SciPy 1.11 之后才有老环境用statsmodels.stats.multitest.multipletests(p, methodfdr_bh)结果一致。sig里额外加了效应量阈值是为了避免大样本下统计显著但差异微小的情况——连接值的绝对差小于 0.05 时我一般不当成有意义的发现。4.3 HTML 报告的结构与导出报告要解决的是「别人不用装 Python 也能看懂」。结构上分三块每个频带的连接矩阵热图、显著边的列表、被试级别的 QC 信息坏道数、有效 epoch 数。热图用 Plotly单文件 HTML 可以直接发邮件。import plotly.graph_objects as go def heatmap(mat, band_idx, roi_names, title): fig go.Figure(go.Heatmap( zmat[:, :, band_idx], xroi_names, yroi_names, colorscaleRdBu_r, zmid0, hovertemplate%{y} → %{x}br值%{z:.3f}extra/extra)) fig.update_layout(titletitle, width520, height500, margindict(l80, r20, t50, b80)) return fig figs [heatmap(mat, i, names, f{b} 频带 wPLI) for i, b in enumerate(BANDS)] html [htmlheadmeta charsetutf-8, titleROI 连接分析报告/title/headbody, h1ROI 连接分析报告/h1] for i, fig in enumerate(figs): # 只有第一张内嵌 plotly.js其余复用避免文件膨胀 html.append(fig.to_html(full_htmlFalse, include_plotlyjs(i 0))) html.append(/body/html) with open(report.html, w, encodingutf-8) as f: f.write(\n.join(html))include_plotlyjs是关键参数。设成cdn文件小但离线打开图表空白发给不联网的合作方就是灾难设成True会把 plotly.js 内嵌进文件单个约 3 MB所以只在第一张图嵌入、后面全部设False否则每张图都重复内嵌文件能涨到几十兆。还有一种做法是include_plotlyjsdirectory把 js 单独存成文件放在同目录适合报告页数多、需要拆分的场景。报告里最好再加一张表格列出经过 FDR 校正后显著的边、对应频带、效应量和校正后 p 值。热图看整体模式表格看具体结论两者配合才不容易被误读。5. 避坑排查五个让连接值失真的常见问题这一章是我自己踩过、也见过别人反复踩的几处。每一条都按「现象 → 原因 → 解决」写遇到对应症状可以直接对照。5.1 所有 ROI 两两连接都偏高且数值接近现象热图上几乎每一格都是暖色没有明显的区块结构左右半球、前后区域之间的差异很小。原因平均参考把公共信号注入了每个通道。每个电极都含有一份相同的参考信号任意两个电极之间就多了一条零相位滞后的通路wPLI 虽然对零相位不敏感但参考信号经过滤波和加窗之后会产生非零的相位扰动仍然会抬高基线。解决换成双侧乳突参考或者做 CSD然后重算。判断方法很简单把参考方式改一遍如果连接矩阵的整体水平明显下降、区块结构浮出来那原来的结果就是参考污染。报告里要注明用了哪种参考不同参考的结果不能混在一张图里比较。5.2 同侧相邻 ROI 的连接远高于远距离 ROI现象左额和左中央、左中央和左顶之间的连接值明显高于左额和右枕距离越近连接越强呈现清晰的随距离衰减。原因容积传导。头皮上的电场会扩散到邻近电极距离越近共享的信号成分越多相干和 PLV 对这种现象尤其敏感。解决先做 CSD再观察衰减是否减弱如果仍然明显改用 wPLI 或 PLI并在报告里明确说明未做源空间重建。想彻底解决只能做源定位但源空间的连接估计又引入逆问题的不确定性属于另一个层级的坑。实操上距离效应不可能完全消除只要它不主导组间差异就可以接受。5.3 delta 频带的连接矩阵在不同试次间几乎不相关现象把同一被试的 epoch 随机分两半分别算 delta 连接矩阵两个矩阵的相关只有 0.2 左右alpha 频带同样操作能到 0.8 以上。原因窗长不足。2 秒窗对 1 Hz 的下限只有两个周期相位估计的方差极大连接值基本是噪声。解决delta 单独用 4 秒窗重算或者把 delta 的下限抬到 2 Hz 并在方法里注明。另一个办法是增加 epoch 数量做平均但连接值不是简单平均能降方差的样本量带来的改善有限。如果研究问题必须覆盖 1 Hz就得接受 delta 结论的置信度低于其他频带别把它和 alpha 的结果放在同等权重上解读。5.4 某个被试的单个 ROI 与其他被试差异巨大现象组水平分析里某一个被试的某个 ROI 连接值比均值高出一大截去掉这个被试结论就变了。原因该 ROI 里有坏道没有插值或者插值后没有重新检查。ROI 平均是把几个电极的时序直接相加除以个数一个幅值异常大的坏道会把整个 ROI 的波形带偏连接值随之失真。解决回到第 3 章的坏道检测把阈值调紧一点重跑确认插值后的通道波形正常。批量处理时把每个被试的坏道数和插值通道写进日志组水平分析前先扫一遍日志坏道数明显偏多的被试单独标记。敏感性分析里做一次留一法看结论是否依赖单个被试这比事后争论要不要剔除更靠谱。5.5 HTML 报告发出去图表空白或加载极慢现象本地打开报告一切正常发给合作方后热图区域一片白或者浏览器卡住十几秒。原因include_plotlyjs设成了cdn对方环境访问不到外部资源或者反过来每张图都设了True文件体积过大。解决统一成「第一张内嵌、其余复用」的写法文件控制在 5 MB 以内。如果报告要归档另存一份纯静态图版本PNG 或 SVG作为附件避免几年后库版本变化导致老报告打不开。导出前用浏览器的无痕模式打开一次能模拟对方没有缓存的环境。6. 进阶验证用替代数据和分半信度给结果上保险连接分析的结论最怕两种情况一是把噪声当成了组间差异二是结果换一批数据就消失。前者靠置换检验和 FDR 挡了一部分后者得靠内部一致性来验证。我习惯在出报告之前加一步分半信度把每个被试的 epoch 随机分成两半各自算连接矩阵然后看两半矩阵在边层面上的相关。alpha 和 beta 一般能到 0.7 以上delta 低于 0.4 就要警惕。这个数字不写进最终结论但会决定我敢不敢把某个频带的结果放进报告正文。更严格一点的做法是相位随机化替代数据。对每个 ROI 的时间序列做傅里叶变换保留幅值谱、把相位随机打乱逆变换回时域再算一遍 wPLI。重复几百次得到一个零分布真实连接值只有明显超出这个分布的分位数才算脱离了「信号自身的频谱结构就能解释」的范围。import numpy as np def phase_randomize(x, rng): x: (n_times,)返回相位随机化后的序列 X np.fft.rfft(x) phase rng.uniform(0, 2 * np.pi, X.shape[0]) phase[0] 0.0 # 直流分量相位固定 if X.shape[0] % 2 0: phase[-1] 0.0 # 奈奎斯特分量同样固定 return np.fft.irfft(np.abs(X) * np.exp(1j * phase), nlen(x))这段代码只处理单通道批量做的时候对每个 ROI 独立随机化再用同一套连接函数重算循环 200–500 次。phase[0]和奈奎斯特分量必须固定否则逆变换出来的序列均值会漂、长度会对不上。替代数据法的计算量是原分析的几百倍一般只对关键频带和关键边做不要全量跑。最后说个习惯。我现在每次跑完连接分析都会强制走一遍「换参考重算一次、分半信度看一眼、报告用无痕模式打开一次」这三步哪怕时间紧也至少做第一步。有一次图省事跳过了参考对比结果整组结论建立在平均参考的伪同步上返工时所有统计重来。那之后这套检查就成了固定动作宁可晚交一天也不交一份自己不敢复现的报告。希望帮到你。本文还有配套的精品资源点击获取