CorTec系统神经科学计算 Python 语言(SyNCoPy):用于大规模电生理数据分析的 软件包

 一、引言 

在神经科学领域，脑电图（EEG）、脑磁图（MEG）、皮层脑电图（ECoG）和微电极记录等方法被用于测量源自大脑活动的电磁信号。这些技术的高时间分辨率使得能够在大频率范围内分析大脑活动，这对于理解系统神经科学中大脑区域的功能互联至关重要。该领域的研究人员通常希望识别与特定实验条件相关的大脑活动，例如向受试者呈现刺激的开始。因此，实验任务会重复多次，随后对产生的试次进行平均以减少噪声和方差。试次重复与现代实验设置相结合，使用越来越多的记录位点（通道）和高采样率，可能产生非常大（>10 GB）的数据集。使用这些数据集，通道之间的全对全连接计算等标准算法在内存有限的笔记本电脑或台式计算机上可能无法执行，需要工作站或高性能计算（HPC）系统，而这些系统使用起来可能很复杂。此外，最近人们对使用科学Python技术栈作为数据分析的开源环境产生了浓厚的兴趣。

在本研究中，我们推出了系统神经科学Python计算（SyNCoPy），这是一个用于大规模电生理数据分析的Python包，它将易于使用的、类似FieldTrip的应用程序编程接口（API）与HPC系统上分布式工作流的内置支持相结合。

 二、相关软件包 

用于神经电磁数据分析的科学软件包包括用于MATLAB的FieldTrip、EEGLAB、NUTMEG和Brainstorm，以及用于Python的MNE Python和Elephant。SPM软件还包括MEG/EEG分析功能。

FieldTrip是一个MATLAB工具箱，于2011年首次发布，此后一直在积极发展。其功能包括预处理、多元时间序列和连接性分析以及源定位。它带有数据浏览器、交互式数据可视化和广泛的文档。功能API由强大的主函数（例如ft_preprocessing、ft_freqanalysis和ft_connectivityanalysis）和许多较小的辅助函数组成。大多数函数可以用输入数据和配置结构作为输入参数调用，并返回包含配置副本的输出数据结构，作为应用于数据的操作历史记录，以及将分析重新应用于不同输入数据的方式。

EEGLAB自2004年以来一直在开发，是一个用于处理连续和事件相关EEG、MEG和其他电生理数据的交互式MATLAB工具箱。它包括图形用户界面（GUI）和API，并通过插件接口支持用户贡献的代码。功能包括交互式可视化、伪迹去除、独立成分分析（ICA）、时频分析和源建模。

Brainstorm软件包是用MATLAB和Java编写的，但可以作为独立应用程序运行，无需MATLAB许可证。它专注于复杂的GUI，并提供一些批处理功能。

Elephant是一个用于电生理数据分析的Python库，专注于电极尖峰序列数据和时序记录的通用分析功能。

MNE Python是一个支持数据预处理、源定位、统计分析和分布式脑区域之间功能连接估计的Python包。它建立在科学Python生态系统之上，有许多贡献者，并与使用Neuromag FIF文件格式的其他应用程序良好集成。MNE具有广泛的绘图功能和文档，包括公开可用的示例数据集和教程。它通过Python的joblib模块支持单台机器上多核的并行化，但目前不支持HPC系统的直接并行化。API是细粒度函数和直接定义在数据对象上的方法的组合。MNE专注于EEG和MEG数据以及局部场电位（LFP）的分析，支持伪迹去除、时/频分析和源建模。

我们开发SyNCoPy是为了补充MNE和Elephant的一些功能，并提供易于使用的、类似FieldTrip的API、对时间离散尖峰数据集的支持以及HPC系统的内置并行化。

 三、SyNCoPy架构 

上述软件解决方案已得到充分验证，并与SyNCoPy共享不同的功能。然而，它们都不是为处理超大数据集和在HPC系统上进行分布式计算而设计的。SyNCoPy通过支持试次并行核外计算的架构来支持这一用例。SyNCoPy的核心数据结构由元数据和多维数据数组组成，但数据数组默认不加载到内存中。相反，当请求计算时，数据以试次为单位从硬盘上存储的分层数据格式5（HDF5）容器中流式传输，结果以类似方式写回磁盘。元数据以JavaScript对象表示法（JSON）格式存储。这种方法允许对具有大量试次的超大数据集进行内存高效处理，并便于基于试次的并行化。并行化通过采用成熟的Dask库实现，每个Dask作业一次处理一个试次。在标准计算机上，如果有足够的内存，试次可以顺序或并行使用多个核心处理。在HPC或基于云的系统上，Dask调度器通常将计算作业分布在多个节点上以实现并行化。这意味着可以使用当今的HPC系统并行处理大量试次。

在笔记本电脑上启动SyNCoPy（左）以处理多试次数据集。当在Jupyter Notebook中执行高级SyNCoPy API函数时，基于NumPy和SciPy的SyNCoPy算法被包装在计算例程中，该例程通过Dask连接到高性能计算集群（或笔记本电脑上的本地集群），并自动将逐试次计算分发到可用资源。作业并行运行（中），每个工作进程一次处理一个作业，并将试次的结果写入磁盘上单个HDF5容器的适当槽位。当所有工作进程完成其分配的作业后，磁盘上的结果即完成，可以从笔记本电脑上的SyNCoPy会话访问（右）。然后可以使用基于matplotlib的SyNCoPy绘图API对结果进行可视化，导出为神经数据无边界（NWB）格式，或可以直接提取NumPy数组用于使用标准科学Python技术栈进行自定义后处理。

SyNCoPy的内部架构和在大型数据集上运行并行计算的推荐设置如图1所示。用户连接到远程JupyterHub实例，例如由机构高性能计算（HPC）集群提供的实例。在创建全局Dask客户端后，运行SyNCoPy分析将使用可用的计算资源。输入数据应驻留在可从集群访问的快速存储上，通常是文件服务器。当用户启动并行计算时，SyNCoPy自动检测并使用Dask集群，并将工作分发到HPC集群节点。节点将结果写入磁盘，SyNCoPy API函数返回的SyNCoPy数据结构指向磁盘上的数据。请注意，结果数据永远不会直接通过网络传输，也永远不会完全加载到内存中。对于后处理，API提供了与matplotlib（用于绘图）以及NumPy和pyNWB（用于数据导出）的接口。

图1. SyNCoPy架构和并行处理的典型设置。SyNCoPy计算函数（作为ComputationalRoutine运行）可以附加到任何正在运行的Dask客户端，因此可以利用Dask生态系统的全部灵活性，例如轻松部署到云资源。

SyNCoPy提供专门的数据结构和在其上实现并行核外计算的通用方法，即ComputationalRoutine。用户暴露的函数（高级SyNCoPy API，例如syncopy.connectivityanalysis）在内部评估用户指定的配置，然后使用ComputationalRoutine机制执行通常处理单个试次数据的代码。根据全局Python环境，ComputationalRoutine顺序或通过Dask并行执行每试次代码（另见esi-acme），以与并行化后端交互，例如在HPC集群上运行的Slurm作业调度器。SyNCoPy分析脚本与硬件环境无关，这意味着分析可以在单台机器（如笔记本电脑）上本地开发和运行，相同的代码以后可以部署在分布式计算资源上。

 四、功能概述 

SyNCoPy的当前功能可分为数据处理、预处理、时域锁定分析、频域分析和基于连接性的分析等广泛类别。

4.1 数据结构和数据处理

数据处理类别包括使用SyNCoPy内部数据格式加载和保存数据的函数，以及一些转换数据的函数，即导入数据并将其导出为其他文件格式。SyNCoPy的核心数据结构通常包含多维数据数组和元数据。在磁盘上，数据表示为HDF5文件，当数据加载到内存中时，它们作为NumPy数组可用。SyNCoPy不直接读取电生理记录系统生成的文件；它目前支持从NWB、HDF5或NumPy格式的文件导入数据。数据结构可分为连续数据的数据类型和离散数据的数据类型。AnalogData类通常用于存储原始电生理数据，即多通道、定期采样的模拟数据，包含一个或多个试次。如果数据源中没有可用的试次信息，用户通常会创建试次定义来定义试次。对于许多分析类型，应用延迟选择以确保数据时间锁定到特定事件（如刺激开始），这会产生TimeLockData实例。输出实数或复数频谱数据的算法将这些结果存储在SpectralData类的实例中，那些产生通道-通道交互信息（连接性度量）的算法返回CrossSpectralData类的实例。离散数据类SpikeData和EventData分别用于存储尖峰和事件。SpikeData类可以存储在外部尖峰排序软件（如SpyKING CIRCUS）中识别的尖峰，包括每个尖峰周围的原始波形。EventData类用于存储事件时间，通常与其他数据类结合使用。

所有数据类都可以从NumPy数组和特定数据类型的元数据（例如AnalogData实例的采样频率）初始化。为了在初始化期间也便于内存安全的数据处理，产生单试次NumPy数组的Python生成器可以直接输入到相应的SyNCoPy数据类构造函数中。为了提高与其他软件包的互操作性，提供了在MNE Python和SyNCoPy的数据结构之间转换的函数。我们还提供了以NWB格式保存和加载数据的函数。

4.2 预处理

SyNCoPy的预处理函数作用于AnalogData实例，支持信号的去趋势、归一化和滤波，包括低通、高通、带通和带阻滤波器。还支持时间序列数据的重采样和降采样。

4.3 时频分析

SyNCoPy提供对AnalogData类型输入进行频域分析和时频分析的函数。（多）锥形傅里叶变换（MTMFFT）算法使用单个锥形窗口或基于离散长椭球序列（DPSS）的多个锥形对时间序列数据进行频谱分析。有效频率平滑宽度可以直接以赫兹为单位用tapsmofrq参数控制，如同FieldTrip中一样。SyNCoPy中可用的单个锥形从SciPy的信号模块导入。得到的频谱可以使用FOOOF方法（拟合振荡和一过f）进行后处理。滑动窗口短时傅里叶变换也可用，以及Welch方法，该方法基于短时间、修改周期图的时域平均来估计功率谱。非正交连续小波变换和超小波在SyNCoPy中也可用于时频分析，后者可以以高时间和频率分辨率揭示快速瞬态振荡。

4.4 连接性分析

连接性分析模块揭示通道之间的功能连接性。它提供互谱密度估计（CSD）、相干性、成对相位一致性（PPC）、非参数格兰杰因果性和互相关算法。运行连接性分析需要SpectralData输入。如果传入AnalogData实例，将使用默认参数隐式运行MTFFT分析以获得SpectralData实例。

4.5 统计

SyNCoPy提供计算其数据类任意轴上的均值、中位数、标准差和方差的函数。可以计算SpectralData类型输入的试次间相干性。还实现了刀切法，可用于计算相干性或格兰杰因果性结果的置信区间。可以计算SpikeData实例的 peri刺激时间直方图（PSTH）。

4.6 绘图和实用函数

我们为各种SyNCoPy数据类型提供绘图函数，包括AnalogData、SpectralData和SpikeData。SyNCoPy绘图函数旨在让科学家在数据分析管道开发和项目演示期间快速轻松地概览其数据，但不提供可直接发表的图形。这些函数在内部使用matplotlib，用户可以根据需要对生成的图形进行后处理。

