Diffusion-weighted in MRI 学习笔记（二）

最近准备写一篇纤维束的文章，于是重新回顾了我之前做的笔记。
其实上一份关于弥散像的学习笔记已经很详细了，基本科普了大致的流程。这个笔记（二）主要是记录一些补充性的知识，以及我的个人体会。所有这些笔记或多或少都会有些知识的遗漏，毕竟我也不是专业弄这个方向的。我建议要着手这方面的同学，去看看官方文档，获取更详细的知识。

MRtrix3官方文档地址：https://mrtrix.readthedocs.io/en/latest/index.html

MRtrix的主要重点功能包括5个：（1）约束球形反卷积估计纤维方向 constrained spherical deconvolution（2）概率流线算法 a probabilisitic streamlines algorithm（3）基于fixels的表观纤维密度分析和纤维截面 fixel-based analysis of apparent fibre density and fibre cross-section（4）定量结构连通性分析 quantitative structural connectivity analysis（5）纤维方向分布图像的非线性空间配准 non-linear spatial registration of fibre orientation distribution images

重要概念：

（1）fiber 还是fibre。 fiber和fibre的区别：fiber是美国拼法，fibre是英国拼法。fiber的意思是（动植物的）纤维，纤维质。fibre的意思是（食物中的）纤维素，（织物或绳等的）纤维制品，（人或动物身体组织及天然物质的）纤维。我的神经解剖学老师教学时，他写的是：有髓神经纤维 myelinated nerve fiber. 而mrtrix文档里面写的都是fibre. 我觉得两种写法都是ok的，只是习惯不一样，就如果grey matter 和 gray matter类似。
（2）fixels，dixels和voxels。voxels是大家做功能像和结构像最常用的，就是3维体素，例如VBM。不过voxels只能代表体积，所以mrtrix作者提出了fixels和dixels，这两种在弥散像分析上比较常用。
fixels是为描述fibre bundle element建立的，dixels是为directional element建立的。提到这些，首先要说到一个名词，就是FOD（fibre orientation distribution）纤维方向分布。在之前的笔记中，我们说到过一个前提假设，就是一条纤维束的最大扩散率取值是主要纤维束上大部分体素的主要方向，因此FOD是连续的，用于拟合纤维束的结果。而如果以体素水平将FOD切割开，那么每一个体素就可以称之为fixels，实际就是voxels叠加了方向变量（4维数组）。因此fixels这个名词的创建就是如此来的。由此可以知道，voxel-based morphology 对应的就是fixel-based morphology，VBM&FBM。所以也就有一些FBM的组统计分析。这种方法是为了更精确的寻找纤维束的某一局部的具体变化效应。fod2fixels函数可以实现该切割步骤。
dixels解释起来比较复杂，它是沿某个预定方向进行的函数样本。这个我就跳过了，因为fixels用的很广泛，而dixels用的少。感兴趣的看官方文档。
（3）全局强度归一化。为什么要进行归一化？原因在于，柏芝内会存在一些compartment(隔间)，一个体素的compartment会影响邻近的compartment的估算，由此产生一些误差。因此作者推荐使用全局强度归一化来减少这种估算误差。但是这个时候又出现新的问题，就是由于核磁的磁场不均性（这个是每个机器的固有误差），如何定义参考区域来执行强度归一化。
所以，作者推荐的方法是，首先进行偏置场校正，这个可以用ANTs的偏置场校正函数dwibiascorrect执行。
参考强度值：理想状态是用每个受试者的脑脊液b=0的位置的强度中位数作为参考值进行归一化，因为脑脊液的位置一般不太会有病理改变。但是由于dwi像的低分辨率，很难获取这个值，所以作者推荐使用受试者全脑白质的中值进行归一化，注意是中位数（median），不是均值。这个前提条件是患者没有广泛的全脑解剖性病理改变，因此中值被认为是可以接受的。
由此，mrtrix的流程就是，先进行偏置场校正，然后中值归一化。

正文

（1）预处理
代码都在我上一份笔记中，我在这里只简单回顾。首先对dwi数据进行降噪，去伪影，头动校正，失真校正，其中头动和失真校正在mrtrix中统称为几何畸变校正。然后紧接着就是偏置场校正，mask估算。这个mask估算步骤中，可以使用默认值，不过需要自己检查下，mask是否符合，如果太小或太大，建议调整下默认值。
（2）纤维方向分布FOD
上一篇提到过，在FOD过程中，对于直直的纤维束估算是很简单的，但是当遇到交叉纤维束或者混合纤维束，就会极大的不准。所以mrtrix提出的是CSD，约束球面反卷积的方法来估算这类纤维束的值。整个过程就是先计算响应函数，然后利用响应函数进行约束球面反卷积，从而计算出FOD。当然，说起来复杂，其实就是一个函数dwi2fod全都实现了。具体的数学公式医学生不需要了解，只要大概了解流程就行。
在FOD过程中，多壳多组织估计是最好的，MSMT，单壳会差一些。但是这个其实受限于数据，对于研究所来说，在收集数据阶段，就会设置好这些参数，但是对于一般临床医院，可能有时候只能用单壳去计算。
（3）强度归一化
这个已经说过了。

AFD apparent fibre density 表观纤维密度

mrtrix介绍了如何实现AFD的过程。AFD是和单壳单组织联用的。我学习以后很受用，这里简单记录下。

预处理

降噪、去伪影、头动和畸变校正、mask估算。

AFD过程

偏置场校正，被试间的强度归一化，upsampling。
（1）偏置场校正上面说过了。
（2）被试间的强度归一化：
首先使用 dwinormalise group 命令，进行组水平的强度归一化。这会生成校正后的图像，和 组水平的FA模板，以及白质mask。
如果后期这个组单独加入了几个额外的被试，可以使用 dwi2tensor命令，如下。用之前的组FA模板和白质mask对这个新被试进行强度归一化。

dwi2tensor new_subject/dwi_denoised_unringed_preproc_unbiased.mif -mask new_subject/dwi_temp_mask.mif - | tensor2metric - -fa - | mrregister -force ../dwinormalise/fa_template.mif - -mask2 new_subject/dwi_temp_mask.mif -nl_scale 0.5,0.75,1.0 -nl_niter 5,5,15 -nl_warp - /tmp/dummy_file.mif | mrtransform ../dwinormalise/fa_template_wm_mask.mif -template new_subject/dwi_denoised_unringed_preproc_unbiased.mif -warp - - | dwinormalise new_subject/dwi_denoised_unringed_preproc_unbiased.mif - ../dwinormalise/dwi_output/new_subject.mif

fixel-based analysis：
如果要使用fixel-based analysis，一个重要的点在于只能使用单壳单组织响应函数，由此所有的FOD纤维束取向分布都会以这个函数进行表示，这样就可以用于表示单个纤维密度的单位。
（3）然后进行 upsampling步骤。这一步可以增加解剖对比度，从而改善图像和后续的所有分析步骤质量。命令如下。

mrgrid IN/dwi_denoised_unringed_preproc_unbiased_normalised.mif regrid -vox 1.25 IN/dwi_denoised_unringed_preproc_unbiased_normalised_upsampled.mif

这里是把图像采样到了1.25mm体素水平。如果dwi图像的原始分辨率本来就小于1.25，就可以跳过。一般临床医院的设备很难达到，所以得执行这一步。
（4）计算brain mask
前面是计算的白质mask，这里是计算brain mask。

dwi2mask IN/dwi_denoised_unringed_preproc_unbiased_normalised_upsampled.mif IN/dwi_mask_upsampled.mif

（5）FOD和CSD
前面说过了，dwiextract 和dwi2fod
（6）生成特定的FOD template
这一部分是根据组内被试的FOD图像，生成组平均的 template，这样用template能代表总体所有。使用的是population_template函数
（7）将组内所有被试注册到FOD template
由于不同被试具有差异，所以需要把组内所有被试对齐到模板上。用mrtransform函数
（8）FOD模板切割成单个fixel
把FOD模板切割成单个fixel，这样就后面可以进行fixel-based analysis。用的是fod2fixel命令。生成fixel_mask
（9）将组内所有被试FOD对齐到 FOD template
这里作者推荐在对齐的时候，不要进行FOD重定向，而是在fixel分割后再进行，会更精确些。

mrtransform IN/wmfod.mif -warp IN/subject2template_warp.mif -reorient_fod no IN/fod_in_template_space_NOT_REORIENTED.mif

