QiitaにfMRIデータの解析事例が...ない...!!
先日、「qiitaっていろんな記事あるけど...pythonとか機械学習使ってfMRIデータを解析した例ってあまりないんじゃないか？」と思いました。そこで、実際にqiitaにfMRIのfMRIデータを解析した事例がないか検索してみました。
あれ？まじでない...？！
ツールや神経科学分野の解説記事はあるんですが、実物のfMRIデータで解析しているのは見当たりませんでした。
そこで今回は、fMRIデータを使って機械学習の解析をします。脳画像といっても、所詮は三次元の画像です。脳を真っ二つにスライスした二次元画像の断面図が何枚も重なって形成された形になっていて、一つ一つの画素が「ボクセル」と呼ばれる単位で情報を持っているだけです。4Dの脳画像になれば、三次元の各ボクセル情報に(fMRIでスキャンしていた時の)時系列情報が含まれているだけの話です。
脳画像はNIfTY形式と呼ばれる特殊なフォーマットで公開されていることが多く、これを読み込む必要があります。今回は学生や研究者がよく使っていると言われているMatlabではなく、pythonのNilearnというライブラリを使って脳画像を読み込んで解析を行います(誰でも解析できることを念頭に置くために、有料であるMatlabはやめました)。機械学習の解析自体は、機械学習をしたことがある人は皆使ったことや聞いたことがあるであろうscikit-learnのSVMで学習して予測を行います。
また、機械学習を使ったアプローチですが、今回はある特定の脳領域をマスキングしたfMRIデータを用いて、「脳活動からその時何をfMRIでみていたのか」を予測する解析を行います。脳解析の業界では脳情報デコーディングという名称で呼ぶみたいですが、ただの機械学習です(業界の人たちごめんなさい)。従来は何かしらの実験タスクで行われた行動を説明変数として、fMRIで行った時に統計的有意に反応した脳領域を独立変数として検出するなどといったエンコーディングのようなアプローチを取っていたのですが、それとは真逆の方策を取っているのでデコーディングと呼称されたのが由来だったと思います(つまり、タスクによって統計的に有意に励起された脳領域を見るのではなく、ある脳領域の活動からタスクを予測する)。
実際に使用したデータについて
今回は、fMRIデータの読み込みなどに利用する Nilearn:Machine learning for Neuro-Imaging in Pythonのデコーディング解析のチュートリアルを参考にして解析していますが、そこで使用されているThe Haxby 2001 experimentで公開されている脳画像を使用します。
A introduction tutorial to fMRI decoding
上記のチュートリアル通りに進めると、以下のコードで簡単にThe Haxby 2001 experimentの被験者のデータセットを集めることができると書かれていました。
from nilearn import datasets
# デフォルトでは、二番目の被験者のデータが取得される
haxby_dataset = datasets.fetch_haxby()
しかし、データを提供しているNIRTCのデータベース上の問題らしいのですがダウンロードすることができなかったので、直接ダウンロードしてローカルに保存する手段をとりました。こちらのURLから二番目の被験者のデータが入っているsubj2-2010.01.14.tar.gz をダウンロードできるので、データを解析したい方はこちらからダウンロードしてみてください(別に他の被験者番号のデータでも構いません)。
この実験では、被験者がfMRIスキャナの中で様々なカテゴリ(ハサミや人の顔、猫など)の視覚刺激が提示されます。このチュートリアルでは、腹側皮質視覚路(色や形の表象に携わっている)の領域内で記録されたfMRI脳活動から被験者がどのカテゴリをfMRIでみているのかを予測します。
解析手順
※ ソースコード上では、ダウンロードしたフォルダの名前をsubj2に変更しています。
画像を読み込んで可視化する
まずは、4DのfMRIデータを読み込んでみましょう。こちらは前述しましたが、ボクセルデータ(三次元データ)である脳画像に実験タスクを行なっていた時の時系列の情報が入っている形式になっています。
しかし、4DのfMRIデータなので、nilearn.plotting.plot_epiなどの3DのfMRIデータ(時系列情報なし)に対応した処理では可視化できません。そこで、可視化した4DのfMRIデータを平均化して一つの3DfMRIデータに変換させます。これは、from nilearn.image import mean_imgをインポートして実現できます。
fmri_filename = "~/Desktop/nilearn/subj2/bold.nii.gz"
from nilearn import plotting
from nilearn.image import mean_img
plotting.view_img(mean_img(fmri_filename), threshold=None)
結果
まずは、脳画像を可視化するところまでできました。たった数行で可視化できるのは素晴らしいことですね。
fMRIデータの特徴量抽出を行い、numpy形式のデータに変換する
numpyという単語が出てきてほっとした方もいるのではないでしょうか。ここでは、nilearn.input_data.NiftiMaskerという機能を使って腹側皮質視覚路のfMRIデータをマスキングして、numpyの形式に変換します。
まずは、腹側皮質視覚路だけのfMRIデータを可視化してみます。
mask_filename = "~/Desktop/nilearn/subj2/mask4_vt.nii.gz"
anat_filename = "~/Desktop/nilearn/subj2/anat.nii.gz"
# 被験者の解剖画像を背景にしてマスクの範囲を表示させてみる
plotting.plot_roi(mask_filename, bg_img=anat_filename,
cmap='Paired')
結果
次に、マスクを作成してみましょう。機械学習の精度を向上させるための方法としてデータを標準化させたいと思います。マスクは、nilearnによる機械学習の準備ができたら2D配列の形で抽出を行います。
from nilearn.input_data import NiftiMasker
masker = NiftiMasker(mask_img=mask_filename, standardize=True) # マスクを標準化する
fmri_masked = masker.fit_transform(fmri_filename) # 標準化させたマスクに対してBOLD信号を当てはめる
ちなみに、masker.generate_report()を実行することによって、マスクのNIfTYイメージを示して様々なパラメータなどを確認することができます。
これで、fmri_maskedという形でマスキングされたデータ行列をnumpy形式として抽出することができました。実際にみてみます。
# fmri_masked はnumpy配列の形式になっています。
print(fmri_masked)
# このデータのサイズは、ボクセル数と時系列点(スキャン数)の総数で構成されています。
# shape[0]はスキャン数, shape[1]はボクセル数
print(fmri_masked.shape)
出力結果
[[ 7.6757896e-01 2.3108697e+00 -2.0519443e-01 ... -1.0261143e+00
8.7993503e-02 2.0720518e+00]
[ 5.5640817e-01 1.6833434e+00 -2.4644937e-01 ... -7.0238107e-01
-3.4570047e-01 2.0341001e+00]
[ 7.6757896e-01 1.9186659e+00 1.0802225e-03 ... -9.9374104e-01
-2.7630943e-01 2.1479552e+00]
...
[-4.2905563e-01 -1.6896105e+00 -7.4150854e-01 ... -1.5440876e+00
1.8054217e+00 -1.6709718e-01]
[-1.4749455e-01 -1.8072717e+00 -2.4644937e-01 ... -1.7707009e+00
1.5452052e+00 7.8169477e-01]
[-2.1788482e-01 -1.4542881e+00 1.0802225e-03 ... -1.6412076e+00
1.2676411e+00 8.9554977e-01]]
(1452, 464)
fmri_masked.shapeで確認した通り、464にはマスキングされた脳領域のボクセル数、1452は時系列情報が格納されています(fMRIはTRという単位の時間分解能でfMRIデータが取得されており、通常は1~2.5秒くらいの間隔で設定されていることが多いですが、特定の脳領域の情報を使って解析するデコーディングやリアルタイムで脳活動を測定して被験者を訓練させるようなニューロフィードバックと呼ばれる手法では、かなり短めのTRが設定されていることが多かった気がします)。
試しに、最初の三つのボクセルの時系列情報を取り出すと以下のようになります。
# 行列から選ばれた最初の三つのボクセルを表示させる。
import matplotlib.pyplot as plt
plt.plot(fmri_masked[5:150, :3])
plt.title('Voxel Time Series')
plt.xlabel('Scan number')
plt.ylabel('Normalized signal')
plt.tight_layout()
ちなみに鋭い方は、fmri_masked[5:150, :3]の部分をどうしてfmri_masked[:150, :3]の形にしてないのか疑問に思うかもしれません。これはダミースキャンといって、fMRIを撮像し始めた最初のスキャンは高い信号値を持っているため実際の解析では使用しないことが多いからです。信号値が高くなるのは、T1効果が安定した状態にならないためです(詳しくは他のサイトでググってください笑)。
実験タスク時の行動ラベルを取り込む
腹側皮質視覚路のfMRIデータをnumpy形式として取り出すことができたので、このデータを使って機械学習を行なっていきます。今回は教師あり学習を行うので、データ以外にラベルが必要です。そこで、ダウンロードされたデータに入っている実験時で被験者がみていたカテゴリの情報をラベルとして取り出す作業をします。
import pandas as pd
# 行動結果の情報を読み込む。スキャン数だけラベルがある
behavioral = pd.read_csv('~/Desktop/nilearn/subj2/labels.txt', delimiter=' ')
出力結果
labels chunks
0 rest 0
1 rest 0
2 rest 0
3 rest 0
4 rest 0
5 rest 0
6 scissors 0
7 scissors 0
8 scissors 0
9 scissors 0
10 scissors 0
11 scissors 0
12 scissors 0
13 scissors 0
14 scissors 0
... ... ...
1426 cat 11
1427 cat 11
1428 cat 11
1429 cat 11
1430 cat 11
1431 cat 11
1432 rest 11
1433 rest 11
1434 rest 11
1435 rest 11
1436 rest 11
1437 scissors 11
1438 scissors 11
1439 scissors 11
1440 scissors 11
1441 scissors 11
1442 scissors 11
1443 scissors 11
1444 scissors 11
1445 scissors 11
1446 rest 11
1447 rest 11
1448 rest 11
1449 rest 11
1450 rest 11
1451 rest 11
1452 rows × 2 columns
このタスクは先述したように視覚認識タスクなので、ラベルは実験条件(被験者に見せた対象物のカテゴリ)が記されています。chunksはセッション数を表しています(被験者にずっと永遠にタスクをやらせると身体への負荷が半端ないので、通常はセッション区切りに体調を確認したりしながら複数回に渡って類似したタスクを行っていきます)。
conditions = behavioral['labels']
出力結果
0 rest
1 rest
2 rest
3 rest
4 rest
5 rest
6 scissors
7 scissors
8 scissors
9 scissors
10 scissors
11 scissors
12 scissors
13 scissors
14 scissors
15 rest
16 rest
17 rest
18 rest
19 rest
20 rest
21 face
22 face
...
1431 cat
1432 rest
1433 rest
1434 rest
1435 rest
1436 rest
1437 scissors
1438 scissors
1439 scissors
1440 scissors
1441 scissors
1442 scissors
1443 scissors
1444 scissors
1445 scissors
1446 rest
1447 rest
1448 rest
1449 rest
1450 rest
1451 rest
Name: labels, Length: 1452, dtype: object
今回はチュートリアルに倣って、猫と顔の実験条件だけを使用します。他にもたくさんのラベルがあるので、興味がある方は他のカテゴリで検証したりしてみてください。
# マスキングした脳画像に対して顔と猫を提示した実験条件時のスキャン数に対応するボクセルを取ってくる
condition_mask = conditions.isin(['face','cat'])
fmri_masked = fmri_masked[condition_mask]
print(fmri_masked.shape)
# ラベルに対しても同じ条件を適応させる
conditions = conditions[condition_mask]
print(conditions.shape)
SVMでデコーディングを行う
今回は、取り出したfmri_masked dataに対してscikit-learnの機械学習ツールボックスを使用する。 デコーダとして、SVMの線形カーネルを使用しています。scikit-learnなので、数行で学習から予測までできます
from sklearn.svm import SVC
svc = SVC(kernel='linear')
# 学習して予測を行う
svc.fit(fmri_masked, conditions)
prediction = svc.predict(fmri_masked)
しかし、これだと解析としては全然ダメなので、交差検証をします(当たり前ですが)
交差検証を行なって予測スコアを算出する
KFold交差検証を行う
from sklearn.model_selection import KFold
import numpy as np
cv = KFold(n_splits=5)
scores = []
for train, test in cv.split(X=fmri_masked):
svc.fit(fmri_masked[train], conditions.values[train])
prediction = svc.predict(fmri_masked[test])
scores.append((prediction == conditions.values[test]).sum()
/ float(len(conditions.values[test])))
# 予測性能の平均
print(np.mean(scores))
出力結果:
0.7628964059196617
また、cross_val_score関数を使って、サブセットごとの評価の計算処理をマシン上に乗っているCPUを複数用いて分散処理を行うことも可能です。(n_jobs=nの形で指定でき、-1と設定すればマシン上で利用可能なすべてのCPUを用いて並列処理を行うことができる)。
from sklearn.model_selection import cross_val_score
cv_score = cross_val_score(svc, fmri_masked, conditions)
print(cv_score)
出力結果：
[0.86363636 0.76744186 0.74418605 0.69767442 0.81395349]
しかし、これでもまだ検証としては不十分です。理由は、セッションごとの影響を考慮した性能評価を行なってないからです。
セッションごとに予測精度の評価を行う
fMRIデータはセッション単位で取得されており、ノイズは特定のセッションごとに自己相関されている形になっています。したがって、これらの影響を考慮した上で交差検証を行う際には、セッションをまたがって予測する必要性があります。
from sklearn.model_selection import LeaveOneGroupOut
cv = LeaveOneGroupOut()
cv_score = cross_val_score(svc,
fmri_masked,
conditions,
cv=cv,
groups=session_label,
)
# セッションごとの検証データでの予測精度を取得する
print(cv_score)
出力結果：
[0.55555556 1. 0.66666667 0.66666667 0.77777778 0.72222222
0.88888889 0.38888889 0.66666667 0.5 0.77777778 0.66666667]
学習したモデルの重みを推定する
最後に、モデルの重みを推定して表示させます。そのために、まずは重みをNIfTY画像に変換する処理を行います。
coef_ = svc.coef_
# numpy arrayの形式になっており、ボクセルごとに一つの係数(重み)を持っている
print(coef_.shape)
出力結果：
(1, 464)
coef_img = masker.inverse_transform(coef_)
print(coef_img)
出力結果：
<class 'nibabel.nifti1.Nifti1Image'>
data shape (40, 64, 64, 1)
affine:
[[ -3.5 0. 0. 68.25 ]
[ 0. 3.75 0. -118.125]
[ 0. 0. 3.75 -118.125]
[ 0. 0. 0. 1. ]]
metadata:
<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'
sizeof_hdr : 348
data_type : b''
db_name : b''
extents : 0
session_error : 0
regular : b''
dim_info : 0
dim : [ 4 40 64 64 1 1 1 1]
intent_p1 : 0.0
intent_p2 : 0.0
intent_p3 : 0.0
intent_code : none
datatype : float64
bitpix : 64
slice_start : 0
pixdim : [-1. 3.5 3.75 3.75 1. 1. 1. 1. ]
vox_offset : 0.0
scl_slope : nan
scl_inter : nan
slice_end : 0
slice_code : unknown
xyzt_units : 0
cal_max : 0.0
cal_min : 0.0
slice_duration : 0.0
toffset : 0.0
glmax : 0
glmin : 0
descrip : b''
aux_file : b''
qform_code : unknown
sform_code : aligned
quatern_b : 0.0
quatern_c : 1.0
quatern_d : 0.0
qoffset_x : 68.25
qoffset_y : -118.125
qoffset_z : -118.125
srow_x : [-3.5 0. 0. 68.25]
srow_y : [ 0. 3.75 0. -118.125]
srow_z : [ 0. 0. 3.75 -118.125]
intent_name : b''
magic : b'n+1'
作成したNIfTY形式の画像をnii.gz形式で保存します
# nii.gzファイル形式で保存する
coef_img.to_filename('haxby_svc_weights.nii.gz')
最後に、重み係数を脳画像として可視化してみます。
# 実際にSVMの重みをプロットしてみる
from nilearn.plotting import plot_stat_map, show
plot_stat_map(coef_img, bg_img=anat_filename,
title="SVM weights", display_mode='yx')
show()
以上で機械学習を用いた一通りのfMRIデータ解析を行いました。
まとめ
今回の解析は個人ブログで投稿する予定だったのですが、QiitaでfMRIデータを解析していた人が見当たらなかったのでこちらで投稿しました(個別ブログはプロフィールに載せています)。
fMRIデータの画像...聞くだけで難しそうと思うかもしれませんが、結局はノイズが一般的な画像より多いことや次元が一つ多いこと、特殊なフォーマットを扱う部分を取り除けばただの画像解析です。どんな人でもマニアックなデータを解析できることを一つの目標にして個人ブログやQiitaを通じて発信しているので、この活動を通してあまり触れることのないデータの解析に少しでも興味を持ってくれたら幸いです。
よかったらLG...なんとかをよろしくお願いします！笑
追伸という名の注意書き
今回は生のfMRI時系列データを使って解析を行いますが、本来はHRF関数と呼ばれるものを畳み込みして課題提示(神経活動)から血流動態(BOLD反応)へ変換する必要があります。さらに言えば、1試行ごと(trial-by-trial)の活動マップ(beta map)を取得するために、一般線形モデルを用いて1st-levelでフィッティングさせる必要があります。このようにすることで、トライアルごとに関連した条件間で予測を行うことが可能になります。本当は様々な処理過程が実際のfMRIデータ解析で必要になってくるのですが、今回はfMRIデータに機械学習を適用させることを目的とした入門なので、先述した処理は行いことにしました(実際にこれだけの手順を踏むとなると、一つに記事に収まりきらないです汗)。
ご了承ください(機会があれば本格的に何か解析したいと思います)