基于Python+MNE的偏侧化差异波(N2pc/Pd/CDA )成分分析方法

之前公众号推送了以下两篇关于偏侧化差异波成分的计算：[1] 基于Matlab、EEGLab和ERPlab的偏侧化差异波(N2pc/Pd/CDA )成分分析方法、[2] 代码篇 || 基于Matlab、EEGLab和ERPlab的偏侧化差异波成分分析。但上述内容主要是基于Matlab, 近期随着Python、MNE使用的增多，如何使用Python+MNE对EEG数据进行偏侧化差异波分析，也逐渐引起研究者的关注。

在本推文中，我们主要基于 Dirk van Moorselaar 等人开发的DvM 工具箱或软件包，进行偏侧化成分分析中核心过程的实现的讲解。

第一部分 前期准备

1 分析环境部署

Python、MNE等分析环境的部署查阅：https://mne.tools/stable/index.html

DvM程序包的下载和部署查阅：https://github.com/dvanmoorselaar/DvM

2 示例数据下载

DvM示例数据下载地址：https://osf.io/j86zp/

第二部分 数据分析核心过程

在核心过程主要分为逐步（Step by Step）和DvM 工具箱对特定的ERP成分进行分析这两个方面的内容，下面我们将依次介绍。

本教程将介绍如何从EEG数据的预处理到创建ERP。要运行此脚本，首先需要确保下载了示例数据。为了能够导入数据，请确保相应的文件均放在DvM工具箱的tutorial文件夹的相应位置，具体如下：

eeg data ('.fif'): .../tutorials/processed/

behavioral data ('.pickle'): .../tutorials/beh/processed

2.1 数据导入

数据保存后，我们首先导入所需的模块：

# import python packages
import sys
import numpy as np
import pandas as pd
import seaborn as sns

import mne.filter

import matplotlib.pyplot as plt

# import DvM toolbox
sys.path.append('/Users/dvm/DvM') # 将路径改为 DvM toolbox 存放的位置
from support.FolderStructure import *
from eeg_analyses.ERP import *

第一步，我们首先导入1个被试的示例数据[1]。数据来源于Dirk van Moorselaar 等人（2021）年发表的文章，所有被试的原始数据都可以在这里（ https://osf.io/n35xa/）下载。在进行数据分析时，我们可以采用如下代码导入数据：

# read in data 
sj = 'erp'
beh, eeg = FolderStructure().loadData(sj = sj)
eeg.baseline = None # reset baseling procedure used during preprocessing

执行代码后，会出现如下内容：

Reading /Users/dvm/DvM/tutorials/processed/subject-erp_all-epo.fif ...
Found the data of interest:
t = -1500.00 ... 1099.61 ms
0 CTF compensation matrices available

Not setting metadata

Not setting metadata

1771 matching events found

Applying baseline correction (mode: mean)

0 projection items activated

2.2 逐步进行ERP差异波成分分析

Condition specific ERPs（step by step）

现在我们已经导入了数据，这样我们实现了建立ERP数据的第一步。在使用DvM计算N2pc之前，我们需要对一些EEG数据预处理的核心步骤进行简要回顾。

2.2.1 选取感兴趣的数据

• 目前，我们的分析只纳入无空间偏向的Block数据（见试验设计）

• 由于这是偏侧化的成分，需要确保我们只选择目标位于水平中线，分心位于垂直中线或没有的试次

# mask data of interest on the basis of behavioral parameters
mask = (beh.target_loc.isin([2,6])) & (beh.dist_loc.isin(['None', '0','4'])) & (beh.spatial_bias == 'no bias')
beh_ = beh[mask]
eeg_ = eeg[mask]

2.2.2 定义偏侧化

在前面的步骤中，我们确保只选择分心物不存在（label='None'）或出现在垂直中线（label='0'或'4'）的试验。而且我们只选择了偏侧目标。也就是说，我们只选择了目标位于右侧（label=6）或左侧（label=2）的试次。接下来，我们需要选择数据，以便我们可以创建同侧（即，与显示的刺激属于同侧视野的电极）和对侧波形（即，与显示的刺激属于相反的对侧视野的电极）。

