pycrostates:脑电微状态分析的Python工具包
导读
本文转载自流浪心球公众号。
微状态分析在本公众号上已有好几篇的学习介绍,亦可参考以下链接:
微状态分析是一种允许调查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)
scipy
scikit-learn
matplotlib
pooch
decorator
jinja2
要求您使用 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 mne
from mne.io import read_raw_eeglab
from pycrostates.datasets import lemon
raw_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 resample
resamples = 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 提取全局场功率峰值
import mne
from mne.io import read_raw_eeglab
from pycrostates.datasets import lemon
raw_fname = lemon.data_path(subject_id='010004', condition='EC')
raw = read_raw_eeglab(raw_fname, preload=True)
raw.pick('eeg')
raw.set_eeg_reference('average')
然后我们可以使用该
extract_gfp_peaks()
函数来提取具有最高全局场功率的样本。min_peak_distance 允许选择 2 个选定峰之间的最小样本数。from pycrostates.preprocessing import extract_gfp_peaks
gfp_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 resample
resample = resample(gfp_data, n_resamples=1, n_samples=100) #extract 1 resample of 100 random high gfp samples.
resample[0]
或者直接拟合聚类算法
from pycrostates.cluster import ModKMeans
n_clusters = 5
ModK = 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_eeglab
from pycrostates.cluster import ModKMeans
from pycrostates.datasets import lemon
raw_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 = 5
ModK = 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_eeglab
from pycrostates.cluster import ModKMeans
from pycrostates.datasets import lemon
from pycrostates.preprocessing import extract_gfp_peaks
from pycrostates.io import ChData
condition = "EO"
subject_ids = ["010020", "010021", "010022", "010023", "010024"]
在本例中,我们首先从GFP峰值(被试水平分析)计算单个拓扑图。然后,我们将各个拓扑图连接到单个数据集中,并将其提交给第二轮聚类(组水平分析)。
import numpy as np
ModK = ModKMeans(n_clusters=5, random_state=42)
n_jobs = 2
individual_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 clustering
ModK.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 pd
ms_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_eeglab
from pycrostates.cluster import ModKMeans
from pycrostates.datasets import lemon
from pycrostates.preprocessing import extract_gfp_peaks, resample
from pycrostates.io import ChData
condition = "EO"
subject_ids = ["010020", "010021", "010022", "010023", "010024"]
在本例中,我们首先提取单个GFP峰。然后,我们将它们连接到单个数据集中,以便将其提交给聚类(组水平分析)。
import numpy as np
ModK = ModKMeans(n_clusters=5, random_state=42)
n_jobs = 2
individual_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 clustering
ModK.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()
我们现在可以将组水平映射回溯到每个单独的记录并提取微状态参数。
import pandas as pd
ms_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.4 指标:选择聚类中心的数量 k
在本教程中,我们将看到如何获得关于要使用的最佳聚类中心数 k 的指示。重要的是要注意,这种选择没有单一的解决方案。它总是在聚类质量、根据文献分析微观状态的可能性和分解表达的方差之间进行权衡。
import matplotlib.pyplot as plt
from mne.io import read_raw_eeglab
from pycrostates.cluster import ModKMeans
from pycrostates.datasets import lemon
raw_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,并在以上指标中进行相应的权衡或妥协
谢谢大家观看,如有帮助,来个喜欢或者关注吧!
本文仅供学习参考,有任何疑问及建议,扫描以下公众号二维码添加交流:
更多学习内容,仅在知识星球发布:
最新评论
推荐文章
作者最新文章
你可能感兴趣的文章
Copyright Disclaimer: The copyright of contents (including texts, images, videos and audios) posted above belong to the User who shared or the third-party website which the User shared from. If you found your copyright have been infringed, please send a DMCA takedown notice to [email protected]. For more detail of the source, please click on the button "Read Original Post" below. For other communications, please send to [email protected].
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。