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

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

【脳MRI画像解析ケーススタディ】施設間差の評価 ~拡散定量値の観点から~

f:id:AIProgrammer:20210224192531p:plain

動かしながら学ぶ PyTorchプログラミング入門

動かしながら学ぶ PyTorchプログラミング入門

  • 作者:斎藤勇哉
  • 発売日: 2020/11/30
  • メディア: Kindle版

目次

目的

  • 多施設MRIデータに含まれる施設間差(MRI装置・撮像条件)の評価
  • 施設間差補正後の評価
  • 拡散定量値を用いて、施設間差を評価

背景

多施設で取得されたMRIデータは、施設ごとに異なる撮像装置および撮像条件で取得されたものであり、施設間差(site effect)を含む。そのため、この施設間差と本来観測したい個人間の差(ex. 疾患・成長・老化による影響)を分離することができず、単に多施設データをするだけでは、感度・特異度に劣る。特に、精神疾患のような脳神経変性がわずかな場合、施設間差による変動に埋もれて、かすかな病変を検出できない。そのため、多施設共同研究において施設間差を除去することは、重要な課題である。

参考資料: mriportal.umin.jp

準備

拡散MRIから計算可能な定量値を用いて、施設間差を評価する。必要なデータは、以下。

  • DWI : 拡散MRI画像
  • DWI mask : 脳領域のバイナリ―マスク画像
  • b-values : b値
  • b-vector : bベクトル

以上の被験者データを、2施設(SiteA, SiteB)で撮像したとする(i.e., 標本バイアスはゼロで、測定バイアスのみを含むデータセット)。施設間差を含むSiteBをSiteB_before、施設間差を除去してSiteAに近づけたSiteBをSiteB_afterとする。

解析手法

1. JHU atlasをSiteAデータに位置合わせ

SiteAにおける各個人脳に、JHU atlasを位置合わせする。

具体的には、標準空間のFAと個人空間のFAの線形・非線形位置合わせをして、変換行列を取得し、標準空間にあるJHU atlasに適応して、JHU atlasを個人空間に移動させる。

準備

次のSiteAデータを使用する。

 .
 ├── drdwi.nii.gz  # Preprocessed DWI
 ├── drdwi_mask.nii.gz  # DWI mask
 └── fa.nii.gz  # FA map

ソースコード

以上のファイルが準備できたら、次のソースコードをShellScriptで実行する。

実行内容は、次の通り。

  1. DWI maskでFAをマスキング
  2. 標準空間のFAと個人空間のFAの線形・非線形位置合わせ
  3. 標準空間ー個人空間の変換行列を取得
  4. 標準空間にあるJHU atlasに変換行列を適応
  5. JHU atlasを個人空間に移動
# Register JHU atlas to subject FA space
# .
# ├── drdwi.nii.gz  # Preprocessed DWI
# ├── drdwi_mask.nii.gz  # DWI mask
# └── fa.nii.gz  # FA map

function masking_map() {
    MASK=drdwi_mask
    OUT_DIR=result
    mkdir -p ${OUT_DIR}
    fslmaths fa -mas ${MASK} ${OUT_DIR}/fa
}

function reg_FMRIB-Atlas_to_FA () {
    OUT_DIR=result
    # mkdir -p ${OUT_DIR}
    flirt -in ${OUT_DIR}/fa -ref ${FSLDIR}/data/standard/FMRIB58_FA_1mm -omat ${OUT_DIR}/FA_to_FMRIB.mat
    fnirt --in=fa --aff=${OUT_DIR}/FA_to_FMRIB.mat --config=FA_2_FMRIB58_1mm --cout=${OUT_DIR}/warp_FA_to_FMRIB
    invwarp -w ${OUT_DIR}/warp_FA_to_FMRIB -o ${OUT_DIR}/warp_FMRIB_to_FA -r fa

    applywarp --in=${FSLDIR}/data/atlases/JHU/JHU-ICBM-labels-1mm --ref=fa --warp=${OUT_DIR}/warp_FMRIB_to_FA --interp=nn --out=${OUT_DIR}/JHU-ICBM-labels_to_FA
    applywarp --in=${FSLDIR}/data/atlases/JHU/JHU-ICBM-tracts-maxprob-thr25-1mm --ref=fa --warp=${OUT_DIR}/warp_FMRIB_to_FA --interp=nn --out=${OUT_DIR}/JHU-ICBM-tracts-maxprob-thr25_to_FA
}

