SPMの結果から、自動で各座標の解剖学的名称を抽出するスクリプト

SPMの結果から領域名を抽出するのに苦労している方は多いと思います。
SPM12からは、GUIを用いて簡単に同定することはできるようになりましたが、
(この方法を知りたい方は、こちらの記事を参照してください)
それでも何十箇所もある場合、マウスで一つ一つ確認するのは骨が折れます。

先日、SPMのMLで、Guillaume Flandinがこのようなメールを流していたことに気づきました。

From within SPM, you can query the Neuromorphometrics atlas given a list
of coordinates:

XYZmm = [
0.5619 -9.1031 -4.7629
1.0137 -24.1197 9.3815
-44.0259 -22.1028 8.9601]’;

for i=1:size(XYZmm,2)
spm_atlas(‘query’,’neuromorphometrics’,XYZmm(:,i))
end

ans =

‘3rd Ventricle’
‘CSF’
‘Left Cerebral White Matter’

任意の座標を指定するならば、SPMに搭載されているNeuromorphometricsのアトラスを用いて、spm_atlasで領域名を求められるよと教えてくれています。

それならば、「SPMの結果から自動で得られるスクリプトを作りたい!」と思いました。

以下のようなCSVが自動で入手できたら嬉しいと思ったわけです。

とりあえずスクリプトだけ欲しい人は一番下までスクロールしてください。

少し勉強したい方は以下を読んでみてください。

例として、以下にあげるような3つの領域の座標を自動で求めたいと思います。

SPMで結果を表示した際、その座標がどこに保存されているかを調べてみたところ、

TabDat

という構造体の中の TabDat.dat に入っていました。

TabDat.dat

ans =

  3×12 の cell 配列

  1 列から 6 列

    {[1.2686e-09]}    {[       1]}    {[1.2686e-09]}    {[4.2618e-09]}    {[     443]}    {[8.6977e-11]}
    {0×0 double  }    {0×0 double}    {0×0 double  }    {0×0 double  }    {0×0 double}    {0×0 double  }
    {0×0 double  }    {0×0 double}    {0×0 double  }    {0×0 double  }    {0×0 double}    {0×0 double  }

  7 列から 12 列

    {[0.0209]}    {[0.0847]}    {[5.0788]}    {[4.9589]}    {[3.5453e-07]}    {3×1 double}
    {[0.0447]}    {[0.0877]}    {[4.8931]}    {[4.7851]}    {[8.5450e-07]}    {3×1 double}
    {[0.0953]}    {[0.1015]}    {[4.6963]}    {[4.6001]}    {[2.1109e-06]}    {3×1 double}

ここに3行12列のセルがあります。3行は、今、求めたい座標の3行です。12列は、SPMの結果の表が12列になっています。座標は第12列に入っているということです。

12列目だけ表示してみます。

TabDat.dat(:,12)

ans =

  3×1 の cell 配列

    {3×1 double}
    {3×1 double}
    {3×1 double}

こうなりますが、cell2matというコマンドをいれるとデータが見えてきます。

 A=TabDat.dat(:,12);
 B=cell2mat(A)

B =

    54
   -15
    18
    39
    33
     6
    30
    21
     3

今、Bは3つの座標が縦に並んでしまいました。
これだと使い勝手が悪いので、Aを転置させてみます。