• 作为感兴趣的电极对，我们选择P7/P8和PO7/PO8（即N2pc分析的典型电极点）

# get electrodes of interest
left_elec, right_elec = ['P7','PO7'], ['P8','PO8']
left_elec_idx = [eeg.ch_names.index(elec) for elec in left_elec]
right_elec_idx = [eeg.ch_names.index(elec) for elec in right_elec]

# split trials with targets left and right from fixation
left_mask = beh_.target_loc == 2
right_mask = beh_.target_loc == 6

# based on left and right trials create contra and ipsi waveforms
contra = np.vstack((eeg_._data[left_mask][:, right_elec_idx], eeg_._data[right_mask][:, left_elec_idx]))
ipsi = np.vstack((eeg_._data[left_mask][:, left_elec_idx], eeg_._data[right_mask][:, right_elec_idx]))

2.2.3 对试次和电极点进行平均

ERP分析的基本假设是，通过对多个试次进行平均，随机的大脑活动将被平均，因此相关的波形仍然存在。在上述细胞中产生的对侧和同侧波形仍然包含所有个体试次。因此，为了创建ERP，我们需要对试次进行平均（即，数组中的第一个维度；epochs X electrodes X times）。此外，由于我们对感兴趣的电极的响应感兴趣，所以我们需要在电极之间平均（即阵列中的第二维度）

# average across trials and electrodes
contra = contra.mean(axis = (0,1))
ipsi = ipsi.mean(axis = (0,1))

2.2.4 基线校正

然而，在我们检查产生的波形之前，我们需要应用基线校正。由于脑电是一种时间分辨信号，它往往含有与实验问题无关的时间漂移。各种内部和外部来源可能会导致时间漂移，这可能会随着时间和电极的变化而变化。为了减少这种漂移的影响，人们习惯于进行所谓的基线校正。基本上，这包括在基线期间（即，在我们的案例中，搜索显示的开始）使用EEG活动来校正刺激后间隔期间的活动（即，在案例中视觉搜索本身）。尽管有各种各样的方法用于基线校正，ERP研究中最常用的方法是从基线和刺激后间隔的每个时间点减去基线周期的平均值。换句话说，在一个时间间隔内计算每个电极的平均电压值，然后从该信号的时间间隔中减去该平均值。

如这里（https://erpinfo.org/order-of-steps）所解释的，基线校正是一个线性操作。这意味着您可以在对所有试验求平均值之前或之后执行基线校正，结果将是相同的。

在我们的例子中，搜索开始于时间点0，我们将在（-200~0 ms）的预刺激窗口内进行基线校正。

# get indices of baseline window
base_time = (-0.2,0)
base_idx = slice(*[np.argmin(abs(eeg.times - t)) + i for t,i in zip(base_time, (0,1))]) # control for slice stop at end - 1

# perform baseline correction
contra -= contra[base_idx].mean()
ipsi -= ipsi[base_idx].mean()

2.2.5 检查ERP

为了检查产生的ERP，我们检查了对侧和同侧波形以及从对侧波形中减去同侧波形所得的差异波的波形。在这一点上，一个常见的步骤是过滤掉产生的波形中存在的高频噪声。在下面的单元格中，我们将首先进行30 Hz 低通滤波，然后仅绘制相关的感兴趣窗口。请注意，数据实际上包含的时间点比我们实际感兴趣的要多，因此可以“安全地”过滤数据，而不必过于担心由滤波导致的边缘伪迹。

运行下面的代码，看看产生的波形是否如预期的那样？此外，还可以使用一些不同的设置来查看最终结果的变化。

# apply a 30-Hz low pass filter to the data 
contra = mne.filter.filter_data(contra, sfreq = 512, l_freq = None, h_freq = 30)
ipsi = mne.filter.filter_data(ipsi, sfreq = 512, l_freq = None, h_freq = 30)

