raw.fifのリストを入力して(1)周波数フィルタ (2) epochにして保存 (3) mneでstcに変換して保存する作業をまとめた関数. 下の方に説明をちょっぴりとサンプルスクリプト追加してます

func_raw_list2stc.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

def filter_raw_list(list_fname, bp_freq):
	u"""
	filter raw files. save file name with "l_freq"_"h_freq".fif
	Arguments
		list_fname: list of raw.fif
		bp_freq: band pass frequency. e.g (0,40)
	return
		list_fname_new: list of filtered raw.fif names
		save raw.fif with with "l_freq"_"h_freq".fif
	"""
	import os
	import mne
	# いろいろ設定確認
	# ファイルリスト

	# filter
	l_freq=bp_freq[0]
	h_freq=bp_freq[1]
	list_fname_new = []
	for ii, fname in enumerate(list_fname):
		raw = mne.io.RawFIF(fname, preload=True, proj=False)
		raw.filter(l_freq,h_freq)
		fname_filt=fname[:-4]+'_'+str(l_freq)+'_'+str(h_freq)+'.fif'
		raw.save(fname_filt)
		list_fname_new.append(fname_filt)

	return(list_fname_new)

def raw_list2epochs(list_fname, tmin, tmax, baseline, event_id,reject):
	u"""
	make raw_list into epochs-epo.fif

	Arguments
		tmin: tmin(sec)
		tmax: tmax(sec)
		baseline: baseline .e.g.(None,0)
		reject: reject level. e.g. dict(grad=8000e-13, mag=8e-12)
		event_id: event_id. e.g. dict(fear_lt =1, fear_rt=2)
	Return
		list_fname_epochs: list of file name of -epo.fif
		save -epo.fif as above according to condition name
	"""
	import os
	import mne
	import numpy as np
	# いろいろ設定確認


	# rawからepochをとりだす
	for ii, fname in enumerate(list_fname):
		print(ii)
		raw = mne.io.RawFIF(fname, preload=True, proj=False)
		picks = mne.pick_types(raw.info, meg=True, eeg=False, eog=False, stim=False,
			exclude='bads')

		# トリガー情報取り出し
		min_dur_trig = 0.003 # (sec) epoch取り出し時にひどく短い持続のものはアーチファクトとして扱う
		events = mne.find_events(raw,stim_channel='STI101', min_duration = min_dur_trig) # min_duration以下の持続のeventは除外

		# epochs取り出し
		epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True, picks=picks, baseline = baseline,
			preload=True, reject = reject)
		if ii == 0:
			epochs_concat = epochs.copy()
		else:
			epochs.times = epochs_concat.times
			epochs_concat = mne.epochs.concatenate_epochs([epochs_concat, epochs])

	# save epochs
	dir_raw = os.path.dirname(list_fname[0])
	os.chdir(dir_raw)
	
	list_fname_epochs = []
	for cond in event_id.keys():
		fname_epochs = cond+'-epo.fif'
		epochs_concat[cond].save(fname_epochs)
		list_fname_epochs.append(fname_epochs)

	# finish message
	import time
	print('\nSaving to ...')
	time.sleep(1)
	print([os.path.basename(fpath) for fpath in list_fname_epochs])
	time.sleep(1)
	print('\nin '+os.path.dirname(list_fname_epochs[0]))
	time.sleep(1)
	print('\nfinished saving -epo.fif')

	return list_fname_epochs

