pythonでfMRIデータ解析をやってみた(脳情報デコーディング入門)

## QiitaにfMRIデータの解析事例が...ない...!!

先日、「qiitaっていろんな記事あるけど...pythonとか機械学習使って**fMRIデータ**を解析した例ってあまりないんじゃないか？」と思いました。そこで、実際にqiitaにfMRIのfMRIデータを解析した事例がないか検索してみました。

![スクリーンショット 2020-04-23 15.18.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580487/685dd352-e246-c620-3570-6af0845b3ecc.png)

あれ？まじでない...？！
ツールや神経科学分野の解説記事はあるんですが、実物の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](https://nilearn.github.io/auto_examples/plot_decoding_tutorial.html#sphx-glr-auto-examples-plot-decoding-tutorial-py)


上記のチュートリアル通りに進めると、以下のコードで簡単にThe Haxby 2001 experimentの被験者のデータセットを集めることができると書かれていました。

```python
from nilearn import datasets
# デフォルトでは、二番目の被験者のデータが取得される
haxby_dataset = datasets.fetch_haxby()
```

しかし、データを提供しているNIRTCのデータベース上の問題らしいのですがダウンロードすることができなかったので、直接ダウンロードしてローカルに保存する手段をとりました。こちらのURLから二番目の被験者のデータが入っているsubj2-2010.01.14.tar.gz をダウンロードできるので、データを解析したい方はこちらからダウンロードしてみてください(別に他の被験者番号のデータでも構いません)。

http://data.pymvpa.org/datasets/haxby2001/

この実験では、被験者がfMRIスキャナの中で様々なカテゴリ(ハサミや人の顔、猫など)の視覚刺激が提示されます。このチュートリアルでは、腹側皮質視覚路(色や形の表象に携わっている)の領域内で記録されたfMRI脳活動から被験者がどのカテゴリをfMRIでみているのかを予測します。


## 解析手順
※ ソースコード上では、ダウンロードしたフォルダの名前をsubj2に変更しています。

### 画像を読み込んで可視化する

まずは、4DのfMRIデータを読み込んでみましょう。こちらは前述しましたが、ボクセルデータ(三次元データ)である脳画像に実験タスクを行なっていた時の時系列の情報が入っている形式になっています。

しかし、4DのfMRIデータなので、nilearn.plotting.plot_epiなどの3DのfMRIデータ(時系列情報なし)に対応した処理では可視化できません。そこで、可視化した4DのfMRIデータを平均化して一つの3DfMRIデータに変換させます。これは、`from nilearn.image import mean_img`をインポートして実現できます。

```python
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)
```

結果
![スクリーンショット 2020-04-23 22.10.45.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580487/ef053a6a-5a2c-03cf-c56a-f72f69ae58a7.png)

まずは、脳画像を可視化するところまでできました。たった数行で可視化できるのは素晴らしいことですね。

### fMRIデータの特徴量抽出を行い、numpy形式のデータに変換する

numpyという単語が出てきてほっとした方もいるのではないでしょうか。ここでは、`nilearn.input_data.NiftiMasker`という機能を使って腹側皮質視覚路のfMRIデータをマスキングして、numpyの形式に変換します。

まずは、腹側皮質視覚路だけのfMRIデータを可視化してみます。

```python
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')
```

結果
![スクリーンショット 2020-04-23 22.25.35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580487/0b416b53-6504-0e68-f3cb-dce3258f4435.png)

次に、マスクを作成してみましょう。機械学習の精度を向上させるための方法としてデータを標準化させたいと思います。マスクは、nilearnによる機械学習の準備ができたら2D配列の形で抽出を行います。

```python
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形式として抽出することができました。実際にみてみます。

```python
# fmri_masked はnumpy配列の形式になっています。
print(fmri_masked)
# このデータのサイズは、ボクセル数と時系列点(スキャン数)の総数で構成されています。
# shape[0]はスキャン数, shape[1]はボクセル数
print(fmri_masked.shape)
```

出力結果

```python
[[ 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が設定されていることが多かった気がします)。

試しに、最初の三つのボクセルの時系列情報を取り出すと以下のようになります。

```python
# 行列から選ばれた最初の三つのボクセルを表示させる。
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()
```

![スクリーンショット 2020-04-23 22.44.52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580487/3984bbd3-7f9e-e120-2ede-bc46a10579cf.png)

ちなみに鋭い方は、`fmri_masked[5:150, :3]`の部分をどうして`fmri_masked[:150, :3]`の形にしてないのか疑問に思うかもしれません。これは**ダミースキャン**といって、fMRIを撮像し始めた最初のスキャンは高い信号値を持っているため実際の解析では使用しないことが多いからです。信号値が高くなるのは、T1効果が安定した状態にならないためです(詳しくは他のサイトでググってください笑)。

### 実験タスク時の行動ラベルを取り込む
腹側皮質視覚路のfMRIデータをnumpy形式として取り出すことができたので、このデータを使って機械学習を行なっていきます。**今回は教師あり学習を行うので、データ以外にラベルが必要です**。そこで、ダウンロードされたデータに入っている実験時で被験者がみていたカテゴリの情報をラベルとして取り出す作業をします。

```python
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はセッション数を表しています(被験者にずっと永遠にタスクをやらせると身体への負荷が半端ないので、通常はセッション区切りに体調を確認したりしながら複数回に渡って類似したタスクを行っていきます)。

```python
conditions = behavioral['labels']
```

出力結果

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

今回はチュートリアルに倣って、猫と顔の実験条件だけを使用します。他にもたくさんのラベルがあるので、興味がある方は他のカテゴリで検証したりしてみてください。

```python
# マスキングした脳画像に対して顔と猫を提示した実験条件時のスキャン数に対応するボクセルを取ってくる
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なので、数行で学習から予測までできます

```python
from sklearn.svm import SVC
svc = SVC(kernel='linear')

# 学習して予測を行う
svc.fit(fmri_masked, conditions)
prediction = svc.predict(fmri_masked)
```

しかし、これだと解析としては全然ダメなので、交差検証をします(当たり前ですが)

### 交差検証を行なって予測スコアを算出する

#### KFold交差検証を行う

```python
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を用いて並列処理を行うことができる)。

```python
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データはセッション単位で取得されており、ノイズは特定のセッションごとに自己相関されている形になっています。したがって、これらの影響を考慮した上で交差検証を行う際には、セッションをまたがって予測する必要性があります。

```python
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画像に変換する処理を行います。

```python
coef_ = svc.coef_

# numpy arrayの形式になっており、ボクセルごとに一つの係数(重み)を持っている
print(coef_.shape)
```

```
出力結果：
(1, 464)
```

```python
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形式で保存します

```python
# nii.gzファイル形式で保存する
coef_img.to_filename('haxby_svc_weights.nii.gz')
```

最後に、重み係数を脳画像として可視化してみます。

```python
# 実際に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()
```

![スクリーンショット 2020-04-24 0.12.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580487/97130651-7524-74f8-926b-05ac47939bce.png)

以上で機械学習を用いた一通りのfMRIデータ解析を行いました。

### まとめ

今回の解析は個人ブログで投稿する予定だったのですが、QiitaでfMRIデータを解析していた人が見当たらなかったのでこちらで投稿しました(個別ブログはプロフィールに載せています)。

fMRIデータの画像...聞くだけで難しそうと思うかもしれませんが、結局はノイズが一般的な画像より多いことや次元が一つ多いこと、特殊なフォーマットを扱う部分を取り除けばただの画像解析です。どんな人でもマニアックなデータを解析できることを一つの目標にして個人ブログやQiitaを通じて発信しているので、この活動を通してあまり触れることのないデータの解析に少しでも興味を持ってくれたら幸いです。

よかったらLG...なんとかをよろしくお願いします！笑


### 追伸という名の注意書き

今回は生のfMRI時系列データを使って解析を行いますが、本来はHRF関数と呼ばれるものを畳み込みして課題提示(神経活動)から血流動態(BOLD反応)へ変換する必要があります。さらに言えば、1試行ごと(trial-by-trial)の活動マップ(beta map)を取得するために、一般線形モデルを用いて1st-levelでフィッティングさせる必要があります。このようにすることで、トライアルごとに関連した条件間で予測を行うことが可能になります。本当は様々な処理過程が実際のfMRIデータ解析で必要になってくるのですが、今回はfMRIデータに機械学習を適用させることを目的とした入門なので、先述した処理は行いことにしました(実際にこれだけの手順を踏むとなると、一つに記事に収まりきらないです汗)。

ご了承ください(機会があれば本格的に何か解析したいと思います)