Pythonで脳波解析：Python MNEのチュートリアル


[Python MNE](https://mne.tools/stable/index.html)はオープンソースの脳磁図（MEG），脳波(EEG)の解析や可視化のツールです．多くのデバイスのデータフォーマットに適用できるため，汎用性が高いと言えるでしょう．この記事では，[最もベーシックなチュートリアル](https://mne.tools/stable/auto_tutorials/intro/plot_10_overview.html#sphx-glr-auto-tutorials-intro-plot-10-overview-py)に沿って，MEGとEEGのMNEによる解析手順を説明します．

# 実行環境
* Mac OS 10.15.3
* Python 3.6.5
* mne 0.20.4

# インストール
[Anaconda](https://anaconda.org/conda-forge/mne)でも[pip](https://pypi.org/project/mne/)でもインストールできます．[公式ページ](https://mne.tools/stable/install/mne_python.html)ではAnacondaの方がおすすめだそうです．

* Anacondaの場合

```bash

conda install -c conda-forge mne
```

* pipの場合

```bash

pip install -U mne
```

また，[GPUの設定](https://mne.tools/stable/install/advanced.html#using-mne-python-with-cuda-nvidia-gpu-acceleration)もできるようなので，フィルタリングやサンプリングで時間がかかるようであれば設定すると良いかもしれません．

```bash

$ conda install cupy
$ MNE_USE_CUDA=true python -c "import mne; mne.cuda.init_cuda(verbose=True)"
Enabling CUDA with 1.55 GB available memory

>>> mne.utils.set_config('MNE_USE_CUDA', 'true')  
# mne.io.Raw.filter()などにn_jobs='cuda'を渡す
```

# MEG/EEG解析のチュートリアル
[こちらのチュートリアル](https://mne.tools/stable/auto_tutorials/intro/plot_10_overview.html#sphx-glr-auto-tutorials-intro-plot-10-overview-py)に沿って，脳波解析の典型的な処理とMNEの基本的なデータ構造について説明します．適当に簡潔にかいつまんで説明するので，英語と脳波解析の知識がある方は元のページを読んでもらったほうがよろしいかと思います．

## 解析手順
* [データのロード](#データのロード)
* [前処理](#前処理)
* [イベント検出](#イベント検出)
* [エポックに切り出す](#エポックに切り出す)
* [時間周波数解析](#時間周波数解析)
* [誘発反応の推定](#誘発反応の推定)

モジュールのインポート

```python
import os
import numpy as np
import mne
```

## データのロード
mneにサンプルとして用意されているデータをロードします（[全データセットの一覧](https://mne.tools/stable/overview/datasets_index.html#datasets)）．ちなみに，MEGやEEGにはデバイスごとにいくつかファイルフォーマットがあるのですが，それらのデータをロードする関数が用意されています（[サポートされているデータフォーマット](https://mne.tools/stable/overview/implementation.html#data-formats)）

```python
sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',
                                    'sample_audvis_filt-0-40_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file)
```
今回は`sample_audvis_filt-0-40_raw.fif`という40Hzでローパスフィルタ済みのデータを使用しています．データ（1.54GB）のダウンロードには1分ほど要しました．<img width="598" alt="downloading.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/a725134d-ad97-9494-0b9b-b1c569f34a0d.png">

`mne.io.read_raw_fif()`という関数により，fifフォーマットのファイルからRawオブジェクトと呼ばれるMNEのデータ構造のデータを読み込んでいます．

ちなみにこのデータの実験では，被験者が視覚刺激と聴覚刺激を左右ランダムに750msごとに呈示されます．また時々画面に笑った顔が表示され，被験者はそれを見たときにボタンを押すというものだそうです．これらの刺激の呈示のことをイベントと呼びます．

ロードしたオブジェクトについて表示します．

```python
print(raw)
print(raw.info)

# output
<Raw | sample_audvis_filt-0-40_raw.fif, 376 x 41700 (277.7 s), ~3.6 MB, data not loaded>
<Info | 15 non-empty values
 bads: 2 items (MEG 2443, EEG 053)
 ch_names: MEG 0113, MEG 0112, MEG 0111, MEG 0122, MEG 0123, MEG 0121, MEG ...
 chs: 204 GRAD, 102 MAG, 9 STIM, 60 EEG, 1 EOG
 custom_ref_applied: False
 dev_head_t: MEG device -> head transform
 dig: 146 items (3 Cardinal, 4 HPI, 61 EEG, 78 Extra)
 file_id: 4 items (dict)
 highpass: 0.1 Hz
 hpi_meas: 1 item (list)
 hpi_results: 1 item (list)
 lowpass: 40.0 Hz
 meas_date: 2002-12-03 19:01:10 UTC
 meas_id: 4 items (dict)
 nchan: 376
 projs: PCA-v1: off, PCA-v2: off, PCA-v3: off, Average EEG reference: off
 sfreq: 150.2 Hz
>
```
MNEには，[Raw](https://mne.tools/stable/generated/mne.io.Raw.html#mne.io.Raw), [Epochs](https://mne.tools/stable/generated/mne.Epochs.html#mne.Epochs), [Evoked](https://mne.tools/stable/generated/mne.Evoked.html#mne.Evoked)オブジェクトと呼ばれるオブジェクトによってデータが管理されています．データの時間的な区切られ方によってこれらは区別されており，Rawでは脳波の収録開始から停止までのすべての時間のデータを含んでおり，イベントにより区切られる前の状態だと考えてもらうと良いと思います．`raw.info`では`raw`のチャンネル数や名前，サンプリングレートなどの情報が保持されています．

このRawオブジェクトのpower spectral density(PSD)とセンサーごとの波形を表示します．`plot()`はインタラクティブシェルではスクロールやスケーリングが可能です．

```python
raw.plot_psd(fmax=50)
raw.plot(duration=5, n_channels=30)
```
<img width="467" alt="psd_and_plot.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/447d8d5e-bc48-0077-afc6-546802540f4a.png">

## 前処理
EEG/MEG解析のノイズ処理などでは定番の独立成分分析（ICA）を行います．簡単に言うと，多チャンネルの信号を独立した成分に分解することで，電源ノイズや瞬きなど脳活動に独立なものを取り除く処理を行います．

ここでは，364チャンネルの信号に対して，20個の独立成分に分解した後，`[1, 2]`と取り除く成分のインデックスを指定しています（実際の除去は，次の`ica.apply(raw)`によって行われます）．`ica_properties()`で，独立成分の情報を表示します．


```python 
# set up and fit the ICA
ica = mne.preprocessing.ICA(n_components=20, random_state=97, max_iter=800)
ica.fit(raw)
ica.exclude = [1, 2]  # details on how we picked these are omitted here
ica.plot_properties(raw, picks=ica.exclude)

# output
Fitting ICA to data using 364 channels (please be patient, this may take a while)
Inferring max_pca_components from picks
Selecting by number: 20 components
Fitting ICA took 2.4s.
    Using multitaper spectrum estimation with 7 DPSS windows
138 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
138 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
```

<img width="523" alt="ica_1.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/347e08bb-e241-7c28-804e-d382dcb55278.png">
<img width="521" alt="ica_2.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/156c89e5-2968-4429-bc63-8d5dac90e7da.png">

指定した独立成分を除去した信号を手に入れるため，先ほどのICAオブジェクトをRawオブジェクトに適用します．

```python
orig_raw = raw.copy()
raw.load_data()
ica.apply(raw)
```

前頭のチャンネルを表示してみると次の2つの図のようにノイズが除去されています．

```python
# show some frontal channels to clearly illustrate the artifact removal
chs = ['MEG 0111', 'MEG 0121', 'MEG 0131', 'MEG 0211', 'MEG 0221', 'MEG 0231',
       'MEG 0311', 'MEG 0321', 'MEG 0331', 'MEG 1511', 'MEG 1521', 'MEG 1531',
       'EEG 001', 'EEG 002', 'EEG 003', 'EEG 004', 'EEG 005', 'EEG 006',
       'EEG 007', 'EEG 008']
chan_idxs = [raw.ch_names.index(ch) for ch in chs]
orig_raw.plot(order=chan_idxs, start=12, duration=4)
raw.plot(order=chan_idxs, start=12, duration=4)

# output
Reading 0 ... 41699  =      0.000 ...   277.709 secs...
Transforming to ICA space (20 components)
Zeroing out 2 ICA components
```
<img width="440" alt="orig_raw.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/da5a8614-4606-bc1a-f980-dccf637f1b0f.png">
<img width="445" alt="prune_raw.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/4356e585-5638-7908-c69c-b22a063b5664.png">

## イベント検出
このサンプルデータにはSTIMチャンネルと呼ばれるものが含まれており，そのチャンネルにイベント種類とその時刻のデータがあります．今回は`STI 014`という名前のチャンネルが該当するため，それを指定して実験のイベントを取り出します．

```python
events = mne.find_events(raw, stim_channel='STI 014')
print(events[:5])  # show the first 5

# output
319 events found
Event IDs: [ 1  2  3  4  5 32]

[[6994    0    2]
 [7086    0    3]
 [7192    0    1]
 [7304    0    4]
 [7413    0    2]]
```

`events`はnumpy arrayであり，最初の列がサンプリングの時刻，最後の列がイベントの種類（ID）を表しています（真ん中の0は無視して良いです）．ここでは，IDが`[ 1  2  3  4  5 32]`のイベントがあるそうです．

| Event ID | Condition |
|:-:|:-:|
| 1  | auditory stimulus (tone) to the left ear  |
| 2  | auditory stimulus (tone) to the right ear  |
| 3  | visual stimulus (checkerboard) to the left visual field  |
| 4  | visual stimulus (checkerboard) to the right visual field  |
| 5  | smiley face (catch trial)  |
| 32  | subject button press  |

Event IDとそのイベントの説明を辞書型で格納します．ちなみにここでのkeyは任意の文字列で構いません．

```python
event_dict = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
              'visual/right': 4, 'smiley': 5, 'buttonpress': 32}
```

時間軸上にそれぞれのイベントの呈示タイミングをプロットしてみます．`event_id=evend_dict`によって先ほど定義したIDとイベントの名前がプロット上で対応づいています．

```python
fig = mne.viz.plot_events(events, event_id=event_dict, sfreq=raw.info['sfreq'],
                          first_samp=raw.first_samp)
```
<img width="457" alt="events.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/b59192dc-f773-244e-ad06-7d69f0a2975a.png">

## エポックに切り出す
MEG/EEG解析で，特にある刺激に対する反応について調べる際に，信号をエポックと呼ばれる単位に分けて，一つのデータとして扱います．刺激に対して前後の一定の区間を切り出すことで，それを一つのデータ単位とします．この切り出し処理のことをエポッキングと呼びます．Rawオブジェクトは収録中の全時間がつながっているものですが，これらを刺激に対する反応の部分のみ切り出したものが，Epochsオブジェクトというデータ型で保持されます．

エポッキングにおいて，しばしばスレッショルドリジェクションと呼ばれる処理が行われます．これは該当の時間中に大きすぎる信号が含まれている場合，その部分は大きなノイズが入り込んでしまったであろうということで，データから除外するというものです．その除外の基準を決めます．

```python
reject_criteria = dict(mag=4000e-15,     # 4000 fT
                       grad=4000e-13,    # 4000 fT/cm
                       eeg=150e-6,       # 150 µV
                       eog=250e-6)       # 250 µV
```

次に，RawオブジェクトからEpochsオブジェクトを作成します．ここでの`mne.Epochs`の引数は以下の通りです．

* `raw`：信号のデータを含んだRawオブジェクト
* `events`：時間とイベントIDの情報を持ったnumpy array
* `event_id`：イベントIDとその名称を対応付ける辞書型
* `tmin`：刺激呈示の何秒からのデータを切り出すか（刺激呈示前の場合マイナスを付ける）
* `tmax`：刺激呈示後何秒までのデータを切り出すか
* `reject`：データ除外の基準

```python
epochs = mne.Epochs(raw, events, event_id=event_dict, tmin=-0.2, tmax=0.5,
                    reject=reject_criteria, preload=True)

# output
319 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
Created an SSP operator (subspace dimension = 4)
4 projection items activated
Loading data for 319 events and 106 original time points ...
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on MAG : ['MEG 1711']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on MAG : ['MEG 1711']
    Rejecting  epoch based on EEG : ['EEG 008']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
10 bad epochs dropped
```

今回は10個のエポックが除外されたようです．

被験者のボタンを押したタイミングは今回扱わないので，視聴覚刺激イベントのみを対象にします．また，条件間でのデータ数に偏りが生じることを防ぐため，`epochs.equalize_event_counts()`で，数が多い条件のイベントを消すことで，各条件の数を揃えています．

```python
conds_we_care_about = ['auditory/left', 'auditory/right',
                       'visual/left', 'visual/right']
epochs.equalize_event_counts(conds_we_care_about)  # this operates in-place
```

`event_dict`で`'auditory/left'`と`auditory/right`のように指定していましたが，Epochsオブジェクトに渡すことで，`auditory`を指定するとこれら両方を指すことになります．ここで聴覚刺激に対するデータと視覚刺激データに分けます．

```python
aud_epochs = epochs['auditory']
vis_epochs = epochs['visual']
del raw, epochs  # free up memory
```

聴覚データでチャンネルを指定して，刺激呈示時間を揃えた加算平均波形を表示します．

```python
aud_epochs.plot_image(picks=['MEG 1332', 'EEG 021'])

# output
136 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
136 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
```
<img width="460" alt="meg_aud.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/db08ef19-7fdd-a68c-4616-37eddf655ea9.png">
<img width="470" alt="meg_viz.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/54f19534-f2ae-454c-7ae8-6e7a59de9f19.png">

補足ですが，Epochsオブジェクトの`get_data()`メソッドによって，numpy arrayの形で時系列波形データを取得できるので，条件ごとにarrayを取ることで，機械学習のモデルに入力することなどが可能です．

```python
aud_epochs.get_data().shape

(136, 376, 106) # (エポック，チャンネル，時刻)
```
以下の周波数処理などもnumpy arrayを返すものも同様です．

## 時間周波数解析
MEG/EEG波形における，刺激呈示前後の時間と周波数の関係について調べるために，時間周波数解析を行います．ここでは，聴覚刺激に対する反応を，wavelet解析によって計算し，それをプロットします．カラーバーに表示されているように，赤が濃くなるほど，その時間のその周波数の振幅が大きいということを示しています．（[もっと詳しいExample](https://mne.tools/stable/auto_tutorials/time-freq/plot_sensors_time_frequency.html#inter-trial-coherence)）


```python
frequencies = np.arange(7, 30, 3)
power = mne.time_frequency.tfr_morlet(aud_epochs, n_cycles=2, return_itc=False,
                                      freqs=frequencies, decim=3)
power.plot(['MEG 1332'])

# output
No baseline correction applied
```
<img width="455" alt="tf.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/f8019161-8a32-0ab8-a6cb-0ff50f86f667.png">

## 誘発反応の推定
MNEのデータ型として，Epochsオブジェクトの加算平均を行ったものがEvokedオブジェクトになります．MEG/EEGデータは脳活動以外のノイズが多く，加算平均を行うことで条件間の差を調べることが多くあります．

ここでは，聴覚データと視覚データの加算波形の違いをプロットします．

```python
aud_evoked = aud_epochs.average()
vis_evoked = vis_epochs.average()

mne.viz.plot_compare_evokeds(dict(auditory=aud_evoked, visual=vis_evoked),
                             legend='upper left', show_sensors='upper right')
```

<img width="562" alt="meg_evoked.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/89aed283-ae1d-57b3-f36a-20d858f85648.png">

<img width="563" alt="eeg_evoked.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/d757ccbd-d425-b1d8-3a04-71ed8ad366f1.png">

<img width="559" alt="grad_evoked.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/f2edf807-85d6-05ce-ca84-e9b00e41bcd9.png">

Evokedオブジェクトはさらに詳しい解析を行うメソッドを備えています．各チャンネルごとの波形を表示すると，前頭のチャンネルで100msで負になり（N100），そこから200msあたりで正のピークを持つ振れ（P200）を確認できます．

```python
aud_evoked.plot_joint(picks='eeg')
```

<img width="554" alt="aud_eeg_evoked.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/e9dea3ea-65e7-1951-a560-efb4864f1d4d.png">

また，トポグラフィーと呼ばれる頭皮上の分布を表示することもできます．ここでは，赤いほど正の，青いほど負の電位を持った脳波が以下のように頭皮上で観測されていることが分かります．

```python
aud_evoked.plot_topomap(times=[0., 0.08, 0.1, 0.12, 0.2], ch_type='eeg')
```
<img width="560" alt="eeg_aud_topo.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/4bd1115d-abec-5f49-fd68-51f0eb64a942.png">

また，条件間の差を見るために`mne.combine_evoked()`によって，2つの条件差のデータを持ったEvokedオブジェクトが生成されます（このとき片方のオブジェクトにマイナス符号をつけます）．このオブジェクトの波形を表示すると，2つの条件間の差の波形を見ることができます．

```python
evoked_diff = mne.combine_evoked([aud_evoked, -vis_evoked], weights='equal')
evoked_diff.pick_types('mag').plot_topo(color='r', legend=False)
```
<img width="425" alt="evoked_diff.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/638747/ac0d9a3a-62c6-6a23-a20d-5430189a2ebe.png">

# おわりに
ここでは基礎的な機能のみ紹介しましたが，MNE公式ページには多くの[Tutorial](https://mne.tools/stable/auto_tutorials/index.html)と[Example](https://mne.tools/stable/auto_examples/index.html)が用意されているので，さらに詳しい分析を行うための最初のステップとして参考にしていただけたらと思います．

また，Pythonには統計や機械学習のライブラリが多くあるので，MEG/EEGの機械学習などによる分析が容易になると思います．