（10）分割对齐后的FOD到fixel水平，得到fibre density， FD。
使用的是fod2fixel
（11）重新定位固定点
把对齐后的分割后的FD，再重新定位固定点。是对第9步的补充。用的fixelreorient函数。
（12）将被试的fixel和模板的fixel进行匹配
在上一步重定位以后，可能会出现一些小误差，例如某个被试的fixel有值，但是在模板上这个fixel是不存在的。主要是为了解决这个问题，所以把两者进行匹配，然后模板fixel为0的地方会赋值给这个被试，让他的fixel也变为0。用的fixelcorrespondence函数
（13）计算FC指标，fibre cross-section，纤维横截面
在功能像里，FC代表的是functional connection，功能连接。在弥散像里，这里FC代表的是纤维横截面。第10步11步12步计算出了每个被试的FD，纤维密度。这里计算的是纤维横截面尺寸。注意的是，密度和横截面尺寸是不同的两个指标。横截面变小一般解释为这个fixel有白质萎缩，但是可能这个fixel的纤维密度和其他地方没区别。所以这两个指标需要独立看待。用的warp2metric函数

warp2metric IN/subject2template_warp.mif -fc ../template/fixel_mask ../template/fc IN.mif

这里作者建议使用$lo{g}^{FC}$来代替FC，因为经过log转换，可以使FC变成以零为中心的正态分布。

mkdir ../template/log_fc
cp ../template/fc/index.mif ../template/fc/directions.mif ../template/log_fc
for_each * : mrcalc ../template/fc/IN.mif -log ../template/log_fc/IN.mif

值得注意的是，这个FC和logFC，是对模板很敏感的，因为他们是基于模板计算出来的，所以在统计时，需要在同一个模板下进行统计，也就是说，在我们的研究中，我们要把健康人组和病人组放在一起，计算出一个template，然后在这个template下生成每个被试的FD和FC，然后再去进行组间比较。这个是我觉得很重要的细节。
（14）组成FDC指标
FDC指标实际是FD和FC的组合测量值，mrcalc …/template/fd/IN.mif …/template/fc/IN.mif -mult …/template/fdc/IN.mif。
（15）在FOD template模板上进行概率性纤维追踪
在FOD template进行概率性纤维追踪后，生成的全脑纤维束图，可以作为全脑连通性的信息，这样后续分析中可以基于此计算局部联通性。
用的tckgen函数
（16）SIFT减少纤维束数量
这个在之前的笔记中有提过，这里不说了，用的tcksift函数。
（17）fixel-fixel连接矩阵
用第16步生成的纤维束计算fixel水平的连接矩阵，这个还是模板数据。用的fixelconnectivity函数
（18）被试的fixel平滑
根据上一步的fixel-fixel连接矩阵，对每个被试的fixel进行数据平滑。

fixelfilter fd smooth fd_smooth -matrix matrix/
fixelfilter log_fc smooth log_fc_smooth -matrix matrix/
fixelfilter fdc smooth fdc_smooth -matrix matrix/

（19）统计分析，FD, FC, FDC
做完上述所有步骤，就可以进行统计分析了，统计分析的指标主要就是FD, FC 和 FDC。

fixelcfestats fd_smooth/ files.txt design_matrix.txt contrast_matrix.txt matrix/ stats_fd/
fixelcfestats log_fc_smooth/ files.txt design_matrix.txt contrast_matrix.txt matrix/ stats_log_fc/
fixelcfestats fdc_smooth/ files.txt design_matrix.txt contrast_matrix.txt matrix/ stats_fdc/

最后生成的stats结果，p值是1-p，所以图像只要在0.95以上，就是p值显著了。别忘了多重校正，作者没提多重校正，我觉得是不是需要TFCE之类的进行cluster的校正，毕竟fixel水平的话，太多了，跟voxel一样。

一个重要的注意点：由于FD和FDC是相对度量，因此他们是没有单位的。作者建议在患者-健康人的对比中，一种呈现结果的方法是用两组平均值的百分比变化来代表参数检验中的统计效应值意义。这个具体的可以看官网文档。

至此就结束了AFD的学习之旅。虽然目前看起来有点麻烦，但是其实很简单，过程都很明了，每一步的作用都很清楚。把这个写成脚本以后，就可以批量使用了。

structural connection 结构连接

上面的FD, FC, FDC是fixel-based analysis的一种方式，当然常用的还有另一种方式，就是不关注纤维束局部的某一段，而是ROI节点间的连接变化，后者这种一般会结合图论方法，去计算拓扑属性指标，着重于纤维束结构的整体变化。
我个人比较喜爱的是后者这种图论的方法。这个方法我在之前的笔记中已经对其如何实现阐述的比较明确了。
虽然官方文档没有统计这一章节，但是作者在commands list里面有提到，就是connectomestats函数，可以实现NBS, TFNBS, none三种反方法的网络统计。
不过我建议在得到矩阵连接后，用一些其他的软件进行图论分析，推荐使用braingraph，可以进行GLM, 置换检验，NBS，SCN等，比较丰富。