def epochs2stc(path_epo, subjects_dir, subject, mri_trans, src, bem):
	u"""
	make epochs-epo.fif into stc file of fsaverage.
	need mri data.
	Arguments
		path_epo: path to -epo.fif
		subjects_dir: path to "subjects_dir"= 2 levels upper of "bem" directory.
		subject: "subject". it may not be necessary.
		mri_trans: usually -trans.fif in -epo.fif
		src: path to "-oct-6-src.fif" usually in "bem" directory.
		bem: path to "-5120-bem-sol.fif" usually in "bem" directory.
	return
		None. save .stc of subject and fsaverage (end with _fsave-lh(rh).stc)
	"""
	import os
	import mne
	import numpy as np
	import matplotlib.pyplot as plt
	import sys

	# epochs: file name(-epo.fif) check
	fname_epo = os.path.basename(path_epo)
	if fname_epo[-8:]!='-epo.fif':
		WARNING = '\033[91m'
		err_message='WARNING: .fif of epochs should be end with "-epo.fif". it will exit...'
		print(WARNING+err_message)
		return
		#sys.exit('hoge')
	# forward solution：エポックで計算
	epochs = mne.read_epochs(path_epo)
	fwd = mne.make_forward_solution(epochs.info, trans=mri_trans, src=src, bem=bem,
		fname=None, meg=True, eeg=False, mindist=5.0, overwrite=False)
	fwd = mne.convert_forward_solution(fwd, surf_ori=True) # need to avoid error below
	#ValueError: You need a free-orientation, surface-oriented forward solution to do depth weighting even when calculating a fixed-orientation inverse.

	# cov: epochsの刺激前区間から計算
	noise_cov = mne.compute_covariance(epochs, tmin=None, tmax=0)

	# inverse operatorの計算
	inv = mne.minimum_norm.make_inverse_operator(epochs.info, fwd, noise_cov, loose=0.2, depth=None)
	# mne.minimum_norm.write_inverse_oeprator(path_epo[:-8]+'-inv.fif',inv) # file save

	# apply inverse : 時系列の信号源データに変換
	#snr_epochs = 1.0 # user lower SNR for single epochs. c.f. http://goo.gl/xOyo0D
	snr = 3.0
	lambda2 = 1.0 / snr ** 2
	method = "dSPM"

	stc = mne.minimum_norm.apply_inverse(epochs.average(), inv, lambda2, method='dSPM',
	 	pick_ori=None) # pick_ori='normal'にすると±ができる
	fname_stc = os.path.basename(path_epo[:-8])
	dir_raw = os.path.dirname(path_epo)
	dir_stc = os.path.join(dir_raw, 'stc')
	if not os.path.exists(dir_stc):
		os.mkdir('stc')
	path_stc = os.path.join(dir_stc, fname_stc)
	stc.save(path_stc) # file save

	# convert to fsaverage
	stc_fsave = stc.morph('fsaverage',subjects_dir=subjects_dir, subject_from=subject)
	path_stc_fsave = os.path.join(dir_stc, fname_stc+'_fsave')
	stc_fsave.save(path_stc_fsave)

	# finish message
	import time
	print('\nfinished saving -lh(rh).stc and _fsave-lh(rh).stc')
	time.sleep(1)
	print('\nin '+ dir_stc)

def epochs_list2stc(list_path_epo, subjects_dir, subject, mri_trans, src, bem):
	u"""
	make epochs-epo.fif into stc file of fsaverage.
	need mri data.
	Arguments
		list_path_epo: list of path to -epo.fif
		subjects_dir: path to "subjects_dir"= 2 levels upper of "bem" directory.
		subject: "subject". it may not be necessary.
		mri_trans: usually -trans.fif in -epo.fif
		src: path to "-oct-6-src.fif" usually in "bem" directory.
		bem: path to "-5120-bem-sol.fif" usually in "bem" directory.
	return
		None. save .stc of subject and fsaverage (end with _fsave-lh(rh).stc)
	"""
	import os
	import mne
	import numpy as np
	import matplotlib.pyplot as plt
	import sys
	
	epochs = mne.read_epochs(list_path_epo[0])
	fwd = mne.make_forward_solution(epochs.info, trans=mri_trans, src=src, bem=bem,
		fname=None, meg=True, eeg=False, mindist=5.0, overwrite=False)
	fwd = mne.convert_forward_solution(fwd, surf_ori=True) # need to avoid error below
	#ValueError: You need a free-orientation, surface-oriented forward solution to do depth weighting even when calculating a fixed-orientation inverse.


	# apply inverse : 時系列の信号源データに変換
	#snr_epochs = 1.0 # user lower SNR for single epochs. c.f. http://goo.gl/xOyo0D
	snr = 3.0
	lambda2 = 1.0 / snr ** 2
	method = "dSPM"

	for path_epo in list_path_epo:
		# epochs: file name(-epo.fif) check
		fname_epo = os.path.basename(path_epo)
		if fname_epo[-8:]!='-epo.fif':
			WARNING = '\033[91m'
			err_message='WARNING: .fif of epochs should be end with "-epo.fif". it will exit...'
			print(WARNING+err_message)
			return
			#sys.exit('hoge')
		# forward solution：エポックで計算
		epochs = mne.read_epochs(path_epo)
		# cov: epochsの刺激前区間から計算
		noise_cov = mne.compute_covariance(epochs, tmin=None, tmax=0)
	
		# inverse operatorの計算
		inv = mne.minimum_norm.make_inverse_operator(epochs.info, fwd, noise_cov, loose=0.2, depth=None)
		# mne.minimum_norm.write_inverse_oeprator(path_epo[:-8]+'-inv.fif',inv) # file save

		stc = mne.minimum_norm.apply_inverse(epochs.average(), inv, lambda2, method='dSPM',
		 	pick_ori=None) # pick_ori='normal'にすると±ができる
		fname_stc = os.path.basename(path_epo[:-8])
		dir_raw = os.path.dirname(path_epo)
		dir_stc = os.path.join(dir_raw, 'stc')
		if not os.path.exists(dir_stc):
			os.mkdir('stc')
		path_stc = os.path.join(dir_stc, fname_stc)
		stc.save(path_stc) # file save

		# convert to fsaverage
		stc_fsave = stc.morph('fsaverage',subjects_dir=subjects_dir, subject_from=subject)
		path_stc_fsave = os.path.join(dir_stc, fname_stc+'_fsave')
		stc_fsave.save(path_stc_fsave)

	# finish message
	import time
	print('\nfinished saving -lh(rh).stc and _fsave-lh(rh).stc')
	time.sleep(1)
	print('\nin '+ dir_stc)