SyNCoPy中的synthdata模块包含创建合成数据集的实用函数，这对培训目的以及测试自定义算法和评估其性能很有用。除了标准过程（如白噪声或泊松散粒噪声以模拟尖峰数据）外，我们还提供红噪声（AR(1)过程）和相位扩散算法来模拟实验LFP信号。

支持（并已并行化）所有SyNCoPy数据类和NumPy数组的基本代数运算（如加法和乘法），允许灵活的合成数据构建和标准操作（如基线校正）。

 五、真实电生理数据集的分步分析管道示例 

以下，我们展示了一个分步分析管道示例，以演示如何使用SyNCoPy分析细胞外电生理数据。作为比较，相同的分析在MATLAB中使用FieldTrip进行。SyNCoPy版本和FieldTrip版本的源代码可在线获取。

图2展示了分析管道和用于处理样本脑信号的SyNCoPy函数。分析中使用的数据集是公开可用的，来自Allen Institute视觉编码—Neuropixels项目，此前已有描述。简而言之，通过高密度Neuropixels细胞外电生理探针同时记录LFP和尖峰活动。这些记录涵盖了小鼠大脑在处理视觉刺激时的各个区域。LFP数据使用Open Ephys记录，尖峰数据使用Kilosort提取。在实验期间，向小鼠呈现不同的视觉刺激。在本研究中，持续250毫秒的全视野闪光刺激被视为刺激时段，而刺激开始前250毫秒的时期被用作基线（图2B）。为了评估连接性分析，从一个样本会话中选择了两个视觉区域（区域A，在Allen数据集中称为VISl区域，对应于初级视觉区外侧部分；区域B，在Allen数据集中称为VISrl区域，对应于初级视觉区嘴侧部分）。在预处理数据以将数据与刺激开始对齐后，对上述时域和频域分析进行了测试。图2C显示了区域A不同试次和通道平均的LFP响应。它表明在视觉刺激呈现后具有短潜伏期的诱发反应。图2D和2E分别显示了样本神经元150次试次的尖峰序列的时域锁定光栅图和peri刺激时间直方图（PSTH）。随后，我们计算了区域A的LFP功率谱、区域A和区域B之间LFP的相干性谱、区域A和区域B之间LFP的成对相位一致性（PPC）谱，以及区域A和区域B之间LFP的非参数格兰杰因果性（GC）谱，作为四种常见的频域分析。

图2. SyNCoPy对电生理数据集示例的分析。（A）使用SyNCoPy函数处理电生理数据的示例分析管道。不同的管道产生图D-I中所示的图，如输入最终绘图程序的箭头上方所示。（B）在持续250毫秒的全视野闪光刺激呈现期间，从清醒小鼠的不同大脑区域记录LFP和尖峰活动。（C）区域A试次和通道平均的LFP响应，时间锁定到刺激开始。（D,E）样本神经元150次试次尖峰活动的时域锁定光栅图（D）和peri刺激时间直方图（E）。（F）区域A试次和通道平均的1-95赫兹频率范围内刺激与基线期间LFP功率比的频谱。黑线反映FieldTrip结果，红色阴影线对应SyNCoPy结果。（G-I）与F相同，但分别为区域A和区域B之间LFP的相干性（G）、区域A和区域B之间LFP的成对相位一致性（H）、区域A和区域B之间LFP的格兰杰因果性（I）。黑线是FieldTrip结果，红色阴影线是SyNCoPy结果。实线是前馈方向，虚线是反馈方向（I）。

为演示目的并与FieldTrip进行直接比较，在SyNCoPy和FieldTrip中都计算了这些分析。为此，首先对数据进行零填充。然后，基于MTMFFT方法并使用Hann窗口，计算每个试次和记录通道在刺激期间和基线期间的功率谱。MTMFFT对时间序列试次数据进行频域分析，要么使用单个锥形（如Hann），要么使用源自离散长椭球序列（DPSS）的多个锥形。对每个记录通道分别单独平均刺激期间和基线期间的功率谱，并计算刺激与基线功率的比值。随后，对通道平均功率比频谱（图2F）。

