mne溯源后的数据初步处理方法

一些溯源后的数据后续的一些读取，处理和可视化的一些记录。
使用jupyter notebook。

导入库

import mne
import numpy as np
import os
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import os.path as op
from mne.datasets import fetch_fsaverage
from mne.datasets import sample

上节说过，stc的数据保存后，

from mne import SourceEstimate
stc.save('1')
#1-lh.stc
#1-rh.stc
#是这两个文件

stc数据的合成

from mne import read_source_estimate
stc_reloaded = read_source_estimate('1', subject='fsaverage')
#通过这个语句就可以重新读取这个源数据

一些绘图可能用到src，把fwd也计算加载进来

pathdir=r"E:\data\301数据汇总\EEG\jxt_analysis\health\1"
file=r"processed_data.set"
path=os.path.join(pathdir, file)
raw_data = mne.io.read_raw_eeglab(path, preload=True)
trans = "fsaverage"
fs_dir = fetch_fsaverage(verbose=True)
subjects_dir = op.dirname(fs_dir)
raw=raw_data.copy()
bem = op.join(fs_dir, "bem", "fsaverage-5120-5120-5120-bem-sol.fif")
src = op.join(fs_dir, "bem", "fsaverage-ico-4-src.fif")
fwd = mne.make_forward_solution(raw.info, trans=trans, src=src, bem=bem, eeg=True, mindist=5.0, n_jobs=None
)

使用这个标记的测试一下，两个大脑都有，成功导入了。

