统计分析

生成好study-specific模板, 下来要做统计分析组间差异了. 如果用fsl提供的vbm方案, 这一步对应的就是fslvbm_3_proc这个命令了.

fslvbm_3_proc

这一步和实验设计和研究问题本身很有关系. 这里讲讲fslvbm_3_proc命令的内部过程吧.

这个文件保存在 $FSLDIR/bin中, 可以用vim查看vim $FSLDIR/bin/fslvbm_3_proc. shell脚本的主体如下(省略帮助和说明文字)

#!/bin/sh
export LC_ALL=Cecho [`date`] [`hostname`] [`uname -a`] [`pwd`] [$0 $@] >> .fslvbmlogmkdir -p stats
cd strucecho "Now running the preprocessing steps and the pre-analyses"/bin/rm -f fslvbm3a
for g in `$FSLDIR/bin/imglob *_struc.*` ; doecho $gecho "${FSLDIR}/bin/fsl_reg ${g}_GM template_GM ${g}_GM_to_template_GM -fnirt \"--config=GM_2_MNI152GM_2mm.cnf --jout=${g}_JAC_nl\"; \$FSLDIR/bin/fslmaths ${g}_GM_to_template_GM -mul ${g}_JAC_nl ${g}_GM_to_template_GM_mod -odt float" >> fslvbm3a
done
chmod a+x fslvbm3a
fslvbm3a_id=`${FSLDIR}/bin/fsl_sub -T 40 -N fslvbm3a -t ./fslvbm3a`
echo Doing registrations: ID=$fslvbm3a_idcd ../statscat <<stage_preproc2 > fslvbm3b
#!/bin/sh

\$FSLDIR/bin/imcp ../struc/template_GM template_GM\$FSLDIR/bin/fslmerge -t GM_merg     \`\${FSLDIR}/bin/imglob ../struc/*_GM_to_template_GM.*\`
\$FSLDIR/bin/fslmerge -t GM_mod_merg \`\${FSLDIR}/bin/imglob ../struc/*_GM_to_template_GM_mod.*\`\$FSLDIR/bin/fslmaths GM_merg -Tmean -thr 0.01 -bin GM_mask -odt char/bin/cp ../design.* .for i in GM_mod_merg ; dofor j in 2 3 4 ; do\$FSLDIR/bin/fslmaths \$i -s \$j \${i}_s\${j} \$FSLDIR/bin/randomise -i \${i}_s\${j} -o \${i}_s\${j} -m GM_mask -d design.mat -t design.con -Vdone
donestage_preproc2chmod a+x fslvbm3bfslvbm3b_id=`${FSLDIR}/bin/fsl_sub -T 15 -N fslvbm3b -j $fslvbm3a_id ./fslvbm3b`echo Doing subject concatenation and initial randomise: ID=$fslvbm3b_idecho "Once this has finished, run randomise with 5000 permutations on the 'best' smoothed 4D GM_mod_merg. We recommend using the -T (TFCE) option. For example:"
echo "randomise -i GM_mod_merg_s3 -o GM_mod_merg_s3 -m GM_mask -d design.mat -t design.con -n 5000 -T -V"

这个shell脚本中, 采用cat命令配合EOF符号生成了另两个脚本: fslvbm3a 和fslvbm3b . 值得一提的是, 这里采用fsl_sub命令调用并行计算进行计算加速. 但是这个命令是基于Sun grid cluster的, 所以对于普通的台式机, 或者工作站, 建议将这一步去掉, 直接运行fslvbm3a或者fslvbm3b.

我们再从头看看这个脚本. 首先来说说生成脚本fslvbm3a的部分.

fslvbm3a

这一步是对统计分析前的准备工作.

使用fsl_reg将所有灰质图像和study-specific的模板配准, 生成${g}_GM_to_template_GM, 之后乘以Jac_nl, 生成${g}_GM_to_template_GM_mod ${g}是指被试样本的名称.

fslvbm3b

1) 进入states 文件夹, 拷贝template_GM到该文件夹.

2) 用fslmerge将(../struc中的所有_GM_to_template_GM)图像沿着”时间”轴拼接concatenate到一起, 生成GM_merge
同时, 将(../struc中所有_GM_to_template_GM_mod)图像沿着”时间”轴concatenate到一起, 生成GM_mod_merg.

当然, 这里的”时间”轴, 只是一个序号而已. 并不没有时间属性. 以下是fslmerge命令的用法:

Usage: fslmerge <-x/y/z/t/a/tr> <output> <file1 file2 .......> [tr value in seconds]-t : concatenate images in time-x : concatenate images in the x direction-y : concatenate images in the y direction-z : concatenate images in the z direction-a : auto-choose: single slices -> volume, volumes -> 4D (time series)-tr : concatenate images in time and set the output image tr to the final option value

3) 对GM_merg 所有图像中信号强度是否0.01做binary(非0即1)的voxel沿着”时间”轴做平均, 设置输出的结果为字符char.
fslmaths GM_merg -Tmean -thr 0.01 -bin GM_mask -odt char

4) 将design.*文件全部拷贝到./stat中. 之后对GM_mod_merg中的每个图片,
对每个图片做高斯平滑, 分别选择平滑核宽度为2mm 3mm 和4mm. 然后用randomise命令做非参数统计. 将刚才进行过高斯平滑的图像作为输入, 输出文件名称为${g}_s${j} ${g}为被试文件名称, ${j}是高斯核宽度. GM_mask 是前一步生成的mask文件, -d后面是设计文件design.mat, -t后面紧跟着的是t检验设置文件,design.con 最后的-V是指在t检验时使用方差平滑.

fslmaths \$i -s \$j \${i}_s\${j}
\$FSLDIR/bin/randomise -i \${i}_s\${j} -o \${i}_s\${j} -m GM_mask -d design.mat -t design.con -V

附录:

1) fslmaths的基本调用格式如下:

fslmaths [-dt <datatype>] <first_input> [operations and inputs] <output> [-odt <datatype>]

fslmaths Usage: fslmaths [-dt <datatype>] <first_input> [operations and inputs] <output> [-odt <datatype>]Datatype information:-dt sets the datatype used internally for calculations (default float for all except double images)-odt sets the output datatype ( default is float )Possible datatypes are: char short int float double input"input" will set the datatype to that of the original imageBinary operations:(some inputs can be either an image or a number)-add   : add following input to current image-sub   : subtract following input from current image-mul   : multiply current image by following input-div   : divide current image by following input-rem   : modulus remainder - divide current image by following input and take remainder-mas   : use (following image>0) to mask current image-thr   : use following number to threshold current image (zero anything below the number)-thrp  : use following percentage (0-100) of ROBUST RANGE to threshold current image (zero anything below the number)-thrP  : use following percentage (0-100) of ROBUST RANGE of non-zero voxels and threshold below-uthr  : use following number to upper-threshold current image (zero anything above the number)-uthrp : use following percentage (0-100) of ROBUST RANGE to upper-threshold current image (zero anything above the number)-uthrP : use following percentage (0-100) of ROBUST RANGE of non-zero voxels and threshold above-max   : take maximum of following input and current image-min   : take minimum of following input and current image-seed  : seed random number generator with following number-restart : replace the current image with input for future processing operations-save : save the current working image to the input filenameBasic unary operations:-exp   : exponential-log   : natural logarithm-sin   : sine function-cos   : cosine function-tan   : tangent function-asin  : arc sine function-acos  : arc cosine function-atan  : arc tangent function-sqr   : square-sqrt  : square root-recip : reciprocal (1/current image)-abs   : absolute value-bin   : use (current image>0) to binarise-binv  : binarise and invert (binarisation and logical inversion)-fillh : fill holes in a binary mask (holes are internal - i.e. do not touch the edge of the FOV)-fillh26 : fill holes using 26 connectivity-index : replace each nonzero voxel with a unique (subject to wrapping) index number-grid <value> <spacing> : add a 3D grid of intensity <value> with grid spacing <spacing>-edge  : edge strength-tfce <H> <E> <connectivity>: enhance with TFCE, e.g. -tfce 2 0.5 6 (maybe change 6 to 26 for skeletons)-tfceS <H> <E> <connectivity> <X> <Y> <Z> <tfce_thresh>: show support area for voxel (X,Y,Z)-nan   : replace NaNs (improper numbers) with 0-nanm  : make NaN (improper number) mask with 1 for NaN voxels, 0 otherwise-rand  : add uniform noise (range 0:1)-randn : add Gaussian noise (mean=0 sigma=1)-inm <mean> :  (-i i ip.c) intensity normalisation (per 3D volume mean)-ing <mean> :  (-I i ip.c) intensity normalisation, global 4D mean)-range : set the output calmin/max to full data rangeMatrix operations:-tensor_decomp : convert a 4D (6-timepoint )tensor image into L1,2,3,FA,MD,MO,V1,2,3 (remaining image in pipeline is FA)Kernel operations (set BEFORE filtering operation if desired):-kernel 3D : 3x3x3 box centered on target voxel (set as default kernel)-kernel 2D : 3x3x1 box centered on target voxel-kernel box    <size>     : all voxels in a cube of width <size> mm centered on target voxel-kernel boxv   <size>     : all voxels in a cube of width <size> voxels centered on target voxel, CAUTION: size should be an odd number-kernel boxv3  <X> <Y> <Z>: all voxels in a cuboid of dimensions X x Y x Z centered on target voxel, CAUTION: size should be an odd number-kernel gauss  <sigma>    : gaussian kernel (sigma in mm, not voxels)-kernel sphere <size>     : all voxels in a sphere of radius <size> mm centered on target voxel-kernel file   <filename> : use external file as kernelSpatial Filtering operations: N.B. all options apart from -s use the default kernel or that _previously_ specified by -kernel-dilM    : Mean Dilation of non-zero voxels-dilD    : Modal Dilation of non-zero voxels-dilF    : Maximum filtering of all voxels-dilall  : Apply -dilM repeatedly until the entire FOV is covered-ero     : Erode by zeroing non-zero voxels when zero voxels found in kernel-eroF    : Minimum filtering of all voxels-fmedian : Median Filtering -fmean   : Mean filtering, kernel weighted (conventionally used with gauss kernel)-fmeanu  : Mean filtering, kernel weighted, un-normalised (gives edge effects)-s <sigma> : create a gauss kernel of sigma mm and perform mean filtering-subsamp2  : downsamples image by a factor of 2 (keeping new voxels centred on old)-subsamp2offc  : downsamples image by a factor of 2 (non-centred)Dimensionality reduction operations:(the "T" can be replaced by X, Y or Z to collapse across a different dimension)-Tmean   : mean across time-Tstd    : standard deviation across time-Tmax    : max across time-Tmaxn   : time index of max across time-Tmin    : min across time-Tmedian : median across time-Tperc <percentage> : nth percentile (0-100) of FULL RANGE across time-Tar1    : temporal AR(1) coefficient (use -odt float and probably demean first)Basic statistical operations:-pval    : Nonparametric uncorrected P-value, assuming timepoints are the permutations; first timepoint is actual (unpermuted) stats image-pval0   : Same as -pval, but treat zeros as missing data-cpval   : Same as -pval, but gives FWE corrected P-values-ztop    : Convert Z-stat to (uncorrected) P-ptoz    : Convert (uncorrected) P to Z-rank    : Convert data to ranks (over T dim)-ranknorm: Transform to Normal dist via ranksMulti-argument operations:-roi <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize> : zero outside roi (using voxel coordinates). Inputting -1 for a size will set it to the full image extent for that dimension.-bptf  <hp_sigma> <lp_sigma> : (-t in ip.c) Bandpass temporal filtering; nonlinear highpass and Gaussian linear lowpass (with sigmas in volumes, not seconds); set either sigma<0 to skip that filter-roc <AROC-thresh> <outfile> [4Dnoiseonly] <truth> : take (normally binary) truth and test current image in ROC analysis against truth. <AROC-thresh> is usually 0.05 and is limit of Area-under-ROC measure FP axis. <outfile> is a text file of the ROC curve (triplets of values: FP TP threshold). If the truth image contains negative voxels these get excluded from all calculations. If <AROC-thresh> is positive then the [4Dnoiseonly] option needs to be set, and the FP rate is determined from this noise-only data, and is set to be the fraction of timepoints where any FP (anywhere) is seen, as found in the noise-only 4d-dataset. This is then controlling the FWE rate. If <AROC-thresh> is negative the FP rate is calculated from the zero-value parts of the <truth> image, this time averaging voxelwise FP rate over all timepoints. In both cases the TP rate is the average fraction of truth=positive voxels correctly found.Combining 4D and 3D images:If you apply a Binary operation (one that takes the current image and a new image together), when one is 3D and the other is 4D,the 3D image is cloned temporally to match the temporal dimensions of the 4D image.e.g. fslmaths inputVolume -add inputVolume2 output_volumefslmaths inputVolume -add 2.5 output_volumefslmaths inputVolume -add 2.5 -mul inputVolume2 output_volumefslmaths 4D_inputVolume -Tmean -mul -1 -add 4D_inputVolume demeaned_4D_inputVolume

2) randomise 非参数统计命令

randomise的详细使用帮助: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Randomise

randomisePart of FSL (ID: 5.0.10)
randomise v2.9Usage:
randomise -i <input> -o <output> -d <design.mat> -t <design.con> [options]Compulsory arguments (You MUST set one or more of):-i <input>  4D input image
    -o <out_root>   output file-rootname

Optional arguments (You may optionally specify one or more of):-D      demean data temporally before model fitting ( demean model as well if required )
    -1      perform 1-sample group-mean test instead of generic permutation test
    -m <mask>   mask image
    -d <design.mat> design matrix file
    -t <design.con> t contrasts file
    -f <design.fts> f contrasts file
    -e <design.grp> exchangeability block labels file
    --effective_design <design2.mat>    alternative design for determining valid permutations
    -q      print out how many unique permutations would be generated and exit
    -Q      print out information required for parallel mode and exit
    -n <n_perm> number of permutations (default 5000, set to 0 for exhaustive)
    -x      output voxelwise corrected p-value images
    --fonly     calculate f-statistics only
    -T      carry out Threshold-Free Cluster Enhancement
    --T2        carry out Threshold-Free Cluster Enhancement with 2D optimisation (e.g. for TBSS data); H=2, E=1, C=26
    -c <thresh> carry out cluster-based thresholding
    -C <thresh> carry out cluster-mass-based thresholding
    -F <thresh> carry out f cluster thresholding
    -S <thresh> carry out f cluster-mass thresholding
    -v <std>    use variance smoothing for t-stats (std is in mm)
    -h,--help   display this message
    --quiet     switch off diagnostic messages
    --twopass   carry out cluster normalisation thresholding
    -R      output raw ( unpermuted ) statistic images
    --uncorrp   output uncorrected p-value images
    -P      output permutation vector text file
    -N      output null distribution text files
    --norcmask  don't remove constant voxels from mask
    --seed <seed>   specific integer seed for random number generator
    --tfce_H <H>    TFCE height parameter (default=2)
    --tfce_D <H>    TFCE delta parameter overide
    --tfce_E <E>    TFCE extent parameter (default=0.5)
    --tfce_C <C>    TFCE connectivity (6 or 26; default=6)
    --vxl       list of numbers indicating voxelwise EVs position in the design matrix (list order corresponds to files in vxf option). caution BETA option.
    --vxf       list of 4D images containing voxelwise EVs (list order corresponds to numbers in vxl option). caution BETA option.
    --permuteBlocks permute exchangeability blocks. Caution BETA option
    --glm_output    output glm information for t-statistics ( unpermuted case only )
    --film  output stats to simulate the output of film

磁共振影像分析之: 基于FSL的VBM分析(2)相关推荐

【文本分析】基于公众需求文本分析的深圳自然博物馆发展策略研究
此文是2021年暑期在深圳大学文化产业研究院黄老师的指导下,与文化艺术管理专业的龚学姐合作完成的,且于2022年发表于<文化深圳>杂志第一期.此文偏向报告形式,没有遵循严格的学术规范,今年 ...
matlab分析数字信号,基于matlab的数字信号分析.doc
基于matlab的数字信号分析.doc 题目二.已知有一频率为50赫兹幅度为1V的方波信号被频率为1600赫兹幅度为2V的正弦波干扰,请设计方案将方波信号提取出来,要求被衰减的信号要衰减60DB以上, ...
java 情感分析_Flink基于Alink中文情感分析示例（Java版本）
Alink是基于Flink的机器学习算法平台,欢迎访问Alink的github获取下载链接及更多信息. alibaba/Alinkgithub.com 情感分析是对带有情感色彩(褒义贬义/正向负向)的 ...
论文浅尝 - ICLR2020 | 知道什么、如何以及为什么：基于方面的情感分析的近乎完整的解决方案...
论文笔记整理:余海阳,浙江大学硕士,研究方向为知识图谱.自然语言处理. 链接:https://arxiv.org/abs/1911.01616 动机基于目标的情感分析或基于方面的情感分析(ABSA) ...
【文本分析】基于粤港澳大湾区博物馆访客评价的文本分析
此文是2021年初在深圳大学文化产业研究院黄老师的指导下,与吴同学合作完成的,且于2021年4月发表于<文化深圳>杂志第一期. 文章目录 1 背景 2 研究思路 3 数据来源 4 数据基本 ...
基于matlab的频率特性测试仪,基于Matlab控制系统频率特性分析法
基于Matlab控制系统频率特性分析法基于Matlab控制系统频率特性分析法本文主要介绍了基于Matlab控制系统的频率特性分析方法.频域稳定性判据以及开环频域性能分析,并获得频率响应曲线等.通过 ...
文献记录(part79)--光学影像序列中基于多视角聚类的群组行为分析
学习笔记,仅供参考,有错必纠光学影像序列中基于多视角聚类的群组行为分析摘要群组行为分析是光学影像序列分析中的一项重要课题, 在近年来引起了人工智能领域研究人员的广泛关注. 与行人个体相比, 群组 ...
使用 FSL 和 TrackVis 分析 DTI 数据
转载原文使用 FSL 和 TrackVis 分析 DTI 数据 Alex / 2018-05-21 / free-learner@163.com 弥散加权成像(Diffusion Weighted ...
基于联合独立成分分析(jICA)的EEG与fMRI耦合与解耦
文章来源于微信公众号(茗创科技),欢迎有兴趣的朋友搜索关注. 导读 fMRI和EEG的有意义整合需要了解这些测量是否反映了相同神经源的活动,即估计神经成像模式之间的耦合与解耦程度.本文提出了一种基于联 ...
ArcGIS土地利用数据制备、分析及基于FLUS模型土地利用预测实践技术
土地资源对人类至关重要,土地是人类赖以生存和发展的物质基础,是一切生产和一起存在的源泉.利用现代化的技术手段及时.准确地获取土地利用现状,以及充分认识土地利用和土地覆盖变化规律,能极大地提高制定土地利 ...

磁共振影像分析之: 基于FSL的VBM分析(2)