function split_JHU-ICBM-labels () {
    INPUT=result/JHU-ICBM-labels_to_FA
    OUT_DIR=result/mask/JHU-ICBM-labels
    mkdir -p ${OUT_DIR}

    fslmaths ${INPUT} -thr  1 -uthr  1 ${OUT_DIR}/01_Middle-cerebellar-peduncle
    fslmaths ${INPUT} -thr  2 -uthr  2 ${OUT_DIR}/02_Pontine-crossing-tract
    fslmaths ${INPUT} -thr  3 -uthr  3 ${OUT_DIR}/03_Genu-of-corpus-callosum
    fslmaths ${INPUT} -thr  4 -uthr  4 ${OUT_DIR}/04_Body-of-corpus-callosum
    fslmaths ${INPUT} -thr  5 -uthr  5 ${OUT_DIR}/05_Splenium-of-corpus-callosum
    fslmaths ${INPUT} -thr  6 -uthr  6 ${OUT_DIR}/06_Fornix
    fslmaths ${INPUT} -thr  7 -uthr  7 ${OUT_DIR}/07_Corticospinal-tract_R
    fslmaths ${INPUT} -thr  8 -uthr  8 ${OUT_DIR}/08_Corticospinal-tract_L
    fslmaths ${INPUT} -thr  9 -uthr  9 ${OUT_DIR}/09_Medial-lemniscus_R
    fslmaths ${INPUT} -thr 10 -uthr 10 ${OUT_DIR}/10_Medial-lemniscus_L
    fslmaths ${INPUT} -thr 11 -uthr 11 ${OUT_DIR}/11_Inferior-cerebellar-peduncle_R
    fslmaths ${INPUT} -thr 12 -uthr 12 ${OUT_DIR}/12_Inferior-cerebellar-peduncle_L
    fslmaths ${INPUT} -thr 13 -uthr 13 ${OUT_DIR}/13_Superior-cerebellar-peduncle_R
    fslmaths ${INPUT} -thr 14 -uthr 14 ${OUT_DIR}/14_Superior-cerebellar-peduncle_L
    fslmaths ${INPUT} -thr 15 -uthr 15 ${OUT_DIR}/15_Cerebral-peduncle_R
    fslmaths ${INPUT} -thr 16 -uthr 16 ${OUT_DIR}/16_Cerebral-peduncle_L
    fslmaths ${INPUT} -thr 17 -uthr 17 ${OUT_DIR}/17_Anterior-limb-of-internal-capsule_R
    fslmaths ${INPUT} -thr 18 -uthr 18 ${OUT_DIR}/18_Anterior-limb-of-internal-capsule_L
    fslmaths ${INPUT} -thr 19 -uthr 19 ${OUT_DIR}/19_Posterior-limb-of-internal-capsule_R
    fslmaths ${INPUT} -thr 20 -uthr 20 ${OUT_DIR}/20_Posterior-limb-of-internal-capsule_L
    fslmaths ${INPUT} -thr 21 -uthr 21 ${OUT_DIR}/21_Retrolenticular-part-of-internal-capsule_R
    fslmaths ${INPUT} -thr 22 -uthr 22 ${OUT_DIR}/22_Retrolenticular-part-of-internal-capsule_L
    fslmaths ${INPUT} -thr 23 -uthr 23 ${OUT_DIR}/23_Anterior-corona-radiata_R
    fslmaths ${INPUT} -thr 24 -uthr 24 ${OUT_DIR}/24_Anterior-corona-radiata_L
    fslmaths ${INPUT} -thr 25 -uthr 25 ${OUT_DIR}/25_Superior-corona-radiata_R
    fslmaths ${INPUT} -thr 26 -uthr 26 ${OUT_DIR}/26_Superior-corona-radiata_L
    fslmaths ${INPUT} -thr 27 -uthr 27 ${OUT_DIR}/27_Posterior-corona-radiata_R
    fslmaths ${INPUT} -thr 28 -uthr 28 ${OUT_DIR}/28_Posterior-corona-radiata_L
    fslmaths ${INPUT} -thr 29 -uthr 29 ${OUT_DIR}/29_Posterior-thalamic-radiation_R
    fslmaths ${INPUT} -thr 30 -uthr 30 ${OUT_DIR}/30_Posterior-thalamic-radiation_L
    fslmaths ${INPUT} -thr 31 -uthr 31 ${OUT_DIR}/31_Sagittal-stratum_R
    fslmaths ${INPUT} -thr 32 -uthr 32 ${OUT_DIR}/32_Sagittal-stratum_L
    fslmaths ${INPUT} -thr 33 -uthr 33 ${OUT_DIR}/33_External-capsule_R
    fslmaths ${INPUT} -thr 34 -uthr 34 ${OUT_DIR}/34_External-capsule_L
    fslmaths ${INPUT} -thr 35 -uthr 35 ${OUT_DIR}/35_Cingulum-cingulate-gyrus_R
    fslmaths ${INPUT} -thr 36 -uthr 36 ${OUT_DIR}/36_Cingulum-cingulate-gyrus_L
    fslmaths ${INPUT} -thr 37 -uthr 37 ${OUT_DIR}/37_Cingulum-hippocampus_R
    fslmaths ${INPUT} -thr 38 -uthr 38 ${OUT_DIR}/38_Cingulum-hippocampus_L
    fslmaths ${INPUT} -thr 39 -uthr 39 ${OUT_DIR}/39_Fornix-Stria-terminalis_R
    fslmaths ${INPUT} -thr 40 -uthr 40 ${OUT_DIR}/40_Fornix-Stria-terminalis_L
    fslmaths ${INPUT} -thr 41 -uthr 41 ${OUT_DIR}/41_Superior-longitudinal-fasciculus_R
    fslmaths ${INPUT} -thr 42 -uthr 42 ${OUT_DIR}/42_Superior-longitudinal-fasciculus_L
    fslmaths ${INPUT} -thr 43 -uthr 43 ${OUT_DIR}/43_Superior-fronto-occipital-fasciculus_R
    fslmaths ${INPUT} -thr 44 -uthr 44 ${OUT_DIR}/44_Superior-fronto-occipital-fasciculus_L
    fslmaths ${INPUT} -thr 45 -uthr 45 ${OUT_DIR}/45_Uncinate-fasciculus_R
    fslmaths ${INPUT} -thr 46 -uthr 46 ${OUT_DIR}/46_Uncinate-fasciculus_L
    fslmaths ${INPUT} -thr 47 -uthr 47 ${OUT_DIR}/47_Tapetum_R
    fslmaths ${INPUT} -thr 48 -uthr 48 ${OUT_DIR}/48_Tapetum_L

    fsladd ${OUT_DIR}/51_Corticospinal-tract ${OUT_DIR}/07_Corticospinal-tract_R ${OUT_DIR}/08_Corticospinal-tract_L > /dev/null
    fsladd ${OUT_DIR}/52_Medial-lemniscus ${OUT_DIR}/09_Medial-lemniscus_R ${OUT_DIR}/10_Medial-lemniscus_L > /dev/null
    fsladd ${OUT_DIR}/53_Inferior-cerebellar-peduncle ${OUT_DIR}/11_Inferior-cerebellar-peduncle_R ${OUT_DIR}/12_Inferior-cerebellar-peduncle_L > /dev/null
    fsladd ${OUT_DIR}/54_Superior-cerebellar-peduncle ${OUT_DIR}/13_Superior-cerebellar-peduncle_R ${OUT_DIR}/14_Superior-cerebellar-peduncle_L > /dev/null
    fsladd ${OUT_DIR}/55_Cerebral-peduncle ${OUT_DIR}/15_Cerebral-peduncle_R ${OUT_DIR}/16_Cerebral-peduncle_L > /dev/null
    fsladd ${OUT_DIR}/56_Anterior-limb-of-internal-capsule ${OUT_DIR}/17_Anterior-limb-of-internal-capsule_R ${OUT_DIR}/18_Anterior-limb-of-internal-capsule_L > /dev/null
    fsladd ${OUT_DIR}/57_Posterior-limb-of-internal-capsule ${OUT_DIR}/19_Posterior-limb-of-internal-capsule_R ${OUT_DIR}/20_Posterior-limb-of-internal-capsule_L > /dev/null
    fsladd ${OUT_DIR}/58_Retrolenticular-part-of-internal-capsule ${OUT_DIR}/21_Retrolenticular-part-of-internal-capsule_R ${OUT_DIR}/22_Retrolenticular-part-of-internal-capsule_L > /dev/null
    fsladd ${OUT_DIR}/59_Anterior-corona-radiata ${OUT_DIR}/23_Anterior-corona-radiata_R ${OUT_DIR}/24_Anterior-corona-radiata_L > /dev/null
    fsladd ${OUT_DIR}/60_Superior-corona-radiata ${OUT_DIR}/25_Superior-corona-radiata_R ${OUT_DIR}/26_Superior-corona-radiata_L > /dev/null
    fsladd ${OUT_DIR}/61_Posterior-corona-radiata ${OUT_DIR}/27_Posterior-corona-radiata_R ${OUT_DIR}/28_Posterior-corona-radiata_L > /dev/null
    fsladd ${OUT_DIR}/62_Posterior-thalamic-radiation ${OUT_DIR}/29_Posterior-thalamic-radiation_R ${OUT_DIR}/30_Posterior-thalamic-radiation_L > /dev/null
    fsladd ${OUT_DIR}/63_Sagittal-stratum ${OUT_DIR}/31_Sagittal-stratum_R ${OUT_DIR}/32_Sagittal-stratum_L > /dev/null
    fsladd ${OUT_DIR}/64_External-capsule ${OUT_DIR}/33_External-capsule_R ${OUT_DIR}/34_External-capsule_L > /dev/null
    fsladd ${OUT_DIR}/65_Cingulum-cingulate-gyrus ${OUT_DIR}/35_Cingulum-cingulate-gyrus_R ${OUT_DIR}/36_Cingulum-cingulate-gyrus_L > /dev/null
    fsladd ${OUT_DIR}/66_Cingulum-hippocampus ${OUT_DIR}/37_Cingulum-hippocampus_R ${OUT_DIR}/38_Cingulum-hippocampus_L > /dev/null
    fsladd ${OUT_DIR}/67_Fornix-Stria-terminalis ${OUT_DIR}/39_Fornix-Stria-terminalis_R ${OUT_DIR}/40_Fornix-Stria-terminalis_L > /dev/null
    fsladd ${OUT_DIR}/68_Superior-longitudinal-fasciculus ${OUT_DIR}/41_Superior-longitudinal-fasciculus_R ${OUT_DIR}/42_Superior-longitudinal-fasciculus_L > /dev/null
    fsladd ${OUT_DIR}/69_Superior-fronto-occipital-fasciculus ${OUT_DIR}/43_Superior-fronto-occipital-fasciculus_R ${OUT_DIR}/44_Superior-fronto-occipital-fasciculus_L > /dev/null
    fsladd ${OUT_DIR}/70_Uncinate-fasciculus ${OUT_DIR}/45_Uncinate-fasciculus_R ${OUT_DIR}/46_Uncinate-fasciculus_L > /dev/null
    fsladd ${OUT_DIR}/71_Tapetum ${OUT_DIR}/47_Tapetum_R ${OUT_DIR}/48_Tapetum_L > /dev/null

    fsladd ${OUT_DIR}/81_Corpus-callosum ${OUT_DIR}/03_Genu-of-corpus-callosum ${OUT_DIR}/04_Body-of-corpus-callosum ${OUT_DIR}/05_Splenium-of-corpus-callosum > /dev/null
    fsladd ${OUT_DIR}/82_Internal-capsule ${OUT_DIR}/56_Anterior-limb-of-internal-capsule ${OUT_DIR}/57_Posterior-limb-of-internal-capsule ${OUT_DIR}/58_Retrolenticular-part-of-internal-capsule > /dev/null
    fsladd ${OUT_DIR}/83_Corona-radiata ${OUT_DIR}/59_Anterior-corona-radiata ${OUT_DIR}/60_Superior-corona-radiata ${OUT_DIR}/61_Posterior-corona-radiata > /dev/null
}