执行代码后，会出现如下内容：

Setting up low-pass filter at 30 Hz

FIR filter parameters

---------------------

Designing a one-pass, zero-phase, non-causal lowpass filter:

- Windowed time-domain design (firwin) method

- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation

- Upper passband edge: 30.00 Hz

- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)

- Filter length: 227 samples (0.443 sec)

Setting up low-pass filter at 30 Hz

FIR filter parameters

---------------------

Designing a one-pass, zero-phase, non-causal lowpass filter:

- Windowed time-domain design (firwin) method

- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation

- Upper passband edge: 30.00 Hz

- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)

- Filter length: 227 samples (0.443 sec)

# set up plotting parameters
times_oi = (-0.2, 0.55)
times_oi_idx = slice(*[np.argmin(abs(eeg.times - t)) + i for t,i in zip(times_oi, (0,1))])

# plot contralateral and ipsilateral waveforms
ax = plt.subplot(1,2,1, xlabel = 'Time', ylabel = 'micro V', title = 'lateralized waveforms')
plt.plot(eeg.times[times_oi_idx], contra[times_oi_idx], label = 'contra', color = 'grey')
plt.plot(eeg.times[times_oi_idx], ipsi[times_oi_idx], label = 'ipsi', color = 'red')
plt.axvline(0, ls = '--', color = 'black')
plt.axhline(0, ls = '--', color = 'black')
plt.legend(loc = 'best')
sns.despine(offset = 10)

# plot difference waveform
ax = plt.subplot(1,2,2, xlabel = 'Time', ylabel = 'micro V', title = 'difference waveform')
diff = contra - ipsi
plt.plot(eeg.times[times_oi_idx], diff[times_oi_idx], color = 'black')
plt.axvline(0, ls = '--', color = 'black')
plt.axhline(0, ls = '--', color = 'black')
sns.despine(offset = 10)

plt.tight_layout()

2.3 使用DvM工具箱进行ERP差异波成分分析

Condition specific ERPs（DvM）

在上面，我们已经在四个单独的步骤中创建了一个差异波波形（1. 选择数据；2.确定偏侧化；3.平均；4.基线校正）。下面的环节将演示了使用DvM工具箱时，如何使用几行代码轻松地组合上述步骤。

2.3.1  实例化ERP对象

为了能够计算ERP，作为第一步，我们需要实例化ERP对象。除了数据（即eeg和行为数据）之外，我们还需要传递至少两个参数，即header和baseline。Header应该是一个字符串，该字符串引用beh数据帧中的列键，其中定义了有关感兴趣的刺激的信息（即，在示例中，具有目标位置的列）。基线是一个胞体，指定将在分析中使用的基线窗口。

# to make use of the ERP class, first instantiate a ERP object
erp = ERP(eeg, beh, header = 'target_loc', baseline = (-0.2,0))

2.3.2  电极翻转

在逐步ERP创建的步骤2中，我们通过指定左电极和右电极来确定侧化。在ERP类中，偏侧化是通过将数据从一个半脑区翻转到另一个半脑区来处理的，所有的试验中，感兴趣的刺激都是在注视的左侧呈现的。换句话说，当计算ERPs时，假设所有刺激都是从注视开始呈现的。

翻转数据之前，我们首先指定感兴趣的窗口，并应用低通滤波器（可选），以便在分段相关时间点之前进行滤波。

erp.selectERPData(time = [-0.2, 0.55], h_filter = 30)
erp.topoFlip(left = [2], header = 'target_loc')

执行代码后，会出现如下内容：

Setting up low-pass filter at 30 Hz

FIR filter parameters

---------------------

Designing a one-pass, zero-phase, non-causal lowpass filter:

- Windowed time-domain design (firwin) method

- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation

- Upper passband edge: 30.00 Hz

- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)

- Filter length: 227 samples (0.443 sec)