def find_label_in_annot(keyword, labels):
	u"""
		find label in annot (e.g. aparc.a2009s)
		Arguments
		:keyword e.g.'fusiform', 'insula'. keyword must be included in label.name
		:labels. usually obtained by mne.read_labels_from_annot
		Return
		:idx index of the labels including keyword in label.name
	"""
	idx = []
	for ii,label in enumerate(labels):
		if label.name.find(keyword)>-1:
			#print(ii)
			idx.append(ii)
	print(idx)
	return(idx)

Readmeの代わり

raw.fifが分割された状態から、mneをかけてstcファイルを保存するまでの関数をまとめてみました

1. filter_raw_list(list_fname, bp_freq)： raw.fifのリストを入力してフィルタをかけていって、raw_0_40.fifみたいな形で保存してファイル名(リスト)を返す
2. raw_list2epochs(list_fname, tmin, tmax, baseline, event_id,reject): raw.fifのリストからevent毎にまとめて-epo.fifを保存、ファイル名(リスト)を返す
3. epochs2stc(path_epo, subjects_dir, subject, mri_trans, src, bem): -epo.fifからMNEをかけて.stcで保存する
4. epochs_list2stc(list_path_epo, subjects_dir, subject, mri_trans, src, bem): 上と同様だが-epo.fifのリストを入力して一度に処理する
5. find_label_in_annot(keyword, labels): aparc.annotなどから欲しいエリア(keywordで指定)のlabelが何番目にあるか探す.

raw.fifのリストから.stcを作るまでの作業を行うサンプルを作ってみました

step1: raw.fifのリストにフィルタ、step2: step1のフィルタをかけたファイルリストからepochのfif(-epo.fif)を作る、step3: -epo.fifからmneをかけてstcを保存に分けてます

#!/usr/bin/env python
# -*-coding: utf-8 -*-

## import modules
import os,sys, urllib, tempfile
import numpy as np
import matplotlib.pyplot as plt
# import custom functions
from git import Repo
url = 'https://gist.github.com/cosacog/68af94e1ac5a5112dc25' # gistのアドレス
fname_func = 'func_raw_list2stc.py' # 一時保存するときのファイル名
tmp_dir = tempfile.mkdtemp()
Repo.clone_from(url, tmp_dir)
sys.path.append(tmp_dir)
import func_raw_list2stc as func
sys.path.remove(tmp_dir)

## common settings
dir_raw = '/net/dir/to/raw_fif'

########## step1. filter raw.fif ####
list_fname = []
list_fname.append(os.path.join(dir_raw,'raw1_trans_tsss.fif'))
list_fname.append(os.path.join(dir_raw,'raw2_trans_tsss.fif'))
list_fname.append(os.path.join(dir_raw,'raw3_trans_tsss.fif'))

bp_freq=(0,40) # Hz

list_fname_filtered = func.filter_raw_list(list_fname, bp_freq)

###########  step2. raw.fif into epochs(-epo.fif) ####
# 注：ここから始める時は、上のlist_fnameを参考にlist_fname_filteredを作る
# settings
# エポックで取り出す区間(sec)
tmin = -0.2
tmax = 0.5
baseline = (None, 0)

# reject条件
reject = dict(grad=8000e-13, mag=8e-12) # default grad=4000e-13, mag=4e-12

# トリガーに名前を割り当て: 注.以下名前で指示するようになるので、ここの対応で間違えるとリカバリできない
event_id = dict(fear_lt = 1,
                fear_rt = 2,
                fear_ct = 4,
                neut_lt = 8,

list_fname_epochs = func.raw_list2epochs(list_fname_filtered, tmin, tmax, baseline, event_id,reject)

########## step3. epochs into stc ###############
# settings
## mne
subject = 'sample_subj'
subjects_dir = '/net/dir/to/subjects_dir/'
#subjects_dir = 'D:/Neurophysiology2/face_hemifield_inamizu/temp_data160120/'
mri_trans = os.path.join(dir_raw,'meg_mri_coregistration-trans.fif') # -trans.fif: MRI->head coordinate info
src = os.path.join(subjects_dir, subject ,'bem', 'sample_subj-oct-6-src.fif') # source space ?
bem = os.path.join(subjects_dir, subject ,'bem', 'sample_subj-5120-bem-sol.fif')# bem model

#func.epochs2stc(path_epo, subjects_dir, subject, mri_trans, src, bem)
func.epochs_list2stc(list_fname_epochs, subjects_dir, subject, mri_trans, src, bem)