目的
- ラット・マウスのMRIデータを前処理
- 拡散MRIから得られる定量値計算
- Atlasを用いた脳領域ごとの定量値算出
予備知識
- Diffusion Kurtosis Imaging（DKI）の基礎と技術開発
準備
準備されたデータは次の通り。
$ tree . . ├── HRT2.nii # 3D-T2WI ├── MToff.nii # MT saturation imaging (OFF) ├── MTon.nii # MT saturation imaging (ON) ├── b1000.bval # b-value ├── b1000.bvec # b-vector ├── b2000.bval # b-value ├── b2000.bvec # b-vector ├── b4000.bval # b-value ├── b4000.bvec # b-vector ├── dwi1000.nii # DWI ├── dwi2000.nii # DWI └── dwi4000.nii # DWI
1. 前処理
前処理の内容は次の通り。
- b valuesの丸め込み
- b valueの値が低い順に並ぶように画像ボリュームを入れ替える
- Eddy currentとボリューム間の位置ずれを補正 (eddy_correct)
- DTI, DKI, NODDIの準備
- Atlas位置合わせに必要なファイルの準備
次のファイルを用意する。
$ tree . . ├── HRT2.nii # 3D-T2WI ├── b1000.bval # b-value ├── b1000.bvec # b-vector ├── b2000.bval # b-value ├── b2000.bvec # b-vector ├── b4000.bval # b-value ├── b4000.bvec # b-vector ├── dwi1000.nii # DWI ├── dwi2000.nii # DWI └── dwi4000.nii # DWI
b valueの丸め込みには、次のスクリプトを用いる。
rat_mround.py
#!/usr/bin/env python # -*- coding: utf-8 -*- import sys, math, argparse __version__ = '1.0 (2019/01/21)' def mround(n, m): c, f = math.ceil(n / m) * m, math.floor(n / m) * m if abs(n - c) < abs(n - f): return float(c) else: return float(f) if __name__ == '__main__': parser = argparse.ArgumentParser( description='''round numbers to the desired multiple (default= 100). if no FILE is given, read standard input.''', formatter_class=argparse.RawDescriptionHelpFormatter, epilog='''examples: $ echo "98 99 102 790 805 1198" | mround.py 100 100 100 800 800 1200 $ cat x.txt 98 99 102 790 805 1198 $ mround.py x.txt 100 100 100 800 800 1200''') parser.add_argument('file', metavar='FILE', nargs='?', type=argparse.FileType('r'), default=sys.stdin, help='read in data from FILE') parser.add_argument('-m', '--mult', type=float, default=100.0, help='multiple (default= 100)') parser.add_argument('-V', '--version', action='version', version=__version__) err = 0 try: args = parser.parse_args() for line in args.file: results = [mround(n, args.mult) for n in map(float, line.split())] print(' '.join(map(str, results))) except Exception as e: print("mround.py: error: %s" % str(e)) err = 1 sys.exit(err)
データを準備したディレクトリで以下のコードを実行。
rat_preprocessing.sh
#!/bin/bash set -u << commentout Execute this program in the rawdata folder. Through PATH for "rat_mround.py" #output files b100020004000.bval #fixed bval(0 or 1000 or 2000 or 4000) b100020004000.bvec #fixed bvec ecb100020004000.nii #image after eddy_corrected DKE.nii # for DKE commentout gzip *nii raw_img=(`ls *1000.nii.gz` `ls *2000.nii.gz` `ls *4000.nii.gz`) bvals=("b1000.bval" "b2000.bval" "b4000.bval") bvecs=("b1000.bvec" "b2000.bvec" "b4000.bvec") # fix bvals (float > int) echo "fix bvals" ## bval for b in ${bvals[0]} ${bvals[1]} ${bvals[2]};do cat $b | rat_mround.py -m 1000 | cut -d " " -f -2 > fix_b0_${b} #100の位で丸めb=0をまとめる cat $b | rat_mround.py -m 1000 | cut -d " " -f 3- > fix_non_b0_${b} #100の位で丸めb=0以外をまとめる done paste -d " " fix_b0_${bvals[0]} fix_b0_${bvals[1]} fix_b0_${bvals[2]} fix_non_b0_${bvals[0]} fix_non_b0_${bvals[1]} fix_non_b0_${bvals[2]} > b100020004000.bval ## bvec for b in ${bvecs[0]} ${bvecs[1]} ${bvecs[2]};do cat $b | cut -d " " -f -2 > fix_b0_${b} cat $b | cut -d " " -f 3- > fix_non_b0_${b} done paste -d " " fix_b0_${bvecs[0]} fix_b0_${bvecs[1]} fix_b0_${bvecs[2]} fix_non_b0_${bvecs[0]} fix_non_b0_${bvecs[1]} fix_non_b0_${bvecs[2]} > b100020004000.bvec rm fix_* # split b0 volume from each bval image and combine these echo "split, combine and make b100020004000.nii.gz" ## split ### b0 for i in `seq 0 2`;do fslroi ${raw_img[i]} b0_${raw_img[i]} 0 2 done fslmerge -t only_b0 b0_${raw_img[0]} b0_${raw_img[1]} b0_${raw_img[2]} ### b1000 fslroi ${raw_img[0]} only_${raw_img[0]} 2 30 ### b2000 fslroi ${raw_img[1]} only_${raw_img[1]} 2 60 ### b4000 fslroi ${raw_img[2]} only_${raw_img[2]} 2 60 ## merge fslmerge -t b100020004000 only_b0 only_${raw_img[0]} only_${raw_img[1]} only_${raw_img[2]} # eddy_correct echo "eddy_correct" eddy_correct b100020004000 ecb100020004000 0 # make img for DKE echo "making img for DKE" ## mean b0 volumes fslmaths only_b0 -Tmean mean_only_b0 ## combine each bval volume fslmerge -t DKE mean_only_b0 only_${raw_img[0]} only_${raw_img[1]} only_${raw_img[2]} # convert .nii.gz into nii gunzip DKE.nii.gz ecb100020004000.nii.gz rm only_* b0_${raw_img[0]}.gz b0_${raw_img[1]}.gz b0_${raw_img[2]}.gz # preparing for AMICO and DKE echo "preparing for AMICO..." mkdir -p ../AMICO/ecb100020004000 cp ecb100020004000.nii b100020004000.bval b100020004000.bvec ../AMICO/ecb100020004000 mv ../AMICO/ecb100020004000/b100020004000.bval ../AMICO/ecb100020004000/ecb100020004000.bval mv ../AMICO/ecb100020004000/b100020004000.bvec ../AMICO/ecb100020004000/ecb100020004000.bvec echo "preparing for DKE..." mkdir -p ../DKE/ecb100020004000 cp DKE.nii b100020004000.bval b100020004000.bvec ../DKE/ecb100020004000 echo " preparing for reg..." #rename cp on.nii MTon.nii cp off.nii MToff.nii cp dwi1000.nii b1000.nii cp dwi2000.nii b2000.nii cp dwi4000.nii b4000.nii mkdir -p ../reg/ue mkdir ../reg/shita cp HRT2.nii ../reg for k in ue shita;do cp b*.nii MTon.nii MToff.nii b100020004000.bval b100020004000.bvec ../reg/$k mv ../reg/$k/b100020004000.bval ../reg/$k/ecb100020004000.bval mv ../reg/$k/b100020004000.bvec ../reg/$k/ecb100020004000.bvec done echo "-----------succeeded... (˘ω˘)-----------" exit
以上の処理が終了すると次のようなディレクトリ構造になる。
$ tree .
.
├── AMICO
│ └── ecb100020004000
│ ├── ecb100020004000.bval
│ ├── ecb100020004000.bvec
│ └── ecb100020004000.nii
├── DKE
│ └── ecb100020004000
│ ├── DKE.nii
│ ├── b100020004000.bval
│ └── b100020004000.bvec
├── preprocessing
│ ├── DKE.nii
│ ├── HRT2.nii
│ ├── MToff.nii
│ ├── MTon.nii
│ ├── b1000.bval
│ ├── b1000.bvec
│ ├── b1000.nii
│ ├── b100020004000.bval
│ ├── b100020004000.bvec
│ ├── b100020004000.nii.gz
│ ├── b2000.bval
│ ├── b2000.bvec
│ ├── b2000.nii
│ ├── b4000.bval
│ ├── b4000.bvec
│ ├── b4000.nii
│ ├── dwi1000.nii
│ ├── dwi2000.nii
│ ├── dwi4000.nii
│ ├── ecb100020004000.ecclog
│ ├── ecb100020004000.nii
│ └── mean_only_b0.nii.gz
└── reg
├── HRT2.nii
├── shita
│ ├── MToff.nii
│ ├── MTon.nii
│ ├── b1000.nii
│ ├── b2000.nii
│ ├── b4000.nii
│ ├── ecb100020004000.bval
│ └── ecb100020004000.bvec
└── ue
├── MToff.nii
├── MTon.nii
├── b1000.nii
├── b2000.nii
├── b4000.nii
├── ecb100020004000.bval
└── ecb100020004000.bvec
8 directories, 43 files
2. DTIおよびNODDI
必要なファイルは次の通り。
$ tree AMICO/
AMICO/
└── ecb100020004000
├── ecb100020004000.bval
├── ecb100020004000.bvec
└── ecb100020004000.nii
MATLABを用いてDTIとNODDIの定量値を計算。 NODDIに関しては、AMICO (Accelerated Microstructure Imaging via Convex Optimization) を用いて計算。
3. DKI
DKIの計算には、DKE (Diffusional Kurtosis Estimator)を用いる。
DKEの実行のために、次のファイルを準備する。
$ tree DKE
DKE
└── ecb100020004000
├── DKE.nii # DWI for DKE
├── b1000_30dir.txt # b-vector for b=1000
├── b2000_60dir.txt # b-vector for b=2000
├── b4000_60dir.txt # b-vector for b=4000
├── DKEParameters.dat # config for DKE
└── run_DKE.bat # batch file for running DKE
DKEParameters.dat
DKEを実行する前に、種々のパラメータを「DKEParameters.dat」に記載する。
studydirに、解析するディレクトリを指定preprocess_options.fn_niiに、解析する画像名を指定preprocess_options.formatに、解析する画像のフォーマットを指定bvalにb-valueを指定ndirに、書くb値の軸数を指定fwhm_img = 1.25*[0.15 0.15 0.15];に関しては、[ ]にDWIの解像度を指定する。
% Wed Mar 02 05:19:26 PM studydir = '.'; subject_list = {''}; preprocess_options.format = 'nifti'; preprocess_options.fn_nii = 'DKE.nii'; fn_img_prefix = 'rdki'; bval = [0 1000 2000 4000]; ndir = [30 60 60]; idx_1st_img = 1; Kmin = 0; NKmax = 3; Kmin_final = 0; Kmax_final = 10; T = 0.1; find_brain_mask_flag = 1; dki_method.no_tensor = 0; dki_method.linear_weighting = 1; dki_method.linear_constrained = 1; dki_method.nonlinear = 0; dki_method.linear_violations = 0; dki_method.robust_option = 0; dki_method.noise_tolerance = 0.09; dti_method.dti_flag = 0; dti_method.dti_only = 0; dti_method.no_tensor = 0; dti_method.linear_weighting = 1; dti_method.b_value = 1e+003; dti_method.directions = {1:30}; dti_method.robust_option = 0; dti_method.noise_tolerance = 0.09; fn_noise = ''; fwhm_img = 1.25*[0.15 0.15 0.15]; fwhm_noise = [0 0 0]; median_filter_method = 0; map_interpolation_method.flag = 0; map_interpolation_method.order = 1; map_interpolation_method.resolution =1; fn_gradients{1} = ['b1000_30dir.txt']; fn_gradients{2} = ['b2000_60dir.txt']; fn_gradients{3} = ['b4000_60dir.txt']; idx_gradients{1} = [1:30]; idx_gradients{2} = [1:60]; idx_gradients{3} = [1:60];
b-vector files
各b値に対するb-vectorの情報をテキストファイル(.txt)に記載する。以下は、「b1000_30dir.txt」のvector情報で「行がdirection」、「列がx, y, z」に対応している。
b1000_30dir.txt
-0.05461659 0.010584959 -0.229761969 0.060003118 0.042641105 -0.224648475 0.027946284 -0.074149337 -0.222725069 -0.079490414 -0.089615515 -0.203803604 -0.057520458 0.109125676 -0.201664472 0.044170229 0.139016049 -0.186035224 0.124634561 -0.0398343 -0.19688837 0.109647411 -0.139868058 -0.155884403 -0.003084178 -0.168457552 -0.165825561 -0.16729917 -0.051926771 -0.158745626 -0.155822964 0.052347851 -0.169895503 -0.024147971 0.199569334 -0.124396672 0.132223386 0.141829461 -0.13522927 0.171615852 0.050897423 -0.154411876 0.197873053 -0.055053012 -0.117051255 -0.01566304 -0.222914558 -0.077131994 -0.119744904 -0.158047999 -0.128715047 -0.223046732 -0.011430705 -0.077492386 -0.202044884 0.10069729 -0.070167629 -0.124438131 0.155175916 -0.127754304 0.070886614 0.214646278 -0.069192765 0.172450697 0.156826558 -0.039392424 0.221773326 0.062258868 -0.053160142 0.185578655 -0.13744913 -0.050535877 0.092913542 -0.203748845 -0.075756566 -0.10940962 -0.206740857 -0.034253882 -0.198599614 -0.114511209 -0.057714398 -0.232647192 0.03795272 0.017901128 -0.13445886 0.192040139 -0.030446267 -0.030555619 0.233444756 -0.02134187
run_DKE.bat
必要な設定とファイルが用意できたら、次のバッチファイルを実行してDKEを実行する。
@echo off cd /d "%~dp0" dke DKEParameters.dat pause
4. 定量値マップ算出結果
拡散MRIから作成した定量値マップは次の通り。
