pycrostates：脑电微状态分析的Python工具包

导读

本文转载自流浪心球公众号。

微状态分析在本公众号上已有好几篇的学习介绍，亦可参考以下链接：

EEG数据的微状态分析教程

EEG微状态资源整合

EEG微状态分析学习资料Mark

微状态分析是一种允许调查EEG记录的时空特征的方法。它包括将多通道 EEG 信号分解为一系列准稳定状态，每个状态的特征在于其头皮电位的空间分布，也称为微状态图或微状态地形。

这种分解基于两个连续步骤：允许定义最能代表所研究数据的拓扑结构的聚类，以及将先前提取的拓扑结构之一分配给一个或多个 EEG 记录的每个时间点的后拟合。

这些方法依赖于将时间点分配给最相似的微状态图，这就是为什么定义两个拓扑之间的距离如何计算很重要的原因。对于微观状态分析，空间相关性的绝对值的倒数被用作距离的度量来执行所有计算。使用绝对值是为了忽略地形极性。

Pycrostates工具包能更便捷地完成脑电微状态的分析，pycrostates.cluster.ModKMeans 通过pycrostates.cluster.ModKMeans.fit() 方法执行 clustering，并通过pycrostates.cluster.ModKMeans.predict() 方法进行backfitting。它还实现了其他方法，以便于分析和显示结果。

第一部分安装

1.1 运行依赖环境

mne(>=1.0)numpy(>=1.20)scipyscikit-learnmatplotlibpoochdecoratorjinja2

要求您使用 Python 3.7 或更高版本。您可以选择通过 conda 或 pip 安装pycrostates。pycrostates 与最新的稳定版 MNE-Python 配合使用效果最好。为确保 MNE-Python 是最新的，请参阅他们的安装说明：https://mne.tools/stable/install/index.html。

1.2 通过 Conda 安装

要通过conda安装最新的稳定版本，请在终端中使用：

$ conda install -c conda-forge pycrostates

1.3 通过 Conda 安装

要通过pip安装最新的稳定版本，请在终端中使用：

$ pip install pycrostates

要安装最新的开发版本，请运行：

$ pip install git+https://github.com/vferat/pycrostates

第二部分使用示例

详细的使用示例教程，请查阅：https://pycrostates.readthedocs.io/en/latest/auto_tutorials/index.html

2.1 预处理

预处理流程中主要有：重采样和提取全局场功率峰值。

2.1.1 重采样

首先导入数据，并进行重参考。

import mnefrom mne.io import read_raw_eeglabfrom pycrostates.datasets import lemonraw_fname = lemon.data_path(subject_id='010004', condition='EC')raw = read_raw_eeglab(raw_fname, preload=True)raw.crop(0, 30)raw.pick('eeg')raw.set_eeg_reference('average')

我们现在可以使用resample（）函数为我们的数据绘制n_samples的n_resamples，其中n_sample定义每个epochs中包含的样本数，n_resamples定义要绘制的epochs数。

from pycrostates.preprocessing import resampleresamples = resample(raw, n_resamples=10, n_samples=150, random_state=40)resamples

我们还可以使用 'coverage' 参数根据我们想要覆盖的原始数据量自动计算前两个参数之一。例如通过设置 n_resamples和 coverage：

resamples = resample(raw, n_resamples=10, coverage=0.5, random_state=40)resamples

resamples = resample(raw, n_samples=150, coverage=0.5, random_state=40)resamples

最后，我们也可以使用这个函数进行重采样Epochs()

epochs = mne.make_fixed_length_epochs(raw, duration=2, preload=True)resamples = resample(epochs, n_samples=150, coverage=0.5, random_state=40)resamples

2.1.2 提取全局场功率峰值

然后我们可以使用该extract_gfp_peaks() 函数来提取具有最高全局场功率的样本。min_peak_distance 允许选择 2 个选定峰之间的最小样本数。

from pycrostates.preprocessing import extract_gfp_peaksgfp_data = extract_gfp_peaks(raw, min_peak_distance=3)gfp_data

这个功能也可以用在Epochs()

epochs = mne.make_fixed_length_epochs(raw, duration=2, preload=True)gfp_data = extract_gfp_peaks(epochs, min_peak_distance=3)gfp_data

gfp_data 可用于其他预处理步骤，例如resample()

from pycrostates.preprocessing import resampleresample = resample(gfp_data, n_resamples=1, n_samples=100) #extract 1 resample of 100 random high gfp samples.resample[0]

或者直接拟合聚类算法

from pycrostates.cluster import ModKMeansn_clusters = 5ModK = ModKMeans(n_clusters=n_clusters, random_state=42)ModK.fit(gfp_data, n_jobs=5)ModK.plot()

2.2 聚类

from mne.io import read_raw_eeglabfrom pycrostates.cluster import ModKMeansfrom pycrostates.datasets import lemonraw_fname = lemon.data_path(subject_id='010017', condition='EC')raw = read_raw_eeglab(raw_fname, preload=True)raw.crop(0, 30)raw.pick('eeg')raw.set_eeg_reference('average')