function split_JHU-ICBM-tracts-maxprob-thr25 () {
    INPUT=result/JHU-ICBM-tracts-maxprob-thr25_to_FA
    OUT_DIR=result/mask/JHU-ICBM-tracts-maxprob-thr25
    mkdir -p ${OUT_DIR}

    fslmaths ${INPUT} -thr  1 -uthr  1 ${OUT_DIR}/01_Anterior-thalamic_Radiation_L
    fslmaths ${INPUT} -thr  2 -uthr  2 ${OUT_DIR}/02_Anterior-thalamic_Radiation_R
    fslmaths ${INPUT} -thr  3 -uthr  3 ${OUT_DIR}/03_Corticospinal-tract_L
    fslmaths ${INPUT} -thr  4 -uthr  4 ${OUT_DIR}/04_Corticospinal-tract_R
    fslmaths ${INPUT} -thr  5 -uthr  5 ${OUT_DIR}/05_Cingulum-cingulate-gyrus_L
    fslmaths ${INPUT} -thr  6 -uthr  6 ${OUT_DIR}/06_Cingulum-cingulate-gyrus_R
    fslmaths ${INPUT} -thr  7 -uthr  7 ${OUT_DIR}/07_Cingulum-hippocampus_L
    fslmaths ${INPUT} -thr  8 -uthr  8 ${OUT_DIR}/08_Cingulum-hippocampus_R
    fslmaths ${INPUT} -thr  9 -uthr  9 ${OUT_DIR}/09_Forceps-major
    fslmaths ${INPUT} -thr 10 -uthr 10 ${OUT_DIR}/10_Forceps-minor
    fslmaths ${INPUT} -thr 11 -uthr 11 ${OUT_DIR}/11_Inferior-fronto-occipital-fasciculus_L
    fslmaths ${INPUT} -thr 12 -uthr 12 ${OUT_DIR}/12_Inferior-fronto-occipital-fasciculus_R
    fslmaths ${INPUT} -thr 13 -uthr 13 ${OUT_DIR}/13_Inferior_Longitudinal-fasciculus_L
    fslmaths ${INPUT} -thr 14 -uthr 14 ${OUT_DIR}/14_Inferior_Longitudinal-fasciculus_R
    fslmaths ${INPUT} -thr 15 -uthr 15 ${OUT_DIR}/15_Superior_Longitudinal-fasciculus_L
    fslmaths ${INPUT} -thr 16 -uthr 16 ${OUT_DIR}/16_Superior_Longitudinal-fasciculus_R
    fslmaths ${INPUT} -thr 17 -uthr 17 ${OUT_DIR}/17_Uncinate-fasciculus_L
    fslmaths ${INPUT} -thr 18 -uthr 18 ${OUT_DIR}/18_Uncinate-fasciculus_R
    fslmaths ${INPUT} -thr 19 -uthr 19 ${OUT_DIR}/19_Superior_Longitudinal-fasciculus-temporal-part_L
    fslmaths ${INPUT} -thr 20 -uthr 20 ${OUT_DIR}/20_Superior_Longitudinal-fasciculus-temporal-part_R

    fsladd ${OUT_DIR}/21_Anterior-thalamic_Radiation ${OUT_DIR}/01_Anterior-thalamic_Radiation_L ${OUT_DIR}/02_Anterior-thalamic_Radiation_R > /dev/null
    fsladd ${OUT_DIR}/22_Corticospinal-tract ${OUT_DIR}/03_Corticospinal-tract_L ${OUT_DIR}/04_Corticospinal-tract_R > /dev/null
    fsladd ${OUT_DIR}/23_Cingulum-cingulate-gyrus ${OUT_DIR}/05_Cingulum-cingulate-gyrus_L ${OUT_DIR}/06_Cingulum-cingulate-gyrus_R > /dev/null
    fsladd ${OUT_DIR}/24_Cingulum-hippocampus ${OUT_DIR}/07_Cingulum-hippocampus_L ${OUT_DIR}/08_Cingulum-hippocampus_R > /dev/null
    fsladd ${OUT_DIR}/25_Forceps ${OUT_DIR}/09_Forceps-major ${OUT_DIR}/10_Forceps-minor > /dev/null
    fsladd ${OUT_DIR}/26_Inferior-fronto-occipital-fasciculus ${OUT_DIR}/11_Inferior-fronto-occipital-fasciculus_L ${OUT_DIR}/12_Inferior-fronto-occipital-fasciculus_R > /dev/null
    fsladd ${OUT_DIR}/27_Inferior_Longitudinal-fasciculus ${OUT_DIR}/13_Inferior_Longitudinal-fasciculus_L ${OUT_DIR}/14_Inferior_Longitudinal-fasciculus_R > /dev/null
    fsladd ${OUT_DIR}/28_Superior_Longitudinal-fasciculus ${OUT_DIR}/15_Superior_Longitudinal-fasciculus_L ${OUT_DIR}/16_Superior_Longitudinal-fasciculus_R > /dev/null
    fsladd ${OUT_DIR}/29_Uncinate-fasciculus ${OUT_DIR}/17_Uncinate-fasciculus_L ${OUT_DIR}/18_Uncinate-fasciculus_R > /dev/null
    fsladd ${OUT_DIR}/30_Superior_Longitudinal-fasciculus-temporal-part ${OUT_DIR}/19_Superior_Longitudinal-fasciculus-temporal-part_L ${OUT_DIR}/20_Superior_Longitudinal-fasciculus-temporal-part_R > /dev/null
}

# 1. Mask FA map
masking_map
# 2. Register JHU atlas to subject FA space
reg_FMRIB-Atlas_to_FA
# 3. Split JHU atlas
split_JHU-ICBM-labels
split_JHU-ICBM-tracts-maxprob-thr25

結果

上記のコードを実行すると、次のような「result」フォルダが生成され、その「mask」フォルダに個人空間に位置合わせされたJHU atlasが格納される。

ディレクトリ構造は、次の通りになる。