B=cell2mat(A')

B =

    54    39    30
   -15    33    21
    18     6     3

今度は、座標らしくなりました。今、この座標は、(54, -15, 18) (39, 33, 6) (30, 21, 3) が縦に並んでいることに注意してください。縦ベクトルの集合です。

spm_atlasは、SPMである座標の解剖学的名称を得たい時に使用します。まさに今回のような目的に使用するわけですね。
使い方は、

spm_atlas(‘query’,’アトラス名’,座標の縦ベクトル)

となっています。

で、この際、座標の行列は縦ベクトルでなければいけません。なので、先程の形式は非常に都合がよいことになります。

早速、(54 -15 18)の名称を調べてみましょう。
このベクトルは、B(:,1)で表すことができます。

spm_atlas('query','Neuromorphometrics',B(:,1))

ans =

    'Right CO central operculum'

求められました。

それでは一気に求めたいと思います。Forループを使えばできます。
今、Bの列数は size(B,2)で求められます。

なので、以下のようになります。

for i=1:size(B,2)
    spm_atlas('query','Neuromorphometrics',B(:,i))
end

ans =
    'Right CO central operculum'

ans =
    'Right Cerebral White Matter'

ans =
    'Right Cerebral White Matter'

求められました!

これをベースに、以下のスクリプトを作成してみました。
SPMで結果を表示した後に、get_spm_names.m を実行すれば、region_names_タイムスタンプ.csv というCSVファイルが作成されます。

%%%% get_spm_names.m
%%%% A script to obtain region names from the SPM result
%%%% K. Nemoto 06 Sep 2019

% mni coordinates are stored in the 12th column of TabDat.dat
XYZmm = cell2mat(TabDat.dat(:,12)');

% enter region names in a cell
region = {};
for i = 1:size(XYZmm,2)
    region{i,1} = spm_atlas('query','neuromorphometrics',XYZmm(:,i));
end

% convert the cell which includes region names to a table
region_T = cell2table(region);

% generate a table with coordinates
coord=XYZmm';
x = coord(:,1);
y = coord(:,2);
z = coord(:,3);
coord_T = table(x,y,z);

% combine tables
T = [coord_T region_T];

% generate a filename
timestamp = datestr(now,'yyyy-mm-dd');
fname = ['region_name_' timestamp '.csv'];

% write CSV files
writetable(T,fname)

関心のある方は、こちらのリンクからどうぞ。(右クリックして「リンク先を保存」で保存してください)

便利さという点では、かなり便利なスクリプトになったのではないかと思います。

SPMの結果から、自動で各座標の解剖学的名称を抽出するスクリプト” へのコメント

  1. 初めまして、千葉大の佐原と申します。
    タスクfMRIの解析をしておりまして、こちらで共有していただいたコマンドを0から作ろうとしていたところでしたので、大変助かりましたし、matlabの勉強にもなりました。

    お聞きしたい点として、こちらで指定されているatlasを他のatlasに変更(例えばconnで用いられるnetworkやYeoのatlasに変更)することは可能でしょうか。

    • ご質問に対する回答としては、「適切なファイルが準備できたら可能です」となります。

      neuromorphometrics アトラスは

      labels_Neuromorhometrics.nii
      labels_Neuromorphometrics.xml

      から構成されています。

      このxmlがポイントかなと思います。

      こちらを参考にされながら、他のアトラスもxmlを作って、spm12/tpm の中にいれたらできるのではないかなと思います。(試したことがないのでそこまでの確証はないですが)

      • ありがとうございます。
        xmlファイルを新たに作成することが出来るか、試してみようと思います。

  2. スクリプトをコピペした上でSPMで結果を表示した後に、get_spm_names.m を実行したところ

    変数 “get_spm_names” またはクラス “get_spm_names.m” は未定義です。

    と表示されました。どのように対処すればよろしいでしょうか。

    • スクリプトを get_spm_names.m として、SPMのフォルダの中に保存していただけませんか?
      それで解決するのではないかと思います。

      • SPMとはDocuments\MATLAB\spm12のspm12の中に保存すればよいのですか?
        確実にspm12の中に保存しているのですが、同じメッセージが出てきてしまいます。

        get_spm_names.m を実行するとはコマンドウィンドウにget_spm_names.mと打てばいいんですか?

        • Matlabからは

          get_spm_names

          とタイプしてください。
          コマンドとしてうつ時、.m は不要となります。

  3. 根本先生
    お忙しいところありがとうございます。
    解析を見直してみます。
    今後ともどうかよろしくお願いいたします。

  4. 根本先生
    ご指導ありがとうございます
    rest-sham の画像を計算したものを one sample t-test にかけた結果となります。
    薄字の部分についての解釈というのは、どうなりますでしょうか。
    論文に結果を書くときには、どう書くのがいいのでしょうか。
    最初の黒字と薄字では左右が異なってることもあります。
    お忙しいところ誠に申し訳ありませんが、ご指導のほどよろしくお願いいたします。

    • rest-shamの画像の計算は正しいでしょうか?
      ここまで違う理由はどこにあるのでしょうか。

      黒字と薄字で左右が異なるのは、クラスターがものすごく大きいからだと思います。
      正しい結果を出せているのかを批判的に吟味することからはじめられた方がいいかなと感じました。

    • 拝見しました。私の思っているものと同じでした。

      確率の%であるということはクラスターサイズの量的な数値ではないとの解釈でよろしいでしょうか。

      はい、そのとおりです。今の場合でしたら、-42 28 -8 という座標はおそらく左下前頭回眼窩面なのですが、
      36.2%の可能性で白質
      27.7%の可能性で左下前頭回眼窩面
      9.6%の可能性で左前頭葉弁蓋部
      といった感じです。

      クラスターサイズとは関係ありません。

      ただ、気になったこととして、タイトルに open rest > open sham となっていますが、
      Design matrixとしては、one sample t-test になっています。

      one sample t-test の帰無仮説は、「信号値は0である」ですから、
      こういう結果になるのは当然です。

      rest と sham の時の画像を 2 sample t-test を行うのが普通だと思うのですが、

      rest-sham の画像を計算したものを one sample t-test にかけたのですか?

  5. 根本先生、お忙しいところありがとうございます。やってみます。ご指導よろしくお願いいたします。

  6. 根本先生 お忙しいところありがとうございます。FWEで結果を出してみます。確率の%であるということはクラスターサイズの量的な数値ではないとの解釈でよろしいでしょうか。

    • すみません、私の思っている%と違うようです。そこの場所のスクリーンショットを見れませんか?このコメントに直接画像は貼れませんが、たとえばDropboxの共有リンクなどを共有していただければ画像を見ることができます。

  7. 根本先生
    お忙しいところ返信誠に、ありがとうございました
    peak-levelのp値は p<0.05 FDR に設定しております
    例えば、cluster-levelでKeが73987というような数字になって、p FDRが0.012 最初の行が黒字で、そのあとに薄い字の行が3行ほどあって、それぞれの領域が異なってでますが、どの領域を選ぶというのは、どうやって決めたらよいのでしょうか。
    いくつかの領域が、%で示されますが、この%には、意味がありますでしょうか。
    基本的、常識的なことなのでしょうが、ご指導のほどよろしくお願いいたします

    • その場合は、p<0.05 FWE にされてはいかがでしょうか。
      %で示されるのは、その座標が、ある解剖学的名称である確率は何%だということを意味します。

  8. 根本先生
    いつもご指導誠にありがとうございます
    日本大学松戸歯学部 有床義歯補綴学 石井智浩 です

    初歩的なことで、申し訳ありません。
    SPMで脳の活動領域の結果を出す際に、活動領域が大きくなり、非常に大きなクラスターサイズの領域ができてしまうのですが、そのような時は、どのような対応をするのが、正しいのでしょうか。
    お忙しいところ誠に申し訳ありませんが、ご教授いただけると幸いです。
    どうかよろしくお願いいたします

    • SPMの結果は、統計の閾値によって決まります。peak-levelのp値をどのように設定されていらっしゃいますか?

  9. 根本先生
    ありがとうございます、small volume correctionの条件も知ることができ、よく理解できないまま進めていた解析が前に進みそうです!形になりましたらぜひご報告させてください。
    自身の勉強不足で試行錯誤しながら行き詰っていたところ、お力を貸していただきありがとうございます。
    神藤

  10. 根本先生
    本当にありがとうございます。small volume correctionの用い方まで知ることができ、データ解析が前に進みそうです!形になったらご報告します。
    自身の勉強不足で行き詰っていたところをお力を貸していただき、ありがとうございました。
    神藤

    • 神藤先生

      無事に進むように願っています。またわからなかったら質問してください。

      根本清貴

  11. 初めまして、東京慈恵会医大大学院の神藤と申します。座標の解剖学的名称についてご質問させてください。ブロックデザインで2群間で出た結果をSPMのアトラスでみると、白質領域を示してしまうことがあったり、SPMのラベルでは灰白質領域を示されるのですが、アトラス「Atlas of the Human Brain, Fourth Edition」で座標を確認すると白質、それもど真ん中だったりします。
    SPMの結果の評価をする際に、SPM上で結果を灰白質のみで抽出することは可能なのでしょうか
    SPMのメーリングリストJISCMailのアーカイブも検索しましたが、うまく見つけられず、こちらで質問させていただきました。よろしくお願いいたします。

    • 神藤先生

      おそらくタスクfMRIの結果をみていらっしゃると想定して返事します。

      まず、SPMでは、厳密に灰白質と白質を区別することはできないです。
      構造画像をsegmentationした際に、灰白質画像は、modulationする前は、「灰白質がある確率」を画像にしたものとなります。
      従って、あるボクセルの信号値が0.8だとすると、80%灰白質なわけですが、20%は灰白質でない(白質や脳脊髄液など)の可能性があるわけですね。

      そのlimitationを前提にしたうえでの方策ですが、

      ひとつの方法としては、Resultで結果を表示する際に、灰白質マスク画像を準備しておいて、そのマスク画像を指定することがあると思います。
      Masking toolbox でマスク画像を作ることができますので、そちらを試されたらいかがでしょうか?

      根本清貴

      • 根本先生

        早速のお返事ありがとうございます。はい、2グループ間でコントラストをとったタスクfMRI結果です。
        灰白質の確率が画像化されていることは理解しておりませんでした。勉強になります。
        Resultの結果をマスクをかけて表示を試みてみました。
        SPMの中に入っている標準化?された灰白質ファイル(SPMフォルダの中のcanonicalフォルダ内のファイル)を用いた場合と
        個々人のデータ全体の灰白質データを平均化して作成したマスクで
        表示結果が若干変わってしまいました。何か間違っているのか、それでいいのか、不安があります。
        もう少し頑張ってみますが、アドバイスいただけますと嬉しいです。

        神藤

        • 神藤先生

          canonicalフォルダ内のファイルは灰白質画像としてはいまいちかと思います。

          以下で灰白質マスク画像を作成してみました。

          spm12/tpm/TPM.nii を使用

          ImCalc で、TPM.nii,1 の 値が0.2以上のものをGM_mask.nii として保存
          Expressionは i1>0.2

          作ったものを以下に置いてみました。

          https://www.dropbox.com/s/aorfrhd9zkysj6n/GM_mask.nii?dl=0

          こちらで試してみてもらえたらと思います。

          根本清貴

          • ありがとうございます、早速ダウンロードさせていただいたものと、自分でもう一度マスクを作ってみました。
            結果表示の時にやはりアトラスでラベルをみると白質が出て来てしまいました。

            https://www.dropbox.com/s/rtuuk1c90zbx4xu/%E3%82%B9%E3%82%AF%E3%83%AA%E3%83%BC%E3%83%B3%E3%82%B7%E3%83%A7%E3%83%83%E3%83%88%202020-03-29%2015.33.38.png?dl=0

            マスクの適応の仕方が間違っているのかも、という不安が襲って来ました。

            Result で「small volume」>「Search volume..」の「image」で作ったマスクファイルを選択する方法で合っていたでしょうか。

            勉強不足で申し訳ございません。

            神藤

          • ありがとうございます、早速ダウンロードさせていただいたものと、自分でもう一度マスクを作ってみました。
            結果表示の時にやはりアトラスでラベルをみると白質が出て来てしまいました。

            https://www.dropbox.com/s/rtuuk1c90zbx4xu/%E3%82%B9%E3%82%AF%E3%83%AA%E3%83%BC%E3%83%B3%E3%82%B7%E3%83%A7%E3%83%83%E3%83%88%202020-03-29%2015.33.38.png?dl=0

            マスクの適応の仕方が間違っているのかも、という不安が襲って来ました。

            Result で「small volume」>「Search volume..」の「image」で作ったマスクファイルを選択する方法で合っていたでしょうか。

            勉強不足で申し訳ございません。

            神藤

          • 神藤先生

            Resultで、コントラストを選んだ後に、
            apply masking
            が出てくると思います。
            そこで、image
            を選択して、
            inclusive
            として進めていただくのがいいかと思います。

            small volume correctionは、その名の通り、もっと限局した関心領域を指定して、
            多重比較補正をかけ直すものです。仮説があることが前提になります。
            search volumeが減るので、多重比較補正に強くなります。そのために使うものです。

            根本清貴

コメントを残す

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください