超详细TMS-EEG数据处理教程（上）

文章来源于微信公众号（茗创科技），欢迎有兴趣的朋友搜索关注。

本教程示例了被试收缩或放松对侧手时，如何处理对初级运动皮层（M1）进行经颅磁刺激（TMS）的脑电图（EEG）。去除TMS伪影最好的办法是要先进行预处理。当TMS伪影去除步骤完成后，就可以继续进行EEG分析了。在本教程中，所要解决的研究问题是手部的预收缩是否会影响脑电图中的TMS诱发电位(TEP)。

EEG信号分析需要干净的数据。但这并不简单，因为同时应用TMS和EEG的难点是TMS线圈电磁场会产生较大的伪影。所以在联合TMS- EEG实验中，数据中不可避免地会出现一些TMS伪影。可以使用适当的设备（高采样率(如> 5 kHz) EEG放大器、为TMS应用而设计的EEG电极等）和事后离线分析中减少伪影，这样有助于数据的进一步分析。

数据信息

测量：使用两个32通道与TMS兼容的BrainAmp直流放大器(BrainProducts)连接到61通道与TMS兼容的EEG脑电帽(EasyCap)记录数据。采样率为5kHz，截止频率为1kHz，位分辨率为0.1微伏。按下图查看等距61通道的排列图。

TMS：使用连接到MagPro X100（Magventure）刺激器的C-B60蝶形线圈，将单个脉冲施加到最大刺激器输出的60%。总共施加了300个脉冲，脉冲间隔为3秒（时间抖动20%）。

线圈的位置是通过电机热点来确定的。使用磁共振（MRI）导航定位保持刺激位置不变。本教程为了凸显伪影，调整了线圈的方向和倾斜度，以诱导产生较强的伪影。

事件markers：在每次TMS脉冲开始时放置一个事件来表示条件。“S 1”表示“放松”条件，“S 3”表示“收缩”条件。

**内存问题**

TMS-EEG研究由于高采样率和较长的记录时间，通常获得的数据都是数GB（很大）。这个时候，你需要注意自己的操作系统是否合适，建议在64位Windows操作系统上运行本教程，运行64位版本的MATLAB，且至少有8GB RAM（运行内存）。关于电脑问题，如果大家需要高性能工作站，请联系我们工程师17373158786（微信同号）~

脉冲和记录

在开始数据分析之前，先从两个方面来考虑实验：1）TMS协议和2）大脑状态的实验操纵（即被试的任务，或缺少任务）。

本教程在编写时考虑了基于试次的分析，所以不能直接应用于连续数据。此外，本教程使用单脉冲颅磁刺激协议期间记录的EEG数据。虽然本教程是以单脉冲研究进行编写的，但其中大部分也适用于多脉冲数据。

伪影类型

下图显示了应用颅磁刺激时通道17的EEG（可以参见以上channel排列图中channel 17的位置）。

接下来，先给大家呈现一些数据中可能存在TMS伪影。虽然根据脑电图和TMS设备的特点和实验设计，并非所有类型的伪影都能在你的数据中观察到，但是理解这些类型的伪影也是非常重要的（有备无患~）。

（1）脉冲伪影：脉冲期间的数据可以认为是丢失的。

（2）铃声/阶跃响应伪影：取决于放大器的范围，一旦电势落在放大器的范围内，可以观察到一个突出的滤波器“振铃”，持续时间大约为7ms，这是由于TMS脉冲的高梯度引起的阶跃响应。下图中，红色的信号反映了脉冲和铃声/阶跃响应的组合。

（3）颅肌伪影：经颅磁刺激脉冲可能导致颅(头皮)肌肉抽搐。这些抽搐不能与由于运动皮层刺激而引起的反应混淆，而纯粹是由于头皮肌肉刺激而引起的抽搐。通常会持续10ms左右。

（4）再充电伪影：根据TMS机器，有可能在数据中观察到机器电容器充电时产生的一个峰值。有些刺激器允许设置给电容器充电的确切时间。在这个数据集中，使用MagPro X100（Magventure）刺激器，在刺激开始后，将充电延迟设置为500ms，会看到如下图的伪影。

（5）衰减伪影：在某些通道中存在类似指数衰减的伪影。这个伪影是很难进行预测的。这是由于TMS在电极引线中引起电流的磁电感应、电极-电解质-界面极化、电极移动和头/颈部/面部肌肉抽搐之间的相互作用而产生的。最差情况下，此伪影可以持续1s，但更常见的是持续50-150ms。