2.3.3  创建同侧和对侧波

我们可以使用ipsiContra函数创建ERP。此功能有多种选项，可在下面进行检查。目前，我们忽略了其中的几个选项，只尝试模拟我们逐步过程中计算的波形。

# inspect ipsiContra docstring
erp.ipsiContra?

请记住，我们以前使用的相关设置如下。

• 仅包括来自无空间偏移块的数据

• 仅包括侧向目标（label=2（左）和 6（右））

• 只包括在中线使用分心器或不使用分心器的试验

# calculate ERP
erp.ipsiContra(sj = sj, left = [2], right = [6], l_elec = ['P7','PO7'], 
               r_elec = ['P8','PO8'], midline = {'dist_loc': ['0','4','None']},
               conditions = ['no bias'], cnd_header = 'spatial_bias',
               erp_name = 'target_N2pc', )

假设所有的刺激都是正确的，则进行基线校正（模式：平均值）。

2.3.4  检查结果

此时，您可能会惊讶地发现没有明确的输出。这是因为，对于ipsiContra函数的每次调用，都将使用函数调用中指定的相关信息创建一个输出文件。下面的代码显示了如何读入这些数据，这些数据存储为mne evoked object。显示由mne实现的所有打印选项超出了本教程的范围。在这里，我们重点讨论如何使用Object对象绘制差异波形（并将结果与2.2 逐步分析的结果进行比较）。

# read in evoked object
erp = mne.Evoked(FolderStructure().FolderTracker(['erp','target_loc'], 'sj_erp-target_N2pc-no bias-ave.fif'))

   执行代码后，会出现如下内容：

Found the data of interest:        t =    -199.22 ...     550.78 ms (0.11 * 20 + 0.12 * 40 + 0.05 * 21 + 0.06 * 23 + 0.07 * 41 + 0.05 * 43 + 0.14 * 120 + 0.13 * 140 + 0.07 * 121 + 0.07 * 123 + 0.07 * 141 + 0.07 * 143)        0 CTF compensation matrices available        nave = 346 - aspect type = 100 No projector specified for this dataset. Please consider the method self.add_proj.

# specify ipsi and contralateral waveforms
contra_elec, ipsi_elec = ['P7','PO7'], ['P8','PO8']
contra_elec_idx = [eeg.ch_names.index(elec) for elec in contra_elec]
ipsi_elec_idx = [eeg.ch_names.index(elec) for elec in ipsi_elec]

# extract ipsi and contralateral data
contra = erp._data[contra_elec_idx].mean(axis = 0)
ipsi = erp._data[ipsi_elec_idx].mean(axis = 0)

# plot contralateral and ipsilateral waveforms
ax = plt.subplot(1,2,1, xlabel = 'Time', ylabel = 'micro V', title = 'lateralized waveforms')
plt.plot(erp.times, contra, label = 'contra', color = 'grey')
plt.plot(erp.times, ipsi, label = 'ipsi', color = 'red')
plt.axvline(0, ls = '--', color = 'black')
plt.axhline(0, ls = '--', color = 'black')
plt.legend(loc = 'best')
sns.despine(offset = 10)

# plot difference waveform
ax = plt.subplot(1,2,2, xlabel = 'Time', ylabel = 'micro V', title = 'difference waveform')
diff = contra - ipsi
plt.plot(erp.times, diff, color = 'black')
plt.axvline(0, ls = '--', color = 'black')
plt.axhline(0, ls = '--', color = 'black')
sns.despine(offset = 10)

plt.tight_layout()

2.3.5  其他

ERP class还提供了几个额外的功能。例如，下面的代码片段显示了如何分别针对空间和无空间偏向Block计算存在干扰物的ERP，仅使用正确的反应和反应时间分割（即快速反应和慢速反应）。

# read in data 
beh, eeg = FolderStructure().loadData(sj = sj)
eeg.baseline = None # reset baseling procedure used during preprocessing

