1. 目的
2. GBSSとは
2.1. データの準備
2.2. 解析パラメータの設定
2.3. gbss_1_T1wpreproc：T1WIの前処理
2.4. gbss_2_DWIwpreproc：拡散定量値の前処理
2.5. gbss_3_skelpreproc：拡散定量値をスケルトンに投影
2.6. randomise：スケルトンに投影された定量値画像を入力したGLMと並べ替え検定（permutation test）
3. おまけ
1. 目的
- Gray matter-Based Spatial Statistics: GBSS
2. GBSSとは
Gray matter-Based Spatial Statistics（GBSS）は、灰白質の統計解析をするための手法。TBSSの灰白質版。
灰白質の中心線（skeleton）に定量値を投影する。通常の脳画像の統計解析では、脳構造の個人差を除外するために空間的「平滑化」を用いる。しかし、平滑化の程度に原則がなく、平滑化をかけては情報があいまいになり、MRIの高空間分解能を生かせないという問題がある。一方、GBSSでは、灰白質の中心線と思われるところにskeletonを生成し、そこに個人ごとの定量値を投影するという手法をとる。これにより、平滑化せずに群間比較をすることができるため、平滑化による問題を回避できるという利点がある。
GBSS解析では、次のような処理をする。
- データの準備
- 解析パラメータの設定
gbss_1_T1wpreproc：T1WIの前処理gbss_2_DWIwpreproc：拡散定量値の前処理gbss_3_skelpreproc：拡散定量値をスケルトンに投影randomise：スケルトンに投影された定量値画像を入力したGLMと並べ替え検定（permutation test）
2.1. データの準備
まず、次のデータを準備する。
- b0フォルダ：b0 (DWI, b=0 s/mm^2)
- T1wフォルダ：3D-T1WI
- 拡散定量値：FA, MD, FW, etc.
フォルダ構造は次のようにする。
ここでは、健常者10名（ID: Con0001~Con0010）と患者10名（ID: Pat0001~Pat0010）いることを想定する。各フォルダに入れるファイル名は同じにする必要がある（例：FA/Subj001.nii.gz, T1w/Subj001.nii.gzは同じファイル名）。
.
├── FA
│ ├── Con0001.nii.gz
│ ├── Con0002.nii.gz
│ ├── ...
│ └── Pat0010.nii.gz
├── FW
│ ├── Con0001.nii.gz
│ ├── Con0002.nii.gz
│ ├── ...
│ └── Pat0010.nii.gz
├── MD
│ ├── Con0001.nii.gz
│ ├── Con0002.nii.gz
│ ├── ...
│ └── Pat0010.nii.gz
├── T1w
│ ├── Con0001.nii.gz
│ ├── Con0002.nii.gz
│ ├── ...
│ └── Pat0010.nii.gz
└── b0
├── Con0001.nii.gz
├── Con0002.nii.gz
├── ...
└── Pat0010.nii.gz
2.2. 解析パラメータの設定
解析に用いるパラメータを設定する。
PEDIR、ECHOSPACINGはepi_regコマンドのパラメータであり、それぞれ位相エンコード方向とecho spacing timeを指定する。SKELETON_THRはスケルトン化する際のしきい値である。
MAP_LIST=$(ls | grep -v b0 | grep -v T1w | grep -v stats) SUBJ_LIST=$(ls T1w | cut -d . -f1) PEDIR='-y' # Phase encoding direction: "epi_reg" config ECHOSPACING='0.0380544' # Echo spacing time (s): "epi_reg" config SKELETON_THR=0.2 # Skeleton threshold
2.3. gbss_1_T1wpreproc：T1WIの前処理
T1WIの前処理内容は、次の通り。
- 脳頭蓋除去
- 灰白質をセグメント
- 灰白質を標準脳（MNI152）に位置合わせ
- 標準空間上の各被験者の灰白質を平均化
以上の処理内容を実行するために、関数gbss_1_T1wpreprocを定義。
function gbss_1_T1wpreproc() {
echo "T1w preprocessing..."
# T1w preprocessing
for SUBJ in ${SUBJ_LIST}; do
## Segment GM
### Skull-stripping
bet T1w/${SUBJ} T1w/${SUBJ}_skull_stripped -f 0.3 -R -S -B
### Segment T1w into CSF/GM/WM
fast -t 1 -g -B -b -p -o T1w/${SUBJ}_fast T1w/${SUBJ}_skull_stripped
## Register GM to MNI
### Register T1w to MNI
flirt -ref ${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
-in T1w/${SUBJ}_fast_restore \
-omat T1w/${SUBJ}_indiv2MNI.mat
fnirt --in=T1w/${SUBJ}_fast_restore \
--aff=T1w/${SUBJ}_indiv2MNI.mat \
--cout=T1w/${SUBJ}_indiv2MNI_warp \
--config=T1_2_MNI152_2mm
### applywarp to GM and move it into MNI
applywarp --ref=${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
--in=T1w/${SUBJ}_fast_pve_1 \
--warp=T1w/${SUBJ}_indiv2MNI_warp \
--out=T1w/${SUBJ}_fast_pve_1_inMNI &
done
## Create 4D all_GM image
### GM_pve1 in MNI are merged into one volume
mkdir stats
fslmerge -t stats/all_GM $(imglob T1w/*_fast_pve_1_inMNI.nii.gz)
### Create merged GM_pve1 (all_GM) mask and apply it to all_GM
fslmaths stats/all_GM -max 0 -Tmin -bin stats/mean_GM_mask -odt char
fslmaths stats/all_GM -Tmean stats/mean_GM
}
関数が定義できたら、実行する。
gbss_1_T1wpreproc
処理が完了すると、statsフォルダにすべての被験者の灰白質画像がまとめられた「all_GM.nii.gz」とそれを平均化した「mean_GM.nii.gz」が生成される。
2.4. gbss_2_DWIwpreproc：拡散定量値の前処理
拡散定量値の前処理の過程は、次の通り。
- b0をT1WIに位置合わせ
- b0 to T1WIとT1WI to 標準脳のWarp Fieldを使って、拡散定量値マップを標準空間に移動
- 標準空間上にあるすべての被験者の拡散定量値マップをひとつの4D volume dataにまとめる
- 平均灰白質マスク
mean_GM_maskで、3の拡散定量値マップをマスク処理
function gbss_2_DWIwpreproc() {
echo "Diffusion MAP preprocessing..."
# Diffusion MAP preprocessing
for SUBJ in ${SUBJ_LIST}; do
## Register b0 to MNI
### Register b0 to T1WI
epi_reg --epi=b0/${SUBJ} \
--t1=T1w/${SUBJ} --t1brain=T1w/${SUBJ}_skull_stripped \
--echospacing=${ECHOSPACING} --pedir=${PEDIR} \
--wmseg=T1w//${SUBJ}_fast_seg_2.nii.gz \
--out=b0/${SUBJ}_BBR
for MAP in ${MAP_LIST}; do
### applywarp to diffusion map and move it to T1w based on epi_reg
flirt -in ${MAP}/${SUBJ} \
-out ${MAP}/${SUBJ}_inT1 \
-ref T1w/${SUBJ}_skull_stripped \
-init b0/${SUBJ}_BBR.mat -applyxfm
### applywarp to dMAP in T1w and move it to MNI
applywarp --ref=${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
--in=${MAP}/${SUBJ}_inT1 \
--warp=T1w/${SUBJ}_indiv2MNI_warp \
--out=${MAP}/${SUBJ}_inMNI
done
done
## 4D all_${dMAP}
for MAP in ${MAP_LIST}; do
### Merge dMAP in MNI into one volume, all_${MAP}
fslmerge -t stats/all_${MAP} $(imglob ${MAP}/*_inMNI.nii.gz)
### Masking dMAP in MNI using mean GM mask
fslmaths stats/all_${MAP} -mas stats/mean_GM_mask stats/all_${MAP}
done
}
関数が定義できたら、実行する。
gbss_2_DWIwpreproc
処理が完了すると、拡散定量値の種類に応じてstatsフォルダに「all_
以下は、標準空間上における灰白質領域のFA画像である。
2.5. gbss_3_skelpreproc：拡散定量値をスケルトンに投影
拡散定量値をスケルトンに投影するまでの手順は、次の通り。
- 平均灰白質画像からスケルトンを生成
- 灰白質スケルトンをしきい値処理をしてスケルトンマスクを生成
- スケルトンの距離マップ（distance map）を生成
- 各拡散定量値マップを、スケルトンに投影
function gbss_3_skelpreproc() {
echo "Projecting MAP into skeleton..."
# Skeletonize dMAP
## Create GM skeleton
tbss_skeleton -i stats/mean_GM -o stats/mean_GM_skeleton
## threshold GM skeleton
fslmaths stats/mean_GM_skeleton -thr ${SKELETON_THR} -bin stats/mean_GM_skeleton_mask
## Create skeleton distancemap
fslmaths stats/mean_GM_mask -mul -1 -add 1 -add stats/mean_GM_skeleton_mask stats/mean_GM_skeleton_mask_dst
distancemap -i stats/mean_GM_skeleton_mask_dst -o stats/mean_GM_skeleton_mask_dst
for MAP in ${MAP_LIST}; do
## Project masked dMAP in MNI into skeleton
tbss_skeleton -i stats/mean_GM \
-p ${SKELETON_THR} \
stats/mean_GM_skeleton_mask_dst \
${FSLDIR}/data/standard/LowerCingulum_1mm \
stats/all_${MAP} \
stats/all_${MAP}_skeletonised
done
}
処理が完了すると、statsフォルダに灰白質のスケルトンmean_GM_skeleton.nii.gz、灰白質のスケルトンの距離マップmean_GM_skeleton_mask_dst.nii.gz、拡散定量値が投影されたスケルトンall_<MAP>_skeletonised.nii.gzが保存される。
2.6. randomise：スケルトンに投影された定量値画像を入力したGLMと並べ替え検定（permutation test）
まず、GLMのデザインマトリックス（計画行列）とコントラストを設定する。
今回は、健常者10名（ID: Con0001~Con0010）と患者10名（ID: Pat0001~Pat0010）のデータがある。これらのファイルを昇順に並び替えると、先に健常者10名の画像、次に患者10名の画像の順に並ぶ。
Con0001.nii.gz
Con0002.nii.gz
Con0003.nii.gz
...
Pat0010.nii.gz
次に、GLMのデザインマトリックス（計画行列）とコントラストを決める設定ファイルを生成する。design_ttest2 <出力ファイル> <健常者数> <患者数>でコマンドを実行。
design_ttest2 stats/design 10 10
statsフォルダに、デザインマトリックス（design.mat）とコントラスト（design.con）が生成される。
デザインマトリックス（design.mat）の中身を確認。
/Matrixの一列目は健常者データであるかどうか、二列目は患者データであるかを0, 1で表している。行の順番は、被験者ファイルの順番（昇順）に対応する。したがって、これらは対応があるようにしておかなければならない。
cat stats/design.mat
/NumWaves 2
/NumPoints 20
/PPheights 1 1
/Matrix
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
コントラスト（design.con）の中身を確認してみる。
/Matrix一列目は健常者の偏回帰係数、二列目は患者の偏回帰係数に対するもので、行は別々のコントラストである。この場合、一行目は健常者>患者の検定、二行目は健常者<患者の検定に相当する。
cat stats/design.con
/NumWaves 2
/NumContrasts 2
/PPheights 1 1
/Matrix
1 -1
-1 1
デザインマトリックス（計画行列）とコントラストの確認ができたら、randomiseコマンド使ってGLMと並べ替え検定（permutation test）を実行する。
randomiseコマンドの各オプションは、次の通り。
- -i：入力画像
- -m：マスク画像
- -o ：出力画像
- -o ：デザインマトリックス
- -o ：デザインコントラスト
- -n：並べ替え検定の数
- –T2：2D最適化を用いたTFCE
- -x：voxel-wiseのcorrected P値マップ
- –uncorrp：un-corrected P値マップ
- -R：統計値マップ
for MAP in ${MAP_LIST}; do
randomise -i stats/all_${MAP}_skeletonised \
-o stats/gbss_${MAP} \
-m stats/mean_GM_mask \
-d stats/design.mat \
-t stats/design.con \
-n 10000 --T2 -V -x --uncorrp -R
done
次のようなファイルが、生成される。
TFCEを用いた”健常群＞患者群”の検定で、FWE補正をされたP値マップ（gbss_FA_tfce_corrp_tstat1.nii.gz）を確認する。ここで、得られたP値マップは1-P値のマップであることに注意する。つまり、P<.05を有意とするのであれば、P値マップで0.95-1.00の値を見ればよい。
fsleyes ${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
stats/mean_GM_skeleton -cm Green \
stats/gbss_FA_tfce_p_tstat1 -cm Red-Yellow -dr 0.95 1
スケルトンは細いため有意差が見づらい場合がある。そのような時、tbss_fillコマンドが役に立つ。
tbss_fillコマンドの基本的な使い方は、以下の通り。
tbss_fill <P値マップ> <しきい値> <平均FA画像> <出力画像>
TFCEを用いた”健常群＞患者群”の検定で、FWE補正をされたP値マップ（gbss_FA_tfce_corrp_tstat1.nii.gz）の有意差があった領域のみを0.95でしきい値処理をして抽出し、その領域を膨張させる。
tbss_fill stats/gbss_FA_tfce_p_tstat1 0.95 stats/mean_GM stats/gbss_FA_tfce_p_tstat1_fill
赤く表示されている領域は、健常群が患者群よりも有意（FWE-corrected）にFA値が大きいことを示している。
3. おまけ
大量に定量値があり、それらをすべて検定する場合、有意差があったかどうかをすべて確認するのは大変である。そこで、各定量値画像のP値マップが0.95以上の値を持つかどうかを判定し、有意差があった場合のみ、tbss_fillコマンドを実行する。
for PMAP in $(ls stats/ | grep tfce_corrp); do
PMAX=$(fslstats stats/${PMAP} -R | cut -d " " -f2)
echo ${PMAP} >>stats/tmp1.txt
echo ${PMAX} >>stats/tmp2.txt
if [ $(echo "${PMAX} > 0.95" | bc) == 1 ]; then
tbss_fill stats/${PMAP} 0.95 stats/mean_GM stats/${PMAP}_fill
fi
done
paste stats/tmp* >stats/tmp_corrected_P_report.txt
echo -e "$(cat stats/tmp_corrected_P_report.txt)\n\n\n$(cat stats/tmp_corrected_P_report.txt | sort -r -n -k 2)" \
>stats/corrected_P_report.txt
rm stats/tmp*
上のコマンドを実行すると、statsフォルダに各検定とそのP値マップの最大値が記された「corrected_P_report.txt」が出力される。
検定結果を、ファイル名でソート（上段）したものと、P値でソートしたもの（下段）に分けて保存している。
cat stats/corrected_P_report.txt
gbss_FA_tfce_corrp_tstat1.nii.gz 0.992000
gbss_FA_tfce_corrp_tstat2.nii.gz 0.416000
gbss_FW_tfce_corrp_tstat1.nii.gz 0.361839
gbss_FW_tfce_corrp_tstat2.nii.gz 0.997261
gbss_MD_tfce_corrp_tstat1.nii.gz 0.389816
gbss_MD_tfce_corrp_tstat2.nii.gz 0.985748
gbss_FW_tfce_corrp_tstat2.nii.gz 0.997261
gbss_FA_tfce_corrp_tstat1.nii.gz 0.992000
gbss_MD_tfce_corrp_tstat2.nii.gz 0.985748
gbss_FA_tfce_corrp_tstat2.nii.gz 0.416000
gbss_MD_tfce_corrp_tstat1.nii.gz 0.389816
tbss_FW_tfce_corrp_tstat1.nii.gz 0.361839
著者情報： 斎藤 勇哉
順天堂大学医学部 大学院医学研究科 放射線診断学講座所属脳MRI 画像解析が専門であり、テーマは①神経変性疾患の機序解明、②医用人工知能の開発、③多施設データのハーモナイゼーション、④速読が脳に与える影響や学習効果、⑤SNS解析を用いたマーケティング戦略の改善。
医療分野に関わらず、自然言語処理・スクレイピング・データ分析・Web アプリ開発を得意とし、企業や他大学の研究を支援。
主な使用言語は、Python、Shell Script、MATLAB、HTML、CSS
