ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

基于YASA的Python多导睡眠图自动分析与睡眠分期实践

基于YASA的Python多导睡眠图自动分析与睡眠分期实践 简介YASA 是一个面向睡眠研究者与脑电分析开发者的开源 Python 工具箱专用于多导睡眠图PSG数据的自动睡眠分期、睡眠纺锤波/慢波/快速眼动事件检测、伪影剔除以及频谱分析与催眠图统计。资源包共 254 个文件、约 115.78MB包含核心 Python 源码、Jupyter Notebook 示例、HTML 文档页面、PNG 示意图、配置文件及少量测试数据其中 py 文件对应算法实现ipynb 可交互演示html 与 png 方便离线查阅文档和结果可视化整体结构规范适合按需检索与二次开发。当前已有 803 人学习下载适合具备 NumPy、Pandas、MNE 基础并希望用命令行方式处理睡眠脑电数据的中高级 Python 用户。包内不仅提供完整算法实现还附带了文档站点与 notebook 演示能帮助快速上手纺锤波检测、睡眠分期和频谱分析等流程并据此扩展自己的研究管线。1. 多导睡眠图分析为什么需要YASA这类Python包多导睡眠图PSG一晚上会记录十几个通道的脑电和生理信号光是一次完整的睡眠分期就需要标注几百个30秒片段。手动做既慢又容易不一致遇到连续几天记录人力成本几乎不可接受。YASAYet Another Spindle Algorithm就是把这件事打包好的Python工具箱。它覆盖自动睡眠分期、纺锤波/慢波/REM检测、伪影剔除、频谱分析和催眠图统计底层依赖NumPy、Pandas和MNE设计上以Jupyter Lab为主要交互界面。如果你已经能用Python读入EDF数据但不想从零实现特征工程和检测器YASA可以省掉大量样板代码。对于睡眠研究人员、生物医学工程师和算法工程师这是套上手就能开始跑流程的工具。2. 环境准备Python、MNE与YASA的安装和EDF数据加载2.1 为什么YASA选择依赖NumPy、Pandas和MNEYASA没有自己实现一套IO系统而是选择站在MNE的肩上是由于MNE已经解决了EEG数据读取问题。MNE的Raw对象自带通道类型和采样率信息能让YASA以统一形式接收EDF、EDF、BrainVision、BDF等格式数据。NumPy和Pandas则分别负责数值计算与结果组织事件检测返回的DataFrame、分期返回的Series都遵循Pandas接口这样后续筛选、可视化、统计都可以直接用熟悉的Pandas方法与现有分析流程无缝对接。这带来的直接收益是只要数据是MNE可读的YASA就能处理。你不需要为每个通道的物理意义重新编码也不需要维护多个文件读取分支。另外YASA的算法核心依赖SciPy的滤波与Hilbert变换、SciKit-Learn的机器学习模型这些都是活跃维护的库远比自己造轮子稳定。作为开源Python包YASA的源码完全公开你可以在pip和源码仓库看到所有检测函数的具体实现。当内置参数无法满足要求时直接修改源码后重装开发版也是常见做法。这种透明度是闭源商业软件不具备的。2.2 安装YASA的几种方式及常见Python环境问题我一般会先创建一个干净的虚拟环境再开始安装。如果你用的是Anaconda可以这样conda create -n sleep python3.9 conda activate sleep pip install --upgrade pip pip install yasa这段命令创建了名为sleep的conda环境Python版本锁定在3.9然后激活环境升级pip后安装yasa。YASA会带上NumPy、Pandas、MNE、SciPy和SciKit-Learn这些核心依赖。为什么强调Python 3.9而不是最新版本因为部分EDF读取工具或C扩展在发布不久的新Python版本上可能还没有对应的二进制wheel包3.9目前是最稳妥的兼容选择。如果你正在用vscode python环境配置开发环境装完包后记得把Jupyter kernel指向当前虚拟环境否则会看到ModuleNotFoundError。如果你不使用conda只靠系统Python可以这样建venvpython3 -m venv sleep-env source sleep-env/bin/activate # macOS/Linux pip install -U yasa jupyter在Windows上激活命令换成sleep-env\Scripts\activate。如果是linux系统安装python后直接跑这条命令建议先确认pip是20.3以上旧版pip可能无法解析yasa的某些依赖。我习惯在安装完成后验证一下版本python -c import yasa; print(yasa.__version__)如果输出版本号说明环境已经就绪。要注意YASA对Python版本有最低要求我建议在Python 3.8以上低于3.9虽然也能跑但个别功能在旧版本上会有弃用提示。提示不要在你的Anaconda base环境里直接pip install yasa。base环境往往装有大量科学计算包依赖冲突很难排查单独的环境能让你快速回滚。2.3 读取EDF/PSG数据从MNE到YASAYASA本身不提供文件读取入口而是接收MNE的Raw对象或numpy数组。所以第一步是用MNE读取睡眠脑电数据import mne raw mne.io.read_raw_edf(subj01_psg.edf, preloadTrue, verbose0) raw.set_channel_types({C4: eeg, EOG1: eog, EMG1: emg}) print(raw.info)这里preloadTrue表示一次性把数据读进内存。对一个8小时、512Hz采样率的PSG记录来说内存占用通常在200MB以下开发机都能接受。verbose0关闭启动日志。set_channel_types是关键一步YASA需要靠通道类型找到EEG、EOG和EMG如果类型不匹配后续功能会被跳过或报错。raw.info会打印通道名、采样率和通道类型这是排查问题时的第一站。拿到Raw对象后可以提取为数组和采样率import numpy as np data raw.get_data() sf raw.info[sfreq] print(data.shape, sf)data是二维数组第一维是通道第二维是时间序列。YASA的事件检测函数可以直接接收data和sf不需要额外包装。这样设计的好处是你能在传给YASA之前先做预处理比如带通滤波、重采样或剔除坏通道预处理不会成为检测流程的瓶颈。如果EDF文件里只有少量通道实际记录与声明不一致MNE会直接报错。常见的做法是先用raw.rename_channels统一通道名再用set_channel_types设类型。值得一提的是如果你导入的是.fif或者.bdf读取API分别是read_raw_fif和read_raw_bdf其他流程完全一致。如果记录设备的采样率是1000Hz而YASA内部参考频率为128Hz或256Hz可以通过raw.resample(256)降低采样率。对睡眠脑电而言256Hz完全足够捕捉纺锤波频段信息还能减少检测时间。我通常在加载数据之后就调用resample而不是把原始1000Hz数据直接传给YASA。在继续睡眠分期之前把数据加载这一步走通后面所有功能才有一个稳定的起点。3. 自动睡眠分期与事件检测纺锤波、慢波和REM3.1 YASA的分期原理基于特征的分类而不是端到端网络自动睡眠分期是YASA最常用的功能之一。与直接喂原始波形给深度网络的方式不同sleep_staging先从30秒的脑电、眼电和肌电片段中提取若干特征再使用预训练的随机森林分类器输出每一期的类别。特征包括各频带功率、功率谱斜率、眼动和肌电的统计量。这种设计的直接好处是模型决策可以追溯从返回的特征表中能看到哪一段是因为theta功率升高而被分为N1也能看到REM期肌电幅度下降的证据而不是面对一个无法解释的黑盒。这个选择也划出了YASA的适用边界。预训练模型是在公共睡眠数据集上训练的如果你的样本来自婴儿、老年病患或神经退行性疾病人群脑电形态偏移会让分期准确率下降。正确的做法是用你的标注数据对随机森林做迁移学习或者至少用k折交叉验证来估算误差再决定是否直接采用默认权重。若只想做快速原型验证默认模型已经能给出合理的粗粒度分期。3.2 用yasa.sleep_staging得到睡眠阶段序列分期函数的调用非常直接staging yasa.sleep_staging( raw, eeg_nameC4, eog_nameEOG1, emg_nameEMG1, hypnohypno.txt )eeg_name、eog_name和emg_name分别指定脑电、眼电、肌电通道hypno是可选参数传入专家标注的催眠图文件路径后可以立即得到评估指标。返回的staging对象调用predict()即可获得预测结果sleep_stages staging.predict() print(sleep_stages.value_counts())sleep_stages是一个Pandas Series索引是30秒对应的时间戳值按AASM标准编码0Wake1N12N23N34REM。value_counts()能快速查看各期占比这个分布如果明显偏向某一期通常说明通道或标记有问题。不同版本YASA对编码定义是固定的但官方文档中提醒过如果你使用RK标准需要自行把W/N1/N2/N3/REM映射到对应数字。为了确认分期是否存在明显边界跳变我一般会画一个阶梯图import pandas as pd import matplotlib.pyplot as plt stage_df sleep_stages.reset_index() stage_df.columns [time, stage] plt.figure(figsize(12, 3)) plt.step(stage_df[time], stage_df[stage], wherepost) plt.yticks([0, 1, 2, 3, 4], [Wake, N1, N2, N3, REM]) plt.show()这段代码把时间作为横轴、分期整数作为纵轴画阶梯图。wherepost使得每个30秒区间保持恒定值符合睡眠分期是分段常数的性质。绘图时如果发现Wake/N1交替频繁到几秒一跳通常不是生理现象而是EMG通道噪声太大。这种情况下回到MNE端对EMG做带通滤波会明显改善结果。3.3 睡眠纺锤波检测的参数边界纺锤波是N2期最具标志性的脑电事件频率范围通常在11-16Hz持续0.3秒以上。YASA的spindles_detect允许你精确控制定义sp yasa.spindles_detect( data, sf, ch_names[C3, C4], freq_range(11, 16), duration(0.3, 2.0), relative_power_threshold0.2, couplingTrue )ch_names指定参与检测的通道列表多通道时逐通道检测结果会带通道列。freq_range控制带通滤波上下限部分研究者用9-16Hz来捕获慢纺锤波这时可以放宽到(9, 16)。duration是事件持续时间的上下限下限太低容易把高频瞬态噪声当作纺锤波上限太高则可能把两段靠近的纺锤波合并成一次。relative_power_threshold是检测核心参数它比较事件片段与周围背景信号的功率默认0.2表示纺锤波频带功率比背景高20%。阈值越高事件数越少但特异性更高。couplingTrue会额外计算纺锤波与慢波之间的相位幅度耦合。不过需要注意coupling仅在单通道检测时有效多通道传入时这个参数会被忽略。如果只关心事件数量建议保持默认False以节省计算时间。下面是一个简化参数参照表当你想针对不同人群调整检测灵敏度和特异性时可以先从这里入手参数默认值作用与边界条件freq_range(11,16)纺锤波频带迷你纺锤可用16-24Hz但误检会增加duration(0.3,2.0)事件持续时间窗下限小于0.3s会捕获高频伪影relative_power_threshold0.2相对背景功率阈值越高事件越少、越保守couplingFalse是否计算相位幅度耦合仅单通道可用调用sp.summary()即可获取每个事件的起止时间、峰值频率、幅度和持续时间等列。我倾向于先在一小段已知有明确纺锤波的N2数据上跑参数确认检测框与人工标注吻合后再全量运行。这样比盲目批量执行要高效得多。3.4 慢波和REM检测的注意点慢波检测使用yasa.sw_detect参数结构和纺锤波类似但是频带低得多通常在0.5-4.5Hz持续时间0.3到1.5秒。慢波对低频基线漂移非常敏感如果低估信号中的直流偏移漂移会被当作慢波检测出来。我的做法是先做一次1Hz以下的高通滤波再调用sw_detect这样得到的慢波密度和幅度都更可靠。如果多通道同时检测慢波YASA会输出每个通道的结果你可能还需要跨通道聚合事件避免同一个慢波被相邻通道重复计数。快速眼球运动检测使用yasa.rem_detect需要至少一导EOG和一导EEG。EOG里天然混入脑电和肌电活动YASA通过回归扣除EEG成分来提取纯眼动信号。如果EEG和EOG位置很靠近回归后的残差可能包含大量脑电活动导致REM事件过多。这时可以预先对EOG做20-30Hz的带通滤波把肌电和眼动频段分开。REM检测结果一般会以DataFrame形式返回包含事件起始终止时间和幅度。在解读REM事件密度时最好同时观察同时间段的睡眠分期因为N1期也偶尔会出现较小的眼动容易混淆。4. 伪影剔除与频谱分析让手动预处理变成自动化流程4.1 伪影剔除yasa.art_detect的阈值逻辑长程PSG记录中的眨眼、眼球运动、动作伪影会严重干扰事件检测。YASA提供art_detect它按固定窗口遍历数据计算窗口内的峰峰值、绝对幅度、梯度等指标根据统计阈值给出是否标记为伪影的布尔值。art yasa.art_detect( data, sf, window1.0, methodboth, threshold5, filter_bandauto )window是滑动窗口长度单位秒。method可选amplitude、gradient或bothboth会同时考虑幅度和梯度两个判别条件。threshold的单位是绝对偏差倍数默认5表示窗口统计值超过中位数5倍绝对偏差时认为包含伪影。filter_bandauto会让YASA根据信号类型自动选择滤波范围这一步主要是避免肌电高频噪声的干扰。返回的art是一个布尔数组长度等于窗口数量。可以用NumPy把这些窗口转成时间段再在事件检测结果里做排除。下面是把伪影窗口索引打印出来的常见做法import numpy as np bad_windows np.where(art)[0] print(f伪影窗口数: {len(bad_windows)}) bad_times bad_windows * 1.0 # 窗口起始秒拿到时间点后最简单的方式是用MNE的raw.annotations集合后续所有下游步骤都能直接使用这些注释标记。注意art_detect面对的是连续窗口而不是具体事件段所以在做事件级分析时你还需要把事件起止时间与伪影时间段做重叠判断。使用Pandas的IntervalIndex.overlaps是最顺手的方式events sp.summary()[[Start, Duration]].copy() events[End] events[Start] pd.to_timedelta(events[Duration], units) clean ~pd.IntervalIndex.from_arrays( events[Start], events[End] ).overlaps( pd.IntervalIndex.from_arrays( pd.to_timedelta(bad_times, units), pd.to_timedelta(bad_times 1.0, units) ) )这里第一行把纺锤波事件的起始与持续时间转为时间戳格式并计算结束时间第二行用IntervalIndex构造事件区间和伪影区间overlaps返回是否存在重叠~取反后就得到干净事件掩码。有了clean掩码后续统计可以直接过滤事件不需要重跑检测。这个方法同样适用于慢波和REM事件。4.2 频谱分析与1/f斜率数据清洗完成后最常用的评价睡眠深度的指标是频带功率。YASA的bandpower函数一次性返回多个频带比手动写Welch谱要方便bp yasa.bandpower( data, sf, ch_names[C3, C4], bands[ (delta, 0.5, 4), (theta, 4, 8), (alpha, 8, 12), (sigma, 12, 16), (beta, 16, 30) ] )bands中的每个元组包含频带名称、下限与上限。返回的DataFrame每一行对应通道和频段的组合列包含绝对功率与相对功率。相对功率在个体间横向比较时更稳定因为绝对功率会受皮肤阻抗和放大器增益影响。假如你想研究睡眠剥夺后的delta功率变化相对功率比绝对功率更稳定因为个体间头皮电阻差异会直接乘以绝对功率。YASA返回的DataFrame中Rel列就是相对功率可以直接用。如果你关心sigma频带与纺锤波时间的一致性可以把各通道的bandpower保存下来再与前面的事件检测结果做时间对齐。1/f斜率也是近年睡眠研究中常被提到的一个指标。YASA通过yasa.compute_slope计算功率谱高段默认20-40Hz上的对数斜率。这段区间没有明显的振荡峰能更干净地反映背景活动。需要注意的是1/f估计对数据长度比较敏感至少要用5分钟以上的连续数据短片段回归出来的斜率方差大得难以解释。4.3 相位幅度耦合PAC怎么算慢波相位与纺锤波幅度之间的耦合是记忆巩固研究的核心指标之一。YASA提供yasa.pac做通用计算pac yasa.pac( data[0, :], sf, f_phase(0.5, 4.5), f_amp(11, 16), methodozkurt )f_phase和f_amp分别是相位频率范围和幅度频率范围methodozkurt是推荐方法之一它对噪声的稳健性较好。method参数可选项还包括plv、corr等。plv对相位同步更敏感但受包络幅度影响小ozkurt更适合调制幅度变化明显的数据。如果你不确定选哪个可以先跑一个小片段对比两者方差选择分离度高的。计算时YASA先对信号进行带通滤波再用Hilbert变换提取相位和瞬时包络最后在给定频段上计算耦合强度。因为Hilbert变换在信号首尾处有边界效应建议传入的数据去掉首尾0.5秒或者提前裁剪否则耦合值会被边界伪影抬高。如果你需要分阶段统计可以先把数据切成30秒片段循环调用pac再按分期结果分组平均。5. 用催眠图统计和批量汇总技巧把分析链路收敛到一张表5.1 从hypnogram到睡眠统计yasa.hypno_analysis睡眠分期完成后我们通常需要得到总睡眠时间、睡眠效率、入睡后清醒时间WASO和阶段转换次数等统计量。YASA的hypno_analysis直接接受分期序列hypno yasa.hypno_analysis(sleep_stages, sf_hypno1/30) print(hypno)sleep_stages就是第3章得到的按30秒排列的阶段序列sf_hypno设为1/30表示每个样本代表30秒。返回的DataFrame包含总记录时长、总睡眠时间、睡眠效率、各期占比、觉醒次数等关键字段。尤其值得注意的是YASA对“睡眠开始”的定义是从入睡后的第一个连续睡眠段算起不包括入睡前的清醒片段。如果你习惯WASO从熄灯到起床的全过程计算需要手动截断关灯前和起床后的记录或者用start_time和end_time参数框定分析区间。5.2 进阶技巧用Pandas对齐多个通道的事件结果最后一个技巧用于批量处理多导记录。当你有多个通道的纺锤波或慢波检测结果时YASA的summary()返回一个DataFrame里面已经包含了通道名、事件开始时间、持续时间和峰值频率。为了按睡眠阶段汇总我们需要把每个事件对齐到它所在的30秒分期区间summary sp.summary() stage_at_event sleep_stages.reindex( summary[Start].dt.round(30s), methodffill ) summary[Stage] stage_at_event.values这里dt.round(30s)把事件的开始时间舍入到最近的30秒边界然后通过reindex找到对应的分期。methodffill用前一个有效阶段填补舍入后可能出现的空值这个近似在分析整夜数据时基本无损。要更精确的话可以用searchsorted把事件起始时间映射到分期序列索引上但多导数据量不大时round更直观。有了带阶段的summary后即可快速生成汇总表grouped summary.groupby([Channel, Stage]).agg( count(Start, count), duration_ms(Duration, mean), frequency_hz(Frequency, mean) )这个分组聚合会按通道和睡眠阶段生成事件数、平均持续时间、平均峰值频率。例如C4通道在N2期内出现200次纺锤波平均持续时间0.8秒峰值频率13.1Hz这些信息可以直接作为论文或报告中的event table。对于多条记录只需要把每个受试者的summary拼接成长表all_summary pd.concat([s1.summary(), s2.summary()], ignore_indexTrue)之后再按受试者ID分组执行同样的聚合就能把十几条记录浓缩成一张表。这个技巧同样适用于慢波和REM检测因为它们的summary结构一致。使用这套流程整夜从加载数据到输出统计表只需要几步而且每一步的结果都是DataFrame随时可以保存为CSV或画图。本文还有配套的精品资源点击获取
返回列表