.
├── drdwi.nii.gz
├── drdwi_mask.nii.gz
├── fa.nii.gz
├── fa_to_FMRIB58_FA_1mm.log
└── result
    ├── FA_to_FMRIB.mat
    ├── JHU-ICBM-labels_to_FA.nii.gz
    ├── JHU-ICBM-tracts-maxprob-thr25_to_FA.nii.gz
    ├── fa.nii.gz
    ├── mask
    │   ├── JHU-ICBM-labels
    │   │   ├── 01_Middle-cerebellar-peduncle.nii.gz
    │   │   ├── 02_Pontine-crossing-tract.nii.gz
    │   │   ├── 03_Genu-of-corpus-callosum.nii.gz
    │   │   ├── 04_Body-of-corpus-callosum.nii.gz
    │   │   ├── 05_Splenium-of-corpus-callosum.nii.gz
    │   │   ├── 06_Fornix.nii.gz
    │   │   ├── 07_Corticospinal-tract_R.nii.gz
    │   │   ├── 08_Corticospinal-tract_L.nii.gz
    │   │   ├── 09_Medial-lemniscus_R.nii.gz
    │   │   ├── 10_Medial-lemniscus_L.nii.gz
    │   │   ├── 11_Inferior-cerebellar-peduncle_R.nii.gz
    │   │   ├── 12_Inferior-cerebellar-peduncle_L.nii.gz
    │   │   ├── 13_Superior-cerebellar-peduncle_R.nii.gz
    │   │   ├── 14_Superior-cerebellar-peduncle_L.nii.gz
    │   │   ├── 15_Cerebral-peduncle_R.nii.gz
    │   │   ├── 16_Cerebral-peduncle_L.nii.gz
    │   │   ├── 17_Anterior-limb-of-internal-capsule_R.nii.gz
    │   │   ├── 18_Anterior-limb-of-internal-capsule_L.nii.gz
    │   │   ├── 19_Posterior-limb-of-internal-capsule_R.nii.gz
    │   │   ├── 20_Posterior-limb-of-internal-capsule_L.nii.gz
    │   │   ├── 21_Retrolenticular-part-of-internal-capsule_R.nii.gz
    │   │   ├── 22_Retrolenticular-part-of-internal-capsule_L.nii.gz
    │   │   ├── 23_Anterior-corona-radiata_R.nii.gz
    │   │   ├── 24_Anterior-corona-radiata_L.nii.gz
    │   │   ├── 25_Superior-corona-radiata_R.nii.gz
    │   │   ├── 26_Superior-corona-radiata_L.nii.gz
    │   │   ├── 27_Posterior-corona-radiata_R.nii.gz
    │   │   ├── 28_Posterior-corona-radiata_L.nii.gz
    │   │   ├── 29_Posterior-thalamic-radiation_R.nii.gz
    │   │   ├── 30_Posterior-thalamic-radiation_L.nii.gz
    │   │   ├── 31_Sagittal-stratum_R.nii.gz
    │   │   ├── 32_Sagittal-stratum_L.nii.gz
    │   │   ├── 33_External-capsule_R.nii.gz
    │   │   ├── 34_External-capsule_L.nii.gz
    │   │   ├── 35_Cingulum-cingulate-gyrus_R.nii.gz
    │   │   ├── 36_Cingulum-cingulate-gyrus_L.nii.gz
    │   │   ├── 37_Cingulum-hippocampus_R.nii.gz
    │   │   ├── 38_Cingulum-hippocampus_L.nii.gz
    │   │   ├── 39_Fornix-Stria-terminalis_R.nii.gz
    │   │   ├── 40_Fornix-Stria-terminalis_L.nii.gz
    │   │   ├── 41_Superior-longitudinal-fasciculus_R.nii.gz
    │   │   ├── 42_Superior-longitudinal-fasciculus_L.nii.gz
    │   │   ├── 43_Superior-fronto-occipital-fasciculus_R.nii.gz
    │   │   ├── 44_Superior-fronto-occipital-fasciculus_L.nii.gz
    │   │   ├── 45_Uncinate-fasciculus_R.nii.gz
    │   │   ├── 46_Uncinate-fasciculus_L.nii.gz
    │   │   ├── 47_Tapetum_R.nii.gz
    │   │   ├── 48_Tapetum_L.nii.gz
    │   │   ├── 51_Corticospinal-tract.nii.gz
    │   │   ├── 52_Medial-lemniscus.nii.gz
    │   │   ├── 53_Inferior-cerebellar-peduncle.nii.gz
    │   │   ├── 54_Superior-cerebellar-peduncle.nii.gz
    │   │   ├── 55_Cerebral-peduncle.nii.gz
    │   │   ├── 56_Anterior-limb-of-internal-capsule.nii.gz
    │   │   ├── 57_Posterior-limb-of-internal-capsule.nii.gz
    │   │   ├── 58_Retrolenticular-part-of-internal-capsule.nii.gz
    │   │   ├── 59_Anterior-corona-radiata.nii.gz
    │   │   ├── 60_Superior-corona-radiata.nii.gz
    │   │   ├── 61_Posterior-corona-radiata.nii.gz
    │   │   ├── 62_Posterior-thalamic-radiation.nii.gz
    │   │   ├── 63_Sagittal-stratum.nii.gz
    │   │   ├── 64_External-capsule.nii.gz
    │   │   ├── 65_Cingulum-cingulate-gyrus.nii.gz
    │   │   ├── 66_Cingulum-hippocampus.nii.gz
    │   │   ├── 67_Fornix-Stria-terminalis.nii.gz
    │   │   ├── 68_Superior-longitudinal-fasciculus.nii.gz
    │   │   ├── 69_Superior-fronto-occipital-fasciculus.nii.gz
    │   │   ├── 70_Uncinate-fasciculus.nii.gz
    │   │   ├── 71_Tapetum.nii.gz
    │   │   ├── 81_Corpus-callosum.nii.gz
    │   │   ├── 82_Internal-capsule.nii.gz
    │   │   └── 83_Corona-radiata.nii.gz
    │   └── JHU-ICBM-tracts-maxprob-thr25
    │       ├── 01_Anterior-thalamic_Radiation_L.nii.gz
    │       ├── 02_Anterior-thalamic_Radiation_R.nii.gz
    │       ├── 03_Corticospinal-tract_L.nii.gz
    │       ├── 04_Corticospinal-tract_R.nii.gz
    │       ├── 05_Cingulum-cingulate-gyrus_L.nii.gz
    │       ├── 06_Cingulum-cingulate-gyrus_R.nii.gz
    │       ├── 07_Cingulum-hippocampus_L.nii.gz
    │       ├── 08_Cingulum-hippocampus_R.nii.gz
    │       ├── 09_Forceps-major.nii.gz
    │       ├── 10_Forceps-minor.nii.gz
    │       ├── 11_Inferior-fronto-occipital-fasciculus_L.nii.gz
    │       ├── 12_Inferior-fronto-occipital-fasciculus_R.nii.gz
    │       ├── 13_Inferior_Longitudinal-fasciculus_L.nii.gz
    │       ├── 14_Inferior_Longitudinal-fasciculus_R.nii.gz
    │       ├── 15_Superior_Longitudinal-fasciculus_L.nii.gz
    │       ├── 16_Superior_Longitudinal-fasciculus_R.nii.gz
    │       ├── 17_Uncinate-fasciculus_L.nii.gz
    │       ├── 18_Uncinate-fasciculus_R.nii.gz
    │       ├── 19_Superior_Longitudinal-fasciculus-temporal-part_L.nii.gz
    │       ├── 20_Superior_Longitudinal-fasciculus-temporal-part_R.nii.gz
    │       ├── 21_Anterior-thalamic_Radiation.nii.gz
    │       ├── 22_Corticospinal-tract.nii.gz
    │       ├── 23_Cingulum-cingulate-gyrus.nii.gz
    │       ├── 24_Cingulum-hippocampus.nii.gz
    │       ├── 25_Forceps.nii.gz
    │       ├── 26_Inferior-fronto-occipital-fasciculus.nii.gz
    │       ├── 27_Inferior_Longitudinal-fasciculus.nii.gz
    │       ├── 28_Superior_Longitudinal-fasciculus.nii.gz
    │       ├── 29_Uncinate-fasciculus.nii.gz
    │       └── 30_Superior_Longitudinal-fasciculus-temporal-part.nii.gz
    ├── warp_FA_to_FMRIB.nii.gz
    └── warp_FMRIB_to_FA.nii.gz

2. 絶対差マップの生成とJHU atlasを用いたサンプリング

MRIデータに含まれる施設間差を可視化するため、絶対差マップを生成する。視覚評価による定性的評価だけでなく定量評価もできるように、JHU atlasを用いて絶対差マップをサンプリングする。

準備

SiteA, SiteBのデータを次のように準備する。「data/mask」フォルダは、「1. JHU atlasをSiteAデータに位置合わせ」で生成した「result/mask」フォルダにあたる。