预处理步骤

处理TMS-EEG数据的步骤如下：（欢迎文末留言评论，在公众号后台回复“TMS-EEG处理脚本”，脚本搬运工竭诚为您服务~）

（1）使用ft_definetrial创建一个trial-structure；

（2）用ft_artifact_tms表示TMS脉冲的开始；

（3）使用ft_preprocessing、ft_timelockanalysis和ft_databrowser直观地确定有哪些伪影；

（4）使用ft_rejectartifact和ft_artifact_tms从trial-structure中排除铃声/阶跃响应和再充电伪影；

（5）使用ft_preprocessing去除之前被拒绝的伪影；

（6）使用ft_componentanalysis和ft_rejectcomponent进行独立成分分析，去除指数衰减和颅肌伪影。在这个步骤中，与其他伪影相关的独立成分（例如线路噪声，眨眼/眼跳）也可以被移除；

（7）使用ft_redefinetrial重新创建trial-structure;

（8）用ft_interpolatenan在先前被铃声/阶跃响应和再充电伪影占用的段中插入空白；

（9）使用ft_preprocessing进行进一步的处理，如基线校正、去趋势和滤波。

数据可视化检查

首先使用ft_databrowser浏览试次，ft_databrowser比较方便，可以直接从磁盘或内存中浏览数据。使用ft_preprocessing从磁盘中读取感兴趣的片段（即试次）。为此，首先需要创建一个试次矩阵，它指定磁盘上数据的哪些部分将表示为试次。这个矩阵有三列（或更多）和有多少试次就有多少行。使用ft_definetrial创建试次矩阵。

输出的cfg变量包含cfg.trl中的trial-structure。在本教程的后面，我们还需要将这个trial-structure复制到另一个MATLAB变量中。

从ft_definetrial获得的cfg结构包含了足够的信息以供ft_preprocessing将我们的数据从磁盘读取到trials中。由于读取数据需要很长时间（5-10分钟），因此可以先找到处理过的数据。比如：

如果是从磁盘上的原始数据文件中读取试次，请使用以下方法：

如果是从磁盘上的数据文件中读取试次，请使用下面的代码来保存处理过的数据结构以备后续使用。

现在，使用ft_databrowser直观地检查数据。在刺激前进行基线校正。

因为与实际的脑电图信号相比，经颅磁刺激脉冲的振幅是巨大的。所以使用“垂直”和“水平”按钮旁边的“+和-”（用蓝色方框标记部分），可以在两个轴上调整缩放。尝试使用试次旁边的箭头“<>”按钮进行浏览。

请注意，如果调整水平缩放，试次按钮将更改为段按钮，表明正在浏览一个试次中的段。通过点击“channel”按钮（用红色方框标记部分）来选择想要绘制的channel。如果想知道哪个通道代表图中的哪条线，单击“identify”按钮并单击图中的一条线。相应的channel将显示在线的上方。浏览这些试次，会发现有很多噪音。可以通过单击和拖动试次中的任何区域并单击选择来标记想要拒绝的数据段。也可以稍后拒绝特定段，删除整个试次，或用nans替换它。现在，跳过此步，尝试用独立成分分析（ICA）去除生理伪影。

为了检查TMS伪影，创建数据锁时平均值，使得在锁时平均值中更容易发现伪影。

为了腾出系统内存，可用以下语句清除data_tms_raw结构。

想要绘制的平均数据在data_tms_avg结构中显示。看下图：

时间用time字段表示，振幅存储在avg中。当检查这个avg时，可以看到它的大小是61x10000。所以，在这种情况下，你应该已经猜到，行表示通道，列表示每个通道中的时间点。dimord显示了数据字段中的维度代表什么。这里的dimord是chan_time。

现在在单独的窗口中绘制所有channel的数据。可以使用以下脚本关闭先前所有的图形窗口。

需要花点时间看看所有channel，看是否能够确定存在的伪影。

这里将仔细观察靠近刺激点的channel 17。

如果稍微调整一下绘图范围，我们可以很容易地将铃声/阶跃响应(~0-0.0045)与颅肌(~0.0045 - 0.015)伪影区分开来。不过，你也可以手动使用缩放按钮或以下代码：

在这个channel中，我们可以发现铃声/阶跃反应、颅肌、指数衰减和再充电伪影。下面的代码能够突出显示这个channel中的伪影。

