三亩地 三亩地SAN MU DI · CODE DIARY
ARTICLE DETAIL

日记详情

真实记录编程学习的某一天,欢迎挑你感兴趣的翻一翻。

从零开始学fMRI:用Python+SPM12处理你的第一个静息态数据(保姆级避坑指南)

从零开始学fMRI:用Python+SPM12处理你的第一个静息态数据(保姆级避坑指南)

从零开始学fMRI:用Python+SPM12处理你的第一个静息态数据(保姆级避坑指南)

当你第一次接触静息态功能磁共振成像(fMRI)数据时,面对DICOM文件、NIfTI格式、预处理流程这些专业术语,可能会感到无从下手。本文将以OpenNeuro公开数据集为例,带你用Python和SPM12完成从原始数据到功能连接分析的全流程,重点解决实际操作中的常见问题。

1. 实验环境搭建与数据准备

1.1 软件安装与配置

处理fMRI数据需要以下核心工具链:

  • SPM12:最常用的fMRI分析工具包,需配合MATLAB运行
  • dcm2niix:DICOM转NIfTI格式的高效工具
  • Python环境:推荐Anaconda+Jupyter Notebook组合

安装步骤:

# 安装dcm2niix (Linux/macOS) git clone https://github.com/rordenlab/dcm2niix.git cd dcm2niix make && sudo make install # Python依赖库 pip install nibabel numpy scipy matplotlib

注意:SPM12需要MATLAB R2015b或更高版本支持,安装后需在MATLAB中添加路径:

addpath('/path/to/spm12') savepath

1.2 数据获取与检查

从OpenNeuro下载数据集(如ds000247)后,首先验证数据完整性:

import os import nibabel as nib data_dir = '/path/to/ds000247' subj_dir = os.path.join(data_dir, 'sub-01', 'func') # 检查功能像文件 func_file = os.path.join(subj_dir, 'sub-01_task-rest_bold.nii.gz') img = nib.load(func_file) print(f"数据维度:{img.shape}") # 应显示(x,y,z,time)四维

常见问题排查:

  • 若报错"File cannot be read",可能是下载中断导致文件损坏
  • 维度异常可能是扫描参数设置问题

2. 数据预处理全流程

2.1 DICOM到NIfTI格式转换

虽然OpenNeuro提供NIfTI数据,但实际项目中常需处理原始DICOM:

dcm2niix -z y -o output_dir -f %p_%s input_dicom_dir

关键参数说明:

  • -z y:启用Gzip压缩
  • -f %p_%s:命名模式(患者_序列)
  • -o:输出目录

2.2 预处理步骤详解

使用SPM12进行预处理的MATLAB批处理脚本框架:

matlabbatch{1}.spm.spatial.realign.estwrite.data = {'/path/to/func.nii,1'}; matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.quality = 0.9; matlabbatch{1}.spm.spatial.realign.estwrite.roptions.which = [2 1]; spm_jobman('run', matlabbatch);

预处理各阶段注意事项:

步骤关键参数典型值常见问题
头动校正Quality0.9位移>3mm建议剔除
时间层校正Orderascending需确认扫描顺序
空间标准化Voxel size[3 3 3]模板选择影响结果
平滑FWHM[6 6 6]过大会损失细节

提示:预处理顺序应为:时间层校正→头动校正→配准→标准化→平滑

3. 静息态功能连接分析

3.1 基于ROI的分析流程

以默认模式网络(DMN)为例,提取后扣带回(PCC)时间序列:

from nilearn import datasets, input_data # 载入PCC种子点 pcc_coords = [(0, -52, 18)] masker = input_data.NiftiSpheresMasker( seeds=pcc_coords, radius=6, detrend=True, standardize=True ) time_series = masker.fit_transform(func_file)

3.2 全脑功能连接计算

使用Python计算功能连接矩阵:

from nilearn.connectome import ConnectivityMeasure # 使用AAL模板划分90个脑区 atlas = datasets.fetch_atlas_aal() masker = input_data.NiftiLabelsMasker( labels_img=atlas.maps, standardize=True, memory='nilearn_cache' ) # 计算相关矩阵 correlation_measure = ConnectivityMeasure(kind='correlation') correlation_matrix = correlation_measure.fit_transform([time_series])[0]

可视化连接矩阵:

import matplotlib.pyplot as plt plt.figure(figsize=(10, 10)) plt.imshow(correlation_matrix, cmap='coolwarm', vmin=-1, vmax=1) plt.colorbar() plt.title('功能连接矩阵') plt.show()

4. 实战中的典型问题解决方案

4.1 头动校正失败处理

当SPM报错"Image does not match dimensions"时:

  1. 检查NIfTI文件头信息:
    print(nib.load('func.nii').header)
  2. 确认所有图像具有相同维度和体素大小
  3. 尝试重新转换DICOM文件

4.2 空间标准化异常

常见于非标准扫描参数,解决方法:

  1. 手动指定模板:
    matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.template = {'/path/to/EPI.nii'};
  2. 调整正则化参数:
    matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.reg = [0.1 0.1 0.1];

4.3 功能连接结果异常排查

若得到全高/全低相关矩阵:

  1. 检查预处理是否包含去线性趋势步骤
  2. 确认是否进行了全局信号回归
  3. 验证时间层校正顺序是否正确
# 检查时间序列质量 plt.plot(time_series) plt.xlabel('时间点') plt.ylabel('信号强度') plt.title('ROI时间序列') plt.show()

5. 分析结果解读与报告

5.1 关键指标计算

计算网络属性:

from nilearn.connectome import sym_matrix_to_vec from scipy.spatial.distance import squareform # 提取上三角矩阵 vec = sym_matrix_to_vec(correlation_matrix) # 计算平均连接强度 mean_connectivity = np.mean(vec) print(f"平均连接强度:{mean_connectivity:.3f}")

5.2 结果可视化技巧

使用Circos图展示强连接:

from nilearn import plotting # 显示前10%强连接 plotting.plot_connectome( correlation_matrix, node_coords=atlas.coords, edge_threshold='90%', title='功能连接网络' )

5.3 数据保存规范

建议采用BIDS格式组织结果:

results/ ├── sub-01/ │ ├── func/ │ │ ├── sub-01_task-rest_bold_preproc.nii.gz │ │ ├── sub-01_task-rest_bold_confounds.tsv │ │ └── sub-01_task-rest_bold_connectivity.mat └── dataset_description.json
← 返回列表