data/
├── SiteB_after
│   ├── drdwi.nii.gz  # Preprocessed DWI
│   ├── drdwi_mask.nii.gz  # DWI mask
│   └── map  # Diffusion map
│       ├── ad.nii.gz
│       ├── ak.nii.gz
│       ├── fa.nii.gz
│       ├── md.nii.gz
│       ├── mk.nii.gz
│       ├── rd.nii.gz
│       └── rk.nii.gz
├── SiteB_before
│   ├── drdwi.nii.gz
│   ├── drdwi_mask.nii.gz
│   └── map
│       ├── ad.nii.gz
│       ├── ak.nii.gz
│       ├── fa.nii.gz
│       ├── md.nii.gz
│       ├── mk.nii.gz
│       ├── rd.nii.gz
│       └── rk.nii.gz
├── SiteA
│   ├── drdwi.nii.gz
│   ├── drdwi_mask.nii.gz
│   └── map
│       ├── ad.nii.gz
│       ├── ak.nii.gz
│       ├── fa.nii.gz
│       ├── md.nii.gz
│       ├── mk.nii.gz
│       ├── rd.nii.gz
│       └── rk.nii.gz
└── mask       # In SiteA space
    ├── JHU-ICBM-labels
    │   ├── 01_Middle-cerebellar-peduncle.nii.gz
    │   ├── 02_Pontine-crossing-tract.nii.gz
    │   ├── 03_Genu-of-corpus-callosum.nii.gz
    │   ├── 04_Body-of-corpus-callosum.nii.gz
    │   ├── 05_Splenium-of-corpus-callosum.nii.gz
    │   ├── 06_Fornix.nii.gz
    │   ├── 07_Corticospinal-tract_R.nii.gz
    │   ├── 08_Corticospinal-tract_L.nii.gz
    │   ├── 09_Medial-lemniscus_R.nii.gz
    │   ├── 10_Medial-lemniscus_L.nii.gz
    │   ├── 11_Inferior-cerebellar-peduncle_R.nii.gz
    │   ├── 12_Inferior-cerebellar-peduncle_L.nii.gz
    │   ├── 13_Superior-cerebellar-peduncle_R.nii.gz
    │   ├── 14_Superior-cerebellar-peduncle_L.nii.gz
    │   ├── 15_Cerebral-peduncle_R.nii.gz
    │   ├── 16_Cerebral-peduncle_L.nii.gz
    │   ├── 17_Anterior-limb-of-internal-capsule_R.nii.gz
    │   ├── 18_Anterior-limb-of-internal-capsule_L.nii.gz
    │   ├── 19_Posterior-limb-of-internal-capsule_R.nii.gz
    │   ├── 20_Posterior-limb-of-internal-capsule_L.nii.gz
    │   ├── 21_Retrolenticular-part-of-internal-capsule_R.nii.gz
    │   ├── 22_Retrolenticular-part-of-internal-capsule_L.nii.gz
    │   ├── 23_Anterior-corona-radiata_R.nii.gz
    │   ├── 24_Anterior-corona-radiata_L.nii.gz
    │   ├── 25_Superior-corona-radiata_R.nii.gz
    │   ├── 26_Superior-corona-radiata_L.nii.gz
    │   ├── 27_Posterior-corona-radiata_R.nii.gz
    │   ├── 28_Posterior-corona-radiata_L.nii.gz
    │   ├── 29_Posterior-thalamic-radiation_R.nii.gz
    │   ├── 30_Posterior-thalamic-radiation_L.nii.gz
    │   ├── 31_Sagittal-stratum_R.nii.gz
    │   ├── 32_Sagittal-stratum_L.nii.gz
    │   ├── 33_External-capsule_R.nii.gz
    │   ├── 34_External-capsule_L.nii.gz
    │   ├── 35_Cingulum-cingulate-gyrus_R.nii.gz
    │   ├── 36_Cingulum-cingulate-gyrus_L.nii.gz
    │   ├── 37_Cingulum-hippocampus_R.nii.gz
    │   ├── 38_Cingulum-hippocampus_L.nii.gz
    │   ├── 39_Fornix-Stria-terminalis_R.nii.gz
    │   ├── 40_Fornix-Stria-terminalis_L.nii.gz
    │   ├── 41_Superior-longitudinal-fasciculus_R.nii.gz
    │   ├── 42_Superior-longitudinal-fasciculus_L.nii.gz
    │   ├── 43_Superior-fronto-occipital-fasciculus_R.nii.gz
    │   ├── 44_Superior-fronto-occipital-fasciculus_L.nii.gz
    │   ├── 45_Uncinate-fasciculus_R.nii.gz
    │   ├── 46_Uncinate-fasciculus_L.nii.gz
    │   ├── 47_Tapetum_R.nii.gz
    │   ├── 48_Tapetum_L.nii.gz
    │   ├── 51_Corticospinal-tract.nii.gz
    │   ├── 52_Medial-lemniscus.nii.gz
    │   ├── 53_Inferior-cerebellar-peduncle.nii.gz
    │   ├── 54_Superior-cerebellar-peduncle.nii.gz
    │   ├── 55_Cerebral-peduncle.nii.gz
    │   ├── 56_Anterior-limb-of-internal-capsule.nii.gz
    │   ├── 57_Posterior-limb-of-internal-capsule.nii.gz
    │   ├── 58_Retrolenticular-part-of-internal-capsule.nii.gz
    │   ├── 59_Anterior-corona-radiata.nii.gz
    │   ├── 60_Superior-corona-radiata.nii.gz
    │   ├── 61_Posterior-corona-radiata.nii.gz
    │   ├── 62_Posterior-thalamic-radiation.nii.gz
    │   ├── 63_Sagittal-stratum.nii.gz
    │   ├── 64_External-capsule.nii.gz
    │   ├── 65_Cingulum-cingulate-gyrus.nii.gz
    │   ├── 66_Cingulum-hippocampus.nii.gz
    │   ├── 67_Fornix-Stria-terminalis.nii.gz
    │   ├── 68_Superior-longitudinal-fasciculus.nii.gz
    │   ├── 69_Superior-fronto-occipital-fasciculus.nii.gz
    │   ├── 70_Uncinate-fasciculus.nii.gz
    │   ├── 71_Tapetum.nii.gz
    │   ├── 81_Corpus-callosum.nii.gz
    │   ├── 82_Internal-capsule.nii.gz
    │   └── 83_Corona-radiata.nii.gz
    └── JHU-ICBM-tracts-maxprob-thr25
        ├── 01_Anterior-thalamic_Radiation_L.nii.gz
        ├── 02_Anterior-thalamic_Radiation_R.nii.gz
        ├── 03_Corticospinal-tract_L.nii.gz
        ├── 04_Corticospinal-tract_R.nii.gz
        ├── 05_Cingulum-cingulate-gyrus_L.nii.gz
        ├── 06_Cingulum-cingulate-gyrus_R.nii.gz
        ├── 07_Cingulum-hippocampus_L.nii.gz
        ├── 08_Cingulum-hippocampus_R.nii.gz
        ├── 09_Forceps-major.nii.gz
        ├── 10_Forceps-minor.nii.gz
        ├── 11_Inferior-fronto-occipital-fasciculus_L.nii.gz
        ├── 12_Inferior-fronto-occipital-fasciculus_R.nii.gz
        ├── 13_Inferior_Longitudinal-fasciculus_L.nii.gz
        ├── 14_Inferior_Longitudinal-fasciculus_R.nii.gz
        ├── 15_Superior_Longitudinal-fasciculus_L.nii.gz
        ├── 16_Superior_Longitudinal-fasciculus_R.nii.gz
        ├── 17_Uncinate-fasciculus_L.nii.gz
        ├── 18_Uncinate-fasciculus_R.nii.gz
        ├── 19_Superior_Longitudinal-fasciculus-temporal-part_L.nii.gz
        ├── 20_Superior_Longitudinal-fasciculus-temporal-part_R.nii.gz
        ├── 21_Anterior-thalamic_Radiation.nii.gz
        ├── 22_Corticospinal-tract.nii.gz
        ├── 23_Cingulum-cingulate-gyrus.nii.gz
        ├── 24_Cingulum-hippocampus.nii.gz
        ├── 25_Forceps.nii.gz
        ├── 26_Inferior-fronto-occipital-fasciculus.nii.gz
        ├── 27_Inferior_Longitudinal-fasciculus.nii.gz
        ├── 28_Superior_Longitudinal-fasciculus.nii.gz
        ├── 29_Uncinate-fasciculus.nii.gz
        └── 30_Superior_Longitudinal-fasciculus-temporal-part.nii.gz

ソースコード

以上の準備が完了したら、以下のソースコードをShellScriptで実行する。

具体的な、処理は次の通り。

  1. 各施設の拡散定量値マップをDWI maskでマスキング(※ただし、SiteBに関しては、SiteB_beforeのDWI maskのみを使用して、SiteB_before, SiteB_afterのマップをマスキング)
  2. SiteAとSiteB_beforeの位置合わせ(剛体変換)
  3. SiteAーSiteB_beforeの変換行列を計算
  4. SiteBの拡散定量値マップをSiteA空間に移動(※SiteB_afterもSiteAーSiteB_beforeの変換行列を用いて、SiteAに移動)
  5. SiteAとSiteBにおける拡散定量値マップの絶対差を計算
  6. JHU atlasで絶対差マップをサンプリング
  7. 絶対差マップを表示
# Analyze map difference
# .
# └── data
#     ├── SiteB_after  # After harmonization
#     │   └── map
#     ├── SiteB_before  # Before harmonization
#     │   └── map
#     ├── SiteA
#     │   └── map
#     └── mask  # ROI in SiteA space
#         ├── JHU-ICBM-labels
#         └── JHU-ICBM-tracts-maxprob-thr25

function masking_map() {
    for SITE in SiteB_after SiteB_before SiteA; do
        INPUT_DIR=data/${SITE}/map
        OUT_DIR=result/masked_map/${SITE}
        MASK=data/${SITE}/drdwi_mask
        mkdir -p ${OUT_DIR}
        for MAP in $(ls data/SiteA/map | grep nii); do
            if [ ${SITE} != 'SiteB_after' ]; then
                fslmaths ${INPUT_DIR}/${MAP} -mas ${MASK} ${OUT_DIR}/${MAP}
            else
                fslmaths ${INPUT_DIR}/${MAP} -mas data/SiteB_before/drdwi_mask ${OUT_DIR}/${MAP}
            fi
        done
    done
}

