【FSL】Gray matter-Based Spatial Statistics: GBSS

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
}

関数が定義できたら、実行する。

1	`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
}

関数が定義できたら、実行する。

1	`gbss_2_DWIwpreproc`

処理が完了すると、拡散定量値の種類に応じてstatsフォルダに「all_.nii.gz」が生成される。

以下は、標準空間上における灰白質領域の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 <出力ファイル> <健常者数> <患者数>でコマンドを実行。

1	`design_ttest2 stats/design` `10 10`

statsフォルダに、デザインマトリックス（design.mat）とコントラスト（design.con）が生成される。

デザインマトリックス（design.mat）の中身を確認。

/Matrixの一列目は健常者データであるかどうか、二列目は患者データであるかを0, 1で表している。行の順番は、被験者ファイルの順番（昇順）に対応する。したがって、これらは対応があるようにしておかなければならない。

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一列目は健常者の偏回帰係数、二列目は患者の偏回帰係数に対するもので、行は別々のコントラストである。この場合、一行目は健常者>患者の検定、二行目は健常者<患者の検定に相当する。

1	`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コマンドの基本的な使い方は、以下の通り。

1	`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値でソートしたもの（下段）に分けて保存している。

1	`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