labelをいじるfunction、label内の時系列データのうち活動の多い方から10%のvertexの平均をプロットするfunction

mne_stc_misc.py

def half_ydir_label(label_orig):
    '''
    要mne等
    labelをy軸方向に半分にするスクリプト
    params: label_orig. 元のlabel
    return: label_half. y軸方向に半分になったlabel
    '''
    ypos_label = label_orig.pos[:,1]
    idx_half_label = np.where(ypos_label >= np.median(ypos_label)) # example: y positions are anterior

    pos_label_half = label_orig.pos[idx_half_label,:]
    vertices_half = label_orig.vertices[idx_half_label]
    values_half = label_orig.values[idx_half_label]
    name_half = label_orig.name[:-3] + '_half'+label_orig.name[-3:]

    label_half = label_orig.copy()
    label_half.pos = pos_label_half
    label_half.vertices = vertices_half
    label_half.values = values_half
    label_half.name= name_half

    return label_half



def plot_conditions_in_label_10percentile(list_path_stc, list_dict_labels, percentile_vertices, t_range, src, subjects_dir, is_label_smooth):
    '''
    list_path_stcのファイルから読み込んだデータを重ね書き(plot)するfunction
    list_dict_labelsに入ったラベルの領域毎にsubplotでプロット
    適当(0.1-0.2 sec)の最大振幅の時刻を基準に10 percentileのverticesを取得して平均振幅を描画

    <params> list_path_stc: stc fileのlist
    <params> list_dict_labels: dict(label と title)のlist
    <params> percentile_vertices: label内verticesから読み込む数. percentile表示(e.g. 0.1で10%)
    <params> t_range: 最大の時刻, vertexを検出する時間幅. e.g. (0.1, 0.2)
    <params> src: いわゆるsrc. ただし中身は mne.add_source_space_distances 後のsrcにしとかないと計算が20分くらいかかる
    <params> subjects_dir: いわゆるsubjects_dir.
    <params> is_label_smooth: True or False. 通常Trueでよいが、Falseでstcの中のvertexの位置が具体的にわかる

    return: subplotで各labelの中のstcデータを平均してプロット。それぞれlabel内で大きい方から10 percentileの頂点を選ぶ
    <out> array_mean_percentile: 平均波形のarray, labelとファイルの数だけ返って来る
    <out> header_csv: 各列のlabel,ファイル名がわかるようにしたheader
    <out> list_label: 検出したlabelのlist, labelとファイルの順番は上記header_csvに準じる
    <out> list_latency: 検出した頂点潜時のlist, 順番は上に同じ
    '''
    # settings
    time_range4peak_detect = t_range
    colors = ['red','black','blue','green'] # plot colors: ファイルの種類に対応する

    # main
    import numpy as np
    import matplotlib.pyplot as plt
    import mne
    header_csv = 'time'
    idx = 0
    list_label = []
    list_latency = []
    for idx_plot, dict_label_title in enumerate(list_dict_labels):
        label = dict_label_title['label']
        title = dict_label_title['title']
        idx_hemi = 0 if label.hemi == 'lh' else 1

        plt.subplot(int(np.ceil(len(list_dict_labels)/2)), 2, idx_plot+1)
        legend_text = [os.path.basename(path)[:-4] for path in list_path_stc]
        for ii, path_stc in enumerate(list_path_stc):
            stc = mne.read_source_estimate(path_stc)
            stc_in_label = stc.in_label(label)
            n_vertices = stc_in_label.data.shape[0]
            stc_data = stc_in_label.data
            idx_toi = np.where((stc_in_label.times >= time_range4peak_detect[0]) &
                               (stc_in_label.times <= time_range4peak_detect[1]))[0]
            # detect peak latency and vertex
            # idx_vertex_peak, idx_time_peak = np.where(stc_data[:,idx_toi] == np.max(stc_data[:,idx_toi]))
            stc_mean_data_toi = np.mean(stc_data[:,idx_toi], axis=0)
            idx_time_peak = np.where(stc_mean_data_toi == np.max(stc_mean_data_toi))[0] # 全vertex平均の最大時刻のindex
            idx_vertices_peak = np.argsort(stc_data[:,idx_time_peak[0]],axis=0)[-int(np.ceil(n_vertices * percentile_vertices)):]
            mean_stc_data_percentile = np.mean(stc_data[idx_vertices_peak,],axis=0) #

            # detected latency
            latency_detected = stc_in_label.times[idx_toi[0] + idx_time_peak][0]
            list_latency.append(latency_detected)

            # plot
            plt.plot(stc.times, mean_stc_data_percentile, color = colors[ii])
            plt.xlim(stc.times[0],stc.times[-1])
            plt.title(title)

            # data concatenate: 三項演算子
            array_mean_percentile = mean_stc_data_percentile if idx == 0 else np.r_['1,2,0', array_mean_percentile, mean_stc_data_percentile]

            # csv header
            header_csv += ',' + title + '-' + os.path.basename(path_stc)[:-7]

            # label
            stc_in_label.data[:,:] = 0
            stc_in_label.data[idx_vertices_peak,0]=1 #ここだけlabelになる
            label_peak_both = mne.stc_to_label(stc_in_label, src=src, smooth=is_label_smooth, connected=False, subjects_dir=subjects_dir)
            label_peak_hemi = label_peak_both[idx_hemi] # rh or lh only
            label_peak_hemi.subject = stc.subject
            label_peak_hemi.name = os.path.basename(path_stc)[:-7] + ' ' + title
            list_label.append(label_peak_hemi)

            idx += 1

    plt.legend(legend_text)
    array_mean_percentile = np.r_['1,2,0', stc.times, array_mean_percentile]

    return array_mean_percentile, header_csv, list_label, list_latency

def find_label_in_annot(label_name, labels):
    '''
    mne.read_labels_from_annotで読みだしたlabelsからほしい名前のlabelの番号を取り出す
    <params>
        label_name: ほしいlabelの名前の一部'fusiform'とか
        labels: mne.read_labels_from_annotで読みだしたlabels
    <return>
        labels_sub: 取り出したlabel
        idx_labels: 取り出したlabelの番号.labels[idx_labels]とかするとほしいlabelが出せるはず
    '''
    import numpy as np
    labels_sub = []
    idx_labels = np.array([])
    ii = 0
    for label in labels:
        idx = label.name.find(label_name)
        if idx > -1:
            labels_sub.append(label)
            idx_labels = np.append(idx_labels, ii)
        ii += 1
    return labels_sub, idx_labels
    

plot_conditions_in_label_10percentile

MNEでstcファイルとlabelオブジェクトからlabel内の活動を取り出し、振幅の高いvertex10%分の平均波形をプロットします

注意点: 「活動が高い」の判定法

• 0.1-0.2秒 (t_rangeで設定)の区間で最大振幅のvertexと時刻を検出
• label内の全vertexの平均波形を出して、その最大振幅の時刻を同定
• 上述で検出された時刻における振幅の高いvertex 10%を取得
  という流れです

改変履歴

• 平均波形データとそのデータの各列に対応する条件、10% vertexのラベルを返り値で出すようにしました (160406). ちと引数増えすぎてすいません.
• 返り値に検出した潜時のリストを出すようにしました(160603)
• 活動の高い時刻の判定法を変更しました(注意点参照) (160606)

式

array_mean_percentile, header_csv, list_label, list_latency = plot_conditions_in_label_10percentile(list_path_stc, list_dict_labels,
 percentile_vertices, t_range, src, subjects_dir, is_label_smooth)

引数

• list_path_stc: .stcファイルパスのリスト. トラブル防止のため絶対パスがよいでしょう
• list_dict_labels: dict(辞書オブジェクト)のリスト. label, titleのkeyあり. labelはmne.read_labelとかで取得できるlabelオブジェクト. titleは対応するプロットのタイトル
• percentile_vertices: 0-1の値. vertexのpercentile値. 0.1だと10 %のverticesを取得する
• t_range: 最大の時刻、vertexを検出する解析の時間幅. e.g. (0.1, 0.2)で0.1-0.2秒の意味
• src: いわゆるsrc (c.f. mne.read_source_spaces). mne.add_source_space_distancesで処理したsrc.fifを読み込んだsrc推奨: 計算が+20分余計にかかるようになります
• subjects_dir: いわゆるsubjects_dir
• is_label_smooth: True or False, Trueが通常. Falseにすると信号源推定に使われたvertexが点で表示される. わかりにくいので実行した方が早いです. 一度は試してみましょう

返り値: プロットが出力されますが、同時に返り値があります

• array_mean_percentile: 平均波形の二次元配列データ. 行が時刻、列が各条件
• header_csv: 各列の条件、ラベルを示すheader. csvに出力するときに使うとよいでしょう
• list_label: 検出した頂点のラベルのリスト. brain.add_label(list_label[0])とかするとラベルが表示できます.名前はheader_csvの中身を参照してください
• list_latency (160603追加): 検出した頂点潜時のリスト.順番はlist_labelと同じ

使い方の例 ： コメントにも注意してご利用ください. プロットは2列になります

import os
import mne 
import numpy as np
import matplotlib.pyplot as plt

# 設定: labels and titles
labels = [label_lv1, label_rv1, label_llo, label_rlo, label_lffa, label_rffa] # labelオブジェクトはmne.read_labelとか利用、個数は任意
titles = ['lt V1','rt V1','lt LO','rt LO','lt Fusiform','rt Fusiform'] # labelsに対応するタイトル(プロットで使う)、labelsの数に合わせる
list_dict_labels = [dict(label = label, title = title) for (label,title) in zip(labels, titles)] #これがfunctionの引数になる

# 設定: stcファイル, ディレクトリ
dir_raw = '/directory/to/mne/de/tsukutta/stc'
list_fname_stc = ['fear_lt_fsave-rh.stc','neut_lt_fsave-rh.stc','objt_lt_fsave-rh.stc'] # 当座ファイルは4つ以内(38行目. colorsの数に依存)
list_path_stc = [os.path.join(dir_raw, fname) for fname in list_fname_stc] # これが引数になる：stcファイルの絶対パスのリスト

# 解析時間幅
t_range = (0.1, 0.2) # 秒単位
p_vertices = 0.1 # 10 パーセンタイルだったら0.1

# ラベルの処理に必要な設定: src, subjects_dir, is_label_smooth
src = mne.read_source_spaces('/dir/to/subjects_dir/subject/bem/hogehoge-src.fif')
subjects_dir = '/dir/to/subjects_dir'
is_label_smooth = True

# メインの処理: 下の1行のみ
array_mean_plot_data, header_csv, list_labels, list_latencies = plot_conditions_in_label_10percentile(list_path_stc, list_dict_labels, p_vertices, t_range, src, subjects_dir, is_label_smooth)

以下のようなプロットができたら成功です. csvの出力、ラベルの利用法は次のcommentにします.

plot_conditions_in_label_10percentile続き

上記作業でできた返り値からcsvを出力する流れ、list_labelsの利用法を示します.

csvの出力 (平均波形のデータ)：上のスクリプトに続く形で

fname_csv = 'fuga.csv' # csvのファイル名.自分で後で思い出せるように名前をつけましょう.
dir_csv = os.path.dirname(list_path_stc[0]) # stcと同じディレクトリにcsvを保存する時
path_csv = os.path.join(dir_csv, fname_csv) 
np.savetxt(path_csv, array_mean_plot_data, delimiter=',', header= header_csv, fmt='%.18f')

とするとfuga.csvがdir_csvの場所に保存されます. excel等で開けますので、さらなる解析、データ提示に活用できます.

list_labelsの利用法: 検出された活動の高いvertexからなるラベルを観察するには

list_labels とか返り値の名前をタイプすると中身が見られますが、labelがまとまった感じ(list)になっていると思います。
各ラベルがどのラベル(上述で生じるプロットのタイトルを利用してます)、どのstcファイル由来かわかるように名前をつけています.

brain = stc.plot(云々)とかで表示した状態で

brain.add_label(list_labels[0], color='red', alpha = 0.5) # alphaは透明度

とかすると該当のラベルが表示されます. 色と番号を変えながら重ね書きすると評価しやすいでしょう.
python的数え方になるので1番目のラベルは0になります.

なお

# ラベルを消したいとき
brain.remove_labels(hemi='lh') 

.annotファイルからほしいlabelのindexを探すfunction: find_label_in_annot

要mne, numpy
ほしいエリアにどのような名前のラベルがついてるかはfreeviewあたりで調べましょう.

# *.annotファイルを読みだしてlabelsの中に代入
labels = mne.read_labels_from_annot(subject='fsaverage', parc='aparc.a2009s', hemi='both', surf_name='inflated', subjects_dir= subjects_dir)

# labelsの中から'cuneus'の単語が入っているlabelsのindexを探す
labels_sub, idxs = find_label_in_annot('cuneus', labels)

#　結果確認
print(labels)
print(labels[idxs])