function correct_digit() {
    for SITE in SiteB_after SiteB_before SiteA; do
        INPUT_DIR=result/masked_map/${SITE}
        OUT_DIR=result/masked_map/${SITE}
        for MAP in $(ls data/SiteA/map | grep nii| cut -d . -f1); do
            if [[ ${MAP} == 'md' || ${MAP} == 'ad' || ${MAP} == 'rd' ]]; then
                fslmaths ${INPUT_DIR}/${MAP} -mul 1000 ${OUT_DIR}/${MAP}
            fi
        done
    done
}

function each_reg_to_SiteA() { # Register SiteB to SiteA
    SiteA_DIR=result/masked_map/SiteA
    for SITE in SiteB_after SiteB_before; do
        SiteB_DIR=result/masked_map/${SITE}
        OUT_DIR=result/SiteB-map_in_SiteA/${SITE}
        # INPUT_PREFIX=drdwi
        mkdir -p ${OUT_DIR}
        # Register SiteB to SiteA
        flirt -in ${SiteB_DIR}/fa \
            -ref ${SiteA_DIR}/fa \
            -dof 6 \
            -omat ${OUT_DIR}/${SITE}_to_SiteA.mat
        # Apply regit transformation
        for MAP in $(ls ${SiteA_DIR} | grep nii); do
            flirt -in ${SiteB_DIR}/${MAP} \
                -ref ${SiteA_DIR}/fa \
                -out ${OUT_DIR}/${MAP} \
                -init ${OUT_DIR}/${SITE}_to_SiteA.mat \
                -applyxfm -interp nearestneighbour
        done
    done
}

function reg_to_SiteA() { # Register SiteB to SiteA
    SiteA_DIR=result/masked_map/SiteA
    SiteB_DIR=result/masked_map
    OUT_DIR=result/SiteB-map_in_SiteA
    mkdir -p ${OUT_DIR}/SiteB_before ${OUT_DIR}/SiteB_after
    # Register SiteB to SiteA
    flirt -in ${SiteB_DIR}/SiteB_before/fa \
        -ref ${SiteA_DIR}/fa \
        -dof 6 \
        -omat ${OUT_DIR}/SiteB_before_to_SiteA.mat
    for SITE in SiteB_after SiteB_before; do
        # Apply regit transformation
        for MAP in $(ls ${SiteA_DIR} | grep nii); do
            flirt -in ${SiteB_DIR}/${SITE}/${MAP} \
                -ref ${SiteA_DIR}/fa \
                -out ${OUT_DIR}/${SITE}/${MAP} \
                -init ${OUT_DIR}/SiteB_before_to_SiteA.mat \
                -applyxfm -interp nearestneighbour
        done
    done
}

function subtract_by_SiteA() {
    SiteA_DIR=result/masked_map/SiteA
    for SITE in SiteB_after SiteB_before; do
        SiteB_DIR=result/SiteB-map_in_SiteA/${SITE}
        OUT_DIR=result/subtracted-map_by_SiteA/${SITE}
        mkdir -p ${OUT_DIR}
        for MAP in $(ls ${SiteA_DIR} | grep nii); do
            fslmaths ${SiteB_DIR}/${MAP} -sub ${SiteA_DIR}/${MAP} -abs ${OUT_DIR}/${MAP}
        done
    done
}

function sampling_map() {
    SiteA_DIR=result/masked_map/SiteA
    MASK_DIR=data/mask
    for SITE in SiteB_after SiteB_before; do
        SiteB_DIR=result/subtracted-map_by_SiteA/${SITE}
        OUT_DIR=result/sampling_subtracted-map/${SITE}
        # [[ -d ${OUT_DIR} ]] && \rm -rf ${OUT_DIR}
        mkdir -p ${OUT_DIR}
        for ATLAS in JHU-ICBM-labels JHU-ICBM-tracts-maxprob-thr25; do
            for MAP in $(ls ${SiteA_DIR} | grep nii | cut -d . -f1); do
                for MASK in $(ls ${MASK_DIR}/${ATLAS}); do
                    fslstats ${SiteB_DIR}/${MAP} -k ${MASK_DIR}/${ATLAS}/${MASK} -M >>${OUT_DIR}/${ATLAS}_${MAP}_mean.txt &
                    fslstats ${SiteB_DIR}/${MAP} -k ${MASK_DIR}/${ATLAS}/${MASK} -S >>${OUT_DIR}/${ATLAS}_${MAP}_sd.txt &
                    wait
                done
            done
        done
    done
}

function show_subtracted-map() {
    for SITE in SiteB_after SiteB_before; do
        MAP_DIR=result/subtracted-map_by_SiteA/${SITE}
        fsleyes ${MAP_DIR}/fa.nii.gz -cm brain_colours_nih \
            ${MAP_DIR}/md.nii.gz -cm brain_colours_nih \
            ${MAP_DIR}/ad.nii.gz -cm brain_colours_nih \
            ${MAP_DIR}/rd.nii.gz -cm brain_colours_nih \
            ${MAP_DIR}/mk.nii.gz -cm brain_colours_nih \
            ${MAP_DIR}/ak.nii.gz -cm brain_colours_nih \
            ${MAP_DIR}/rk.nii.gz -cm brain_colours_nih \
            ${MAP_DIR}/icvf.nii.gz -cm brain_colours_nih \
            ${MAP_DIR}/isovf.nii.gz -cm brain_colours_nih \
            ${MAP_DIR}/odi.nii.gz -cm brain_colours_nih &
    done
}

# 1. Mask map
masking_map
# correct_digit  # Diffusivity (MD, AD, RDの桁を他のMAPと合わせたい場合)
# 2. Register SiteB to SiteA using FA map
reg_to_SiteA
# 3. Subtract SiteB by SiteA map
subtract_by_SiteA
# 4. Sampling map using JHU atlas
sampling_map
# 5. Show subtracted map by SiteA map
show_subtracted-map

結果

コードを実行すると、「result」フォルダが生成され、そのフォルダ内に結果が格納される。

  • masked_map : マスキング後の拡散定量値マップ
  • SiteB-map_in_SiteA : SiteA空間に位置合わせした、SiteBの拡散定量値マップ
  • subtracted-map_by_SiteA : SiteAーSiteBの絶対差マップ
  • sampling_subtracted-map : SiteAーSiteBの絶対差マップをJHU atlasでサンプリングした結果
