【初心者向け】基礎&実践プログラミング

初心者がつまづきやすいところ、最短で実力が身につく方法をお伝えします。

【MRI画像解析】ラット・マウスデータに対する前処理から定量値計算まで

目的

  • ラット・マウスのMRIデータを前処理
  • 拡散MRIから得られる定量値計算
  • Atlasを用いた脳領域ごとの定量値算出

予備知識

  • Diffusion Kurtosis Imaging(DKI)の基礎と技術開発

www.innervision.co.jp

準備

準備されたデータは次の通り。

$ 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. 前処理

前処理の内容は次の通り。

  1. b valuesの丸め込み
  2. b valueの値が低い順に並ぶように画像ボリュームを入れ替える
  3. Eddy currentとボリューム間の位置ずれを補正 (eddy_correct)
  4. DTI, DKI, NODDIの準備
  5. 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から作成した定量値マップは次の通り。

f:id:AIProgrammer:20201104161954p:plain