SPMのメーリングリストでSPMの中の人が明らかにしていますが、SPM12の正式版がそろそろ公開されるようです。SPM12が正式に公開されました。（2014年10月1日）
SPM12はいろいろよい点も加わっていますが、まだ、プラグインで対応していないものも多々あります。
このため、しばらくは、SPM8とSPM12が同時に動かせる環境にあった方がよいですよね。
いい方法がありますので、ご紹介します。
カテゴリ別アーカイブ： neuroimaging
RでのROC解析：ROCRパッケージを使ったROC曲線とAUCの求め方
研究でROC解析を行う必要があり、Rでどうやったらできるのか調べてみました。
そうしたところ、ROCRというパッケージが公開されており、比較的簡単にROC解析を行い、グラフを作成できることがわかりました。
- ROCRパッケージのインストール（Ubuntu)
既にRはインストールされているとします。Ubuntuの場合、ROCRパッケージはapt経由で簡単に入手できます。
$ sudo apt-get install r-cran-rocr
これでインストール完了です。
ROC解析に必要なものは、何らかの指標と、それが属するグループの一覧です。具体例を挙げると、以下のようになります。
第1列に指標、第2列に属するグループ（0か1)が記載されています。
これをroc_data.txtとという名前で保存することとします。保存したディレクトリをRのワーキングディレクトリとします。
0.9706 1 0.9572 1 0.4854 1 0.8003 1 0.1419 1 0.4218 1 0.9157 1 0.7922 1 0.9595 1 0.6557 1 0.0357 1 0.8491 1 0.934 1 0.6787 1 0.7577 1 0.7431 1 0.3922 1 0.6555 1 0.1712 1 0.706 1 0.4797 0 0.4551 0 0.0374 0 0.081 0 0.2984 0 0.7597 0 0.1404 0 0.3853 0 0.0238 0 0.5513 0 0.0551 0 0.306 0 0.4991 0 0.6909 0 0.7593 0 0.3472 0 0.0614 0 0.0507 0 0.0575 0 0.6407 0
ROCRはRを立ち上げた後に、library(ROCR)で起動できます。
$ R > library(ROCR)
先ほどのroc_data.txtをrocdataという変数に読み込みます。変数名は何でもいいのですが、ここではそうします。
read.tableという関数で表を読み込めるので、それを使います。
rocdata <- read.table("roc_data.txt")
ROCRのステップ1はpredictionで、値と属するグループを指定します。
今、rocdataは20行2列の行列になっています。1列目はrocdata[,1]で、2列目はrocdata[,2]であらわすことができますので、以下のように記載します。
pred <- prediction(rocdata[,1], rocdata[,2])
ここで、変数predはpredictionの頭文字4文字です。もちろん、別の名前でもかまいません。
私はよくカンマを忘れるので、カンマも忘れないようにしましょう。
ここで、何が行われているかというと、データを大きい順にソートし、真陽性(TP)、偽陽性(FP)、偽陰性(FN)、真陰性(TN) の数を算出します。
ROCRのステップ2はperformanceです。ここでは、感度、すなわち真陽性率 (TP/(TP+FN)で定義)と、1-特異度、すなわち偽陽性率(FP/(FP+TN)で定義)を求めます。以下のようにタイプします。
perf <- performance(pred, "tpr", "fpr")
ここでtprはtrue positive rateを、fprはfalse positive rateを意味します。
それでは、ROC曲線を描きます。非常に簡単です。
plot(perf)
そうすると、次のようなグラフが現れるはずです。
グラフをPNG形式で保存するには、次のように行うことで、roc-curve.pngという名前でワーキングディレクトリに保存されます。
png("roc-curve.png")
plot(perf)
dev.off()
先ほどのperformanceの際に”auc”と指定するとAUCも計算されます。ただ、1クッション入れる必要があります。具体的な方法はこちらのサイトに記載されていましたが、それを転載します。
auc.tmp <- performance(pred,"auc") auc <- as.numeric(auc.tmp@y.values)
performance(pred,”auc”)の結果をauc.tmpという変数に代入し、
auc.tmpの中からy.valuesの値を取り出して、その値を変数aucに代入します。
最後に変数aucを表示させてみます。
auc [1] 0.8
これでAUCが0.8だということがわかります。
ここからさらに一歩踏み込んで、正診率を求めたいと思います。
正診率、感度、特異度は以下で定義されます。
正診率＝(TP+TN)/総数
感度(真陽性率)＝TP/(TP+FN)
特異度＝TN/(FP+TN)
今、総数は、rocdataの行数を求めればよいですから、nrow(rocdata)で求められます
カットオフ値を少しずつずらした時に、TP, FP, FN, TNは変わっていきますので、その一覧を表に出力しましょう。
表の列は以下のようにしたいと思います
Cutoff TP FP FN TN Sensitivity Specificity Accuracy
table <- data.frame(Cutoff=unlist(pred@cutoffs), TP=unlist(pred@tp), FP=unlist(pred@fp), FN=unlist(pred@fn), TN=unlist(pred@tn), Sensitivity=unlist(pred@tp)/(unlist(pred@tp)+unlist(pred@fn)), Specificity=unlist(pred@tn)/(unlist(pred@fp)+unlist(pred@tn)), Accuracy=((unlist(pred@tp)+unlist(pred@tn))/nrow(rocdata)) )
これで、tableを表示させると、以下のように表示されます。
> table Cutoff TP FP FN TN Sensitivity Specificity Accuracy 1 Inf 0 0 20 20 0.00 1.00 0.500 2 0.9706 1 0 19 20 0.05 1.00 0.525 3 0.9595 2 0 18 20 0.10 1.00 0.550 4 0.9572 3 0 17 20 0.15 1.00 0.575 5 0.9340 4 0 16 20 0.20 1.00 0.600 6 0.9157 5 0 15 20 0.25 1.00 0.625 7 0.8491 6 0 14 20 0.30 1.00 0.650 8 0.8003 7 0 13 20 0.35 1.00 0.675 9 0.7922 8 0 12 20 0.40 1.00 0.700 10 0.7597 8 1 12 19 0.40 0.95 0.675 11 0.7593 8 2 12 18 0.40 0.90 0.650 12 0.7577 9 2 11 18 0.45 0.90 0.675 13 0.7431 10 2 10 18 0.50 0.90 0.700 14 0.7060 11 2 9 18 0.55 0.90 0.725 15 0.6909 11 3 9 17 0.55 0.85 0.700 16 0.6787 12 3 8 17 0.60 0.85 0.725 17 0.6557 13 3 7 17 0.65 0.85 0.750 18 0.6555 14 3 6 17 0.70 0.85 0.775 19 0.6407 14 4 6 16 0.70 0.80 0.750 20 0.5513 14 5 6 15 0.70 0.75 0.725 21 0.4991 14 6 6 14 0.70 0.70 0.700 22 0.4854 15 6 5 14 0.75 0.70 0.725 23 0.4797 15 7 5 13 0.75 0.65 0.700 24 0.4551 15 8 5 12 0.75 0.60 0.675 25 0.4218 16 8 4 12 0.80 0.60 0.700 26 0.3922 17 8 3 12 0.85 0.60 0.725 27 0.3853 17 9 3 11 0.85 0.55 0.700 28 0.3472 17 10 3 10 0.85 0.50 0.675 29 0.3060 17 11 3 9 0.85 0.45 0.650 30 0.2984 17 12 3 8 0.85 0.40 0.625 31 0.1712 18 12 2 8 0.90 0.40 0.650 32 0.1419 19 12 1 8 0.95 0.40 0.675 33 0.1404 19 13 1 7 0.95 0.35 0.650 34 0.0810 19 14 1 6 0.95 0.30 0.625 35 0.0614 19 15 1 5 0.95 0.25 0.600 36 0.0575 19 16 1 4 0.95 0.20 0.575 37 0.0551 19 17 1 3 0.95 0.15 0.550 38 0.0507 19 18 1 2 0.95 0.10 0.525 39 0.0374 19 19 1 1 0.95 0.05 0.500 40 0.0357 20 19 0 1 1.00 0.05 0.525 41 0.0238 20 20 0 0 1.00 0.00 0.500
Accuracyがもっとも高いところを見つけるには、
max(table$Accuracy)
とします。そうすると、今は、
> max(table$Accuracy) [1] 0.775
となりますので、該当するところをみると、感度70%、特異度85%、正診率77.5%達成できるということがわかりました。
感度と特異度の曲線も簡単に書けます。predictionまで行った後に、次のようにします。
perf <- performance(pred, "sens", "spec")
png("sens-spec-curve.png")
plot(perf)
dev.off()
これで、下図のような感度、特異度の曲線がsens-spec-curve.pngという名前で保存されます。
比較的簡単に求められるので便利です。
VBMミニ勉強会＠秋葉原 6月28日15:00-18:00
5人の方の申し込みをいただいたので、開催決定です！
岩手医大の山下典生先生も助っ人で参加してくださることになりました！
※場所が確定しました。秋葉原駅徒歩5分のAP秋葉原です。
※事前登録制です。参加希望の方はコメント欄からご連絡をお願いします。
すぐできるVBMが無事に出版に至りましたが、「本は買ったけれどもやる時間がない」という方もいらっしゃるのではないかと思います。また、「ここはどうなの？」と質問をしたいという方もいらっしゃると思います。
そこで、VBMのミニ勉強会を開催させていただきます。
日時：2014年6月28日（土）15:00-18:00
場所：AP秋葉原
参加費：3000円程度（会議室の実費代のみ。）
持参していただくもの：SPMがセットアップされている自分のパソコン
人数：最小5名〜最大20名
テキスト：「すぐできるVBM」（前もってご準備ください）
内容：実習形式
・VBMを実際に動かしてみる
・自分の解析でわからないことの質疑応答
参加希望の方は、下記コメントにお名前（ニックネーム可）とメールアドレスを書きこんでいただけますでしょうか。
VBMに関心のある方はご検討ください。
第1回 IMU脳画像解析セミナー「VBM／ROI解析を極める」@岩手
友人の山下典生先生が、2014年8月31日に『第1回 IMU脳画像解析セミナー「VBM／ROI解析を極める」』を開催します。
SPMでのVBMだけでなく、3D-SlicerやITK-SNAPも利用するとのこと。とても楽しみです。
早速申し込みしました。
関心のある方はこちらをごらんください。
VBMでの全脳容積(Total Brain Volume; TBV)の求め方
「すぐできるVBM」に、全脳容積の求め方を記載していませんでした…。
以下の方法で、求めることができます。
必要なものは以下のとおりです。
- get_totals.m
- 灰白質画像と白質画像
容積を簡便に求めることができるスクリプトは、Ged Ridgway氏によるget_totals.mです。右クリックで「名前をつけて保存」とし、SPMのディレクトリに保存してください。
Segmentで最初に作られる画像であるc1*.niiとc2*.niiを準備します。ここで*はワイルドカードを意味し、「c1からはじまり、そのあとは何でもOKで、最後は.niiで終わる」という意味です。c1画像は厳密には確率画像なのですが、標準化前のこの画像のボクセル値は、ほぼ容積とイコールであると考えられています。
それでは実際にやっていきましょう。
まずは、灰白質の容積を求めます。Matlabのコマンドウィンドウから以下のようにタイプします。
gmv = get_totals
(gmvは変数名なので、gray_matterでも、gmでもかまいません）
そうすると、ファイルを選択するダイアログがあらわれます。ここで、容積を求めたいc1画像を指定してください。ここでは、ためしに5つのファイルを選びました。そうすると、以下のような表示になると思います。（もちろん、数字は選択した画像で変わります）
gmv =
640.9410
664.4302
617.7647
874.1366
670.1864
次に、白質の容積を求めます。Matlabのコマンドウィンドウから以下のようにタイプします。
wmv = get_totals
同様に、ファイルを選択するダイアログがあらわれますので、容積を求めたいc2画像を指定します。重要なのは、c1画像と同じ順番で選択してください。以下のような表示になると思います。
wmv =
431.9188
463.1192
418.5571
572.9889
488.4050
そうしたら、2つを足して全脳容積を求めましょう。
tbv = gmv + wmv
tbv =
1.0e+03 *
1.0729
1.1275
1.0363
1.4471
1.1586
これは、1.0×10^3（すなわち1000）にそれ以下の数字をかけたものということです。つまり、
tbv =
1072.9
1127.5
1036.3
1447.1
1158.6
ということになります。
これを選択してコピーしても使えますが、ファイルに書き出しましょう。
コンマ区切りファイルというものに出力することができます。
出力ファイル名をtotal_brain_volume.csvとします。
csvwrite('total_brain_volume.csv',tbv)
csvwriteの書式は csvwrite(出力ファイル名, 変数名) となります。
変数tbvの内容をtotal_brain_volume.csvに出力したいので、上記のようになります。
そうすると、ワーキングディレクトリにtotal_barin_volume.csvが生成され、Excelのような表計算ソフトでも、テキストエディタでも開くことができます。
最後に少し応用しましょう。灰白質容積(GMV), 白質容積(WMV), 全脳容積(TBV)をすべてCSVファイルに出力したいとしましょう。
その場合は、まず、すべての情報をひとつの変数volumesにまとめます。
Matlabのコマンドウィンドウから次のようにタイプしてください。
volumes = [gmv wmv tbv]
そうすると下記のようになります。
volumes =
1.0e+03 *
0.6409 0.4319 1.0729
0.6644 0.4631 1.1275
0.6178 0.4186 1.0363
0.8741 0.5730 1.4471
0.6702 0.4884 1.1586
左からgmv, wmv, tbvという順番になっています。
これを、先程と同じようにしてCSVファイルに保存します。ファイル名をvolumes.csvとします。
csvwrite('volumes.csv',volumes)
このようにして、脳容積の情報を保存できます。
SPMでの「1対多」の解析
臨床では、1人の症例で健常者と比較して容積低下部位を求めたいことがあります。
SPMでどのようにするかというと、シンプルです。
- DesignはTwo sample t-testを用います
- Group 1に1例のみ選びます
- Group 2に健常者の複数の画像を指定します
- IndependenceはYesを選びます
- （ここがキモ）VarianceはEqualを選びます。もし、Unequalにするとエラーとなります。1例では分散を計算することができないからです。
そうすると、Design matrixを作るための設定画面は下図のようになるはずです。
これを普通にEstimateして、Contrast managerで症例の容積が低下している領域を求めたかったら、[-1 1]とすれば大丈夫です。
Twitterで疑問を投げかけてくださった@silverjet_jpさん、ありがとうございました。
MRIConvertやdcm2niiを用いたDICOM→NIFTI変換の方法
画像解析のためには、DICOM画像をNIFTI画像に変換する必要があります。
主に使われるソフトウェアとして、MRI Convertとdcm2niiがあります。
開発者からスクリーンショットを使ってよいという許可をいただきましたので、
MRI Convertおよびdcm2niiの使用方法をPDFで公開します。
Flexible factorial designを用いたVBMの縦断解析
様々な方々から「2群の縦断解析をしたいんだけれども、どうすればいいの？」と質問されてきて、その都度、きちんとした答えをできずにいました。
ちょうど、今日、ある方から質問をうけたので、これを機にまとめてみることにします。
一言で言うならば「Flexible factorial modelを使うといいよ」となります。しかし、これはそんな簡単ではありません。
以下にその方法をまとめます。
根拠となっている資料として、”Contrast weights in flexible factorial design with multiple groups of subjects“を挙げます。これは、SPMのML上で流れたドキュメントで、Flexible factorial designにおけるコントラストを様々な条件から検討しています。このドキュメントの他にも、SPM-ML上の議論も参考にしています。できるだけ正確を期していますが、間違っている可能性もゼロではありませんので、もし間違いを見つけた方がいらっしゃったらぜひご指摘ください。
例として、以下のような状況を想定します。上述のContrast weights in flexible factorial design with multiple groups of subjectsでとりあげられている状況と同じ状況です。
- 被験者は11人。6人(Subject1(S1)-S6)がグループ1、5人(S7-S11)がグループ2
- それぞれ3回MRIを撮像。従って、11×3=33の画像があることになる
以下、SPM12bのスクリーンショットを用いながら説明します。基本的にSPM8でも同じです。
SPMのメニューから”Basic models”をクリックします。
- Design matrixの作成
- Flexible factorialの指定とFactorの設定
まずは、DesignでFlexible factorialを指定します。その後、Factorを指定します。今、指定すべきFactorは3つです。各々に対してindependenceとvarianceを指定します。
- subject
- group
- time
これは必ず入れないといけないとのことです。independenceはyes, varianceはequalとなります。
次はgroupです。これはindependenceはyes, varianceはunequalとなります。
縦断解析ですので、timeが入ります。これはindependenceはno, varianceはequalとなります。
すぐに見直すことができるように表にしておきます。
factor independence variance subject yes equal group yes unequal time no equal - subjectの指定
- Conditionsの指定
- Main effectsとInteractionの指定
- Flexible factorial designのコントラスト
- groupの主効果
- timeの主効果
- group x timeの交互作用
Specify Subjects or all Scans & Factorsで、”Subjects”を選択します（all Scans & Factorsもいけそうなのですが、Design matrixがやや異なるので、ここは無難にSubjectsとします）。そして、Subjectsに、Subject 1のtime 1, 2, 3の3つのファイルを指定します。
Conditionは行列で指定します。行は画像の順番です。列は指定したFactorの順番です。Factor1のSubjectは今、ひとつひとつ指定することで考慮されていますので、groupとtimeを指定します。被験者1-6(S1-S6)はグループ1, S7-S11はグループ2ですので、S1-S6のconditionは次のようになります。
1 1
1 2
1 3
同様に、S7-S11のconditionは次のようになります。
2 1
2 2
2 3
その後、Main effects(主効果)とInteraction(交互作用)の指定をします。
subject, group, timeの主効果を見るために、Main effectとしてFactor 1, 2, 3をそれぞれ指定します。
さらに、group x timeの交互作用もみたいので、Interactionとして、2 3とします。(Factor2がgroup, Factor 3がtimeだからです)
それ以降の設定(mask, global normalization)は通常のVBMと同様です。
この設定を適当な名前(flex_fact_design_matrix)で保存し、実行して、次のようなDesign matrixが得られたらOKです。
ここで気を付けなくてはいけないのは、Subjectは最初に指定しても、Design matrix上では最後に来ることです。
最初の2列(赤)がgroupの主効果, 次の3列(黄緑)がtimeの主効果, 次の6列(青)がgroup x timeの交互作用, 最後の11列(紫)がsubjectの主効果となっています。
ここまでできたらEstimateしてください。
Design matrixの設定もややこしいのですが、コントラストはもっとややこしいことになっています。ここは上述の”Contrast weights in flexible factorial design with multiple groups of subjects”に従って説明します。
今、Design matrixの各列の要素を列挙すると以下のようになっています。
| G1 | G2 | T1 | T2 | T3 | G1T1 | G1T2 | G1T3 | G2T1 | G2T2 | G2T3 | S1 | S2 | … | S6 | S7 | S8 | … | S11 |
ここでGはgroup, Tはtime, Sはsubjectを意味します。
これを見ると、Interactionの中にgroupとtimeの要素が入っていることがわかります。コントラスト作成の際には、ここのところに気をつける必要があります。nの数によってコントラストに重みづけをしなければならないからです。
| G1 | G2 | T1 | T2 | T3 | G1T1 | G1T2 | G1T3 | G2T1 | G2T2 | G2T3 | S1 | S2 | … | S6 | S7 | S8 | … | S11 |
| 1 | -1 | 0 | 0 | 0 | 1/3 | 1/3 | 1/3 | -1/3 | -1/3 | -1/3 | 1/6 | 1/6 | … | 1/6 | -1/5 | -1/5 | … | -1/5 |
ここでG1T1, G1T2, G1T3は皆等分なので、シンプルに3で割っただけです。G2T1なども同様に-1/3となります。S1-S6, S7-S11は人数にあわせて1/6と-1/5となっています。
ただ、これをひたすらSPMに打つのは非常に苦痛ですね。ここで、Matlabの便利な機能を使います。
Matlabでは、1がn個続く時は、 ones(1,n) とします。
これは応用が効いて、もし、1/3がn個続く時は、 ones(1,n)/3 とします。
-1/4がn個続く時は、-ones(1,n)/4 となります。
そして、ゼロがn個続くときは zeros(1,n) となります。
従って、今の場合は、以下のように書くことができます。私はこれをSPMのcontrast managerに入力しています。
1 -1 zeros(1,3) ones(1,3)/3 -ones(1,3)/3 ones(1,6)/6 -ones(1,5)/5
すると、下図に示すようなコントラストになるはずです。
| G1 | G2 | T1 | T2 | T3 | G1T1 | G1T2 | G1T3 | G2T1 | G2T2 | G2T3 | S1 | S2 | … | S6 | S7 | S8 | … | S11 |
| 0 | 0 | 1 | 0 | -1 | 6/11 | 0 | -6/11 | 5/11 | 0 | -5/11 | 0 | 0 | … | 0 | 0 | 0 | … | 0 |
timeの主効果は、経時的にvolumeが小さくなると仮定して、[1 0 -1]と考えています。G1T1, G1T2, G1T3に関しては、今はT1がG1とG2で不均等になっていますので、その比率を考慮して、G1T1, G1T2, G1T3に関しては、6/11×[1 0 -1], G2T1, G2T2, G2T3に関しては、5/11×[1 0 -1]となっています。そして、subjectのfactorに関してはこの中に3つの時系列をすべて包含していますので、コントラストに影響しませんので、すべて0となります。
これも先ほどと同様にonesとzerosを上手に使って書いてみましょう。以下をcontrast managerに入力します。
zeros(1,2) [1 0 -1] 6/11*[1 0 -1] 5/11*[1 0 -1] zeors(1,11)
すると、コントラストは下図のようになるはずです。
| G1 | G2 | T1 | T2 | T3 | G1T1 | G1T2 | G1T3 | G2T1 | G2T2 | G2T3 | S1 | S2 | … | S6 | S7 | S8 | … | S11 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | -1 | -1 | 0 | 1 | 0 | 0 | … | 0 | 0 | 0 | … | 0 |
交互作用で用いられる項目は、GxTxのみです。これも次のように書きなおすことができます。
zeros(1,2) zeros(1,3) [1 0 -1] [-1 0 1] zeros(1,11)
コントラストは下図のようになります。
これでひと通り検定することができるはずです。
いざまとめてみると、結構すっきりした感じになりました。
縦断データをお持ちの方は、是非トライしてみてください。VBMだけでなく、脳血流SPECTなどにも応用可能です。
すぐできるVBM: 精神・神経疾患の脳画像解析 SPM12対応
ついにできました！VBMの教科書です。
薄い本（150ページ弱です）ですが、VBMについての基礎知識がぎっちり詰まっています。
この本のいいところは、以下のようなところです。
- 日本初の形態MRIの解析法であるVoxel-based morphometry (VBM)の方法を丁寧に紹介した本です。VBMだけを取り上げた本は世界にもないので、ある意味世界初です。
- これから正式版が公開されるSPM12をもとに解説しています。SPMは4年ごとにアップデートされますが、ここ最近は正式版の公開は遅れています。2014年なのにまだSPM12の正式版がリリースされていないわけですから。したがって、今購入していただければ、最低でも5年間は使えます。SPMの使い方の基本も説明していますので、これからSPMを使ってみたいという方にもオススメです。
- サンプルデータがついています。自分のデータセットがなくても、灰白質の抽出(Segmentation)から統計解析まで実際に行なってみることができます。
- スクリーンショットが豊富です。すべて私が自分でスクリーンショットをとっているので、大きな間違いはないはずです。
- 動作検証をWindowsとMacで行なっています（Linuxでも行なっています）。なので、WindowsユーザーでもMacユーザーでもきちんとできます。
- VBMを応用した認知症診断補助ソフトウェアVSRADについても開発者の松田先生自らの解説が書かれています。
構成は以下の通りになっています。
第1章 はじめに
第2章 画像解析の準備
2.1 画像解析に必要なコンピュータのスペック
2.2 画像のファイル形式
2.3 画像解析の一連の流れ
2.4 SPMのインストール
2.5 サンプルデータ
第3章 SPMの動作の基本
3.1 SPMのウィンドウ構成
3.2 作業ディレクトリの設定
3.3 画像の表示
3.4 AC-PC位置合わせ
3.5 MATLABスクリプトを用いた半自動AC-PC補正プログラム
第4章 VBMの概要と前処理
4.1 VBMとは
4.2 分割化
4.3 DARTEL
4.4 バッチ（Batch）処理
第5章 統計モデルと結果表示
5.1 事前準備
5.2 群間比較（two-sample t-test）
5.3 相関解析
5.4 要因の分散分析
第6章 VBMの臨床への応用
6.1 VSRAD®
6.2 J-ADNI
6.3 VSRAD®の臨床応用例
ちなみに、総ページのだいたい半分が第5章に割かれています。VBMの前処理についての情報はそれなりにあるのですが、統計処理のところについての説明は非常に少ないので、統計処理について重点を置いています。
値段が5000円＋消費税8％で5400円となっていますが、ここまで情報が詰め込まれているVBMの解説書はありませんので、形態MRI解析を志していて英語の文献を読むだけの気力がない方にはお手頃だと思います。（実は英語でもVBMの統計の部分の情報は系統的なものは少ないです…）
Amazonでも購入いただけます。リンクは下記。
UbuntuへのConnectome mapperのインストール方法
Connectomeに関心が集まっている昨今、様々なソフトウェアが発表されています。
それらの中にConnectome Mapping Toolkitがあります。Connectome Mapperがその根幹をなすものですが、Connectome mapperを動かすためには、
- Diffusion Toolkit
- Freesurfer
- FSL
- MRTrix
- Camino
- Camino-Trackvis
- MITK
- Nipype
を事前に準備しなくてはならず、そのうえでConnectome mapperを入れないといけません。
結構大変な作業です。
順天堂放射線科の神谷先生と情報交換をしながら、Ubuntuにこれらのソフトをいれ、Connectome mapperを起動させる手順書を作成しました。
悩んでいる人は試してみる価値があると思いますので、試してみてください。
PDFはこちらからダウンロードできます。
なお、Lin4Neuroの最新版(20140313)には、FreesurferとDiffusion Toolkitを除くこれらのソフトのすべてが入っています。めんどくさがりの方は、最新版のLin4Neuroを使用することをおすすめします。
FSLユーザーのためのシェルスクリプトチュートリアル
FSLを開発しているFMRIBグループは、FSL courseのスライドおよびデータセットを公開してくださっています。ここで勉強しているときに、さりげなくシェルスクリプトのチュートリアルもあることが気になっていました。
先日、時間をとってこのチュートリアルを丁寧に見ていったところ、非常に勉強になることが多く、シェルスクリプトの理解が深まりました。FSLユーザーのみならず、脳画像解析でUNIX系のOS(Macも含む)を使う方々にとって有益な情報が多いと思いましたので、日本語にしておきたいと思い、FMRIBのSteve Smith教授の承諾を得て、このチュートリアルを翻訳させていただきました。
オリジナルのサイト自体は、スライド形式になっており、レクチャーと演習部分にわかれています。しかし、翻訳していく段階でこれを完全に模倣するよりは、レクチャーと演習をそれぞれ別の媒体として準備するのも勉強するのには都合がいいかもと考え、それぞれ別のPDFを準備しました。
関心のある方は、以下のリンクからPDFをご覧いただけたらと思います。
FSLViewで画像が正しい向きで表示されない時
FSLViewで画像を表示すると、思ったような向きで表示されないときがあります。
下に例を示します。
ここでは、上にある水平断、冠状断が角度が変になってしまっています。
しかし、よくみると、FSLViewは画像の前後をきちんと認識しています。
画像のまわりにあるアルファベットはそれぞれ
A: anterior
P: posterior
L: left
R: right
I: inferior
S: superior
を意味しますが、正しくなっています。
それでは、「正しく」表示するにはどうしたらよいでしょうか。
続きを読む
SPMで画像の中心を原点に設定するスクリプト
岩手医大の山下典生先生スクリプト第２弾です。
NIfTI画像を扱っていると、ときどき画像の原点(Origin)が変なところに設定されているケースに遭遇することがあります。このスクリプトは、SPMで、NIfTI画像の原点を、画像の中心に設定するスクリプトです。多くの画像の原点が変なところになっているような時に役立つと思います。
使い方は簡単です。
- setorigin_center.mをダウンロード（右クリックで名前をつけてリンク先を保存）し、SPMのフォルダに入れます。
- Matlabのコマンドウィンドウから、以下のようにタイプします。
- 画像選択のウィンドウが開きますので、原点を画像の中心にセットしたい画像を選択してください。
>>setorigin_center
これだけです。
実際に手元にあったデータでやってみました。
左側が設定前です。原点がかなり後ろに設定されています。
右側がスクリプトを走らせた後です。この場合は、幸運にも画像の中心がちょうどACと合致しています。ここから微調整でAC-PCを設定した方がずっと楽ですよね。
SPM8でSPECT/PET画像を一括で解剖学的標準化するスクリプト
SPMでSPECT画像やPET画像を標準化するとき、Normalise (Est & Write)から行います。しかし、この時、Source ImageとImages to Writeにいちいち画像ファイルを指定せねばならず、苦痛を伴います。
私の友人の岩手医大の山下典生先生が、この問題を解決するスクリプトを書いて下さいました。
使い方はいたって簡単です。
- fy_normalise.mをダウンロード（右クリックで名前をつけてリンク先を保存）し、SPMのフォルダの中に保存します。
- Matlabから以下のようにタイプします。
- すると画像を選択するダイアログが出現しますので、標準化したいファイルを選択します。
- 次に、ディレクトリを聞かれますので、設定ファイルを保存したいディレクトリを指定します。
- そうすると、保存先にspm_normalise_20140117_1621112.matのようなファイルが保存されます。
- 次に、SPMのBatchを起動し、ファイルを開くから先ほど保存したmatファイルをしてします。
- そうすると、Subjectが指定された状態になります。
- パラメータを確認し、大丈夫であることを確認したら、実行すれば解剖学的標準化がはじまります。
fy_normalise
非常に便利なので、関心のある方はどうぞお使いください。なお、fy_normalise(”,1)とすると、一度に標準化を走らせることもできますが、設定ミスで後悔しないように、設定を確認してから走らせることをお勧めします。
山下先生、ありがとうございました！
VBM8の前処理のベンチマーク
VBMの前処理にはどのくらいの時間がかかりますか？という質問を受けましたので、様々なマシンで計測してみました。
3人の3D-T1 MRI画像（ボクセルサイズ 1x1x1mm）をVBM8でDARTELでノーマライズしています。
前処理の設定をvbm8_preproc.matという名前で保存して、以下のベンチマークスクリプトを用意しました。とても単純なMATLABのスクリプトで、処理の前後にtic, tocというコマンドを入れることで前処理にかかった時間を測定するというものです。
%vbm8_benchmark.m
tic
spm('defaults','pet');
spm_jobman('initcfg');
load vbm8_preproc.mat;
spm_jobman('run',matlabbatch);
toc
その結果は、以下のようになりました。
| Machine | CPU | Core | Thread | Memory | OS | Matlab | Processed time |
| Workstation BTO | AMD Opteron Processor 2425HE x 2 | 12 | 12 | 24GB | Xubuntu 12.04 | R2013b | 74m31s |
| Desktop BTO | Core i7-960 3.2GHz | 4 | 8 | 24GB | Xubuntu 12.04 | R2011b | 34m47s |
| Desktop BTO | Core i7-960 3.2GHz | 4 | 8 | 24GB | Xubuntu 12.04 | R2013b | 33m20s |
| Mac Book Pro (Late 2011) | Core i7-2640M 2.8GHz | 2 | 4 | 8GB | Mac OSX 10.9 | R2011b | 30m29s |
| Thinkpad T430s | Core i7-3520M 2.9GHz | 2 | 4 | 16GB | Windows 7 | R2013a | 35m0s |
| Thinkpad T430s | Core i7-3520M 2.9GHz | 2 | 4 | 16GB | Xubuntu 12.04 | R2013b | 29m12s |
| Desktop BTO | Core i7-3930K 3.2GHz | 6 | 12 | 32GB | Xubuntu 12.04 | R2013b | 26m52s |
この結果をまとめると、以下のようになります。
- CPUが一番影響します。新しいCPU、かつ1コアあたりの性能が高いCPUだとよい結果となります。
- これはSPM全体に言えることですが、マルチコアは効果はたいしてありません。逆に言えば、SPMを複数立ち上げて処理させることで処理時間を短縮できる可能性があります。
- メモリは8GB以上あれば、さほど変わりありません。
- WindowsとLinuxを比べるとLinuxに軍配があがります。
- Matlabもバージョンが最新の方が処理速度が改善されているようです。
ということで、SPM-VBMという観点からでは、1コアあたりの周波数ができるだけ速いCPUを入手することをお勧めします。
Non-stationary correction for VBM on SPM8
Below is an excerpt from SPM-ML
(https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=SPM;72ea5d13.1109)
There is an option for the newest SPM8 releases (>=4010), which enables non-stationarity
correction on the cluster level. Please set in spm_defaults.m:
defaults.stats.rft.nonstat = 1;However, this option will result in a much longer calculation if you print the results table.
This is maybe also the reason that this option is not set as default.
SPMで複数の群の比較を行いたい時のDesign Matrix
SPM-MLに勉強になる話題が流れていたので、共有します。
出典はこちら。
https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=SPM;a65dd354.1307
コントロールAと疾患B、疾患Cという3群がある時、Design Matrixに全部を入れてしまった方がいいのか、
もしくは、コントロールと疾患A、コントロールと疾患Bというように、別々にDesign Matrixを作った方がいいのかという質問です。
これに対し、Cyril Pernetが非常にわかりやすいたとえを使って説明をしています。
ただ、少しだけtypoがあったのでそこを補足して説明します。
VBMの観点からのSPMの理解（一般線型モデルの理解）
2013年9月に、包括脳MRIチュートリアルで、「VBMの観点からのSPMの理解」なる話をさせていただきました。
VBMに限らず、SPMで用いられる一般線形モデル（GLM）の基本をまとめてみました。
SPMユーザには多少はお役に立てる内容かと思いますので、講義の配布資料を公開させていただきます。
Non-stationary correction
Christian Gaser made the following comment on non-stationary correction in SPM8.
Furthermore, SPM8 (release >= 4010) provides a non-stationarity correction, that is based on smoothed RPV-maps, which should be more reliable than the unsmoothed RPV maps in VBM8 or in the NS-toolbox. To enable this function set:
defaults.stats.rft.nonstat = 1;
in spm_defaults. However, this option will result in a much longer calculation if you print the results table.
This is maybe also the reason that this option is not set as default.