result/
├── SiteB-map_in_SiteA
│   ├── SiteB_after
│   │   ├── ad.nii.gz
│   │   ├── ak.nii.gz
│   │   ├── fa.nii.gz
│   │   ├── md.nii.gz
│   │   ├── mk.nii.gz
│   │   ├── rd.nii.gz
│   │   └── rk.nii.gz
│   ├── SiteB_before
│   │   ├── ad.nii.gz
│   │   ├── ak.nii.gz
│   │   ├── fa.nii.gz
│   │   ├── md.nii.gz
│   │   ├── mk.nii.gz
│   │   ├── rd.nii.gz
│   │   └── rk.nii.gz
│   └── SiteB_before_to_SiteA.mat
├── masked_map
│   ├── SiteB_after
│   │   ├── ad.nii.gz
│   │   ├── ak.nii.gz
│   │   ├── fa.nii.gz
│   │   ├── md.nii.gz
│   │   ├── mk.nii.gz
│   │   ├── rd.nii.gz
│   │   └── rk.nii.gz
│   ├── SiteB_before
│   │   ├── ad.nii.gz
│   │   ├── ak.nii.gz
│   │   ├── fa.nii.gz
│   │   ├── md.nii.gz
│   │   ├── mk.nii.gz
│   │   ├── rd.nii.gz
│   │   └── rk.nii.gz
│   └── SiteA
│       ├── ad.nii.gz
│       ├── ak.nii.gz
│       ├── fa.nii.gz
│       ├── md.nii.gz
│       ├── mk.nii.gz
│       ├── rd.nii.gz
│       └── rk.nii.gz
├── sampling_subtracted-map
│   ├── SiteB_after
│   │   ├── JHU-ICBM-labels_ad_mean.txt
│   │   ├── JHU-ICBM-labels_ad_sd.txt
│   │   ├── JHU-ICBM-labels_ak_mean.txt
│   │   ├── JHU-ICBM-labels_ak_sd.txt
│   │   ├── JHU-ICBM-labels_fa_mean.txt
│   │   ├── JHU-ICBM-labels_fa_sd.txt
│   │   ├── JHU-ICBM-labels_md_mean.txt
│   │   ├── JHU-ICBM-labels_md_sd.txt
│   │   ├── JHU-ICBM-labels_mk_mean.txt
│   │   ├── JHU-ICBM-labels_mk_sd.txt
│   │   ├── JHU-ICBM-labels_rd_mean.txt
│   │   ├── JHU-ICBM-labels_rd_sd.txt
│   │   ├── JHU-ICBM-labels_rk_mean.txt
│   │   ├── JHU-ICBM-labels_rk_sd.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_ad_mean.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_ad_sd.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_ak_mean.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_ak_sd.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_fa_mean.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_fa_sd.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_md_mean.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_md_sd.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_mk_mean.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_mk_sd.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_rd_mean.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_rd_sd.txt
│   │   ├── JHU-ICBM-tracts-maxprob-thr25_rk_mean.txt
│   │   └── JHU-ICBM-tracts-maxprob-thr25_rk_sd.txt
│   └── SiteB_before
│       ├── JHU-ICBM-labels_ad_mean.txt
│       ├── JHU-ICBM-labels_ad_sd.txt
│       ├── JHU-ICBM-labels_ak_mean.txt
│       ├── JHU-ICBM-labels_ak_sd.txt
│       ├── JHU-ICBM-labels_fa_mean.txt
│       ├── JHU-ICBM-labels_fa_sd.txt
│       ├── JHU-ICBM-labels_md_mean.txt
│       ├── JHU-ICBM-labels_md_sd.txt
│       ├── JHU-ICBM-labels_mk_mean.txt
│       ├── JHU-ICBM-labels_mk_sd.txt
│       ├── JHU-ICBM-labels_rd_mean.txt
│       ├── JHU-ICBM-labels_rd_sd.txt
│       ├── JHU-ICBM-labels_rk_mean.txt
│       ├── JHU-ICBM-labels_rk_sd.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_ad_mean.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_ad_sd.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_ak_mean.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_ak_sd.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_fa_mean.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_fa_sd.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_md_mean.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_md_sd.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_mk_mean.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_mk_sd.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_rd_mean.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_rd_sd.txt
│       ├── JHU-ICBM-tracts-maxprob-thr25_rk_mean.txt
│       └── JHU-ICBM-tracts-maxprob-thr25_rk_sd.txt
└── subtracted-map_by_SiteA
    ├── SiteB_after
    │   ├── ad.nii.gz
    │   ├── ak.nii.gz
    │   ├── fa.nii.gz
    │   ├── md.nii.gz
    │   ├── mk.nii.gz
    │   ├── rd.nii.gz
    │   └── rk.nii.gz
    └── SiteB_before
        ├── ad.nii.gz
        ├── ak.nii.gz
        ├── fa.nii.gz
        ├── md.nii.gz
        ├── mk.nii.gz
        ├── rd.nii.gz
        └── rk.nii.gz

以上のソースコードを各被験者データに対して、実行した後、以下のコマンドで結果をまとめるとよい。以下のコマンドをShellScriptで実行すると、「summary」フォルダが生成され、そのフォルダにすべての被験者の結果がまとめられる。

# .
# ├── Subj001
# ├── Subj002
# ├── Subj003
# ├── ...
# └── Subj100

rm -r summary
ls -F | grep / | cut -d / -f1 | grep -v temp_ID | tr "\n" "\t" | sed -e 's/^/Region\t/' | sed -e 's/$/\n/' >temp_ID

for ATLAS in JHU-ICBM-labels JHU-ICBM-tracts-maxprob-thr25; do
    ATLAS_DIR=data/mask/${ATLAS}
    ls $(ls -F | grep / | cut -d / -f1 | head -n 1)/${ATLAS_DIR} >temp_${ATLAS}_name

    for SITE in SiteB_before SiteB_after; do
        OUT_DIR=summary/${SITE}
        mkdir -p ${OUT_DIR}
        for MAP in fa md ad rd mk ak rk icvf isovf od; do
            SAMPLE_DIR=result/sampling_subtracted-map/${SITE}/${ATLAS}_${MAP}_mean.txt
            paste temp_${ATLAS}_name */${SAMPLE_DIR} >${OUT_DIR}/temp_${ATLAS}_${MAP}_mean_value
            cat temp_ID ${OUT_DIR}/temp_${ATLAS}_${MAP}_mean_value >${OUT_DIR}/${ATLAS}_${MAP}_mean_result.csv
        done
    done
done

