load raw and remontage EEG channels。要mne python。説明は下の方にあります。

load_raw_n_remontage.py

import mne
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import os
'''
EEG
21: Th1, 22: lt iliac crest, 23: rt. iliac crest
24.25: lt abductor hallucis, 26.27: rt abductor hallucis
61: lt.apb, 62: rt.apb
'''
#%% function
def rename_chnames(raw, path_chnames):
    '''
    rename channel names
    '''
    import pandas as pd
    import numpy as np
    # load chnames
    ch_names_series = pd.read_csv(path_chnames, header=None)
    ch_names = [ch_name for ch_name in ch_names_series[0]] # to list
    raw_chnames = raw.ch_names
    if not len(raw_chnames) == len(ch_names):
        raise ValueError("Channel names length does not match.")
    for idx, ch_name in enumerate(ch_names):
        dict_rename = {raw_chnames[idx]:ch_name}
        raw.rename_channels(dict_rename)
    return raw

def get_remontaged_eeg(eeg_inst, remontage_info):
    '''
    return remontage eeg data by 2d np.array
    '''
    isEpochs = _judge_inst_class(eeg_inst)
    import numpy as np
    eeg_data_remontage = _get_remontaged_eeg(eeg_inst, remontage_info, isEpochs)
    info = _get_remontage_info(eeg_inst, remontage_info)
    if isEpochs:
        eeg_remontage = mne.EpochsArray(eeg_data_remontage, info, tmin=eeg_inst.tmin)
        eeg_remontage.events = eeg_inst.events
        eeg_remontage.event_id = eeg_inst.event_id
    else:
        eeg_remontage = mne.io.RawArray(eeg_data_remontage, info)
    return eeg_remontage

def remontage_eeg(raw, path_remontage_info):
    '''
    remontage eeg channels. basically bipolar
    remove eeg channels and replace by remontaged channels
    rename channel names
    '''
    import pandas as pd
    import numpy as np
    # load raw-> copy
    raw_meg = raw.copy()
    picks_meg = mne.pick_types(raw_meg.info, meg=True, eeg=False, stim=False)
    raw_meg.pick_channels([raw_meg.ch_names[pick] for pick in picks_meg])
    raw_stim = raw.copy()
    picks_stim = mne.pick_types(raw_stim.info, meg=False, eeg=False, stim=True)
    raw_stim.pick_channels([raw_stim.ch_names[pick] for pick in picks_stim])
    raw_eeg = raw.copy()
    picks_stim = mne.pick_types(raw_eeg.info, meg=False, eeg=True, stim=False)
    raw_eeg.pick_channels([raw_eeg.ch_names[pick] for pick in picks_stim])
    # load remontage info
    remontage_info = pd.read_csv(path_remontage_info)
    raw_eeg_remontaged = get_remontaged_eeg(raw_eeg, remontage_info)
    raw_meg.add_channels([raw_eeg_remontaged, raw_stim])
    return raw_meg

def _get_g2data(eeg_data, eeg_ch_names, g2list, isEpochs):
    '''
    get G2 data. averaged if multiple channels are included
    '''
    if g2list[0]=='None':
        return 0
    g2data = []
    for idx, item_g2 in enumerate(g2list):
        if not item_g2 in eeg_ch_names:
            raise ValueError("G2:{0} is not in channel list.".format(item_g2))
        if isEpochs:
            g2_ch_data = eeg_data[:,eeg_ch_names.index(item_g2),:]
        else:
            g2_ch_data = eeg_data[eeg_ch_names.index(item_g2),:]
        # summation
        if idx==0:
            g2data = g2_ch_data
        else:
            g2data += g2_ch_data
    g2data /= idx+1
    return g2data

def _get_g2list(g2item):
    '''
    get G2 channel list
    '''
    g2list = g2item.split(',')
    g2list[0] = g2list[0].split('[')[-1]
    g2list[-1] = g2list[-1].split(']')[0]
    g2list = [s.strip() for s in g2list]
    g2list = [s.strip() for s in g2list]    
    return g2list

def _get_remontaged_eeg(inst, remontage_info, isEpochs):
    '''
    calculate remontaged eeg
    '''
    eeg_data = inst.get_data()
    eeg_ch_names = inst.ch_names
    axis_concat = 1 if isEpochs else 0
    for idx, item in remontage_info.iterrows():
        print("Processing {0}th channel ...".format(idx))
        g1item = item['G1']
        if not g1item in eeg_ch_names:
            raise ValueError("G1:{0} is not in channel list.".format(g1item))
        idx_g1_ch = eeg_ch_names.index(g1item)
        g2item = item['G2']
        g2list = _get_g2list(g2item)
        g2data = _get_g2data(eeg_data, eeg_ch_names, g2list, isEpochs)
        g1data = eeg_data[:,idx_g1_ch,:] if isEpochs else eeg_data[idx_g1_ch,:]
        eeg_ch_data = g1data - g2data
        if isEpochs:
            epoch_shape = eeg_ch_data.shape
            eeg_ch_data_reshape = eeg_ch_data.reshape((epoch_shape[0], 1, epoch_shape[1]))
        else:
            eeg_ch_data_reshape = eeg_ch_data.reshape((1, len(eeg_ch_data)))
        if idx==0:
            eeg_data_remontage = eeg_ch_data_reshape
        else:
            eeg_data_remontage = np.concatenate([eeg_data_remontage, eeg_ch_data_reshape],axis=axis_concat)
    return eeg_data_remontage

def _get_remontage_info(inst, remontage_info):
    '''
    get inst.info for remontaged data
    '''
    info = inst.info
    info_premitive = mne.create_info(list(remontage_info['new_name']), info['sfreq'], ch_types='eeg')
    info['ch_names'] = info_premitive['ch_names']
    info['chs'] = info_premitive['chs']
    info['nchan'] = info_premitive['nchan']
    return info

def _judge_inst_class(inst):
    '''
    '''
    # judge class of inst
    classes = [mne.epochs.Epochs, mne.io.fiff.raw.Raw]
    if not type(inst) in classes:
        raise ValueError("Data must be epochs or raw. Type(data):{0}, which must be either {1}"
            .format(type(inst), classes))
    isEpochs = type(inst)==classes[0] # epochs
    return isEpochs

load_raw_n_remontage.py: raw.fif を読み込んで脳波チャンネルのリモンタージュ

• raw.fifの脳波チャンネルをリモンタージュします。いわゆるbipolarのモンタージュを組みやすいように考えました。

前提

• MNE python
• gistをインポートするlibrary-> こちら

ファイルの内容

• load_raw_n_remontage.py: スクリプト本体：これは基本変更しません。
• ch_names.csv: チャンネル名変更のリスト.64chの一部をリネームしてます。サンプルはこちら。適宜変更してください。
• remontage_info.csv: リモンタージュの情報.ch_names.csvの名前を基に組み合わせを作っています。サンプルはこちら。適宜変更してください。
• ch_names.csvの基になるチャンネル名リストのcsvは、

import pandas as pd
raw = mne.io.read_raw_fif('/path/to/raw.fif')
pd.Series(raw.ch_names).to_csv('/path/to/ch_names_orig.csv', index=False)

で作成できます。

使い方

• raw.fifの読み込み-> チャンネル名のリネーム -> リモンタージュの流れにしてます。
• チャンネル名のリネームは必須ではないですが、EEG001, EEG002...と言うのは使い勝手が悪いので、リネームした方がよいでしょう。

import os
import mne
from import_gist import *
# info
dir_raw = '/dir/to/raw'
fname_raw = 'raw.fif'
path_raw = os.path.join(dir_raw, fname_raw)
path_chnames = '/path/to/ch_names.csv'
path_remontage_info = '/path/to/remontage_info.csv'
url_gist = 'https://gist.github.com/cosacog/7b593940742233fb49e000940ce5b215' # このgist
# load gist script
mdl = import_gist(url_gist) # ref. https://github.com/cosacog/import_gist
# load raw.fif
raw = mne.io.read_raw_fif(path_raw, preload=True)
# rename channels
raw_renamed = mdl.rename_chnames(raw, path_chnames)
# remontage eegs
raw_remontaged = mdl.remontage_eeg(raw, path_remontage_info)
# raw_remontaged.plot()で確認