脑功能MRI数据处理手册系列教程(二)
cac55 2024-10-08 00:51 73 浏览 0 评论
脑功能MRI数据处理手册
作者：范嘉晨
指导老师：刘加成
单位：东南大学（江苏省分子影像与功能影像重点实验室）
目录：
- 本手册介绍
- 相关软件
- 图像格式
- Bold预处理
- Resting-Bold建模、统计分析、结果可视化
- Task-Bold建模、统计分析、结果可视化
- DTI预处理
- TBSS及DTI结果可视化、FDT
- 结构像预处理
- VBM及结果可视化
- 附录一：task-Bold处理相关参数
- 附录二：VBM处理相关参数
- 附录三：resting-Bold处理相关参数
- 附录四：DTI处理相关参数
- 附录五：行为学范式与磁共振扫描对应关系
四.Bold预处理
声明：1.图像文件命名，尽量用字母+数字，如T001C，代表治疗组第一个人第三次实验数据；2.原始图像存放处尽量与图像结果输出处分开；3.图像文件存放目录绝对不要出现中文字符，以免SPM在读取路径时候出错。
预处理的目的是为分析数据做准备，其中一些操作是为了检测和修复数据中的潜在伪影，这些伪影可能是MR扫描仪本身或由被试者所造成的，其他操作是为了后期处理数据做准备。
Bold预处理步骤（包括静息态与任务态）：
打开MATLAB，输入spm，打开spm窗口，在SPM8的Menu中找到Utils下拉列表，点击CD，选择要处理图像的位置。
1.转换图像格式（针对DICOM的原始文件）
在SPM8的Menu中找到DICOMImport，点击此按钮后选择DICOM files,可以在Fit后空格中输入“^”+”特定字母”筛选符合条件的图像，再右击Select All，选择完图像后单击Done按钮。随后出现Output directory，选择输出目录，点击Done，开始图像格式转换。转换结束后（MATLAB里会显示DONE）可以在输出目录看到以s开头的文件（比如静息态单个受试有7920张图片，则转换完后会出现240对（480个）以s开头的文件，每一对文件（.hdr/.img）包含一个时间点采集的Volume,即33层一个Volume。）。
同样，需要把与Bold同层扫描的T1WI轴位图像进行格式转换。步骤同上，此时会产生一对（.hdr/.img）文件。
当然也可以通过Batch的方式批量转换图像格式：点击SPM Menu里的Batch，点击SPM—— Util —— Dicom Import。鼠标在打开的Dicom Import上右击Replicate Module，复制多个Dicom Import。其中Dicom files选择原始dicom文件，Output directory选择输出目录。（注意复制完后，才在每个Dicom Import里输入相应序列的文件），点击左上角绿色三角Run，MATLAB里出现Done代表运行结束，图像格式转换结束。（注意：我们可以在如下图所示一样，选择转换后图像的格式（.hdr/.img的NIfTI还是.nii的NIfTI），记得统一所有数据的格式即可。）最后，保存此batch档，方便下次直接load。
2. 变形校正
最常见的fMRI采集方法，即回波平面成像（EPI），受到气体和组织交界区的伪影影响，如鼻窦或耳道等。这些都是由气体-组织相邻界面引起主磁场（B0）不均匀所致，它表现为两种形式，即信号丢失和几何变形。信号丢失表现为邻近气体-组织界面的脑区信号减低，如眶额皮质和外侧颞叶。鉴别信号丢失模式的一个有用的方法就是将功能图像与其一致的结构图像进行匹配融合。
除了信号丢失，fMRI图像也可以在同一区域发生空间变形。在MRI图像上，当在编码的空间信息施加梯度时，磁场的不均匀性就会导致该图像组织结构的定位错误。大多数情况，前额叶皮质和眶额叶皮质区域是扭曲变形的。一般沿所用磁共振脉冲序列的相位编码方向发生扭曲变形，通常是Y轴（前后）。这些扭曲变形使得功能磁共振数据与结构图像很难匹配一致。使用具有B0磁场特征的场图能在一定程度上矫正磁场不均匀性的影响。大多数磁共振扫描仪均有使用场图的脉冲序列。它们通常在两个不同的回波时间获取图像进行工作。两幅图像之间的相位差别可以用来计算局部磁场的不均匀性，并可以使用这些值来创建一个图测量每个体素移动的距离，对这个图进行反转，可以确定每个体素数据的原始位置。
如果采用变形校正，则应该对校正后的图像进行检查，并且和校正前的图像进行比较，以确保变形校正操作并未引入任何伪影（如果场图有问题，这就有可能发生。）
变形校正可以使用FSL来进行，但是必须有其需要的B0场图及相关文件。
SPM8目前无法进行变形校正。
3.层间时间校正（Slice timing）
几乎所有的fMRI数据采集都使用二维磁共振成像，扫描1次获得1层数据。在某些情况下，层面的获取按升序或降序。另一种方法也称为间隔数据采集，按顺序每隔1层数据采集1次，这样采集到一半的层面数据（如奇数层面），接着采集另一半（如偶数层面）。使用二维采集，就意味着在所述图像的不同部分的数据在系统的不同时间被获取，这些差异最大可达几秒钟（取决于TR）。不同体素采集时间之间存在差异，这些差异对fMRI数据的分析是不利的。事件的发生时间（如一个任务中的试验）是用来创建一个统计模型，代表任务诱发的预期信号。这个数据在每个时间点的数据进行比较，然而，这种分析是假设图像的所有数据是同时获得的，这样就会导致大脑内模型和变化的数据不匹配。层间时间校正就是为了解决不同层面采集时间的不一致。层间时间校正最常见的方法是选择一个参考层面，然后插值所有其他层面的数据，使之与参考层面的时间相一致。这会导致数据集中每个层面在同一个时间点活动。因为扫描仪和脉冲序列不同，在应用层间时间校正时，必须知道扫描序列的相关参数。实践经验显示，相对较短的重复时间（TR≤2秒）和间隔采集，事件相关分析对层面时间问题的解决更为稳健。
本实验数据（静息态、任务态的TR都为2秒，间隔扫描，每个volume都是33层，静息态与任务态扫描序列的区别就是扫描总次数不相同：240次与265次。）
点击SPMMenu的Slice timing按钮，打开BatchEditor,首先Data里输入Dicom Import后的数据（s开头），去掉前5个图像，共计235个文件输入，Number of Slices 输入层数33，TR输入2秒，TA输入1.9394（或者直接输入TR-TR/33），Slice order要注意是间隔扫描，所以输入（1:2:33 2:2:32）或者（1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32）；Reference Slice 建议选择中间层面 33 。所有信息输入完毕后点击左上角绿色三角，开始运行。跑完后MATLAB里会显示Done。此时可以发现我们保存输出文件的目录下多了470个（235对）.hdr/.img文件。这些文件即为层间时间校正后的图像。
4.头动校正（Realignment）
再配合的研究对象也会伴有头部的运动（如因为吞咽，呼吸所致的头动），这对fMRI数据的影响很大。它主要包括两方面：1.大体积运动（bulk motion）会涉及头部整体的运动，而导致时间序列的后续图像不匹配，标准的头动校正技术就是为其而设计，办法就是对照一个参考图像重新调整时间序列的图像。大体积运动对激活图像的影响较大，多见于图像边缘，因为图像边缘的某一像素本来没有脑组织信号，经运动无中生有，图像信号变化较大。不同的运动可能会表现为不同的模式：如正激活或负激活（垂直运动）；一侧正激活，一侧负激活（左右运动）；或在眶额叶皮质的大片正或负激活（沿左右轴旋转）。这种伪影还常见于脑室边缘。
其次，头部运动可能会干扰磁共振信号本身。头部动作时，质子进入相邻的体素，于是扫描重建信号不能准确反映体素所在组织的性质，这就是所谓的自旋历史效应，它会引起一层或一组图像信号强度的明显变化。如果用交替采集，则会出现明暗相间的条纹。这种运动不能用标准头动校正技术纠正，但用ICA等方法或自旋历史校正可能有效。
头动校正（又称为对齐）目的是减少头部运动引起的fMRI时间序列的图像偏差。简而言之，就是给fMRI时间序列中的每幅图像分配一个参考扫描，然后对原始数据重新校正。
头部校正工具通常假设可以用刚体变换校正头部运动，即位置可以改变（移或沿三个轴旋转），但形状不能改变，因此只能校正头部的整体运动。
头动校正和层间时间校正之间可以相互影响，而影响这些的本质是由它们进行的顺序决定的。如果先进行头动校正，那么在同一时间点采集的数据就可能被移动到不同的层面，这样在层间时间校正算法里名义上的采集时间将和实际的采集时间不同。此外，如果有明显的跨平面的运动或旋转，那么头动校正后插入的体素将混入不同时间点采集的数据，再次因层间时间校正导致潜在错误。
如果层间时间校正先进行，就可能造成运动相关信号强度差异（可能非常大）随着时间增大。此外，任何跨平面运动或旋转都将导致与前述相同的名义上的采集时间和实际不同的错误。使用前瞻性头动校正可以在某种程度上减少这个问题，该算法不需要使用插值减少运动影响，因此不会把不同时间点的数据混合。
如果坚持使用层间时间校正，则在头动校正后使用，因为运动对体素信号的影响可能会非常大。然而，如果数据采集使用相对较短的TR则可以避免层间时间校正。
在SPM8Menu中点击Realign(Estimate and Reslice),Data输入层间时间校正完后以as开头的235个文件，其他参数默认。点击左上角绿色三角开始运行，运行结束后MATLAB里会出现Done。
头动校正结束后会输出470个前缀为ras的（235对.hdr/.img）文件；还会输出以rp开头的头动参数txt文件，里面包含头动六个参数（平移三个、旋转三个）。此时可以在SPMGraphics上看到对应此txt文件的六参数图像。（一般建议平移≤±1.5mm，旋转≤±1。）当然对于病人数据来说，标准可以放宽点，平移不超过一个体素直径。还会输出一对meanas开头的文件（.hdr/.img），这个文件是这些图像的一个平均，用于下一步Coregister中作为Reference Image。
5.Coregister（功能像对齐到共面像）
在SPM8的Menu中选择Coregister(Estimate),点击后在弹出的batch中选择图像输入Reference Image (mean开头的文件） 和Source Image（输入同层扫描T1WI文件，s开头），其他参数默认，点击左上角绿色三角run，结束后MATLAB会显示Done。这一步结束后不会产生任何文件，因为它直接在之前的T1文件上修改了
6.空间标准化（Normalization）
大多数情况，我们希望把脑功能成像从个体推广到更大范围的人群。这就需要采集所有个体的数据，然而每个人大脑的大小和形状是不同的，为了能够和其他人大脑在空间上一一对齐，就首先需要进行空间变换。把数据变换到一个共同空间的处理过程称为被试间配准或空间标准化。
尽管不同个体大脑的大小和形状差异很大，但人类的大脑结构还是表现出了高度的一致性。除了那些具有大脑遗传病的人，每个人的大脑都有胼胝体连接的两个半球，不同个体之间形状差异很小，不同的个体有共同的脑沟标志。更深入的研究发现脑组织的细微结构存在较大差异，即使是同卵双胞胎也可能有结构上的重大差异。空间标准化的目的是为了把每个个体的脑图像变换到标准空间，以便减少个体间的差异从而成功进行有意义的组间分析。
在功能MRI的发展中，空间标准化最常用的模板是加拿大蒙特利尔神经病学研究所发明的MNI模板。这些模板是基于MRI的自动化识别模板，而不是基于解剖标志识别的。最先使用的模板称为MNI305,模板，随后通过将高分辨率的图像配准到MNI305创建了另一个称为ICBM-152的模板，需要注意的是，MNI305和ICBM-152模板有一些细微差异。
（共面像对齐到高分辨率结构像T1.nii，高分辨率结构像配准到标准空间ICBM space template）：在SPM8 Menu中点击Normalise (Estimate &Write),在Source Image（Coregister后文件即s开头的T1文件）、Images toWrite（Coregister后文件即s开头的T1文件）、Template Image中（SPM8自带templates里的T1.nii）输入相应文件。点击左上角的绿色三角Run，在MATLAB里显示done即表示完成，会在原来存放s开头的T1文件处出现一个以sn.mat结尾的文件以及一对ws开头的（.hdr/.img)文件。
（将功能像对齐到共面像sT1，共面像对齐到高分辨率结构像T1.nii，高分辨率结构像配准到标准空间ICBM，这三步串联起来，将本地空间的功能像变换到标准空间。）：再次点击SPM menu下的Normalise(Write),在Parameter File中输入上一步Normalise(Estimate &Write)后产生的文件（sn.mat文件），在Imagesto Write处输入Realign(Estimate & Reslice )后产生的文件（ras开头的235个文件），其他参数默认，点击左上角绿色三角run，MATLAB中显示done就是完成。此时，在保存ras文件处会出现一系列wras文件，共235对（.hdr/.img）。
串联多步变换变为一个单一的变换，而不是在每一步变换中都对功能像进行插值，减少每一步插值过程中累积的误差和模糊。如果需要对一幅图像连续进行多次插值变换，就应该使用高阶插值方法，如正弦插值方法来最小化错误和模糊。
7.空间平滑（Smoothing）
空间平滑需要对图像应用一个滤波器以消除高频信息。对fMRI数据选择应用空间平滑功能有许多原因：1.通过删除高频信息（即小规模的图像变化），使图像平滑增加了较大空间信号的信噪比。因为在fMRI研究中，大多数激活跨越许多像素点，获得较大的特征性信号可能带来的好处要远远大于失去较小特征的坏处。2.较小像素的采集可以帮助减少伪影易感区域的信号丢失，空间平滑可以帮助克服采用小体素发生的噪声增加。3.当数据是由不同的个体组合起来的时候，众所周知，功能区的空间位置存在变异，而不能由空间标准化进行校正，通过模糊空间数据，空间平滑可以帮助减少个体间不匹配，这是以牺牲空间分辨率为代价的。4.有一些分析方法（特别是高斯随机场理论），需要一个特定的空间平滑核。
空间平滑最常见的方法是用一个三维高斯滤波器（或核）的三维图像的卷积。高斯滤波器施加平滑核的数量由分布的宽度决定。在统计中，这就是通常描述的标准偏差，而在图像处理中，分布宽度被描述为半峰值全宽（半高全宽FWHM），这就在最大值一半的地方测量峰的宽度，FWHM越大，图像就越平滑。
平滑图像所应用的图像的平滑核不一定是相同的，平滑核描述了相邻像素点之间的相关性。
图像平滑处理程度：1.如果对图像进行平滑处理是为了减少噪声，那么应该应用不超过想要检测的激活信号大小的一个滤波器。2.如果进行平滑处理时为了减少解剖变异的影响，那么最优平滑宽度取决于要成像研究人群的数量变化及在何种程度上可以通过空间标准化方法减少这种变化。3.如果进行平滑处理是以确保统计分析高斯随机场理论的有效性，那么两倍体素大小的FWHM是合适的。一般来说，建议宁可少做一些而不是做过多平滑处理，因此，我们推荐体素尺寸的两倍作为一个合理的FWHM起点。
在SPM的Menu里，点击Smooth，在Imagesto Smooth里输入上一步Normalise (Write)后产生的wras开头的235个文件。其他参数默认，点击左上绿色三角Run，在MATLAB里出现Done代表smooth跑完了。此时会生成235对（.hdr/.img）以swras开头的文件。
8.批处理
一步一步按上述步骤做，不仅繁琐而且容易出错。可以将上述步骤做成一个batch档，这样既方便，而且对于大量数据的分析也方便。当然利用脚本命令运行磁共振数据分析会更好，但是编码语言的学习不是一蹴而就的。
点击SPM的Menu里Batch，按下图设置即可，所有含<-X的条目代表需要输入相应文件进去，当所有X都消失后，此时点击左上角绿色三角，即可运行。当MATLAB里显示Done即表示batch跑完。（注意：在Named Directory Selector里面的Directory里需输入原始文件的目录。Slice Timing中的Session需要输入s开头的经Dicom Import后的原始图片；Coregister Estimate 中的Source Image需要输入s开头的经Dicom Import后的T1WI共面图像。）我们可以在左上角File里选择Save Batch and Script,方便下次直接点击SPM Menu里的Batch后直接Load进来，避免再一次的输入。
我们可以在左上角File里选择Save Batch and Script,方便下次直接点击SPM Menu里的Batch后直接Load进来（选_job.m的文件），避免再一次的输入。
上面这种方法只是处理了单个受试的一个Bold序列，如果一个患者有多个Bold序列或有多个受试的Bold序列时，如何才能提高预处理效率呢？
首先完成上面的步骤建立名为preprocessing的batch档（preprocessing.m和preprocessing_job.m）；其次，点击SPM Menu里的Batch，点击左上角basicIO,打开Run Batch Jobs，按下图设置。其中Job Files里放入（preprocessing.m和preprocessing_job.m），Directory放入输出目录，其下面第一个NIfTI Images 放入Dicom Import后以s开头的240-5个原始图像文件，第二个NIfTIImages放入s开头的T1WI共面像。在Run Batch Jobs上鼠标右击点击Replicate Module即可复制多个Run Batch Jobs。复制完后，再在里面分别输入每个序列的相关文件（绝对不要在输入完每个序列的相关文件后再复制，这样每个Run Batch Jobs 里就都是同一个序列的文件）。左上角按绿色的三角Run，跑完后MATLAB里会显示Done。通过这种方式，我们就可以一次跑几个甚至几十个序列的预处理，大大减少了预处理出错的可能，提高数据预处理的效率。
预处理注意事项
预处理是正确建模及统计分析的前提，一定要保证所有相同序列预处理的参数是一致的，不要遗漏任何序列的图像；此外尽量选用统一版本的MATLAB及SPM，如所有Stroop任务功能磁共振的数据都使用MATLAB-R2015b及SPM8分析。
相关推荐