修改后的 Kmeans 可以用 n_clusters要计算的集群中心数来实例化。默认情况下，修改后的 Kmeans 仅适用于 EEG 数据，但这可以通过picks参数进行修改。可以在类定义期间定义 random_state 以获得可重现的结果。

n_clusters = 5ModK = ModKMeans(n_clusters=n_clusters, random_state=42)

大多数方法都需要拟合修改后的 Kmeans。这可以通过数据结构mne.io.Raw或mne.epochs.Epochs数据结构来完成。请注意，根据您的设置，您可以进行更改n_jobs=1以使用并行处理并减少计算时间。

ModK.fit(raw, n_jobs=5)

现在我们的算法已经安装好了，我们可以可视化集群中心，也称为微状态图或微状态拓扑图ModK.plot()。请注意，此方法使用Info拟合实例的对象来显示拓扑。

ModK.plot()

cluster_centers_由于以下属性，可以将集群中心作为一个 numpy 数组访问：

ModK.cluster_centers_

集群中心可以使用重新排序ModK.reorder_clusters()：

ModK.reorder_clusters(order=[3, 2, 4, 0, 1])ModK.plot()

并使用重命名ModK.rename_clusters()：

ModK.rename_clusters(new_names=['A', 'B', 'C', 'D', 'F'])ModK.plot()

由于方法，地图极性可以反转ModK.invert_polarity() 。请注意，它只影响可视化，它在反向拟合期间没有影响，因为极性被忽略了。

ModK.invert_polarity([False, False, True, True, False])ModK.plot()

最后，改进的Kmeans可用于预测使用ModK的微状态分割。predict（）方法。默认情况下，注释为“坏”的段不会被标记，但可以使用reject_By_ANTATION参数更改此行为。通过将因子参数设置为>0（默认情况下不平滑因子=0），可以对输出序列执行平滑，而half_window_size参数用于指定平滑时间跨度。最后，reject_edges参数允许不分配每个记录（或每个历元）的第一段和最后一段，因为它们可能不完整。它对raw的影响应该很小，但在处理epochs时可能很重要。

segmentation = ModK.predict(raw, reject_by_annotation=True, factor=10, half_window_size=10, min_segment_length=5, reject_edges=True)segmentation.plot(tmin=1, tmax=5)

2.3 组水平分析

2.3.1 单个聚类簇的组水平分析

from mne.io import read_raw_eeglabfrom pycrostates.cluster import ModKMeansfrom pycrostates.datasets import lemonfrom pycrostates.preprocessing import extract_gfp_peaksfrom pycrostates.io import ChDatacondition = "EO"subject_ids = ["010020", "010021", "010022", "010023", "010024"]

在本例中，我们首先从GFP峰值（被试水平分析）计算单个拓扑图。然后，我们将各个拓扑图连接到单个数据集中，并将其提交给第二轮聚类（组水平分析）。

import numpy as npModK = ModKMeans(n_clusters=5, random_state=42)n_jobs = 2individual_cluster_centers = list()for subject_id in subject_ids: # Load Data raw_fname = lemon.data_path(subject_id=subject_id, condition=condition) raw = read_raw_eeglab(raw_fname, preload=True) raw = lemon.standardize(raw) raw.pick("eeg") # For sake of time, we only use 30s of recording. raw.crop(0, 30) raw.set_eeg_reference("average") # Extract Gfp peaks gfp_peaks = extract_gfp_peaks(raw) # Subject level clustering ModK.fit(gfp_peaks, n_jobs=n_jobs) individual_cluster_centers.append(gfp_peaks.get_data())group_cluster_centers = np.hstack(individual_cluster_centers)group_cluster_centers = ChData(group_cluster_centers, ModK.info)# Group level clusteringModK.fit(group_cluster_centers, n_jobs=n_jobs)ModK.plot()

我们可以根据需要重新整理聚类结果。

ModK.reorder_clusters(order=[2, 3, 4, 1, 0])ModK.rename_clusters(new_names=["A", "B", "C", "D", "F"])ModK.plot()

我们现在可以将组水平映射回溯到每个单独的记录并提取微状态参数。

import pandas as pdms_data = list()for subject_id in subject_ids: # Load Data raw_fname = lemon.data_path(subject_id=subject_id, condition=condition) raw = read_raw_eeglab(raw_fname, preload=True) raw = lemon.standardize(raw) raw.pick("eeg") raw.crop(0, 30) # for sake of time raw.set_eeg_reference("average") segmentation = ModK.predict(raw, factor=10, half_window_size=8) d = segmentation.compute_parameters() d["subject_id"] = subject_id ms_data.append(d)ms_data = pd.DataFrame(ms_data)ms_data

2.3.2 单个GFP值的组水平分析

