脑电分析系列[MNE-Python-16]| 脑电数据的Epoching处理

2022-08-17 10:37:00 浏览数 (1)

代码语言:javascript复制
import os.path as op
import numpy as np
import mne
import matplotlib.pyplot as plt

在MNE中,epochs是指单个试验或一小段时间锁定的原始数据的集合。

首先,读入原始样本数据:

代码语言:javascript复制
data_path = mne.datasets.sample.data_path()
fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw.fif')
raw = mne.io.read_raw_fif(fname)
raw.set_eeg_reference('average', projection=True)  # 设置EEG 平均参考

要创建时间锁定的epochs,首先需要一组包含有关时间信息的事件。

这里使用刺激通道来定义事件。

代码语言:javascript复制
order = np.arange(raw.info['nchan'])
"""
交换两个通道的绘图顺序,以将触发通道显示为第十个通道。
"""
order[9] = 312
order[312] = 9  
raw.plot(n_channels=10, order=order, block=True)
plt.show()

注意底部的STI 014频道。触发通道用于将所有事件合并到单个通道中。

从上图上可以看到在整个记录中有几个振幅不同的脉冲。这些脉冲对应于在采集过程中呈现给受试者的不同刺激。

脉冲的值为1、2、3、4、5和32。要从原始数据创建事件列表,MNE中只需调用一个专门用于此目的的函数。

由于事件列表只是一个numpy数组,所以也可以手动创建一个。

如果是从外部源(如单独的事件文件)创建事件,则应注意将事件与原始数据正确对齐。

代码语言:javascript复制
events = mne.find_events(raw)
print('Found %s events, first five:' % len(events))
print(events[:5])

"""
绘制事件以了解范例
为图例指定颜色和event_id字典。
"""
event_id = {'Auditory/Left': 1, 'Auditory/Right': 2,
            'Visual/Left': 3, 'Visual/Right': 4,
            'smiley': 5, 'button': 32}
color = {1: 'green', 2: 'yellow', 3: 'red', 4: 'c', 5: 'black', 32: 'blue'}

mne.viz.plot_events(events, raw.info['sfreq'], raw.first_samp, color=color,
                    event_id=event_id)
plt.show()
代码语言:javascript复制
320 events found
Event IDs: [ 1  2  3  4  5 32]
Found 320 events, first five:
[[27977     0     2]
 [28345     0     3]
 [28771     0     1]
 [29219     0     4]
 [29652     0     2]]

如上面所示,事件列表包含三列。

第一列对应于样本编号,要将此转换为秒,可以将采样数除以使用的采样频率。

第二列是在转换时保留给触发器通道的旧值,目前没有使用。

第三列是触发ID(脉冲幅度)。

这里说明一下为什么这些样本看起来与绘制的数据不一致。

例如,第一个事件的样本编号为27977,应该转换为大约46.6秒(27977 / 600)。但是查看脉冲时,可以在3.6秒时看到第一个脉冲。这是因为Neuromag记录有一个first_samp属性,它表示系统启动和录制开始之间的偏移量。Neuromag记录数据的first_samp等于25800。这意味着使用raw.plot看到的第一个样本的样本号25800。

一般来说,在使用时不需要担心这个偏移量,因为它在MNE函数中已经被考虑进去了,不过最好要注意这一点。

为了确认一下,我们将事件与原始数据一起绘制。注意垂直线(事件)如何与STI 014上的脉冲很好地对齐。

代码语言:javascript复制
raw.plot(events=events, n_channels=10, order=order)
plt.show()

在本文中,我们只对触发器1、2、3和4感兴趣。这些触发器对应于听觉和视觉刺激。

这里的event_id可以是int、int列表或dict。使用dict可以将这些id分配给不同的类别。当使用int或列表时,这个信息就会丢失。

首先,我们为mne.Epochs构造函数定义一些参数,tmin和tmax指的是与事件相关的偏移量,并使用epoch来封装事件前200毫秒到事件后500毫秒的数据。

代码语言:javascript复制
tmin, tmax = -0.2, 0.5
event_id = {'Auditory/Left': 1, 'Auditory/Right': 2,
            'Visual/Left': 3, 'Visual/Right': 4}

经过上面步骤,构建了epochs所需的参数。为了得到一些有意义的结果,我们还希望将这些epochs作为基线。基线化计算基线期间的平均值并相应地调整数据。epochs构造函数默认使用从tmin到0.0秒的基线周期,作为元组的第一个元素的None指的是时间窗口的开始(本例中为-200 ms)

这里定义了阈值来去除噪声。阈值被定义为epoch时间窗口内的峰到峰的值。定义为T/m表示梯度计,T表示磁强计,V表示EEG和EOG电极。

代码语言:javascript复制
baseline = (None, 0.0)
reject = {'mag': 4e-12, 'eog': 200e-6}
epochs = mne.Epochs(raw, events=events, 
                    event_id=event_id, tmin=tmin,
                    tmax=tmax, baseline=baseline,
                    reject=reject,
                    picks=('meg', 'eog'))  # 这里只包含 MEG 和 EOG
代码语言:javascript复制
289 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
Created an SSP operator (subspace dimension = 3)
4 projection items activated

下面就绘制epochs查看结果。顶部的数字是指ID号,从图中可以看到145个事件中有128个通过了rejection process。

代码语言:javascript复制
epochs.plot(block=True)
plt.show()

通过绘制drop日志,来查看为什么剔除epoch(一般被伪影等污染的epoch数据需要被剔除)。

代码语言:javascript复制
epochs.plot_drop_log()
plt.show()

要获得诱发响应,只需执行epoch.average()。它默认只包含数据通道。为了便于举例,我们还使用pick来包含EOG通道。

代码语言:javascript复制
picks = mne.pick_types(epochs.info, meg=True, eog=True)
evoked_left = epochs['Auditory/Left'].average(picks=picks)
evoked_right = epochs['Auditory/Right'].average(picks=picks)

注意,这里使用了前斜杠('/')来分隔实验条件的各个因素。我们可以使用这些“标签”来选择例如所有左侧试验(包括视觉左侧和听觉右侧)。

代码语言:javascript复制
epochs_left = epochs['Left']

# ... or to select a very specific subset. This is the same as above:
evoked_left_again = epochs['Left/Auditory'].average(picks=picks

下面绘制诱发响应

代码语言:javascript复制
evoked_left.plot(time_unit='s')
evoked_right.plot(time_unit='s')
plt.show()

0 人点赞