1. 目的
2. 準備
3. 必要なpackageのdownload
4. 必要なmoduleのimport
5. CSVファイルの読み込み
6. 2変数の定義
7. MAEの計算
8. RMSEの計算
9. 使用したコードまとめ
著者別アーカイブ： 斎藤 勇哉
【Python】画像の類似度の計算（Dice係数）
1. 目的
2. 準備
3. フォルダ構造
4. 画像類似度の計算（Dice係数）
4.1. calc_dice.py
4.2. result_diceindex.csv
5. FSLでやりたい場合
5.1. fsl_dicecalc.sh
【Python】カッパ（k）係数の計算
1. 目的
2. scikit-learnのインストール
3. データの準備
3.1. label.csv
4. ソースコード
4.1. calc_kappa.py
5. 実行
6. 結果の解釈
7. Neural Network Consoleを使っている場合
7.1. voutput_result.csv
7.2. calc_kappa_nnc.py
【Python】縦断データの可視化（プロット・箱ひげ図・バイオリン図・信頼区間・ヒストグラム）
1. 目的
2. 準備
2.1. open-visualizationsのダウンロード
2.2. ライブラリのインストール
3. チュートリアル1 (plotnineを用いる場合)
3.1. ライブラリの読み込み
3.2. 保存用のフォルダを用意
3.3. データの読み込み
3.4. データの選択
3.5. プロット
3.6. プロットと直線
3.7. グループごとのプロットの位置を微妙に変える
3.8. プロットの色を変更
3.9. 箱ひげ図 (boxplots)
3.10. バイオリン図 (violin plot)
3.11. 信頼区間 (CI bar)
3.12. 各グループの平均を直線で結ぶ
3.13. プロット・箱ひげ図・バイオリン図・信頼区間
4. チュートリアル2 (matplotlibを使う場合)
4.1. ライブラリの読み込み
4.2. 保存用のフォルダを用意
4.3. データの初期化
4.4. プロット
4.5. プロットと直線
4.6. グループごとのプロットの位置を微妙に変える
4.7. the amount of jitter and create a dataframe containing the jittered x-axis values
4.8. 信頼区間 (CI bar)
4.9. バイオリン図 (violin plot)
4.10. 2群のBeforeとAfterをそれぞれプロット
4.11. さらに信頼区間の追加
4.12. プロット・箱ひげ図・バイオリン図・信頼区間
5. 高画質で保存したい場合
5.1. plotnineの場合
5.2. matplotlibの場合
【Python】相関行列 (Correlation Matrix)とヒートマップ (Heat Map)の作り方
1. 目的
2. シンプルなヒートマップ
2.1. ライブラリのインポート
2.2. データの読み込み
2.3. 相関行列の計算
2.4. ヒートマップの作成
3. プロットの大きさを相関係数に応じて変える
3.1. heatmapzのインストール
3.2. ライブラリのインポート
3.3. データの読み込み
3.4. ヒートマップの作成
【Python】箱ひげ図の作り方
1. 目的
2. 1つの図に3群の結果をプロット
2.1. データ準備
2.2. ソースコード
2.3. 結果確認
2.3.1. Boxplot
2.3.2. Boxplot_with_dot
2.3.3. Violinplot
3. 1つの図に3群の結果を各領域ごとにプロット
3.1. データ準備
3.2. ソースコード
3.3. 結果
4. 1つの図に3つの変数に対して4群の結果を3パターンプロット
4.1. データ準備
4.2. ソースコード
4.3. 結果
【FreeSurfer】FreeSurferの結果のまとめ方 ～aparc, aseg, wmparc, brainstem～
1. 目的
2. フォルダ構造
3. aparcとasegのまとめ方
4. wmparcのまとめ方
5. Brain Stemの計測
6. 出力結果
6.1. aseg
6.2. wmparc
6.3. brainstem
【FSL】FLICAを用いたマルチモーダル独立成分分析
1. 目的
2. リポジトリ
3. 必要なデータ
4. 注意点
5. 実行
5.1. run_FLICA.m
6. 結果
6.1. index.html
6.2. コンポーネントごとの結果（All_in_one_page）
6.3. モダリティごとの結果（FA, MD, AD, RD, etc.）
7. GBSSへの適応
7.1. index.html
7.2. コンポーネントごとの結果（All_in_one_page）
7.3. モダリティごとの結果（FA, MD, AD, RD, etc.）
【FSL】fslstatsを使ったROI解析
1. 目的
2. fslstats
2.1. binary maskを使う場合
2.2. index maskを使う場合
3. 使用例
3.1. フォルダ構造
3.2. コード
3.3. 結果
【FSL】MRI画像のレジストレーション（Registration）と関心領域（ROI）解析
1. 目的
2. JHU White-matter labels & tractography atlasでROI解析する方法
2.1. ディレクトリ構造
2.2. 実行方法
2.3. ソースコード
3. Desikan Killiany AtlasでROI解析する方法
3.1. 前提条件
3.2. ディレクトリ構造
3.3. 実行方法
3.4. ソースコード
4. FreeSurferのwm.seg.nii.gzをDWI空間に位置合わせする方法
4.1. 必要なファイル
4.2. ソースコード
5. 標準空間にあるROIを各被験者脳に位置合わせ
6. bertのwmparcをMNI空間に移動後、個人脳（Perfusion image）にレジストレーション
【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
}
関数が定義できたら、実行する。
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
【QSM】QSM解析 ～MEDIの使い方～
1. 目的
2. 準備
2.1. MEDI toolboxにパスを通す
2.2. QSM解析実行のコード準備
2.3. Write_DICOM.mを修正
2.4. QSM_3Dを集める
3. 実行
3.1. 作業フォルダ(analyzeフォルダ)へ移動
3.2. runme.mを実行
3.3. Outputファイルの確認
4. 複数の被験者データをまとめて処理したい場合
4.1. run_sequential.m
1. 目的
Cornell大学が開発したQSM software (MEDI)を使って解析
(Cornell MRI Research Lab. QSMより、URL: http://pre.weill.cornell.edu/mri/pages/qsm.html)
2. 準備
####MEDI toolboxのダウンロード
こちらからMEDI_toolboxをダウンロード
MEDI_toolboxの中身はこのような感じ。
2.1. MEDI toolboxにパスを通す
MATLABを起動し、パスの設定からサブフォルダも追加をクリックしfunctionsフォルダを選択した後、保存。
2.2. QSM解析実行のコード準備
QSM解析を実行するためのコードは以下。
以下をrunme.m(※任意の名前でOK)で保存し、作業フォルダとなるanalyzeフォルダ(※任意の名前でOK)に保存する。
clear all;clc;close all;
% Set path
% MEDI_set_path
% STEP 1: Import data
[iField,voxel_size,matrix_size,CF,delta_TE,TE,B0_dir,files]=Read_DICOM('AXL_QSM');
%%%% Generate the Magnitude image %%%%
iMag = sqrt(sum(abs(iField).^2,4));
%%%%% provide a Mask here if possible %%%%%%
if (~exist('Mask','var'))
% Mask = genMask(iField, voxel_size);
Mask = BET(iMag,matrix_size,voxel_size);
end
%%%%% provide a noise_level here if possible %%%%%%
if (~exist('noise_level','var'))
noise_level = calfieldnoise(iField, Mask);
end
%%%%% normalize signal intensity by noise to get SNR %%%
iField = iField/noise_level;
% STEP 2a: Field Map Estimation
%%%%%Estimate the frequency offset in each of the voxel using a
%%%%%complex fitting %%%%
[iFreq_raw N_std] = Fit_ppm_complex(iField);
% STEP 2b: Spatial phase unwrapping %%%%
iFreq = unwrapPhase(iMag, iFreq_raw, matrix_size);
% STEP 2c: Background Field Removal
%%%% Background field removal %%%%
[RDF shim] = PDF(iFreq, N_std, Mask,matrix_size,voxel_size, B0_dir);
% CSF Mask for zero referencing
R2s=arlo(TE,abs(iField));
Mask_CSF = extract_CSF(R2s, Mask, voxel_size);
% STEP 3: Dipole Inversion
save RDF.mat RDF iFreq iFreq_raw iMag N_std Mask matrix_size...
voxel_size delta_TE CF B0_dir Mask_CSF;
%%%% run MEDI %%%%%
QSM = MEDI_L1('lambda',1000, 'smv', 5, 'merit');
%%% save to DICOM
% ignore warnings...
Write_DICOM(QSM,files,'QSM')
2.3. Write_DICOM.mを修正
解析後のQSMの保存には、Write_DICOM.mを使用する。(runme.mを実行すれば自動で実行される。)
修正前のファイルでは、dicomwriteでエラーが発生する。
そのため、MEDI_toolbox/functions/Write_DICOM.mを以下のコードにそっくりそのまま書き換える。
やっていることは、dicomwriteのCreateModeをCopyからCreateへチェンジ。
function Write_DICOM(M,files,outdir,opts)
%WRITE_DICOM Summary of this function goes here
% Detailed explanation goes here
defopts.SeriesDescription = 'QSM';
defopts.SeriesInstanceUID = [];
defopts.SeriesNumber = [];
defopts.Convert = @(x) convert(x);
defopts.Window = 0.500;
defopts.Level = 0;
defopts.flag3D = 0;
% defopts.EchoNumber = [];
% defopts.EchoTime = 0.0;
if nargin<4; opts=struct; end
deffields=fieldnames(defopts);
for i=1:length(deffields)
if ~isfield(opts, deffields{i})
opts.(deffields{i})=defopts.(deffields{i});
end
end
if isfield(files,'M')
filenames=files.M;
elseif isfield(files,'R')
filenames=files.R;
elseif isfield(files, 'info3D')
opts.flag3D=1;
else
error('No filenames (M nor R) found.');
end
flag_signed=min(M(:))<0;
if ~opts.flag3D
if size(M,3) ~= size(filenames,1)
error([num2str(size(filenames,1)) ' filenames given for ' num2str(size(M,3)) ' slices.']);
end
end
if isempty(opts.SeriesInstanceUID)
opts.SeriesInstanceUID=dicomuid;
end
progress='';
mkdir(outdir);
warning('off','images:dicom_add_attr:invalidAttribChar');
load_flag=1;
insert_flag=~opts.flag3D;
for slice=1:size(M,3)
if load_flag
if opts.flag3D
filename=files.info3D.filename;
else
filename=filenames{slice,end};
end
info = dicominfo(filename);
info.SeriesDescription = opts.SeriesDescription;
info.SeriesInstanceUID = opts.SeriesInstanceUID;
if isempty(opts.SeriesNumber)
opts.SeriesNumber=info.SeriesNumber*100;
end
info.SeriesNumber = opts.SeriesNumber;
info.SOPInstanceUID = dicomuid;
info.InstanceNumber = slice;
if opts.flag3D
load_flag=0;
end
end
if opts.flag3D
item=files.info3D.slice2index{slice};
% info.PerFrameFunctionalGroupsSequence.(item).PlanePositionSequence.Item_1.ImagePositionPatient;
info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.SOPInstanceUID = dicomuid;
end
im = M(:,:,slice);
if (isfield(info, 'SmallestImagePixelValue'))
info.SmallestImagePixelValue=opts.Convert(min(im(:)));
end
if (isfield(info, 'LargestImagePixelValue'))
info.LargestImagePixelValue=opts.Convert(max(im(:)));
end
if (isfield(info, 'RescaleIntercept'))
info.RescaleIntercept=0;
end
if (isfield(info, 'RescaleSlope'))
info.RescaleSlope=1;
end
info.WindowCenter=opts.Convert(opts.Level);
info.WindowWidth=opts.Convert(opts.Window);
% if opts.flag3D
% info.PerFrameFunctionalGroupsSequence.Item_1.PixelValueTransformationSequence.Item_1.RescaleIntercept=0;
% info.PerFrameFunctionalGroupsSequence.Item_1.Private_2005_140f.Item_1.RescaleIntercept=0;
% info.PerFrameFunctionalGroupsSequence.Item_1.PixelValueTransformationSequence.Item_1.RescaleSlope=1;
% info.PerFrameFunctionalGroupsSequence.Item_1.Private_2005_140f.Item_1.RescaleSlope=1;
% info.PerFrameFunctionalGroupsSequence.Item_1.FrameVOILUTSequence.Item_1.WindowCenter=opts.Convert(opts.Level);
% info.PerFrameFunctionalGroupsSequence.Item_1.FrameVOILUTSequence.Item_1.WindowWidth=opts.Convert(opts.Window);
% end
info.SamplesPerPixel=1;
info.BitsAllocated=16;
info.BitsStored=16;
info.HighBit=15;
info.PixelRepresentation=flag_signed;
if size(M,3)==slice
insert_flag=1;
end
if insert_flag
outfile=fullfile(outdir,sprintf('IM%05d.dcm', slice));
print_progress(outfile);
if opts.flag3D
f=fieldnames(info.PerFrameFunctionalGroupsSequence);
f=setdiff(f,files.info3D.slice2index,'stable');
for i=1:length(f)
info.PerFrameFunctionalGroupsSequence=rmfield(info.PerFrameFunctionalGroupsSequence, f{i});
end
for i=1:length(files.info3D.slice2index)
item=files.info3D.slice2index{i};
info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.InstanceNumber=1;
info.PerFrameFunctionalGroupsSequence.(item).PixelValueTransformationSequence.Item_1.RescaleIntercept=0;
info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.RescaleIntercept=0;
info.PerFrameFunctionalGroupsSequence.(item).PixelValueTransformationSequence.Item_1.RescaleSlope=1;
info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.RescaleSlope=1;
info.PerFrameFunctionalGroupsSequence.(item).FrameVOILUTSequence.Item_1.WindowCenter=opts.Convert(opts.Level);
info.PerFrameFunctionalGroupsSequence.(item).FrameVOILUTSequence.Item_1.WindowWidth=opts.Convert(opts.Window);
end
info.NumberOfFrames=length(files.info3D.slice2index);
sz=size(M);
M=reshape(opts.Convert(M), sz(1), sz(2), 1, []);
M=permute(M, [2 1 3 4]);
if isfield(files, 'slices_added')
if files.slices_added
warning('Removing empty slice at bottom of volume');
M=M(:,:,1:end-1);
end
end
%dicomwrite(M,outfile, ...
% 'CreateMode', 'copy', 'WritePrivate', true, info);
dicomwrite(M,outfile, ...
'WritePrivate', true, info);
else
if isfield(files, 'slices_added')
if files.slices_added
warning('Removing empty slice at bottom of volume');
M=M(:,:,1:end-1);
end
end
%dicomwrite(opts.Convert(M(:,:,slice))',outfile, ...
% 'CreateMode', 'copy', 'WritePrivate', true, info);
%
dicomwrite(opts.Convert(M(:,:,slice))',outfile, ...
'WritePrivate', true, info);
end
end
end
fprintf('\n');
function print_progress(arg)
num=length(progress);
num=num-numel(regexp(progress, '\\\\'));
for ii=1:num; fprintf('\b'); end
progress=['Writing file ' arg];
progress=regexprep(progress,'\','\\\\');
fprintf(progress);
end
function y=convert(x)
if flag_signed
y=int16(x*1000);
else
y=uint16(x*1000);
end
end
end
2.4. QSM_3Dを集める
各フォルダには、強度画像と位相画像が入っている。
TEを7回変えて撮像している。
国際脳QSMの撮像は、1回撮像に128 slicesなので強度画像と位相画像はそれぞれ、896枚(=128 slices ×7 phase)
(Cornell MRI Research Lab. QSMより、URL: http://pre.weill.cornell.edu/mri/pages/qsm.html)
QSM_3Dにある強度画像と位相画像をanalyzeフォルダのrawdata(※必ずrawdata)にまとめて保存。
3. 実行
3.1. 作業フォルダ(analyzeフォルダ)へ移動
MATLABを起動し、analyzeフォルダへ移動。
3.2. runme.mを実行
runme.mをMATLABへDrag & Dropし、QSM解析を実行。
3.3. Outputファイルの確認
runme.mを実行後、
- QSM
- results
フォルダが生成される。
QSMフォルダにQSM画像がDICOM形式で保存される。
(Cornell MRI Research Lab. QSMより、URL: http://pre.weill.cornell.edu/mri/pages/qsm.html)
4. 複数の被験者データをまとめて処理したい場合
被験者ごとにフォルダを作成し一つのフォルダにまとめます。
さらに、後で紹介するまとめてMEDIを実行するスクリプトも入れておきます。
各被験者フォルダには強度画像と位相画像のDICOMがまとめて入ったrawdataフォルダがあります。
以下のスクリプトを実行すると、すべての被験者のQSM mapが計算できます。
4.1. run_sequential.m
clear all;clc;close all;
% change direct to study folder
selpath = uigetdir('Select the folder including all subject');
cd(selpath);
% list folder
folderlist = dir(selpath);
folderlist = folderlist(~ismember({folderlist.name}, {'.', '..'}));
folderlist = folderlist([folderlist.isdir]);
% run MEDI
for i = 1:length(folderlist)
disp(['processing ' folderlist(i).name ' ...'])
main(folderlist(i).name, selpath);
end
% define function of MEDI processings
function main(folder, basepath)
% initialize
clearvars -except selpath folderlist folder basepath i
% move to subject folder
cd(folder)
% STEP 1: Import data
[iField,voxel_size,matrix_size,CF,delta_TE,TE,B0_dir,files]=Read_DICOM('rawdata');
%%%% Generate the Magnitude image %%%%
iMag = sqrt(sum(abs(iField).^2,4));
%%%%% provide a Mask here if possible %%%%%%
if (~exist('Mask','var'))
% Mask = genMask(iField, voxel_size);
Mask = BET(iMag,matrix_size,voxel_size);
end
%%%%% provide a noise_level here if possible %%%%%%
if (~exist('noise_level','var'))
noise_level = calfieldnoise(iField, Mask);
end
%%%%% normalize signal intensity by noise to get SNR %%%
iField = iField/noise_level;
% STEP 2a: Field Map Estimation
%%%%%Estimate the frequency offset in each of the voxel using a
%%%%%complex fitting %%%%
[iFreq_raw N_std] = Fit_ppm_complex(iField);
% STEP 2b: Spatial phase unwrapping %%%%
iFreq = unwrapPhase(iMag, iFreq_raw, matrix_size);
% STEP 2c: Background Field Removal
%%%% Background field removal %%%%
[RDF shim] = PDF(iFreq, N_std, Mask,matrix_size,voxel_size, B0_dir);
% CSF Mask for zero referencing
R2s=arlo(TE,abs(iField));
Mask_CSF = extract_CSF(R2s, Mask, voxel_size);
% STEP 3: Dipole Inversion
save RDF.mat RDF iFreq iFreq_raw iMag N_std Mask matrix_size...
voxel_size delta_TE CF B0_dir Mask_CSF;
%%%% run MEDI %%%%%
QSM = MEDI_L1('lambda',1000, 'smv', 5, 'merit');
%%% save to DICOM
% ignore warnings...
Write_DICOM(QSM,files,'QSM')
% back to study folder
cd(basepath)
end
著者情報： 斎藤 勇哉
順天堂大学医学部 大学院医学研究科 放射線診断学講座所属脳MRI 画像解析が専門であり、テーマは①神経変性疾患の機序解明、②医用人工知能の開発、③多施設データのハーモナイゼーション、④速読が脳に与える影響や学習効果、⑤SNS解析を用いたマーケティング戦略の改善。
医療分野に関わらず、自然言語処理・スクレイピング・データ分析・Web アプリ開発を得意とし、企業や他大学の研究を支援。
主な使用言語は、Python、Shell Script、MATLAB、HTML、CSS
【FSL】Voxel-BasedAnalysis: VBA
1. 目的
2. VBAとは
2.1. ファイルの準備
2.2. 脳解剖学的標準化
2.3. 平滑化
2.4. GLMと並べ替え検定（permutation test）
3. おまけ
1. 目的
- Voxel-BasedAnalysis: VBA
2. VBAとは
Voxel-Based Analysis（VBA）は、脳構造解析手法の一つであり、現在では脳科学の分野において幅広く用いられている。VBAは、観測者が関心領域を設定して解析する古典的なマニュアル計測とは異なり、自動処理によって全脳をボクセルごとに解析するため、評価の客観性が高い。
VBAの解析手順は、次の通り。
- ファイルの準備
- 脳解剖学的標準化
- 平滑化
- GLMと並べ替え検定（permutation test）
2.1. ファイルの準備
VBA解析をしたファイルを準備する。
ここでは、拡散MRIから計算したFA画像を準備する。健常者10名（ID: Con0001~Con0010）と患者10名（ID: Pat0001~Pat0010）いることを想定する。
. ├── Con0001_FA.nii.gz ├── Con0002_FA.nii.gz ├── Con0003_FA.nii.gz ... └── Pat0010_FA.nii.gz
2.2. 脳解剖学的標準化
非線形変換を用いて、全ての被験者のFA画像を標準FA画像（FMRIB58_FA）に位置合わせする。
以下のコマンドを実行する。
mkdir preprocessing
for SUBJ in $(ls | grep _FA.nii.gz | cut -d . -f1); do
# 位置合わせ
## Affine transform
flirt -in ${SUBJ}.nii.gz -ref ${FSLDIR}/data/standard/FMRIB58_FA_1mm.nii.gz \
-dof 12 -omat preprocessing/${SUBJ}_indiv2std.mat
## non-linear transform
fnirt --in=${SUBJ}.nii.gz --aff=preprocessing/${SUBJ}_indiv2std.mat --config=FA_2_FMRIB58_1mm.cnf \
--cout=preprocessing/${SUBJ}_warp_indiv2std.nii.gz --iout=preprocessing/${SUBJ}_instd.nii.gz
done
位置合わせ方法の詳細は、こちらの記事を参考に。
https://www.nemotos.net/?p=4513
2.3. 平滑化
非線形変換を用いた位置合わせで、脳解剖学的標準化をしたとしても、ボクセル単位でみると完全に標準化できているわけではない。そこで、そのようなエラーを抑えるために平滑化処理をする。
まず、標準FA画像に位置合わせをした、被験者ごとのFA画像をまとめて4D-FA画像（all_FA.nii.gz）とし、4D-FA画像の平均画像からマスク画像を生成する。
mkdir stats fslmerge -t stats/all_FA.nii.gz preprocessing/*_instd.nii.gz fslmaths stats/all_FA.nii.gz -Tmean -bin mask stats/all_FA_mask.nii.gz
次に、ガウシアンフィルタを用いた平滑化処理をする。
ここでは、よく用いられるσ=3 voxel（FWHM=約7mm）のガウシアンフィルタで平滑化する。
fslmaths stats/all_FA.nii.gz -fmean -s 3 stats/all_FA_smoothed_s3.nii.gz
2.4. GLMと並べ替え検定（permutation test）
ガウシアンフィルタで平滑化した灰白質画像を用いて、健常群と患者群の群間比較をする。
まず、GLMのデザインマトリックス（計画行列）とコントラストを設定する。
今回は、健常者10名（ID: Con0001~Con0010）と患者10名（ID: Pat0001~Pat0010）のデータがある。すべての被験者のFAをlsコマンドでみると、先に健常者10名のFA、次に患者10名のFAが並んでいることが分かる。
ls |grep nii
Con0001_FA.nii.gz
Con0002_FA.nii.gz
Con0003_FA.nii.gz
...
Pat0010_FA.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：統計値マップ
design_ttest2 stats/design 10 10
randomise -i stats/all_FA_smoothed_s3.nii.gz \
-m stats/all_FA_mask.nii.gz \
-o stats/fslvba \
-d stats/design.mat \
-t stats/design.con \
-n 10000 --T2 -V -x --uncorrp -R
次のようなファイルが、生成される。
TFCEを用いた”健常群＞患者群”の検定で、FWE補正をされたP値マップ（fslvba_tfce_corrp_tstat1.nii.gz）を確認する。ここで、得られたP値マップは1-P値のマップであることに注意する。つまり、P<.05を有意とするのであれば、P値マップで0.95-1.00の値を見ればよい。
fsleyes $FSLDIR/data/standard/MNI152_T1_2mm \
stats/fslvba_tfce_corrp_tstat1 \
-cm Red-Yellow -dr 0.95 1
3. おまけ
検定数が多い場合、有意差があったかどうかをすべて確認するのは大変である。そこで、一目で有意差があるかどうかを判断できるように、各検定ごとのP値マップの最大値を自動計測し、テキスト（corrected_P_report.txt）としてまとめる。
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
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
fslvba_tfce_corrp_tstat1.nii.gz 0.982391
fslvba_tfce_corrp_tstat2.nii.gz 0.993181
fslvba_tfce_corrp_tstat2.nii.gz 0.993181
fslvba_tfce_corrp_tstat1.nii.gz 0.982391
著者情報： 斎藤 勇哉
順天堂大学医学部 大学院医学研究科 放射線診断学講座所属脳MRI 画像解析が専門であり、テーマは①神経変性疾患の機序解明、②医用人工知能の開発、③多施設データのハーモナイゼーション、④速読が脳に与える影響や学習効果、⑤SNS解析を用いたマーケティング戦略の改善。
医療分野に関わらず、自然言語処理・スクレイピング・データ分析・Web アプリ開発を得意とし、企業や他大学の研究を支援。
主な使用言語は、Python、Shell Script、MATLAB、HTML、CSS