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続き

上記作業でできた返り値から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])