接手第一批fMRI数据那会儿我压根没把“分割”俩字放心上。当年从扫描仪拷回来的文件夹一堆DICOM文件我满脑子想的都是“赶紧跑预处理”结果第一步就卡了壳——软件不认这些格式文件东一个西一个连哪个序列是功能像、哪个是结构像都分不清。后来才意识到fMRI原始数据分割听起来像是个不起眼的杂活实际上是整个数据分析流程能不能顺利跑通的命门。这篇笔记就是把我这几年处理原始数据时摸索出来的一套分割、整理和校验的方法记录下来给刚入坑fMRI数据处理的研究生以及被多被试、多session数据折磨到崩溃的科研助理做个参考。这里说的“分割”不是指把大脑影像分割成灰质、白质那种组织分割而是数据层面的拆分与组织把一堆混乱的原始文件整理成一个个干净、可追溯、能被下游工具直接读取的标准文件。搞清楚这件事后面所有预处理步骤才有意义。1. 为什么“分割”不只是一个搬运文件的机械动作1.1 扫描仪导出的原始数据到底长什么样MRI扫描仪输出到光盘或服务器上的数据绝大多数是DICOM格式。这个格式本身没有问题它是医学影像的通用标准医院PACS系统全靠它运行。但问题在于DICOM文件到了科研分析环境中非常”笨重“一个功能像序列可能包含几百上千个独立文件一个slice一个文件还夹杂着定位像、校准像、各种名字看不懂的序列。分析工具像FSL、SPM、AFNI虽然也能读DICOM但效率低而且容易把不同序列搞混。还有一种情况有些研究组已经用配套软件把DICOM转成了4D的NIfTI文件一个文件包含整个时间序列比如说bold_4D.nii.gz。那是不是就不需要分割了不一定。后面要做slices到volume的一致性检查、剔除前几个不稳定volume、按条件切分回归量都需要把4D文件拆开处理。“分割”在这条流程里的意思是双重的一是把DICOM按序列拆成NIfTI二是把4D NIfTI按时间点或空间位置拆开。1.2 分割要解决的是三个层面的混乱第一层是序列层面一次扫描会出来几十个序列要能从中挑出T1结构像、BOLD功能像、场图field map、扩散像等并且正确配对不张冠李戴。第二层是时间层面功能像通常是4D数据分割成单个3D volume后才能逐帧检查头动、信号异常多run实验也需要把连续采集的数据切回每个run对应的独立文件。第三层是空间层面个别场景下要单独提取某个slice或某个区域来检查比如确认slice order设置是否对应实际采集顺序。如果这三层混乱没有理清就急着去做头动校正、配准往往会在结果里埋雷而且一旦进到后面步骤再回头排查原始数据问题代价极大。所以我现在的习惯是拿到数据的第一天先不做任何高难度处理专心把数据组织结构搞清楚。这一步踏实了后面全流程都能受益。2. 分割前先花十分钟摸清DICOM的数据结构2.1 看懂Patient/Study/Series/Image四个层级DICOM的数据组织方式类似于医院病历的归档逻辑。最顶层是Patient对应一个受试者往下是Study对应一次检查比如某天在扫描仪上做的一次完整扫描再往下是Series对应一个扫描序列比如T1_MPRAGE、BOLD_rest、fieldmap_phase最底层是Image对应单个slice或单个信号采集对应的文件。要注意的是不同厂家的扫描仪导出的文件夹命名习惯完全不同西门子经常是带一堆序号的全大写目录GE是.IMA后缀居多飞利浦则喜欢用随机字符串命名。靠文件夹名字猜内容很不靠谱一定要靠DICOM头部信息来识别。2.2 用工具快速列出一份“扫描清单”我习惯在分割前先用dcm2niix的查看模式打印一份序列摘要不需要任何图形界面终端敲一行命令就行dcm2niix -v -o /tmp/dump /path/to/raw_data-v会输出每个序列的详细信息包括SeriesNumber、SequenceName、ImageType、TR/TE、矩阵大小、层数等。把这份输出存成一个文本文件作为原始数据的索引档案。如果数据是.IMA或是嵌套目录dcm2niix会自动递归找到所有DICOM文件所以不用自己写find命令去逐个统计。如果只想快速知道有哪些序列、每个序列多少张图还可以用dcmdumpDCMTK工具包配合grep来提取关键标签比如dcmdump $(find /path/to/raw_data -type f | head -n 1) | grep -E SeriesDescription|SeriesNumber|EchoTime|RepetitionTime这一步的核心目的是在动手分割之前先把“这一批数据里到底有什么”这个问题回答清楚。我见过有人拿到数据就盲目批量转换转完才发现漏了场图序列或者把定位像当成功能像转了出来浪费了半天时间。3. 两条主流分割路线的实操对比3.1 路线一dcm2niix 一把梭DICOM转NIfTI这是绝大多数情况下的首选。dcm2niixChris Rorden出品是目前转换DICOM到NIfTI的事实标准跨平台支持几乎所有主流扫描仪的私有格式还能直接按BIDS风格输出JSON文件。我最常用的转换命令长这样dcm2niix -z y -o /data/BIDS/sub-01/func -f sub-01_task-rest_run-1 %d /path/to/SERIES_FOLDER逐个解释参数-z y输出压缩的NIfTI文件.nii.gz能省下大概一半的磁盘空间。-o指定输出目录。我习惯提前按BIDS目录结构建好转换完直接落到正确位置。-f指定输出文件名模板。%d会被替换为SeriesNumber加序列描述避免重名更精细的做法是自己拼名字比如sub-01_task-rest_run-1。最后的参数指向某个序列的文件夹如果直接指向受试者总目录dcm2niix默认会把所有序列都转出来。对新手来说一次全转然后手动挑文件也是个可行策略但前提是你已经对着2.2的清单核对过一遍。dcm2niix在转换时还会自动生成.json文件里面记录了TR、TE、翻转角、slice timing等关键参数。这个文件别删后续做时间层校正、设计矩阵时全靠它。还有一个隐藏好处部分扫描仪比如西门子的VB系列在DICOM头里的slice order信息不完整甚至错误dcm2niix会尝试从SequenceName和私有标签推断虽然不能保证100%正确但至少比一个外行空手查DICOM头靠谱得多。3.2 路线二fslsplit 按时间点把4D序列切成3D文件很多时候你手里的数据已经是4D NIfTI比如别人已经帮你转好了或者你从OpenNeuro这类公开数据库下载的数据。这时候要做的是把4D数据按时间点切开比如检查每个volume的头动、去掉前5个预扫描volume、把静息态数据按连续帧切块做滑动窗口分析。FSL里的fslsplit一行命令搞定fslsplit sub-01_task-rest_bold.nii.gz sub-01_task-rest_vol -t参数-t表示按时间轴第4维切分输出文件会得到sub-01_task-rest_vol0001.nii.gz、sub-01_task-rest_vol0002.nii.gz这样按顺序编号的3D文件。默认编号是四位数如果你的数据超过9999个volume实际上几乎不可能需要加-z参数调整填充位数。比如fslsplit sub-01_task-rest_bold.nii.gz vol_ -t -z如果不希望压缩输出可以加--nocompress但一般没必要。切完之后要验证一下数量对不对总共应该等于原始4D文件的时间点数。用一条命令确认fslhd sub-01_task-rest_bold.nii.gz | grep dim4 # 或者 python -c import nibabel as nib; print(nib.load(sub-01_task-rest_bold.nii.gz).shape)我在实际工作中还经常用到一个反向操作fslmerge把切开的文件再合并回去。这个在特殊场景下很有用比如你手动去掉了一些坏volume之后可以fslmerge -t把剩下的文件重新拼成一个完整的4D序列继续跑后续流程。3.3 如果只想抽特定几个时间点或某个slice有一个需求频率也很高只提取某个时间点或者只看某个位置的slice。比如你想确认第10个volume和第100个volume之间有没有明显的全局信号跳变或者你想把第20层的图像单独导出来画个图。fslroi最顺手# 提取第10个volume索引从0开始 fslroi input_4D.nii.gz vol_10.nii.gz 10 1 # 提取第20层slice保留所有时间点 fslroi input_4D.nii.gz slice_20.nii.gz 0 -1 0 -1 20 1 # 提取第5到第10个time point共6个volume fslroi input_4D.nii.gz sub_bold_05_10.nii.gz 5 6它的参数逻辑是输入 输出 xmin xsize ymin ysize zmin zsize没有写的维度默认全取。写完参数建议自己拿fslhd核对一下输出的维度这个ROI提取虽然不难但三维坐标写反的情况我见过不少次。如果你更习惯用Pythonnibabel的方式同样直观import nibabel as nib import numpy as np img nib.load(input_4D.nii.gz) data img.get_fdata() # 取第10个volume索引为9 vol data[:, :, :, 9] # 取第20层slice所有时间点 sl data[:, :, 19, :] # 保存 nib.save(nib.Nifti1Image(vol, img.affine), vol_10.nii.gz)这个路子适合要顺带做点数据检查的情况比如算一下某个slice的时间序列均值、看有没有全0序列等。4. 分割后的一致性检查参数不能拍脑袋4.1 用fslhd逐项核对关键元数据分割操作完成后最大的风险在于文件路径是对的命名也对但里面的参数和实验设计对不上。比如你计划里写TR2000ms实际扫描TR1500ms或者你打算用隔层采集interleaved实际是顺序采集sequential。这些如果不在源头纠正后面所有套用固定参数的脚本都会乖乖地把错误放大。我的检查清单固定包含这几项全部来自fslhd或nibabel的输出参数检查目的dim1/dim2/dim3空间维度是否和扫描协议一致如64x64x30dim4时间点数是否等于实际采集volume数pixdim1/2/3体素大小功能像常见3x3x3mm或3.4x3.4x3mmpixdim4TR数值注意单位是秒srow_x/y/z头动初始位置信息可以顺带确认方向是否左右颠倒qform_code/sform_code坐标映射是否有效通常应该是1或2有时dcm2niix转出来的文件pixdim4会显示一个很小的小数比如0.000833那是因为某些扫描仪把TR记录为毫秒但文件内部按秒又有误差。遇到这种情况不要急着改数字先回看.json文件里的RepetitionTime字段以那个为准后续预处理时手动传入正确的TR即可。4.2 用BIDS风格命名文件和目录而不是随手取名分割整理出来的数据最终要服务的是下游所有人包括未来的你自己。我在吃过几次亏之后彻底倒向了BIDSBrain Imaging Data Structure命名规范。它不只是一套命名规则更是一套能机器可读的数据组织标准。简单来说功能像放在func/目录文件名按sub-label_task-task_run-run_bold.nii.gz的格式来结构像放anat/场图放fmap/。这个规范的好处是FSL、SPM、fMRIPrep等工具可以直接识别你换一个人来分析数据也不会一头雾水。一个典型的BIDS目录结构长这样sub-01/ ├── anat/ │ └── sub-01_T1w.nii.gz ├── func/ │ ├── sub-01_task-rest_bold.nii.gz │ └── sub-01_task-rest_bold.json ├── fmap/ │ ├── sub-01_phasediff.nii.gz │ └── sub-01_magnitude1.nii.gz └── sub-01_scans.tsv当然如果你只是自己跑通一条分析流程不打算公开发布数据BIDS看起来有点重。但至少文件名里要带上subjID_task_scanDate这类信息别用001.nii.gz这种根本追溯不到出处的名字。4.3 建立一份扫描记录表比任何记忆都可靠分割、命名都做完了还差最后一件事把每一步整理的依据记录下来。我习惯在每个被试的目录下放一个scans_notes.tsv里面包含列名示例值original_series_number8dicom_series_descriptionBOLD_EPI_restoutput_filenamesub-01_task-rest_bold.nii.gznum_volumes240TR_sec2.0TE_ms30slice_order_sourcedcm2niix jsonnotes前5个volume为dummy已计划剔除这份表不只是给别人看的更是给两周后的自己看的。人的记忆在大量数据面前非常不可靠尤其是同时处理二三十个被试的时候记录表就是你的导航系统。5. 分割过程中几个让人血压升高的翻车现场5.1 多回波数据被当成单回波一股脑合并现在很多研究组用multiband多回波序列比如Siemens的CMRR序列每个时间点会采集两个或更多个回波TE不同。这种数据在DICOM层面是同一个序列里的多个Series或者在同一个Series里交替存储。如果直接用dcm2niix默认参数去转它可能会把多个回波合并到一个4D文件里或者生成带_e1、_e2后缀的多个文件具体行为取决于扫描仪和版本。我遇到过的最坑的情况是三个回波被合并成了一个3倍时间点数的文件第一眼看过去还以为是正常数据结果后面做glm时激活图一团糟。排查了好久才发现时间点数不对。处理多回波数据时务必用2.2的清单确认好每个回波的图像数量然后分割成独立的回波文件。dcm2niix有一个-m参数可以设置合并策略但我更推荐先把各回波拆开等预处理阶段再按需要合并或做加权平均这样始终保留原始信息。5.2 slice order方向搞反时间层校正直接废掉时间层校正slice timing correction非常依赖slice order信息。slice order指的是扫描仪采集slice的先后顺序常见的有sequential descending从顶部到底部、sequential ascending从底部到顶部、interleaved隔层采集等。如果分割时没有正确记录slice order或者转换工具推断错了后面的时间层校正就是在错误的方向上插值相当于给每个voxel的时间序列引入系统性的相位错误。怎么核对slice order最靠谱的方式是看扫描仪的DICOM头或扫描参数页。西门子数据可以在.json文件里看SliceTiming字段GE和飞利浦则需要翻序列参数。dcm2niix生成的.json里通常有SliceTiming数组里面有每个slice的采集时间点把这个数组画出来或者输出到文本和扫描协议里的采集方式对照一下基本就能确认。我曾经在分不清SeqDesc和SeqAsc的情况下硬跑了整个流程结果发现事件相关设计的统计结果一片乱麻浪费了整整一周。现在我的习惯是做时间层校正之前先把slice order打印出来贴在工位上。5.3 压缩格式带来的兼容性陷阱dcm2niix默认输出.nii.gz我很喜欢因为省空间。但有一个坑老版本的SPM比如SPM8、SPM12某些配置对.nii.gz的支持不稳定某些插件和脚本会直接报“Unknown file type”错误。遇到这种情况不一定是文件损坏多半是工具不认压缩格式。解法很简单解压即可gunzip sub-01_task-rest_bold.nii.gz或者用fslchfiletype统一转格式fslchfiletype NIFTI input.nii.gz output.nii类似的兼容性坑还有nii和.img/.hdr两种NIfTI封装。fslchfiletype NIFTI转成.img/.hdr对老工具最友好。我的建议是数据存档用.nii.gz给具体工具跑之前先确认输入要求省得在报错上反复折腾。5.4 在原始DICOM目录里直接改文件新手最容易犯的一个错误是在扫描仪导出的原始文件夹里直接做修改重命名DICOM文件、移动子文件夹、甚至为了腾空间删除某些序列。这么做非常危险一来原始DICOM文件是研究对象出了问题很难重新获取二来很多DICOM文件之间的关联依赖文件名和内部标签手动乱改可能破坏序列完整性最后实验记录和伦理审查往往要求保留原始数据不留痕迹地改动会让你后面写方法部分时无据可依。我自己的硬性铁律是任何处理都在副本上进行。拿到原始数据先做整体备份然后在另一个工作目录里操作。分割产生的中间文件可以被随意删除重建但原始文件夹始终保持复制过来后的初始状态。6. 分割结果如何与下游预处理无缝衔接6.1 分割后先做一轮简便QC别急着跑全流程分割整理好的数据在进入fMRIPrep或者FSL全套流程之前值得花几分钟做一次视觉和数值上的QC。我的做法分三步第一看整体信号。用fslview或fsleyes加载一个中间volume检查有没有明显的全脑信号丢失、头部外伪影、明显的死区。第二看时间序列稳定性。用fslmeants提取全脑均值时间序列fslmeants -i sub-01_task-rest_bold.nii.gz -o global_signal.txt然后快速画个图观察有没有突然的大幅跳变或缓慢漂移。如果发现某些volume的信号特别异常记下来后面处理。第三估算头动。我一般先跑一次mcflirt不看最终结果只看它输出的.par文件里平移和旋转参数的峰值头动超过3mm或3度的数据直接标记为高危。这三步做完数据能跑什么级别的分析心里就有底了。6.2 把分割脚本化集成到批处理或BIDS流程里手动处理一个被试还可以接受但二三十个被试还手动搞纯属折磨自己。我强烈建议把分割和整理过程写成脚本。最简单的形式是一段Bash循环处理多个被试的DICOM数据自动创建BIDS目录结构然后调用dcm2niix转换。#!/bin/bash for subj in sub-01 sub-02 sub-03; do mkdir -p /data/$subj/{anat,func,fmap} dcm2niix -z y \ -o /data/$subj/func \ -f ${subj}_task-rest_bold \ /raw/$subj/SERIES_8 done脚本化之后最大的好处是可复现无论谁在什么时间重跑只要原始数据不变得到的分割结果完全一致。这比命令行记在聊天记录里可靠多了。更进一步可以研究一下dcm2bids这个工具它专门用来把DICOM转成BIDS格式通过一个JSON配置文件映射序列和输出文件名整个分割过程只需要两条命令dcm2bids_scaffold -o /data/BIDS dcm2bids -d /raw/sub-01 -p 01 -c dcm2bids_config.json配置文件的写法不复杂核心就是把序列描述和BIDS标签对应起来{ descriptions: [ { dataType: func, modalityLabel: bold, customLabels: task-rest_run-1, criteria: { SeriesDescription: BOLD_rest } } ] }我第一次跑通这套批量转换的时候最大的感慨是之前花在人工整理数据上的时间简直可以攒出一个完整的周末。早点把这些重复劳动交给脚本专业的人做专业的事数据整理同样值得认真对待。从最原始的DICOM一堆散件到规规矩矩的BIDS目录中间经过的就是分割和整理这一步。它不复杂但特别考验耐心和细心。以我自己的经验分割时偷的懒都会在后面的预处理和统计分析阶段加倍还回来。宁可多花一晚上把数据检查得清清楚楚也别在结果一团糟的时候才回头追悔。