# instantiate a ERP object
erp = ERP(eeg, beh, header = 'dist_loc', baseline = (-0.2,0))

# flip topography
erp.selectERPData(time = [-0.2, 0.55], h_filter = 30, excl_factor = dict(correct = [0]))
erp.topoFlip(left = ['2'], header = 'dist_loc')

# calculate ERP
erp.ipsiContra(sj = sj, left = ['2'], right = ['6'], l_elec = ['P7','PO7'], 
               r_elec = ['P8','PO8'], midline = {'target_loc': [0,4]},
               conditions = ['no bias','bias'], cnd_header = 'spatial_bias',
               erp_name = 'dist_tuned', RT_split = True)

执行代码后，会出现如下内容：

Reading /Users/dvm/DvM/tutorials/processed/subject-erp_all-epo.fif ...

Found the data of interest:
t = -1500.00 ... 1099.61 ms
0 CTF compensation matrices available
Not setting metadata
Not setting metadata
1771 matching events found
Applying baseline correction (mode: mean)
0 projection items activated
Dropped 177 epochs: 5, 6, 11, 14, 15, 20, 29, 32, 33, 34, 41, 43, 44, 50, 53, 55, 57, 60, 62, 64, 66, 68, 70, 77, 78, 79, 80, 81, 99, 102, 116, 126, 128, 132, 134, 139, 142, 149, 156, 159, 161, 173, 175, 181, 182, 190, 192, 195, 199, 201, 204, 215, 234, 244, 245, 246, 247, 263, 275, 280, 282, 295, 316, 321, 324, 336, 340, 350, 351, 353, 360, 365, 378, 383, 384, 398, 403, 414, 436, 454, 458, 469, 470, 483, 484, 497, 502, 503, 508, 513, 537, 539, 552, 561, 562, 621, 628, 632, 651, 654, 663, 669, 671, 698, 701, 711, 714, 720, 727, 743, 745, 746, 765, 787, 798, 800, 808, 812, 813, 814, 815, 821, 824, 829, 833, 837, 850, 851, 854, 855, 857, 858, 859, 862, 863, 866, 868, 869, 873, 874, 875, 877, 879, 880, 881, 882, 884, 885, 888, 890, 891, 894, 895, 896, 905, 906, 932, 935, 1004, 1005, 1027, 1028, 1076, 1078, 1234, 1255, 1393, 1408, 1414, 1456, 1466, 1571, 1607, 1623, 1649, 1669, 1704
Dropped 177 trials after specifying excl_factor
NOTE DROPPING IS DONE IN PLACE. PLEASE REREAD DATA IF THAT CONDITION IS NECESSARY AGAIN
Setting up low-pass filter at 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 227 samples (0.443 sec)

It is assumed as if all stimuli are presented right
Applying baseline correction (mode: mean)
Applying baseline correction (mode: mean)
Applying baseline correction (mode: mean)
Applying baseline correction (mode: mean)
Applying baseline correction (mode: mean)
Applying baseline correction (mode: mean)

在使用工具箱时，请记得引用：van Moorselaar, D., Daneshtalab, N., & Slagter, H. A. (2021). Neural mechanisms underlying distractor inhibition on the basis of feature and/or spatial expectations. Cortex, 137, 232-250.

为了方便读者，我们已经将文中提到的核心内容上传至百度云盘，公众号对话框中回复“DvM”即可获取下载链接。或直接直接使用如下链接获取内容。

链接: https://pan.baidu.com/s/13FYb4YwUSwy0s3ePOlsD1w 

提取码: hgja 

注：如遇链接失效，请留言。以便小编及时处理。同时也欢迎大家交流使用心得，以便让更多的小伙伴能不在为学习工具而感到苦恼，尽早解放双手，提升科研效率。

==流浪心球 精品推荐==

ERP数据分析入门指导课程：

ERP数据分析高阶进击：

ERP数据质量检测新指标（SME）网络研讨会笔记

EEG/ERP研究中避坑指南：

科研利器