当前位置： 首页 > news >正文

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

news 来源：原创 2025/5/2 6:12:01

文章目录

      • 导入库
  • Yeo2011_7Networks_N1000
      • 绘制一些圆球来代表区域大小和强度
    • 单网络绘制和扩展的方式
    • AI补充一下背景知识
      • 📚 **背景与研究来源**
      • 🧠 **7 个功能网络的定义**
      • 📂 **标签数据获取**
      • 🔍 **标签文件内容解析**
      • 🛠 **使用注意事项**
      • 📊 **典型应用场景**

一些溯源后的数据后续的一些读取，处理和可视化的一些记录。
使用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 Normalize

intensities = []
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.hemi
    src_idx = 0 if hemi == 'lh' else 1
    vertex_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()

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

相关文章：

  • ubuntu系统/run目录不能执行脚本问题解决
  • 从单任务到多任务：进程与线程如何实现并发？
  • python 标准库之 functools 模块
  • 豪越科技：融合低空经济的消防一体化安全管控解决方案
  • openai agent实践
  • 什么是MCP｜工作原理是什么｜怎么使用MCP｜图解MCP
  • 六十天前端强化训练之第二十七天之Pinia 状态管理全解与购物车实战案例
  • 【Linux】I/O 多路转接：select epoll 技术剖析
  • 安卓 vs iOS 文件系统深度解析：开放自由与封闭安全的终极博弈
  • DeepSeek 助力 Vue3 开发：打造丝滑的表格（Table）之添加导出数据功能示例10，TableView15_10带搜索的导出表格示例
  • [DDD架构]不同数据模型DTO、VO、PO、DAO、DO的含义
  • 自动驾驶系统的车辆动力学建模：自行车模型与汽车模型的对比分析
  • Linux：基础IO---文件描述符
  • JavaSE1.0（实战之图书管理系统）
  • FlowMo: 模式搜索+扩散模型提升图像Token化性能
  • 基于Azure Delta Lake和Databricks的安全数据共享（Delta Sharing）
  • C++异常处理完全指南：从原理到实战
  • 操作系统知识点33
  • 31天Python入门——第10天:深入理解值传递·引用传递以及深浅拷贝问题
  • 计算机网络性能优化相关内容详解
  • 见证历史与未来共舞：上海西岸“蝶变共生”对话讲坛圆满举行
  • 此前显示售罄的火车票“五一”前大量放出来了？12306回应
  • 七部门联合发布《终端设备直连卫星服务管理规定》
  • 网商银行2024年年报发布，客户资产管理规模超过1万亿
  • 小核酸药物企业瑞博生物递表港交所，去年亏损2.81亿元
  • 光明日报：回应辅警“转正”呼声，是一门政民互动公开课