labels = mne.read_labels_from_annot('fsaverage',  # 使用fsaverage模板parc='aparc',  # 解剖学分区subjects_dir=subjects_dir)  # FreeSurfer数据路径
# 这个parc是功能分区是annotation的内容
v1_label = [label for label in labels if label.name == 'temporalpole-lh'][0]# 可视化ROI
brain = stc_reloaded.plot(subject='fsaverage', hemi="both",subjects_dir=subjects_dir)
brain.add_label(v1_label, color='pink', alpha=1)

Yeo2011_7Networks_N1000

通过ai知道有几个分析静息态的模版

#加载这两个模版
yeo_7_labels = mne.read_labels_from_annot('fsaverage', parc='Yeo2011_7Networks_N1000',subjects_dir=subjects_dir)
yeo_17_labels = mne.read_labels_from_annot('fsaverage', parc='Yeo2011_17Networks_N1000',subjects_dir=subjects_dir)
yeo_labels = yeo_7_labels#选择第一个
yeo_labels = [label for label in yeo_labels if not label.name.startswith('Unknown')]#排除一下错误的脑区

[<Label | fsaverage, ‘7Networks_1-lh’, lh : 21632 vertices>, <Label |
fsaverage, ‘7Networks_1-rh’, rh : 22519 vertices>, <Label |
fsaverage, ‘7Networks_2-lh’, lh : 29778 vertices>, <Label |
fsaverage, ‘7Networks_2-rh’, rh : 30207 vertices>, <Label |
fsaverage, ‘7Networks_3-lh’, lh : 17511 vertices>, <Label |
fsaverage, ‘7Networks_3-rh’, rh : 17604 vertices>, <Label |
fsaverage, ‘7Networks_4-lh’, lh : 17392 vertices>, <Label |
fsaverage, ‘7Networks_4-rh’, rh : 19402 vertices>, <Label |
fsaverage, ‘7Networks_5-lh’, lh : 11380 vertices>, <Label |
fsaverage, ‘7Networks_5-rh’, rh : 11614 vertices>, <Label |
fsaverage, ‘7Networks_6-lh’, lh : 16635 vertices>, <Label |
fsaverage, ‘7Networks_6-rh’, rh : 22330 vertices>, <Label |
fsaverage, ‘7Networks_7-lh’, lh : 35314 vertices>, <Label |
fsaverage, ‘7Networks_7-rh’, rh : 26056 vertices>, <Label |
fsaverage, ‘FreeSurfer_Defined_Medial_Wall-lh’, lh : 14200 vertices>,
<Label | fsaverage, ‘FreeSurfer_Defined_Medial_Wall-rh’, rh : 14110
vertices>]

brain = stc_reloaded.plot(subject='fsaverage', subjects_dir=subjects_dir,hemi="both",background='white', alpha=0.5)
colors = ['red', 'green', 'blue', 'orange', 'yellow', 'purple', 'cyan']
for label, color in zip(yeo_labels, colors):brain.add_label(label, color=color, alpha=0.7)

我们标记了分析的脑功能区。
只标记一个进行选择

Default_Mode_network = [label for label in yeo_labels if label.name == '7Networks_7-lh'][0]
brain = stc_reloaded.plot(subject='fsaverage', hemi="both",subjects_dir=subjects_dir)
brain.add_label(Default_Mode_network, color='green', alpha=0.5)

我们可以看到这个功能网络有不同的脑区部分。
通过目测，我们发现有六个非联通区域。
可以使用聚类的方法，来把这几个区域给分开，并获取到对应的数据。
我们根据pos的位置来进行聚类分析。

# 使用坐标计算，可以计算出坐标，但是在对应coords_as_verts时会出现问题。因为会出现坐标和verts无法转换。
coords = mne.vertex_to_mni(label.vertices, hemis=1, subject='fsaverage', subjects_dir=subjects_dir)
center = np.mean(coords, axis=0)

咨询ai

Default_Mode_network = [label for label in yeo_labels if label.name == '7Networks_7-lh'][0]#选择这个的原因看后面的补充说明
brain = stc_reloaded.plot(subject='fsaverage', hemi="both",subjects_dir=subjects_dir)
brain.add_label(Default_Mode_network, color='green', alpha=0.5)
#我们可以看到这个pos是点的坐标，根据连通性，我们可以进行聚类的分析
Default_Mode_network.pos

绘制一些圆球来代表区域大小和强度

import numpy as np
from mne.viz import Brain
# 定义感兴趣的时间范围（例如：100ms到200ms）
tmin = 0.1  # 起始时间（秒）
tmax = 0.2  # 结束时间（秒）
# 获取时间点对应的索引
time_indices = np.where((stc_reloaded.times >= tmin) & (stc_reloaded.times <= tmax))[0]
# Plot the source time course
brain = stc_reloaded.plot(subject='fsaverage', hemi="both", subjects_dir=subjects_dir)
# 遍历每个子区域
for label in sub_regions:time_course = stc_reloaded.extract_label_time_course(label, fwd['src'], mode='mean')mean_intensity = np.mean(time_course[:, time_indices])intensities.append(mean_intensity)# 创建颜色映射
#使用那个
cmap = plt.get_cmap('coolwarm')  # 使用 'viridis' 颜色映射
norm = plt.Normalize(vmin=min(intensities), vmax=max(intensities))  # 归一化强度值# 遍历每个子区域
for label, intensity in zip(sub_regions, intensities):# 获取中心顶点索引vertex_idx = label.center_of_mass('fsaverage', subjects_dir=subjects_dir)# 将强度映射到颜色color = cmap(norm(intensity))  # 获取颜色值# 固定焦点大小，或根据顶点数量调整size = len(label.vertices) / 1500  # 可根据需要调整缩放因子# 添加焦点和标签brain.add_foci(vertex_idx, coords_as_verts=True, hemi=label.hemi,  # 使用标签的半球信息color=color, scale_factor=size, alpha=0.7)brain.add_label(label, borders=True)# 打印调试信息print(f"Label: {label.name}, Mean Intensity: {intensity}, Color: {color}")# Add the source estimate data

单网络绘制和扩展的方式

不看大脑，只看这个功能网络的分割出的区域，其强度如何。可以以在这个上面扩展数值的连接和计算方式。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.cm import ScalarMappable
from matplotlib.colors import Normalizeintensities = []
coords = []# 定义时间范围
tmin, tmax = 0.1, 0.2
time_indices = np.where((stc_reloaded.times >= tmin) & (stc_reloaded.times <= tmax))[0]# 遍历每个标签获取坐标和强度
for label in sub_regions:# 计算平均强度time_course = stc_reloaded.extract_label_time_course(label, fwd['src'], mode='mean')mean_intensity = np.mean(time_course[:, time_indices])intensities.append(mean_intensity)# 获取中心点坐标hemi = label.hemisrc_idx = 0 if hemi == 'lh' else 1vertex_idx = label.center_of_mass('fsaverage', subjects_dir=subjects_dir)coord = fwd['src'][src_idx]['rr'][vertex_idx]coords.append(coord)coords = np.array(coords)
print(coords.shape)
intensities = np.array(intensities)
print(intensities.shape)
print(intensities)
# 创建 3D 画布
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')# 配置颜色映射
cmap = plt.get_cmap('coolwarm')
norm = Normalize(vmin=intensities.min(), vmax=intensities.max())
colors = cmap(norm(intensities))# 绘制散点图（颜色和大小映射到强度）
scatter = ax.scatter(coords[:, 0], coords[:, 1], coords[:, 2],c=colors, s=np.abs(intensities) * 100 + 50,  # 点大小动态调整（根据强度）depthshade=False,          # 关闭深度阴影（可选）alpha=0.8
)# # 示例：按顺序连接相邻点（可自定义连接逻辑）
# for i in range(len(coords) - 1):
#     ax.plot(
#         [coords[i, 0], coords[i+1, 0]],
#         [coords[i, 1], coords[i+1, 1]],
#         [coords[i, 2], coords[i+1, 2]],
#         color='gray',
#         linewidth=1,
#         alpha=0.5
#     )# 添加颜色条
cbar = fig.colorbar(ScalarMappable(norm=norm, cmap=cmap),ax=ax,shrink=0.6,label='Intensity'
)# 调整视角
ax.view_init(elev=20, azim=45)  # 仰角 20°, 方位角 45°# 坐标轴标签
ax.set_xlabel('X (mm)')
ax.set_ylabel('Y (mm)')
ax.set_zlabel('Z (mm)')plt.tight_layout()
plt.show()

AI补充一下背景知识

关于 Yeo2011_7Networks_N1000，这是由 Yeo 等人在 2011 年提出的一种 基于静息态 fMRI 的功能性脑网络分区图谱，被广泛用于脑功能连接和网络分析。以下是详细介绍和关键资源：

---

📚 背景与研究来源

• 文献来源：
  Yeo, B. T. T., et al. (2011). “The organization of the human cerebral cortex estimated by intrinsic functional connectivity”. Journal of Neurophysiology.
  • 核心贡献：通过对 1000 名被试的静息态 fMRI 数据进行聚类分析，提出了 7 个和 17 个功能网络的分区方案。
  • 版本说明：
    • N1000：表示该分区基于 1000 名被试的数据训练得到（高可靠性版本）。
    • 7Networks：将全脑皮层划分为 7 个宏观功能网络（如图）。

---

🧠 7 个功能网络的定义

网络编号	功能名称	核心脑区	功能描述
1	视觉网络 (Visual)	枕叶 (Occipital)	视觉信息处理
2	躯体运动网络 (Somatomotor)	中央前/后回 (Pre/Postcentral gyrus)	运动和体感控制
3	背侧注意网络 (Dorsal Attention)	顶内沟 (IPS), 额眼区 (FEF)	空间注意和眼动控制
4	腹侧注意网络 (Ventral Attention)	颞顶联合区 (TPJ), 前岛叶 (AI)	注意重新定向、突显信息检测
5	边缘网络 (Limbic)	前扣带回 (ACC), 眶额皮层 (OFC)	情绪处理和动机调控
6	额顶控制网络 (Frontoparietal Control)	背外侧前额叶 (DLPFC), 后顶叶 (PPC)	执行控制和工作记忆
7	默认模式网络 (Default Mode)	内侧前额叶 (mPFC), 后扣带回 (PCC)	自我参照思维、静息态活动

---

📂 标签数据获取

Yeo 分区的标签文件通常以 FreeSurfer 兼容格式存储，您可以通过以下方式获取：

1. 官方来源：

   • 从 Yeo 实验室网站下载：Yeo2011 7/17 Networks Atlas
   • 包含 fsaverage 空间的分区文件（lh.Yeo2011_7Networks_N1000.annot, rh.Yeo2011_7Networks_N1000.annot）。

2. 通过 MNE-Python 直接加载：

   from mne.datasets import fetch_fsaverage
   from mne import read_labels_from_annot# 获取 fsaverage 路径
   fs_dir = fetch_fsaverage(verbose=True)# 读取 Yeo 7网络标签（左右半球合并）
   yeo_labels = read_labels_from_annot(subject='fsaverage', parc='Yeo2011_7Networks_N1000', subjects_dir=fs_dir,hemi='both'
   )
   

---

🔍 标签文件内容解析

• 文件结构：

  • 每个半球（lh/rh）的标签文件包含 7 个网络 + 1 个内侧壁（Medial Wall）标签。
  • 标签命名格式：7Networks_{网络编号}-{半球}，如 7Networks_1-lh。
  • 内侧壁标签：FreeSurfer_Defined_Medial_Wall-{半球}（非功能网络区域）。

• 顶点数量：

  • 每个网络的顶点数反映其在大脑皮层的空间范围（如默认模式网络通常最大）。

---

🛠 使用注意事项

1. 空间配准：

   • 确保使用与分区模板相同的标准空间（如 fsaverage）进行源估计和标签提取。

2. 网络合并：

   • 若需合并左右半球的同一网络（如双侧视觉网络），需手动操作：

     vis_network = [label for label in yeo_labels if '7Networks_1' in label.name]
     # 废除，使用上面的联合读取的方案
     combined_vis_label = vis_network[0] + vis_network[1]
     

3. 排除内侧壁：

   • 内侧壁标签通常不参与功能分析，需过滤：

     yeo_networks = [label for label in yeo_labels if 'Networks' in label.name]
     

---

📊 典型应用场景

1. 功能网络时间序列提取：

   # 提取默认模式网络（网络7）的时间序列
   dmn_labels = [label for label in yeo_labels if '7Networks_7' in label.name]
   dmn_time_series = stc.extract_label_time_course(dmn_labels, src, mode='mean')
   

2. 网络间功能连接分析：

   • 计算网络间 Pearson 相关性或相位同步性。

---

# 这个提取的是默认模式网络的子区域，并对时间段的平均激活数据进行可视化network_time_courses = []
network_names = []
for label in sub_regions:# 提取该网络区域的时间序列(平均值)time_course = stc_reloaded.extract_label_time_course(label,fwd['src'], mode='mean')network_time_courses.append(time_course.squeeze())  # 移除多余维度network_names.append(label)  # 存储网络区域名称
plt.figure(figsize=(15, 10))  # 增大画布尺寸
times = stc_reloaded.times# 关键优化点：下采样数据 & 标记重要时间窗
plot_step = 10  # 每10个时间点取1个（根据实际数据量调整）
smooth_window = 5  # 滑动平均窗口大小（可选）
right_xlimit = 0.0
left_xlimit = 2.5
for i, (time_course, name) in enumerate(zip(network_time_courses, network_names)):# 数据平滑和下采样smoothed = np.convolve(time_course, np.ones(smooth_window)/smooth_window, mode='same')  # 滑动平均downsampled_t = times[::plot_step]downsampled_data = smoothed[::plot_step]plt.plot(downsampled_t, downsampled_data, label=name, linewidth=1.5, alpha=0.8)  # 适当降低透明度# 聚焦关键时间段（示例：-0.2~1.0秒）
plt.xlim([right_xlimit, left_xlimit]) # 标记事件时间点
plt.axvline(0, color='k', linestyle='--', linewidth=0.8)  # 刺激呈现时间
plt.axvspan(0.1, 0.3, alpha=0.1, color='red')  # 标记感兴趣的时间窗# 优化坐标标签
plt.xticks(np.arange(right_xlimit, left_xlimit, 0.1),  # 设置刻度间隔为0.1秒rotation=45)  # 旋转标签避免重叠
plt.xlabel('Time (s)', fontsize=12)
plt.ylabel('Normalized Activation', fontsize=12)plt.title('Network Activation Time Courses (Smoothed)', fontsize=14)
plt.legend(ncol=2, loc='upper right', fontsize=10)  # 分两列显示图例
plt.grid(True, linestyle=':', alpha=0.5)
plt.tight_layout()  # 自动调整布局
plt.show()

network_time_courses = np.array(network_time_courses)
#(6, 30051)# 计算关联矩阵
correlation_matrix = np.corrcoef(network_time_courses)# 打印结果
print("关联矩阵形状:", correlation_matrix.shape)
print("关联矩阵:\n", correlation_matrix)
import seaborn as sns
import matplotlib.pyplot as plt# 设置画布
plt.figure(figsize=(10, 8))# 绘制热图
sns.heatmap(correlation_matrix, annot=True,  # 在单元格中显示数值fmt=".2f",  # 数值格式为两位小数cmap='coolwarm',  # 颜色映射vmin=-1, vmax=1,  # 相关系数范围xticklabels=network_names,  # X轴标签yticklabels=network_names,  # Y轴标签linewidths=0.5  # 单元格边框宽度
)# 添加标题和标签
plt.title('Network Correlation Matrix', fontsize=14)
plt.xlabel('Networks', fontsize=12)
plt.ylabel('Networks', fontsize=12)# 显示图形
plt.tight_layout()
plt.show()

应该可以根据这个如果相关性都比较强，或者说，这些区域都是正激活（大部分激活）的时候，判断其是主导功能网络。