伪影排除

脚本ft_rejectartifact可以调整试次结构，以排除包含伪影的段。首先必须知道要排除哪些段，可以使用ft_artifact_tms脚本。这个可以用来检测数据中的TMS脉冲。

现在有了铃声/阶跃响应伪影和再充电伪影的配置结构。ft_artifact_tms函数创建一个Nx2的矩阵（表示N个伪影）。由于创建了两个结构，那么将两个结构中的伪影合并到一个结构中，这样就可以直接一步就实现拒绝。

ft_rejectartifact可以以几种方式对包含伪影的试次进行操作：首先，可以选择拒绝整个试次，但在这里不合适，因为每个试次都包含TMS伪影。其次，可以选择用nans替换伪影数据。第三，使用在本教程中采用的方法，通过将试次分割成更小的部分来删除伪影，并保留其他部分。然后，用前面定义和保存的试次矩阵（trl）重建试次。

接下来，加载剔除伪影的数据段。

如果在上一步中已经下载了data_tms_segments（可以直接加载），则跳过下面的代码。否则，也可以使用以下代码读取。

可以看到去除铃声/阶跃响应和再充电伪影的试次段，试次数增目加了。

使用ft_databrowser，既可以浏览分段数据，也可以浏览原始数据。

可以很清晰地看到伪影段已经被标记出来了。

ICA独立成分分析

TMS刺激诱发EEG的独立成分分析后剔除肌肉伪影。这里将使用手动拒绝伪影的方法。将数据分解为独立成分，并拒绝伪影成分，同时注意不要剔除非伪影数据。根据数据大小，运行ICA可能需要很长的时间，同时运行ICA也需要相当大的内存（这里大概需要4GB的额外系统内存）。建议在这里下载另存：comp_tms.mat。然后可以用下面的代码进行加载。

**内存问题**

如果在跑ICA时出现内存不足问题，其实可以先对数据进行降采样率。注意，在进行降采样时，数据将大约是目标采样率的一半进行低通FIR滤波（例如，如果采样率降到1000Hz，数据将以500Hz进行低通滤波）。使用下面的代码进行降采样率。

对降采样率的数据运行ICA之后，重新加载原始数据，并对原始数据应用空间分解矩阵，以避免试次的边缘伪影。

运行完ICA之后，会得到一个类似于data_tms_segments结构的数据。

现在要看一下成分的时间进程，以确定哪些成分需要拒绝。还可以查看成分的地形图，同时查看伪影的空间分布。所有的TMS-EEG伪影在TMS脉冲开始时都是锁时的，可以通过查看成分的锁时均值来简化视觉检查。

现在，可以以浏览channel数据的方式浏览平均数据。与channel数据一样，你可以使用ft_databrowser浏览，也可以使用MATLAB内置的绘图函数浏览。不过要注意的是，如果使用ft_databrowser浏览平均数据（即从ft_timelockanalysis输出的数据），那么只能通过通道浏览，因为没有试次。

由于ICA原则上是一个空间滤波器，可以检查每个成分是如何在空间上加载原始channel数据的。为此，这里将使用ft_topoplotIC。

运行完ICA后，使用ft_databrowser或MATLAB的绘图函数，以及ft_topoplotIC，应该能够找到一个或两个衰减伪影或颅肌伪影成分。我们来看看成分41和56可以明显地看出衰减伪影。这些成分的地形图与刺激位置重叠。我们还可以看到，几乎所有其他成分都包含或多或少的颅肌伪影。因此，可以得出这样的结论，在该数据中，通过ICA无法完全去除颅肌伪影。

将想要删除的伪影成分列出来：

衰减&颅肌：41，56

眼跳：7

眨眼：33

线噪声：31

再充电/肌肉：1，25，52，37，49，50

然后从数据中删除这些成分。使用ft_rejectcomponent，可以删除成分并将数据转换回为channel表示。

现在可以对成分进行删除。

相应成分现在已经被删除，数据回到了以原来的channel表示。

现在我们可以查看一下数据的当前状态。

插值

我们可以看到channel 17上还是存在肌肉伪影。为了防止这个伪影干扰我们进一步的数据处理，还需要对这些肌肉伪影进行剪切和插值。通过删除某些段而造成的空白先由nans填补。然后，对所有包含nans的段进行插值，包含肌肉伪影的时间段也将被nans替代，一同进行插值。