rm temp* summary/*/temp*

すべての被験者の結果をまとめた「summary」フォルダの中身は、次の通り。

summary/
├── SiteB_after
│   ├── JHU-ICBM-labels_ad_mean_result.csv
│   ├── JHU-ICBM-labels_ak_mean_result.csv
│   ├── JHU-ICBM-labels_fa_mean_result.csv
│   ├── JHU-ICBM-labels_icvf_mean_result.csv
│   ├── JHU-ICBM-labels_isovf_mean_result.csv
│   ├── JHU-ICBM-labels_md_mean_result.csv
│   ├── JHU-ICBM-labels_mk_mean_result.csv
│   ├── JHU-ICBM-labels_od_mean_result.csv
│   ├── JHU-ICBM-labels_rd_mean_result.csv
│   ├── JHU-ICBM-labels_rk_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_ad_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_ak_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_fa_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_icvf_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_isovf_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_md_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_mk_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_od_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_rd_mean_result.csv
│   └── JHU-ICBM-tracts-maxprob-thr25_rk_mean_result.csv
└── SiteB_before
    ├── JHU-ICBM-labels_ad_mean_result.csv
    ├── JHU-ICBM-labels_ak_mean_result.csv
    ├── JHU-ICBM-labels_fa_mean_result.csv
    ├── JHU-ICBM-labels_icvf_mean_result.csv
    ├── JHU-ICBM-labels_isovf_mean_result.csv
    ├── JHU-ICBM-labels_md_mean_result.csv
    ├── JHU-ICBM-labels_mk_mean_result.csv
    ├── JHU-ICBM-labels_od_mean_result.csv
    ├── JHU-ICBM-labels_rd_mean_result.csv
    ├── JHU-ICBM-labels_rk_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_ad_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_ak_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_fa_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_icvf_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_isovf_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_md_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_mk_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_od_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_rd_mean_result.csv
    └── JHU-ICBM-tracts-maxprob-thr25_rk_mean_result.csv

例えば、「JHU-ICBM-tracts-maxprob-thr25_fa_mean_result.csv」の中身は、次の通り。

Region Subj001 Subj002 Subj100
01_Anterior-thalamic_Radiation_L.nii.gz 0.030828 0.031981 0.032008
02_Anterior-thalamic_Radiation_R.nii.gz 0.033499 0.030607 0.030954
03_Corticospinal-tract_L.nii.gz 0.036112 0.038567 0.035203
04_Corticospinal-tract_R.nii.gz 0.033464 0.038637 0.036139
05_Cingulum-cingulate-gyrus_L.nii.gz 0.05438 0.035468 0.039964
06_Cingulum-cingulate-gyrus_R.nii.gz 0.037334 0.033929 0.037894
07_Cingulum-hippocampus_L.nii.gz 0.043861 0.022866 0.028342
08_Cingulum-hippocampus_R.nii.gz 0.04599 0.025942 0.039053
09_Forceps-major.nii.gz 0.038693 0.036278 0.037092
10_Forceps-minor.nii.gz 0.033049 0.03118 0.02924
11_Inferior-fronto-occipital-fasciculus_L.nii.gz 0.024478 0.0306 0.024691
12_Inferior-fronto-occipital-fasciculus_R.nii.gz 0.028521 0.02798 0.024997
13_Inferior_Longitudinal-fasciculus_L.nii.gz 0.023941 0.030016 0.025934
14_Inferior_Longitudinal-fasciculus_R.nii.gz 0.02708 0.02826 0.026991
15_Superior_Longitudinal-fasciculus_L.nii.gz 0.025323 0.020107 0.022455
16_Superior_Longitudinal-fasciculus_R.nii.gz 0.02511 0.021673 0.022644
17_Uncinate-fasciculus_L.nii.gz 0.024351 0.022198 0.024761
18_Uncinate-fasciculus_R.nii.gz 0.026874 0.027458 0.027105
19_Superior_Longitudinal-fasciculus-temporal-part_L.nii.gz 0.02856 0.021351 0.031076
20_Superior_Longitudinal-fasciculus-temporal-part_R.nii.gz 0.039501 0.020662 0.028387
21_Anterior-thalamic_Radiation.nii.gz 0.03213 0.031379 0.031527
22_Corticospinal-tract.nii.gz 0.034761 0.038603 0.035712
23_Cingulum-cingulate-gyrus.nii.gz 0.048843 0.0351 0.039315
24_Cingulum-hippocampus.nii.gz 0.045027 0.024657 0.034999
25_Forceps.nii.gz 0.034184 0.032299 0.030945
26_Inferior-fronto-occipital-fasciculus.nii.gz 0.026715 0.029178 0.024856
27_Inferior_Longitudinal-fasciculus.nii.gz 0.025091 0.02938 0.026332
28_Superior_Longitudinal-fasciculus.nii.gz 0.025231 0.020793 0.022546
29_Uncinate-fasciculus.nii.gz 0.025345 0.024466 0.025865
30_Superior_Longitudinal-fasciculus-temporal-part.nii.gz 0.035917 0.020796 0.028871

SiteAとSiteBにおける拡散定量値マップの絶対差マップは、「subtracted-map_by_SiteA」に保存される。

f:id:AIProgrammer:20210226124135p:plain

3. 脳領域ごとに絶対差を図示

JHU atlasでサンプリングした、SiteAーSiteBの拡散定量値マップの絶対差(施設間差)を、箱ひげ図で図示する。

準備

すべての被験者の結果をまとめた「summary」フォルダを「data」フォルダに変更して、使用する。

data/
├── SiteB_after
│   ├── JHU-ICBM-labels_ad_mean_result.csv
│   ├── JHU-ICBM-labels_ak_mean_result.csv
│   ├── JHU-ICBM-labels_fa_mean_result.csv
│   ├── JHU-ICBM-labels_icvf_mean_result.csv
│   ├── JHU-ICBM-labels_isovf_mean_result.csv
│   ├── JHU-ICBM-labels_md_mean_result.csv
│   ├── JHU-ICBM-labels_mk_mean_result.csv
│   ├── JHU-ICBM-labels_od_mean_result.csv
│   ├── JHU-ICBM-labels_rd_mean_result.csv
│   ├── JHU-ICBM-labels_rk_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_ad_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_ak_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_fa_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_icvf_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_isovf_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_md_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_mk_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_od_mean_result.csv
│   ├── JHU-ICBM-tracts-maxprob-thr25_rd_mean_result.csv
│   └── JHU-ICBM-tracts-maxprob-thr25_rk_mean_result.csv
└── SiteB_before
    ├── JHU-ICBM-labels_ad_mean_result.csv
    ├── JHU-ICBM-labels_ak_mean_result.csv
    ├── JHU-ICBM-labels_fa_mean_result.csv
    ├── JHU-ICBM-labels_icvf_mean_result.csv
    ├── JHU-ICBM-labels_isovf_mean_result.csv
    ├── JHU-ICBM-labels_md_mean_result.csv
    ├── JHU-ICBM-labels_mk_mean_result.csv
    ├── JHU-ICBM-labels_od_mean_result.csv
    ├── JHU-ICBM-labels_rd_mean_result.csv
    ├── JHU-ICBM-labels_rk_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_ad_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_ak_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_fa_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_icvf_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_isovf_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_md_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_mk_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_od_mean_result.csv
    ├── JHU-ICBM-tracts-maxprob-thr25_rd_mean_result.csv
    └── JHU-ICBM-tracts-maxprob-thr25_rk_mean_result.csv

ソースコード

以下のソースコードをPythonで実行する。ここでは、「JHU-ICBM-tracts-maxprob-thr25 atlas」を用いた結果に対して、処理を実行する。

処理内容は、次の通り。

  1. CSVデータの読み込み
  2. 左右平均された脳領域のみを抽出
  3. SiteB_beforeで脳領域ごとに平均を算出し、小さい順に脳領域を並び替え
  4. SiteB_beforeとSiteB_afterを結合
  5. 領域名を略語に置換
  6. 箱ひげ図で図示・保存
# Make figure
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
from scipy import stats
import scikit_posthocs as sp
import statsmodels


def import_data(x1='before', x2='after'):
    # Read CSV
    df_before = pd.read_csv(x1, delimiter='\t').dropna(axis=1).iloc[-10:, :]
    df_after = pd.read_csv(x2, delimiter='\t').dropna(axis=1).iloc[-10:, :]
    # Add label
    df_before['Label'] = 'before'
    df_after['Label'] = 'after'
    # Sort columns by values
    sorted_regions = df_before.set_index('Region').mean(
        axis=1).sort_values().index.values.tolist()
    df_before = df_before.set_index('Region').reindex(
        index=sorted_regions).reset_index()
    # Merge before and after
    df_merge = df_before.merge(df_after, how='outer')
    # Melt
    df_merge_melt = pd.melt(
        df_merge, id_vars=['Region', 'Label'], var_name='ID')
    # Replace region name
    df_merge_melt_rpl = df_merge_melt.replace({
        '21_Anterior-thalamic_Radiation.nii.gz': 'SiteB',
        '22_Corticospinal-tract.nii.gz': 'CST',
        '23_Cingulum-cingulate-gyrus.nii.gz': 'CCG',
        '24_Cingulum-hippocampus.nii.gz': 'CH',
        '25_Forceps.nii.gz': 'Forceps',
        '26_Inferior-fronto-occipital-fasciculus.nii.gz': 'IFOF',
        '27_Inferior_Longitudinal-fasciculus.nii.gz': 'ILF',
        '28_Superior_Longitudinal-fasciculus.nii.gz': 'SLF',
        '29_Uncinate-fasciculus.nii.gz': 'UF',
        '30_Superior_Longitudinal-fasciculus-temporal-part.nii.gz': 'SLFT'
    })
    return df_merge_melt_rpl


def box_plt(df, map):
    # Plot
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    sns.boxplot(x='Region', y='value', data=df, hue='Label',
                showfliers=False, palette="Set1", width=0.8)
    # Config
    # fig.subplots_adjust(bottom=0.3)
    handles, labels = ax.get_legend_handles_labels()
    ax.legend(handles[0:2], labels[0:2])
    fig.autofmt_xdate(rotation=45)
    plt.gca().spines['right'].set_visible(False)
    plt.gca().spines['top'].set_visible(False)
    ax.set(xlabel='', ylabel="Δ{}".format(map.upper()))
    # plt.tick_params(pad=10)
    # plt.show()
    # Save Figure
    plt.savefig("./result/Diff_{}.png".format(map.upper()),
                dpi=300, bbox_inches="tight")
    plt.close()


def main():
    map_list = ['fa', 'md', 'ad', 'rd',
                'mk', 'ak', 'rk',
                'icvf', 'isovf', 'odi']

    if not os.path.isdir('result'):
        os.makedirs('result')

    for map in map_list:
        # atlas='JHU-ICBM-labels'
        atlas = 'JHU-ICBM-tracts-maxprob-thr25'
        x1 = "data/SiteB_before/{}_{}_mean_result.csv".format(atlas, map)
        x2 = "data/SiteB_after/{}_{}_mean_result.csv".format(atlas, map)
        # Import data
        df = import_data(x1, x2)
        # Box plot
        box_plt(df, map)


if __name__ == '__main__':
    main()

結果

以下のソースコードを実行すると、「result」フォルダが生成され、拡散定量値マップごとにJHU atlasでサンプリングした、絶対差の箱ひげ図が格納される。

result/
├── Diff_FA.png
├── Diff_MD.png
├── Diff_AD.png
├── Diff_RD.png
├── Diff_MK.png
├── Diff_AK.png
├── Diff_RK.png
├── Diff_ICVF.png
├── Diff_ISOVF.png
└── Diff_ODI.png

以下は、ΔDTI(ΔFA, ΔMD, ΔAD, ΔRD)の箱ひげ図である。

f:id:AIProgrammer:20210226123916p:plain


頑張れ!喝!!の代わりにB!ブックマークを押していただけるとただただうれしいです(^^)! ↓

動かしながら学ぶ PyTorchプログラミング入門

動かしながら学ぶ PyTorchプログラミング入門

  • 作者:斎藤勇哉
  • 発売日: 2020/11/30
  • メディア: Kindle版