类似地，在对信号进行零填充后，测量所选区域对之间的相干性（图2G）、PPC（图2H）和格兰杰因果性（图2I）。对于功率、相干性和PPC，SyNCoPy和FieldTrip之间的结果基本相同，对于GC也非常相似（图2F-I）。

 六、内存基准测试 

6.1 峰值内存消耗—方法

我们研究了SyNCoPy在典型使用场景中几种算法的峰值内存消耗（PMC），即在HPC集群上进行并行处理期间。具体而言，使用了马克斯·普朗克学会马克斯·普朗克计算和数据设施（MPCDF）Raven集群的"small"队列。为了评估作为数据集大小函数的内存消耗，我们使用SyNCoPy的synthdata模块创建了递增大小的合成数据集，并使用SyNCoPy进行处理。我们评估了：（1）使用巴特沃斯+希尔伯特滤波器进行预处理，（2）MTMFFT，（3）MTMFFT f.t.算法（此处f.t.指定使用固定数量的锥形以便更好地比较，如下文更详细解释），（4）小波，以及（5）相干性。起始数据集大小为10个试次、每试次5,000个样本和50个通道，需要约10 MB空间。我们创建了脚本以使用不同数据集大小运行每个算法。每次调用SyNCoPy API函数后，调用Python垃圾收集器以确保测量有意义。在每次运行期间，使用Python的memory_profiler包监测PMC。PMC是提交过程中测量的最高内存消耗量和一个工作进程在一次运行期间测量的最高内存消耗量。我们对每个数据集大小和算法的唯一组合重复该过程20次以获得稳健结果。我们在图3中报告20次运行的平均值和标准差。

图3. SyNCoPy内存效率。所选算法的峰值内存消耗（PMC）作为输入大小的函数。PMC测量基于合成数据。起始数据集大小为10个试次、5,000个样本和50个通道。每个数据点显示20次独立运行的PMC平均值和标准差。（A）PMC在很大程度上与试次数无关。测试数据集的总大小在几乎三个数量级上变化（10个试次到7,000个试次，约10 MB到7 GB），而单个试次的大小保持在1 MB不变。（B）PMC取决于每试次的样本数和算法。样本数（信号长度）从10变化到40,000。（C）PMC取决于通道数和算法。通道数从2变化到250。

6.2 峰值内存消耗—结果

峰值内存消耗（PMC）测量结果如图3所示。首先，我们研究了试次数对PMC的影响（图3A）。我们将试次数从10逐步增加到40,000，同时保持每试次样本数和通道数不变。在每个数据点，我们使用相应算法进行20次独立运行。对于所有算法，PMC基本保持恒定，与试次数无关。PMC最低的是巴特沃斯滤波器，其次是MTMFFT、相干性和小波。其次，我们展示了增加每试次样本数对内存消耗的影响（图3B）。我们将每试次样本数从10逐步增加到40,000，同时保持试次数和通道数不变。对于小波、使用固定数量锥形的多锥形分析（MTMFFT f.t.）和相干性计算，对PMC有线性影响。对于巴特沃斯滤波器，内存消耗基本恒定，因为对于此方法我们采用SciPy的signal.sosfiltfilt实现，它在输入数据的有限部分上工作。对于完整的多锥形分析（MTMFFT），PCM随样本数呈二次方增加：FFT本身的PMC是样本数的线性函数，实现一致频率平滑所需的锥形数（tapsmofrq参数）也随样本数缩放。最后，我们观察到增加通道数对算法PMC的影响（图3C）。对于小波和MTMFFT，对PMC有线性影响。对于巴特沃斯滤波器，PMC再次几乎恒定。相干性显示PMC随通道数呈二次方缩放，这直接源于组合数学。

 七、讨论 

