freesurfer
freesurfer是一个处理大脑3D结构像数据,进行自动皮层和皮下核团分割的工具,用起来非常方便。freesurfer wiki上的教程也非常详细,但是有一点,freesurfer的命令很复杂,很难准确地记住每个参数该怎么设置。本人比较懒,不愿记,也记不住,每次都需要打开wiki进行对照。由于wiki非常详尽,每次都是在一大篇英文中搜索命令。在这里弄一个简洁版,只把分析流程所用到的命令贴在这里,以便查阅。 一、数据处理
freesurfer分析3D,最好是原始的dcm数据,不要进行数据转格式转换和坐标变化,原始数据就可以。所以把数据放在sub1,sub2,sub3…….,这样我们就可以用循环来做了。我们的计算机中心用的是PBS的系统,进行并行运算: #!/bin/bash
# SUBJECTS_DIR is where you want to put your result SUBJECTS_DIR=your_subjects_path
imge=`ls your_subjects_path/sub1/*.IMA|head -n 1` recon-all -s sub1 -i ${imge} recon-all -all -s $sub1
这就是分析一个被试的代码,然后用前面介绍过的sed命令,把sub1替换成其他被试编号,就生成了其他被试代码。然后在端口敲入命令:sh 代码文件。数据就开始分析了。 二、数据检查
数据检查主要tkregister2、 tkmedit、tksurfer三个命令结合起来。但是在这种视觉检查以前,应该先看前面recon-all.log文件是否报错。可是recon- all.log文件好大了,怎么办呢,用grep命令,查找一下文件里面有没有”error”,并输出含有“error”的行数: grep -n “error” recon-all.log
如果都没有错误,那就OK.另外还有一个命令也很有意思,会自动帮我们察看是否存在top错误:
mris_euler_number sub1/surf/lh.orig mris_euler_number sub1/surf/lh.white mris_euler_number sub1/surf/lh.pial
如果这三个命令生成的数字完全一样,就没有top结构问题,然后我们再进行视觉检查,首先检查register:
tkregister2 --mgz --s sub1 --fstal --surf orig
然后检查tkmedit 和tksurfer检查白质、灰质分割问题。
tkmedit sub1 brainmask.mgz rh.white -aux wm.mgz -aux-surface lh.white tksurfer sub1 rh inflated
检查的时候,一般在c130-170这部分slices问题比较严重,总是有脑膜被看成是灰质了,需要编辑brainmask.mgz,把它删掉。具体检查参看freesurfer wiki. 检查完毕后,根据编辑过的地方重新跑一下数据。根据recon-all的步骤,先register,然后是-autorecon2-cp,然后是-autorecon2-wm, -autorecon2-pial,-autorecon3。当register出现问题时,几乎需要完全重新算数据,其他的就从编辑过的最早步骤开始。例如,一个被试修改过cp、wm,那么就从cp开
始重新算。
recon-all -subjid sub1
recon-all -autorecon2-cp -autorecon3 -subjid sub2 recon-all -aurorecon2-wm -autorecon3 -subjid sub3 recon-all -autorecon2-pial -autorecon3 -subjid sub4
数据重新算后,还要一个一个地进行视觉检查一下,才能进行下一步数据统计分析. 三、数据统计分析
开始进行数据分析了。首先把数据对齐到freesurfer自带的fsaverage空间上。数据跑完后,会自动出现在被试目录中。先把每个被试的数据叠加成一个4D的文件,而叠加顺序事先要设计好,把其他相关的变量也放进去。这里举一个简单例子:三组被试,每组两个被试:AA型、AB型、BB型,另外还有被试年龄数据,那么这个文件内容如下: GroupDescriptorFile 1 Title g3v1 Class AA Class AB Class BB Variables Age Input sub1 AA 22 Input sub2 AA 19 Input sub3 AB 20 Input sub4 AB 15 Input sub5 BB 22 Input sub6 BB 18
当编辑好这个这个文件后,命名为g3v1.fsgd,意思是被试分为三组,一个连续变量。就运行下面的命令:
mris_preproc --fsgd g3v1.fsgd --target fsaverage --hemi lh --meas \\thickness --out lh.thickness.00.mgz
mris_preproc --fsgd g3v1.fsgd --target fsaverage --hemi rh --meas \\thickness --out rh.thickness.00.mgz
上面的命令把被试的皮层厚度数据按照g3v1.fsgd中的被试顺序排列成一个4D的文件,左右半球分别生成一个文件。但是这些数据没有进行 smooth.然后我们进行10mm的smooth。可以用mri_info lh.thickness.00.mgz查看文件信息,如是不是有6个被试等。
mris_surf2surf --hemi lh --s fsaverage --sval lh.thickness.00.mgz --\\fwhm 10 --cortex --tval lh.thickness.10.mgz
mris_surf2surf --hemi rh --s fsaverage --sval rh.thickness.00.mgz --\\fwhm 10 --cortex --tval rh.thickness.10.mgz
生成的文件lh.thickness.10.mgz和rh.thickness.10.mgz是最重要的数据文件.可以用mri_info命令察看数据信息。
下一步就要进行数据统计,在进行数据统计之前,我们还要做一些准备,举上面的例子吧。
我想要检验三组被试皮层厚度是否存在差异,总体被试的皮层厚度与年龄的关系,每个组被试的皮层厚度与年龄的关系,不同组被试年龄与皮层厚度的关系是否有差异。我要做这些统计分析,就需要有design matrix. freesurfer默认design matrix从文件中读取,这些文件叫做mtx文件。准备mtx文件需要首先计算回归子的数据,公式为:Nregressors = Nclasses*(Nvariables+1) = 3*(1+1) = 6.这里我们的有三组,一个连续变量。所以有6各回归子。 考察AA与AB的组间差异则为:1 -1 0 0 0 0 考察AA与BB的组间差异则为:1 0 -1 0 0 0 考察AB与BB的组件差异则为:0 1 -1 0 0 0
所有被试皮层厚度与年龄的关系:0 0 0 0.333 0.333 0.333 AA组内皮层与厚度的关系:0 0 0 1 0 0
AA组与AB组皮层与厚度关系是否存在差异:0 0 0 1 -1 0
AA组与AB和BB组皮层厚度与年龄的关系是否存在差异:0 0 0 1 -0.5 -0.5 现在我们准备好了上面的所有文件,分别命名为 1.mtx、2.mtx、3.mtx....... 7.mtx 终于可以最后统计了,命令如下:
mri_glmfit --y lh.thickness.10.mgz --fsgd g3v1.fsgd --C C1.mtx --C \\C2.mtx --C C3.mtx --C C4.mtx .... --C C 7.mtx --surf fsaverage lh --\\cortex --glmdir g3v1.lh
这样结果就存放在g3v1.lh文件夹中,在这个文件夹中,还会生成7个子文件夹,每个文件夹对一个统计比较,里面最重要的文件时一个叫做sig.mgz的文件 四、数据结果的整理和报告
当统计分析完后,需要对数据进行查看和整理,到底显著没有啊。首先肉眼看一下吧,看一下AA与AB有没有差异,p<0.01:
tksurfer fsaverage lh inflated -annot aparc.annot -fthresh 2 -overlay \\g3v1.lh/C1/sig.mgh 如果有激活,下一步就是进行校正,这里就做Clusterwise Correction。
mri_glmfit-sim --glmdir g3v1.lh/C1 --sim mc-z 5000 4 mc-z.abs --sim-sig \\abs --overwrite 这里用的Monte Carlo检验, 5000次, p<0.0001, Z进行绝对值比较,这样可以看双侧了,如果只想看AA>AB,那么就用“mc-z.pos”和 “--sim-sig pos”。overwright就是要把原始数据覆盖,生成一个新的CSD文件,放在C1文件夹中。这是一个文本文件。这一步需要几个小时的时间,这一步完成后,查看校正后的结果:
tksurfer fsaverage lh inflated -annot g3v1.lh/C1/mc-\\z.abs4.sig.ocn.annot -fthresh 2 -curv -gray 或者
tksurfer fsaverage lh inflated -annot g3v1.lh/C1/mc-\\z.abs4.sig.ocn.annot -fthresh 2 -overlay g3v1.lh/C1/mc-\\z.abs4.sig.cluster.mgh 这就是最后的结果。