我们需要重建一个原来的试次结构。也就是说，这里需要重建一个在本教程最开始时创建的trl矩阵。

重新创建好原来的试次结构并将所有数据块放回原来的位置。铃声和再充电伪影造成的数据空白段已经用NaN，即非数字值来替代。

从这两种数据结构的trial字段可以看到，同样有300个试次。在这个数据结构中，每个试次数据（通道x时间）位于试次字段的单元格中。现在要对每个试次含有肌肉伪影的段用nans进行替代。通过对数据的观察就会发现，大约在刺激开始后用nans替代15ms就足够了。

接下来就可以开始插值了。函数ft_interpolatenan通过试次、通道以及使用MATLAB内置的interp1函数插值包含nans的段。

现在将原始数据与处理后的数据进行比较，首先先加载原始数据的锁时average。

使用以下代码，比较原始数据的平均TEP与处理后数据的平均TEP。TEP（经颅磁诱发电位）

到这一步，我们已经完成了对TMS伪影的预处理啦~~

接下来的滤波、去均值、去趋势和降采样等分析步骤，我们紧接着会出一期“TMS-EEG数据处理教程（下）”，欢迎大家点赞、关注、收藏（不迷路哦~）。

超详细TMS-EEG数据处理教程（上）相关推荐

超详细的gnuplot使用教程【2】
超详细的gnuplot使用教程 1.gnuplot参数介绍及演示 1.1首先来解释一下会用到的各类参数以及其解释 1.2 画图实际测试 1.3 其它参数介绍(约定范围.坐标轴设定) 1.3.1 约束画 ...
超详细的 Wireshark 使用教程
原文超详细的 Wireshark 使用教程一.wireshark是什么? wireshark是非常流行的网络封包分析软件,简称小鲨鱼,功能十分强大.可以截取各种网络封包,显示网络封包的详细信息. ...
【02】2022.11最新超详细Vuforia图片识别教程
[02]2022.11最新超详细Vuforia图片识别教程文章目录 [02]2022.11最新超详细Vuforia图片识别教程 1.Vuforia环境搭建 2.License Key获取及注册 3. ...
Anaconda超详细下载安装配置教程(Windows)
Anaconda最新超详细下载安装配置教程(Windows) 命令总结写在最前面 1.查看conda版本: conda --version 2.进入python交互模式: python 3.退出pyt ...
【排序】什么都能排的C语言qsort排序详解【超详细的宝藏级别教程】深度理解qsort排序
[排序]什么都能排的C语言qsort排序详解[超详细的宝藏级别教程]深度理解qsort排序作者: @小小Programmer 这是我的主页:@小小Programmer 在食用这篇博客之前,博主在这里 ...
超详细的MySQL入门教程（四）
MySQL:简单的增删改查查询数据基本语法介绍打印任意值查询表中全部数据查询表中部分字段限定条件查询例1:查询编号值小于指定值的记录例2:查询地址不等于某值的记录例3:查询一级地址等 ...
【数据的存储】浮点数在内存中的存储详解【超详细的保姆级别教程，让面试官心服口服】手撕浮点数存储使用方式
[数据的存储]浮点数在内存中的存储详解[超详细的保姆级别教程,让面试官对你心服口服]手撕浮点数存储使用方式作者: @小小Programmer 这是我的主页:@小小Programmer 在食用这篇博客 ...
华为C8500S 超详细线刷刷机教程
华为C8500S 超详细线刷刷机教程,菜鸟也可以刷机手机, 华为C8500S, 教程华为C8500S详细线刷刷机教程,开始: 鉴于很多网友是新手,都在询问刷机过程,怎样刷机,下面我就做个详细教程, ...
超详细的cmake入门教程【转载】
这篇文章主要介绍了超详细的cmake入门教程,需要的朋友可以参考下源出处超详细的cmake入门教程什么是cmake 在 linux 平台下使用 CMake 生成 Makefile 并编译的流程 ...
DBeaver安装与使用教程（超详细安装与使用教程）
文章预览: DBeaver安装与使用教程(超详细安装与使用教程) 一.DBeaver安装教程 ①下载地址 ②图文安装教程二.DBeaver使用教程 ①mysql数据库为例 1>填写数据库信息 ...

超详细TMS-EEG数据处理教程（上）

超详细TMS-EEG数据处理教程（上）相关推荐

最新文章

热门文章