SyNCoPy是一个用于电生理数据分析的Python包，专注于细胞外电生理学。它通过能够轻松从笔记本电脑扩展到HPC系统从而支持超大数据集，以及类似FieldTrip的API，从类似软件包中脱颖而出。SyNCoPy对大数据的支持基于其架构，该架构（1）允许轻松使用许多科学机构可用的典型HPC系统，（2）仅在需要时将数据从磁盘流式传输到内存，以及（3）将计算隔离在独立计算所需的最小数据量上。我们通过针对多种算法对峰值内存消耗（PMC）进行基准测试，展示了SyNCoPy的内存效率。结果表明，SyNCoPy的架构确实能够提供基本恒定的PMC，与试次数无关。此外，PMC随单试次大小的增加按各自算法的预期进行缩放。

从功能角度来看，SyNCoPy目前专注于原始数据预处理、时频分析和连接性度量。我们预计神经科学用户可能希望将SyNCoPy与包含互补功能的其他成熟软件包（如MNE Python、Elephant等）结合使用。为了促进这一点，我们提供了对转换MNE Python数据结构的支持，以及导入和导出NWB等标准文件格式的功能。此外，SyNCoPy文件格式基于开放标准HDF5和JSON，因此可以被各种语言可用的标准库读取。

SyNCoPy没有图形用户界面，依赖脚本编写。虽然这可能需要完全不懂编程的用户一定的初始时间投入，但我们相信这种方法提供的标准化和增强的可重复性很快就会得到回报。FieldTrip主要基于相同的方法，并已获得庞大的用户社区。为了帮助新用户，SyNCoPy附带完整的API文档，并包括一组演示典型分析工作流程的文章。问题和问题可以在SyNCoPy GitHub存储库上报告和讨论。

 八、局限性 

首先，重要的是要承认内存效率是一种软件需求，在某些情况下与处理速度方面的性能冲突：对于小数据集，一次性将所有内容加载到内存中比按需流式传输数据块更快。然而，对于大数据集，这种计算策略会阻止处理大于（机器RAM的某个分数的）数据集，因此不可行。

其次，SyNCoPy专注于试次并行处理，从我们的角度来看，这在神经科学中是一种非常常见的情景。然而，在某些情况下或对于某些算法，支持沿不同轴的并行化可能是有益的。虽然SyNCoPy确实对某些算法内置了通道并行化的支持，但它通常不支持沿数据集任意轴的并行化。

第三，通过创建自定义计算例程来扩展SyNCoPy的新算法是可能的，但此过程目前需要对并行计算和某些SyNCoPy内部机制有良好的理解，因此面向更高级的用户。

第四，SyNCoPy的目标受众是需要处理更大数据集的神经科学家。数据集大小的确切限制当然取决于特定算法和使用的设置，但始终成立的是单个试次必须轻松放入机器的RAM中，通常是运行计算的HPC集群节点。重要的是要理解，在加载和保存数据时或在算法本身中使用的某些操作将需要在内存中创建一个或多个试次数据的副本。因此，使用大小几乎等于RAM的数据集在现实中是不可行的。这不是SyNCoPy的限制，而是适用于所有计算机操作，包括SyNCoPy内部用于在单个试次数据上实现或运行算法的标准NumPy和SciPy库。所需内存通常是单试次大小的小倍数。

 九、结论 

SyNCoPy为Python中大型电生理数据集的试次工作流提供无缝扩展。在本研究中，我们通过在试次数、每试次样本数和通道数变化的数据集范围内测量多种算法的峰值内存消耗，展示了其扩展到超大数据集的能力。此外，我们展示了如何在真实世界数据集上使用SyNCoPy，以及与使用成熟的FieldTrip工具箱进行的相同分析的直接比较。

SyNCoPy的构建旨在集成到当前的神经科学工具生态系统中。我们希望它将帮助研究人员以可重复的方式处理大数据集，并降低充分利用神经科学现有HPC资源的障碍。