from mne.io import read_raw_eeglabfrom pycrostates.cluster import ModKMeansfrom pycrostates.datasets import lemonfrom pycrostates.preprocessing import extract_gfp_peaks, resamplefrom pycrostates.io import ChDatacondition = "EO"subject_ids = ["010020", "010021", "010022", "010023", "010024"]

在本例中，我们首先提取单个GFP峰。然后，我们将它们连接到单个数据集中，以便将其提交给聚类（组水平分析）。

import numpy as npModK = ModKMeans(n_clusters=5, random_state=42)n_jobs = 2individual_gfp_peaks = list()for subject_id in subject_ids:# Load Data raw_fname = lemon.data_path(subject_id=subject_id, condition=condition) raw = read_raw_eeglab(raw_fname, preload=True) raw = lemon.standardize(raw) raw.pick("eeg")# For sake of time, we only use 30s of recording. raw.crop(0, 30) raw.set_eeg_reference("average")# Extract GFP peaks gfp_peaks = extract_gfp_peaks(raw)# Equalize peak number across subjects gfp_peaks = resample( gfp_peaks, n_resamples=1, n_samples=880, random_state=42 )[0] individual_gfp_peaks.append(gfp_peaks.get_data())individual_gfp_peaks = np.hstack(individual_gfp_peaks)individual_gfp_peaks = ChData(individual_gfp_peaks, raw.info)# Group level clusteringModK.fit(individual_gfp_peaks, n_jobs=n_jobs)ModK.plot()

我们可以根据需要重新整理聚类结果。

ModK.reorder_clusters(order=[0, 2, 4, 3, 1])ModK.rename_clusters(new_names=["A", "B", "C", "D", "F"])ModK.plot()

我们现在可以将组水平映射回溯到每个单独的记录并提取微状态参数。

2.4 指标：选择聚类中心的数量 k

在本教程中，我们将看到如何获得关于要使用的最佳聚类中心数 k 的指示。重要的是要注意，这种选择没有单一的解决方案。它总是在聚类质量、根据文献分析微观状态的可能性和分解表达的方差之间进行权衡。

import matplotlib.pyplot as pltfrom mne.io import read_raw_eeglabfrom pycrostates.cluster import ModKMeansfrom pycrostates.datasets import lemonraw_fname = lemon.data_path(subject_id="010017", condition="EO")raw = read_raw_eeglab(raw_fname, preload=True)raw.crop(0, 30)raw.pick("eeg")raw.set_eeg_reference("average")

Pycrostates 实施了几个指标来评估聚类解决方案的质量而不知道基本事实：- The silhouette score（越高越好） - The Calinski Harabasz score（越高越好） - The Dunn score（越高越好） - The Davies Bouldin score（越低越好）。

我们可以将这些指标应用于使用不同的 n_clusters 值计算的聚类解决方案，并分析哪些可以提供最佳聚类解决方案。

from pycrostates.metrics import ( silhouette_score, calinski_harabasz_score, dunn_score, davies_bouldin_score,)cluster_numbers = range(2,10)silhouette_scores = []calinski_harabasz_scores = []dunn_scores = []davies_bouldin_scores = []for k in cluster_numbers: ModK = ModKMeans(n_clusters=k, random_state=42) ModK.fit(raw, n_jobs=12) # silhouettee silhouette_scores.append(silhouette_score(ModK)) # calinski and harabasz calinski_harabasz_scores.append(calinski_harabasz_score(ModK)) # dunn dunn_scores.append(dunn_score(ModK)) # davies bouldin davies_bouldin_scores.append(davies_bouldin_score(ModK))

2.4.1 Silhouette score

plt.figure()plt.scatter(cluster_numbers, silhouette_scores)plt.xlabel('n_clusters')plt.ylabel('Silhouette score')plt.show()

2.4.2 Calinski Harabasz score

plt.scatter(cluster_numbers, calinski_harabasz_scores)plt.xlabel('n_clusters')plt.ylabel('Calinski Harabasz score')plt.show()

2.4.3 Dunn score

plt.figure()plt.scatter(cluster_numbers, dunn_scores)plt.xlabel('n_clusters')plt.ylabel('Dunn score')plt.show()

2.4.4 Davies Bouldin

plt.figure()plt.scatter(cluster_numbers, davies_bouldin_scores)plt.xlabel('n_clusters')plt.ylabel('Davies Bouldin score')plt.show()

可以看出，对于选择哪个n_cluster值，没有全球共识。Silhoutettete得分似乎n_clusters=3，Calinski-Harabasz得分似乎n_clusters=4，Dunns core得分似乎n_ clusters=8和Davies-Bouldin得分似乎n_clusters=3或4的。因此，科学家需要根据实际研究来选择合适的n_ clusters，并在以上指标中进行相应的权衡或妥协

谢谢大家观看，如有帮助，来个喜欢或者关注吧！

本文仅供学习参考，有任何疑问及建议，扫描以下公众号二维码添加交流：

更多学习内容，仅在知识星球发布：

继续阅读

阅读原文