Python 2.7です。

In [1]:
import sys
print(sys.version_info)
del sys
sys.version_info(major=2, minor=7, micro=13, releaselevel='final', serial=0)

結合性解析の例

In [2]:
# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>
#
# License: BSD (3-clause)
In [3]:
import numpy as np
import mne
mne.set_log_level('INFO')
from mne import io
from mne.connectivity import spectral_connectivity, seed_target_indices
from mne.datasets import sample
from mne.time_frequency import AverageTFR
print(__doc__) # doc stringというらしい 関数の説明らしいけど 深く考えないことにする。
Automatically created module for IPython interactive environment

パラメータの設定

In [4]:
data_path=sample.data_path()
raw_fname=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw.fif'
event_fname=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
# rawデータとeventの読み込み
raw=io.read_raw_fif(raw_fname)
events=mne.read_events(event_fname)
# 不良チャンネルの選択
raw.info['bads']+=['MEG 2443']
# gradiomterだけ選択
picks=mne.pick_types(
    raw.info, meg='grad', eeg=False, stim=False, eog=True, exclude='bads')
# ecposhクラスの作成
event_id,tmin,tmax=3,-0.2,0.5
epochs=mne.Epochs(
    raw,events,event_id,tmin,tmax,picks=picks,baseline=(None,0),
    reject=dict(grad=4000e-13,eog=150e-6),preload=True)
# MEG 2343をseedとする
seed_ch='MEG 2343'
picks_ch_names=[raw.ch_names[i] for i in picks]
# seed-target indicesの計算
seed=picks_ch_names.index(seed_ch)
targets=np.arange(len(picks))
indices=seed_target_indices(seed,targets) # tuple 2 102と102
# waveletの定義
cwt_frequencies=np.arange(7,30,2)
cwt_n_cycles=cwt_frequencies/7.
# connectivityの計算 並列処理は2
sfreq=raw.info['sfreq']
con,freqs,times,_,_=spectral_connectivity( # conは204x11x106の三次元配列
    epochs,indices=indices,method='wpli2_debiased',mode='cwt_morlet',sfreq=sfreq,
    cwt_frequencies=cwt_frequencies,cwt_n_cycles=cwt_n_cycles,n_jobs=1)
Opening raw data file C:\Users\akira\mne_data\MNE-sample-data\MEG\sample\sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Current compensation grade : 0
73 matching events found
4 projection items activated
Loading data for 73 events and 106 original time points ...
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
6 bad epochs dropped
Connectivity computation...
    computing connectivity for 204 connections
    using t=0.000s..0.699s for estimation (106 points)
    frequencies: 9.0Hz..29.0Hz (11 points)
    using CWT with Morlet wavelets to estimate spectra
    the following metrics will be computed: Debiased WPLI Square
    computing connectivity for epoch 1
    computing connectivity for epoch 2
    computing connectivity for epoch 3
    computing connectivity for epoch 4
    computing connectivity for epoch 5
    computing connectivity for epoch 6
    computing connectivity for epoch 7
    computing connectivity for epoch 8
    computing connectivity for epoch 9
    computing connectivity for epoch 10
    computing connectivity for epoch 11
    computing connectivity for epoch 12
    computing connectivity for epoch 13
    computing connectivity for epoch 14
    computing connectivity for epoch 15
    computing connectivity for epoch 16
    computing connectivity for epoch 17
    computing connectivity for epoch 18
    computing connectivity for epoch 19
    computing connectivity for epoch 20
    computing connectivity for epoch 21
    computing connectivity for epoch 22
    computing connectivity for epoch 23
    computing connectivity for epoch 24
    computing connectivity for epoch 25
    computing connectivity for epoch 26
    computing connectivity for epoch 27
    computing connectivity for epoch 28
    computing connectivity for epoch 29
    computing connectivity for epoch 30
    computing connectivity for epoch 31
    computing connectivity for epoch 32
    computing connectivity for epoch 33
    computing connectivity for epoch 34
    computing connectivity for epoch 35
    computing connectivity for epoch 36
    computing connectivity for epoch 37
    computing connectivity for epoch 38
    computing connectivity for epoch 39
    computing connectivity for epoch 40
    computing connectivity for epoch 41
    computing connectivity for epoch 42
    computing connectivity for epoch 43
    computing connectivity for epoch 44
    computing connectivity for epoch 45
    computing connectivity for epoch 46
    computing connectivity for epoch 47
    computing connectivity for epoch 48
    computing connectivity for epoch 49
    computing connectivity for epoch 50
    computing connectivity for epoch 51
    computing connectivity for epoch 52
    computing connectivity for epoch 53
    computing connectivity for epoch 54
    computing connectivity for epoch 55
    computing connectivity for epoch 56
    computing connectivity for epoch 57
    computing connectivity for epoch 58
    computing connectivity for epoch 59
    computing connectivity for epoch 60
    computing connectivity for epoch 61
    computing connectivity for epoch 62
    computing connectivity for epoch 63
    computing connectivity for epoch 64
    computing connectivity for epoch 65
    computing connectivity for epoch 66
    computing connectivity for epoch 67
[Connectivity computation done]
In [5]:
# クラスepochsをfif形式で保存する場合
filename=data_path+'/MEG/sample/aaa-epo.fif'
epochs.save(filename)
# MATLABのMAT形式で保存する場合
import scipy.io
out_data={}
out_data['eopch']=epochs.get_data()
savename=data_path+'/MEG/sample/aaa.mat'
scipy.io.savemat(savename,out_data)
In [6]:
# 描画
# %matplotlib inline # いらないみたい
con[np.where(indices[1]==seed)]=1.0
title='WPLI2-Visual-Seed %s' % seed_ch
layout=mne.find_layout(epochs.info, 'meg')
tfr=AverageTFR(epochs.info,con, times,freqs, len(epochs))
print(type(tfr)) # クラス名表示
tfr.plot_topo(fig_facecolor='w',font_color='k',border='k');
<class 'mne.time_frequency.tfr.AverageTFR'>
No baseline correction applied

ラベル指定した特定の信号源空間のconnectivity

In [7]:
% reset
Once deleted, variables cannot be recovered. Proceed (y/[n])? y
In [8]:
# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
#
# License: BSD (3-clause)
In [9]:
import matplotlib.pyplot as plt
import mne
mne.set_log_level('INFO')
from mne.datasets import sample
from mne.minimum_norm import apply_inverse_epochs, read_inverse_operator
from mne.connectivity import spectral_connectivity
print(__doc__)
data_path=sample.data_path()
subjects_dir=data_path+'/subjects'
fname_inv=data_path+'/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
fname_raw=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw.fif'
fname_event=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
# データ読み込み 空間フィルタ設計済み
inverse_operator=read_inverse_operator(fname_inv)
print(type(inverse_operator)) # クラスです
print(inverse_operator.keys()) # いろいろなクラスの集まりです
raw= mne.io.read_raw_fif(fname_raw)
events=mne.read_events(fname_event)
# 不良チャンネル
raw.info['bads']+=['MEG 2443']
# 脳磁図チャンネルと眼電図のみ利用
picks=mne.pick_types(raw.info,meg=True,eeg=False,stim=False,eog=True,exclude='bads')
# エポッククラスの作成
event_id,tmin,tmax=1,-0.2,0.5
epochs=mne.Epochs(
    raw,events,event_id,tmin,tmax,picks=picks,baseline=(None,0),
    reject=dict(mag=4e-12,grad=4000e-13,eog=150e-6))
# 信号源空間のエポックごとの計算
snr=1.0
lambda2=1.0/snr**2
method="dSPM"  # MNE or sLORETAでもOK
# return_generator=Trueが大事らしい
stcs=apply_inverse_epochs(
    epochs,inverse_operator,lambda2,method,pick_ori="normal",return_generator=True) 
print(stcs) # generatorだそうです
# Label読み込み
names=['Aud-lh','Aud-rh','Vis-lh','Vis-rh']
labels=[mne.read_label(data_path+'/MEG/sample/labels/%s.label' % name)for name in names]
# 加算平均
src=inverse_operator['src'] # クラスSourceSpaces 要は三角メッシュ
label_ts=mne.extract_label_time_course(
    stcs,labels,src,mode='mean_flip',return_generator=True)
print(label_ts) # generatorだそうです
fmin,fmax= 5.,40.
sfreq=raw.info['sfreq']
con,freqs,times,n_epochs,n_tapers=spectral_connectivity(
    label_ts,method='wpli2_debiased',mode='multitaper',sfreq=sfreq,
    fmin=fmin,fmax=fmax,mt_adaptive=True,n_jobs=1) # n_jobs=2だと遅くなる気がする
print(con.shape) # 4x4x25
n_rows,n_cols=con.shape[:2]
Built-in functions, exceptions, and other objects.

Noteworthy: None is the `nil' object; Ellipsis represents `...' in slices.
Reading inverse operator decomposition from C:\Users\akira\mne_data\MNE-sample-data/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif...
    Reading inverse operator info...
    [done]
    Reading inverse operator decomposition...
    [done]
    305 x 305 full covariance (kind = 1) found.
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Noise covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 2) found.
    Source covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 6) found.
    Orientation priors read.
    22494 x 22494 diagonal covariance (kind = 5) found.
    Depth priors read.
    Did not find the desired covariance matrix (kind = 3)
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    2 source spaces read
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Source spaces transformed to the inverse solution coordinate frame
<class 'mne.minimum_norm.inverse.InverseOperator'>
['nave', 'methods', 'noisenorm', 'eigen_fields', 'units', 'fmri_prior', 'noise_cov', 'orient_prior', 'mri_head_t', 'sing', 'nsource', 'info', 'src', 'nchan', 'source_cov', 'projs', 'source_nn', 'depth_prior', 'whitener', 'source_ori', 'coord_frame', 'proj', 'eigen_leads_weighted', 'eigen_leads', 'reginv']
Opening raw data file C:\Users\akira\mne_data\MNE-sample-data\MEG\sample\sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Current compensation grade : 0
72 matching events found
Created an SSP operator (subspace dimension = 3)
4 projection items activated
<generator object _apply_inverse_epochs_gen at 0x000000000C199990>
<generator object _gen_extract_label_time_course at 0x000000000C199F30>
Connectivity computation...
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 1
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 3)
    Created the whitener using a full noise covariance matrix (3 small eigenvalues omitted)
    Computing noise-normalization factors (dSPM)...
[done]
Picked 305 channels from the data
Computing inverse...
(eigenleads need to be weighted)...
Processing epoch : 1
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for 6 connections
    using t=0.000s..0.699s for estimation (106 points)
fmin corresponds to less than 5 cycles, spectrum estimate will be unreliable
    frequencies: 5.7Hz..39.7Hz (25 points)
    using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: Debiased WPLI Square
    computing connectivity for epoch 1
Processing epoch : 2
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 2
Processing epoch : 3
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 3
Processing epoch : 4
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 4
Processing epoch : 5
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 5
Processing epoch : 6
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 6
Processing epoch : 7
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 7
Processing epoch : 8
Extracting time courses for 4 labels (mode: mean_flip)
<ipython-input-9-e830aec1d221>:48: RuntimeWarning: fmin corresponds to less than 5 cycles, spectrum estimate will be unreliable
  fmin=fmin,fmax=fmax,mt_adaptive=True,n_jobs=1) # n_jobs=2だと遅くなる気がする
    computing connectivity for epoch 8
Processing epoch : 9
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 9
Processing epoch : 10
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 10
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 11
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 11
Processing epoch : 12
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 12
Processing epoch : 13
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 13
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 14
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 14
Processing epoch : 15
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 15
Processing epoch : 16
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 16
Processing epoch : 17
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 17
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 18
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 18
Processing epoch : 19
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 19
Processing epoch : 20
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 20
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 21
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 21
Processing epoch : 22
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 22
Processing epoch : 23
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 23
    Rejecting  epoch based on MAG : [u'MEG 1711']
Processing epoch : 24
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 24
Processing epoch : 25
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 25
Processing epoch : 26
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 26
Processing epoch : 27
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 27
Processing epoch : 28
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 28
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 29
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 29
Processing epoch : 30
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 30
Processing epoch : 31
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 31
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 32
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 32
Processing epoch : 33
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 33
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 34
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 34
Processing epoch : 35
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 35
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 36
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 36
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 37
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 37
Processing epoch : 38
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 38
Processing epoch : 39
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 39
Processing epoch : 40
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 40
Processing epoch : 41
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 41
Processing epoch : 42
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 42
Processing epoch : 43
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 43
Processing epoch : 44
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 44
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 45
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 45
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 46
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 46
Processing epoch : 47
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 47
Processing epoch : 48
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 48
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 49
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 49
Processing epoch : 50
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 50
Processing epoch : 51
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 51
Processing epoch : 52
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 52
Processing epoch : 53
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 53
Processing epoch : 54
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 54
Processing epoch : 55
Extracting time courses for 4 labels (mode: mean_flip)
    computing connectivity for epoch 55
[done]
    assembling connectivity matrix
[Connectivity computation done]
(4L, 4L, 25L)
In [10]:
# 描画
fig,axes=plt.subplots(n_rows,n_cols,sharex=True,sharey=True)
plt.suptitle('Between labels connectivity')
for i in range(n_rows):
    for j in range(i + 1):
        if i==j:
            axes[i,j].set_axis_off()
            continue
        axes[i,j].plot(freqs,con[i,j,:])
        axes[j,i].plot(freqs,con[i,j,:])
        if j==0:
            axes[i,j].set_ylabel(names[i])
            axes[0,i].set_title(names[i])
        if i==(n_rows-1):
            axes[i,j].set_xlabel(names[j])
        axes[i,j].set_xlim([fmin,fmax])
        axes[j,i].set_xlim([fmin,fmax])
        # band limits
        for f in [8,12,18,35]:
            axes[i,j].axvline(f,color='k')
            axes[j,i].axvline(f,color='k')
plt.show()

センサ空間(信号源)のconnectivity

In [11]:
% reset
Once deleted, variables cannot be recovered. Proceed (y/[n])? y
In [12]:
# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>
#
# License: BSD (3-clause)
In [13]:
import numpy as np
from scipy import linalg
import mne
mne.set_log_level('INFO')
from mne import io
from mne.connectivity import spectral_connectivity
from mne.datasets import sample
print(__doc__)
# ファイル名設定
data_path=sample.data_path()
raw_fname=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw.fif'
event_fname=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
# ファイル読み込み
raw=io.read_raw_fif(raw_fname)
events=mne.read_events(event_fname)
# 不良チャンネル
raw.info['bads']+=['MEG 2443']
# gradiometerと眼電図のみ選択
picks=mne.pick_types(raw.info,meg='grad',eeg=False,stim=False,eog=True,exclude='bads')
# エポッククラスの作成
event_id,tmin,tmax=3,-0.2,0.5
epochs=mne.Epochs(
    raw,events,event_id,tmin,tmax,picks=picks,baseline=(None,0),
    reject=dict(grad=4000e-13,eog=150e-6))
# 誘発応答を含んだまま結合性の計算
# 基線時間は省く
fmin,fmax=3.,9.
sfreq=raw.info['sfreq']  
tmin=0.0
con,freqs,times,n_epochs,n_tapers=spectral_connectivity(
    epochs, method='pli',mode='multitaper',sfreq=sfreq,fmin=fmin,fmax=fmax,
    faverage=True,tmin=tmin,mt_adaptive=False,n_jobs=1)
print(con.shape) # 204x204x1
# 眼磁図雑音エポックを除去
ch_names=epochs.ch_names
idx=[ch_names.index(name) for name in ch_names if name.startswith('MEG')]
con=con[idx][:,idx]
# conは三次元配列
con = con[:, :, 0]
Built-in functions, exceptions, and other objects.

Noteworthy: None is the `nil' object; Ellipsis represents `...' in slices.
Opening raw data file C:\Users\akira\mne_data\MNE-sample-data\MEG\sample\sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Current compensation grade : 0
73 matching events found
4 projection items activated
Connectivity computation...
    computing connectivity for 20706 connections
    using t=0.000s..0.699s for estimation (106 points)
fmin corresponds to less than 5 cycles, spectrum estimate will be unreliable
    frequencies: 4.2Hz..8.5Hz (4 points)
    connectivity scores will be averaged for each band
    using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: PLI
    computing connectivity for epoch 1
    computing connectivity for epoch 2
    computing connectivity for epoch 3
    computing connectivity for epoch 4
    computing connectivity for epoch 5
    computing connectivity for epoch 6
<ipython-input-13-7ca924564fea>:32: RuntimeWarning: fmin corresponds to less than 5 cycles, spectrum estimate will be unreliable
  faverage=True,tmin=tmin,mt_adaptive=False,n_jobs=1)
    computing connectivity for epoch 7
    computing connectivity for epoch 8
    computing connectivity for epoch 9
    computing connectivity for epoch 10
    Rejecting  epoch based on EOG : [u'EOG 061']
    computing connectivity for epoch 11
    computing connectivity for epoch 12
    computing connectivity for epoch 13
    computing connectivity for epoch 14
    computing connectivity for epoch 15
    computing connectivity for epoch 16
    computing connectivity for epoch 17
    computing connectivity for epoch 18
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    computing connectivity for epoch 19
    computing connectivity for epoch 20
    computing connectivity for epoch 21
    computing connectivity for epoch 22
    computing connectivity for epoch 23
    computing connectivity for epoch 24
    computing connectivity for epoch 25
    computing connectivity for epoch 26
    computing connectivity for epoch 27
    computing connectivity for epoch 28
    computing connectivity for epoch 29
    computing connectivity for epoch 30
    computing connectivity for epoch 31
    computing connectivity for epoch 32
    computing connectivity for epoch 33
    computing connectivity for epoch 34
    computing connectivity for epoch 35
    computing connectivity for epoch 36
    computing connectivity for epoch 37
    computing connectivity for epoch 38
    Rejecting  epoch based on EOG : [u'EOG 061']
    computing connectivity for epoch 39
    computing connectivity for epoch 40
    computing connectivity for epoch 41
    computing connectivity for epoch 42
    computing connectivity for epoch 43
    computing connectivity for epoch 44
    computing connectivity for epoch 45
    computing connectivity for epoch 46
    computing connectivity for epoch 47
    computing connectivity for epoch 48
    Rejecting  epoch based on EOG : [u'EOG 061']
    computing connectivity for epoch 49
    computing connectivity for epoch 50
    computing connectivity for epoch 51
    computing connectivity for epoch 52
    computing connectivity for epoch 53
    computing connectivity for epoch 54
    computing connectivity for epoch 55
    computing connectivity for epoch 56
    computing connectivity for epoch 57
    computing connectivity for epoch 58
    computing connectivity for epoch 59
    Rejecting  epoch based on EOG : [u'EOG 061']
    computing connectivity for epoch 60
    computing connectivity for epoch 61
    computing connectivity for epoch 62
    computing connectivity for epoch 63
    computing connectivity for epoch 64
    computing connectivity for epoch 65
    computing connectivity for epoch 66
    computing connectivity for epoch 67
    assembling connectivity matrix
[Connectivity computation done]
(204L, 204L, 1L)

Mayaviを使います。Jupyterでは具合が悪いのでpython consoleを使います。

In [14]:
print(type(con))
print(con.shape) 
filename=data_path+'/MEG/sample\\aaa.npz'# .npzはndarrayの拡張子 con.npyはダメみたい
np.savez(filename,con=con,idx=idx)
<type 'numpy.ndarray'>
(203L, 203L)

python consoleでは以下の図となります。

In [ ]:
# SpyderのConsoleで実行
import numpy as np
from scipy import linalg
import mne
from mne import io
from mne.datasets import sample
data_path=mne.datasets.sample.data_path()
filename=data_path+'/MEG/sample/aaa.npz'
x=np.load(filename)
con,idx=x['con'],x['idx']
raw_fname=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw.fif'
raw=io.read_raw_fif(raw_fname)
picks=mne.pick_types(raw.info,meg='grad',eeg=False,stim=False,eog=True,exclude='bads')
from mayavi import mlab
mlab.figure(size=(600,600),bgcolor=(0.5,0.5,0.5));
sens_loc=[raw.info['chs'][picks[i]]['loc'][:3] for i in idx]
sens_loc=np.array(sens_loc)
pts=mlab.points3d(
    sens_loc[:,0],sens_loc[:,1],sens_loc[:,2],color=(1,1,1),opacity=1,scale_factor=0.005);
n_con=20  # 最大20のconnection # connectionの最大値
min_dist=0.05  # 5cm以内のセンサは排除
threshold=np.sort(con,axis=None)[-n_con]
ii,jj=np.where(con>=threshold)
con_nodes=list()
con_val=list()
for i, j in zip(ii,jj):
    if linalg.norm(sens_loc[i]-sens_loc[j])>min_dist:
        con_nodes.append((i,j))
        con_val.append(con[i,j])
con_val=np.array(con_val)
vmax=np.max(con_val)# 描画
vmin=np.min(con_val)
for val,nodes in zip(con_val,con_nodes):
    x1,y1,z1=sens_loc[nodes[0]]
    x2,y2,z2=sens_loc[nodes[1]]
    points=mlab.plot3d(
        [x1,x2],[y1,y2],[z1,z2],[val,val],vmin=vmin,vmax=vmax,
        tube_radius=0.001,colormap='RdBu');
    points.module_manager.scalar_lut_manager.reverse_lut = True
mlab.scalarbar(title='Phase Lag Index (PLI)',nb_labels=4)
nodes_shown=list(set([n[0] for n in con_nodes]+[n[1] for n in con_nodes]))
for node in nodes_shown:
    x,y,z=sens_loc[node]
    mlab.text3d(x,y,z,raw.ch_names[picks[node]],scale=0.005,color=(0,0,0));
view=(-88.7,40.8,0.76, np.array([-3.9e-4,-8.5e-3,-1e-2]))
mlab.view(*view)

視覚刺激時の信号源空間のPhase Slope Index (PSI)の計算

In [15]:
% reset
Once deleted, variables cannot be recovered. Proceed (y/[n])? y
In [16]:
# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>
#
# License: BSD (3-clause)

IPython環境になるようで、ちょっと時間がかかります。

In [17]:
import numpy as np
import mne
from mne.datasets import sample
from mne.minimum_norm import read_inverse_operator, apply_inverse_epochs
from mne.connectivity import seed_target_indices, phase_slope_index
print(__doc__)
Built-in functions, exceptions, and other objects.

Noteworthy: None is the `nil' object; Ellipsis represents `...' in slices.
In [18]:
data_path=sample.data_path()
subjects_dir=data_path+'/subjects'
fname_inv=data_path+'/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif' # 計算済み
fname_raw=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw.fif'
fname_event=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
fname_label=data_path+'/MEG/sample/labels/Vis-lh.label'
event_id,tmin,tmax= 4,-0.2,0.3
method="dSPM"  # use dSPM method (could also be MNE or sLORETA)
# データ読み込み
inverse_operator=read_inverse_operator(fname_inv)
print(type(inverse_operator)) # クラスです
print(inverse_operator.keys()) # いろいろなクラスの集まりです
raw=mne.io.read_raw_fif(fname_raw)
events=mne.read_events(fname_event)
picks=mne.pick_types(raw.info,meg=True,eeg=False,stim=False,eog=True,exclude='bads')
# エポッククラスの作成
epochs=mne.Epochs(
    raw,events,event_id,tmin, tmax,picks=picks,baseline=(None, 0),
    reject=dict(mag=4e-12,grad=4000e-13,eog=150e-6))
# 空間フィルタの作成
snr=1.0 
lambda2=1.0/snr**2
stcs=apply_inverse_epochs( # generatorだそうです return_generator=Falseは数が少ないから?
    epochs,inverse_operator,lambda2,method,pick_ori="normal",return_generator=False)
# seed time seriesの計算
label=mne.read_label(fname_label) #Vis-lh.label
src=inverse_operator['src']  # 信号源空間を利用
seed_ts=mne.extract_label_time_course(stcs,label,src,mode='mean_flip')
# seed time courseと信号源推定
# index 0: Vis-lh.labelの経時変化
# index 1..7499: 信号源空間の経時変化
comb_ts=zip(seed_ts,stcs)
# connectivityの計算
vertices=[src[i]['vertno'] for i in range(2)]
n_signals_tot=1+len(vertices[0])+len(vertices[1])
indices=seed_target_indices([0],np.arange(1,n_signals_tot))
# 8~30HzでPhase Slow Index baseline期間は含まず
fmin,fmax,tmin_con=8.,30.,0.
sfreq=raw.info['sfreq'] 
psi,freqs,times,n_epochs,_=phase_slope_index(
    comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
psi_stc=mne.SourceEstimate(psi,vertices=vertices,tmin=0,tstep=1,subject='sample')
Reading inverse operator decomposition from C:\Users\akira\mne_data\MNE-sample-data/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif...
    Reading inverse operator info...
    [done]
    Reading inverse operator decomposition...
    [done]
    305 x 305 full covariance (kind = 1) found.
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Noise covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 2) found.
    Source covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 6) found.
    Orientation priors read.
    22494 x 22494 diagonal covariance (kind = 5) found.
    Depth priors read.
    Did not find the desired covariance matrix (kind = 3)
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    2 source spaces read
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Source spaces transformed to the inverse solution coordinate frame
<class 'mne.minimum_norm.inverse.InverseOperator'>
['nave', 'methods', 'noisenorm', 'eigen_fields', 'units', 'fmri_prior', 'noise_cov', 'orient_prior', 'mri_head_t', 'sing', 'nsource', 'info', 'src', 'nchan', 'source_cov', 'projs', 'source_nn', 'depth_prior', 'whitener', 'source_ori', 'coord_frame', 'proj', 'eigen_leads_weighted', 'eigen_leads', 'reginv']
Opening raw data file C:\Users\akira\mne_data\MNE-sample-data\MEG\sample\sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Current compensation grade : 0
70 matching events found
Created an SSP operator (subspace dimension = 3)
4 projection items activated
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 1
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 3)
    Created the whitener using a full noise covariance matrix (3 small eigenvalues omitted)
    Computing noise-normalization factors (dSPM)...
[done]
Picked 305 channels from the data
Computing inverse...
(eigenleads need to be weighted)...
Processing epoch : 1
Processing epoch : 2
Processing epoch : 3
Processing epoch : 4
Processing epoch : 5
Processing epoch : 6
Processing epoch : 7
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 8
Processing epoch : 9
Processing epoch : 10
Processing epoch : 11
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 12
Processing epoch : 13
Processing epoch : 14
Processing epoch : 15
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 16
Processing epoch : 17
Processing epoch : 18
Processing epoch : 19
Processing epoch : 20
Processing epoch : 21
Processing epoch : 22
Processing epoch : 23
Processing epoch : 24
Processing epoch : 25
Processing epoch : 26
Processing epoch : 27
Processing epoch : 28
Processing epoch : 29
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 30
Processing epoch : 31
Processing epoch : 32
Processing epoch : 33
Processing epoch : 34
Processing epoch : 35
Processing epoch : 36
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 37
Processing epoch : 38
Processing epoch : 39
Processing epoch : 40
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 41
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 42
Processing epoch : 43
Processing epoch : 44
Processing epoch : 45
Processing epoch : 46
Processing epoch : 47
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 48
Processing epoch : 49
Processing epoch : 50
Processing epoch : 51
Processing epoch : 52
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 53
Processing epoch : 54
Processing epoch : 55
Processing epoch : 56
Processing epoch : 57
Processing epoch : 58
Processing epoch : 59
Processing epoch : 60
[done]
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Extracting time courses for 1 labels (mode: mean_flip)
Estimating phase slope index (PSI)
Connectivity computation...
    computing connectivity for 7498 connections
    using t=0.000s..0.300s for estimation (46 points)
fmin corresponds to less than 5 cycles, spectrum estimate will be unreliable
    frequencies: 9.8Hz..29.4Hz (7 points)
    using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: Coherency
    computing connectivity for epoch 1
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: fmin corresponds to less than 5 cycles, spectrum estimate will be unreliable
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 2
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 3
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 4
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 5
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 6
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 7
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 8
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 9
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 10
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 11
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 12
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 13
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 14
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 15
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 16
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 17
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 18
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 19
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 20
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 21
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 22
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 23
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 24
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 25
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 26
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 27
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 28
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 29
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 30
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 31
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 32
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 33
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 34
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 35
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 36
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 37
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 38
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 39
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 40
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 41
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 42
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 43
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 44
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 45
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 46
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 47
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 48
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 49
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 50
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 51
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 52
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 53
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 54
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 55
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 56
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 57
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 58
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
Performance can be improved by not accessing the data attribute before calling this method.
    computing connectivity for epoch 59
Performance can be improved by not accessing the data attribute before calling this method.
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
    computing connectivity for epoch 60
Performance can be improved by not accessing the data attribute before calling this method.
[Connectivity computation done]
Computing PSI from estimated Coherency
[PSI Estimation Done]
<ipython-input-18-30819b23e864>:41: RuntimeWarning: Performance can be improved by not accessing the data attribute before calling this method.
  comb_ts,mode='multitaper',indices=indices,sfreq=sfreq,fmin=fmin,fmax=fmax,tmin=tmin_con)
In [19]:
# データ保存
fname_psi=fname_raw=data_path+'/MEG/sample/psi.npy'
fname_psi_stc=fname_raw=data_path+'/MEG/sample/psi.stc'
np.save(fname_psi,psi)
psi_stc.save(fname_psi_stc)
Writing STC to disk...
[done]
In [20]:
# python console側で読みだすファイル名
print(fname_psi)
print(fname_psi_stc)
print(subjects_dir)
print(fname_label)
C:\Users\akira\mne_data\MNE-sample-data/MEG/sample/psi.npy
C:\Users\akira\mne_data\MNE-sample-data/MEG/sample/psi.stc
C:\Users\akira\mne_data\MNE-sample-data/subjects
C:\Users\akira\mne_data\MNE-sample-data/MEG/sample/labels/Vis-lh.label

3DなのでPython Consoleを使います。
以下のような図が描画されます。
python consoleでは以下の図となります。

In [ ]:
# SpyderのPython Consoleで実行
import numpy as np
import mne
data_path=mne.datasets.sample.data_path()
psi=np.load(data_path+'/MEG/sample/psi.npy')
psi_stc=mne.read_source_estimate(data_path+'/MEG/sample/psi.stc')
psi_stc.subject='sample' # これは何故か保存されてないみたい
subjects_dir=data_path+'/subjects'
v_max=np.max(np.abs(psi))
brain=psi_stc.plot(
    surface='inflated',hemi='lh',time_label='Phase Slope Index (PSI)',
    subjects_dir=subjects_dir,clim=dict(kind='percent',pos_lims=(95, 97.5, 100)))
brain.show_view('medial')
fname_label=data_path+'/MEG/sample/labels/Vis-lh.label'
brain.add_label(fname_label,color='green',alpha=0.7)

信号源空間でのcoherence

In [21]:
% reset
Once deleted, variables cannot be recovered. Proceed (y/[n])? y
In [22]:
# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>
#
# License: BSD (3-clause)
In [23]:
import numpy as np
import mne
from mne.datasets import sample
from mne.minimum_norm import (apply_inverse, apply_inverse_epochs,read_inverse_operator)
from mne.connectivity import seed_target_indices, spectral_connectivity
print(__doc__)
data_path=sample.data_path()
subjects_dir=data_path + '/subjects'
fname_inv=data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
fname_raw=data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
fname_event=data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
label_name_lh='Aud-lh'
fname_label_lh=data_path + '/MEG/sample/labels/%s.label' % label_name_lh
event_id,tmin,tmax=1,-0.2,0.5
method="dSPM" 
Built-in functions, exceptions, and other objects.

Noteworthy: None is the `nil' object; Ellipsis represents `...' in slices.
In [24]:
inverse_operator=read_inverse_operator(fname_inv)
label_lh=mne.read_label(fname_label_lh) # 左聴覚野をseedにする
raw=mne.io.read_raw_fif(fname_raw)
events=mne.read_events(fname_event)
raw.info['bads']+=['MEG 2443']
picks=mne.pick_types(raw.info,meg=True,eeg=False,stim=False,eog=True,exclude='bads')
epochs = mne.Epochs(
    raw,events,event_id,tmin,tmax,picks=picks,baseline=(None,0),
    reject=dict(mag=4e-12,grad=4000e-13,eog=150e-6))
snr=3.0
lambda2=1.0/snr**2
evoked=epochs.average()
stc=apply_inverse(evoked,inverse_operator,lambda2,method,pick_ori="normal")
stc_label=stc.in_label(label_lh) # 左半球に限定
src_pow=np.sum(stc_label.data**2,axis=1)
seed_vertno=stc_label.vertices[0][np.argmax(src_pow)]
seed_idx=np.searchsorted(stc.vertices[0],seed_vertno)
n_sources=stc.data.shape[0]
indices=seed_target_indices([seed_idx],np.arange(n_sources))

# return_generator=Trueにすると信号源データなしで計算できてメモリ節約できる
snr=1.0  
lambda2=1.0/snr**2
stcs=apply_inverse_epochs(
    epochs,inverse_operator,lambda2,method,pick_ori="normal",return_generator=True)
print(stcs) # generatorだそうです
fmin=(8.,13.) # Hz
fmax=(13.,30.) # Hz
sfreq=raw.info['sfreq']
coh,freqs,times,n_epochs,n_tapers=spectral_connectivity(
    stcs,method='coh',mode='fourier',indices=indices,sfreq=sfreq,fmin=fmin,
    fmax=fmax,faverage=True,n_jobs=1)
print('Frequencies in Hz over which coherence was averaged for alpha: ')
print(freqs[0])
print('Frequencies in Hz over which coherence was averaged for beta: ')
print(freqs[1])

# seed1個で信号源のconnectivityの計算
tmin=np.mean(freqs[0])
tstep=np.mean(freqs[1])-tmin
coh_stc= mne.SourceEstimate(
    coh,vertices=stc.vertices,tmin=1e-3*tmin,tstep=1e-3*tstep,subject='sample')
Reading inverse operator decomposition from C:\Users\akira\mne_data\MNE-sample-data/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif...
    Reading inverse operator info...
    [done]
    Reading inverse operator decomposition...
    [done]
    305 x 305 full covariance (kind = 1) found.
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Noise covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 2) found.
    Source covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 6) found.
    Orientation priors read.
    22494 x 22494 diagonal covariance (kind = 5) found.
    Depth priors read.
    Did not find the desired covariance matrix (kind = 3)
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    2 source spaces read
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Source spaces transformed to the inverse solution coordinate frame
Opening raw data file C:\Users\akira\mne_data\MNE-sample-data\MEG\sample\sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Current compensation grade : 0
72 matching events found
Created an SSP operator (subspace dimension = 3)
4 projection items activated
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on MAG : [u'MEG 1711']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 55
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 3)
    Created the whitener using a full noise covariance matrix (3 small eigenvalues omitted)
    Computing noise-normalization factors (dSPM)...
[done]
Picked 305 channels from the data
Computing inverse...
(eigenleads need to be weighted)...
(dSPM)...
[done]
<generator object _apply_inverse_epochs_gen at 0x000000000D6F9708>
Connectivity computation...
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 1
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 3)
    Created the whitener using a full noise covariance matrix (3 small eigenvalues omitted)
    Computing noise-normalization factors (dSPM)...
[done]
Picked 305 channels from the data
Computing inverse...
(eigenleads need to be weighted)...
Processing epoch : 1
    computing connectivity for 7498 connections
    using t=-0.200s..0.499s for estimation (106 points)
    computing connectivity for the bands:
     band 1: 8.5Hz..12.7Hz (4 points)
     band 2: 14.2Hz..29.7Hz (12 points)
    connectivity scores will be averaged for each band
    using FFT with a Hanning window to estimate spectra
    the following metrics will be computed: Coherence
    computing connectivity for epoch 1
Processing epoch : 2
    computing connectivity for epoch 2
Processing epoch : 3
    computing connectivity for epoch 3
Processing epoch : 4
    computing connectivity for epoch 4
Processing epoch : 5
    computing connectivity for epoch 5
Processing epoch : 6
    computing connectivity for epoch 6
Processing epoch : 7
    computing connectivity for epoch 7
Processing epoch : 8
    computing connectivity for epoch 8
Processing epoch : 9
    computing connectivity for epoch 9
Processing epoch : 10
    computing connectivity for epoch 10
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 11
    computing connectivity for epoch 11
Processing epoch : 12
    computing connectivity for epoch 12
Processing epoch : 13
    computing connectivity for epoch 13
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 14
    computing connectivity for epoch 14
Processing epoch : 15
    computing connectivity for epoch 15
Processing epoch : 16
    computing connectivity for epoch 16
Processing epoch : 17
    computing connectivity for epoch 17
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 18
    computing connectivity for epoch 18
Processing epoch : 19
    computing connectivity for epoch 19
Processing epoch : 20
    computing connectivity for epoch 20
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 21
    computing connectivity for epoch 21
Processing epoch : 22
    computing connectivity for epoch 22
Processing epoch : 23
    computing connectivity for epoch 23
    Rejecting  epoch based on MAG : [u'MEG 1711']
Processing epoch : 24
    computing connectivity for epoch 24
Processing epoch : 25
    computing connectivity for epoch 25
Processing epoch : 26
    computing connectivity for epoch 26
Processing epoch : 27
    computing connectivity for epoch 27
Processing epoch : 28
    computing connectivity for epoch 28
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 29
    computing connectivity for epoch 29
Processing epoch : 30
    computing connectivity for epoch 30
Processing epoch : 31
    computing connectivity for epoch 31
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 32
    computing connectivity for epoch 32
Processing epoch : 33
    computing connectivity for epoch 33
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 34
    computing connectivity for epoch 34
Processing epoch : 35
    computing connectivity for epoch 35
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 36
    computing connectivity for epoch 36
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 37
    computing connectivity for epoch 37
Processing epoch : 38
    computing connectivity for epoch 38
Processing epoch : 39
    computing connectivity for epoch 39
Processing epoch : 40
    computing connectivity for epoch 40
Processing epoch : 41
    computing connectivity for epoch 41
Processing epoch : 42
    computing connectivity for epoch 42
Processing epoch : 43
    computing connectivity for epoch 43
Processing epoch : 44
    computing connectivity for epoch 44
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 45
    computing connectivity for epoch 45
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 46
    computing connectivity for epoch 46
Processing epoch : 47
    computing connectivity for epoch 47
Processing epoch : 48
    computing connectivity for epoch 48
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 49
    computing connectivity for epoch 49
Processing epoch : 50
    computing connectivity for epoch 50
Processing epoch : 51
    computing connectivity for epoch 51
Processing epoch : 52
    computing connectivity for epoch 52
Processing epoch : 53
    computing connectivity for epoch 53
Processing epoch : 54
    computing connectivity for epoch 54
Processing epoch : 55
    computing connectivity for epoch 55
[done]
[Connectivity computation done]
Frequencies in Hz over which coherence was averaged for alpha: 
[  8.49926873   9.91581352  11.33235831  12.74890309]
Frequencies in Hz over which coherence was averaged for beta: 
[ 14.16544788  15.58199267  16.99853746  18.41508225  19.83162704
  21.24817182  22.66471661  24.0812614   25.49780619  26.91435098
  28.33089577  29.74744055]
In [25]:
stcs
Out[25]:
<generator object _apply_inverse_epochs_gen at 0x000000000D6F9708>
In [26]:
# データ保存
fname_coh_stc=fname_raw=data_path+'/MEG/sample/coh.stc'
coh_stc.save(fname_coh_stc)
Writing STC to disk...
[done]

SpyderのConsoleで実行すると以下の図となります。

In [ ]:
# SpyderのPython Consoleで実行
import mne
data_path=mne.datasets.sample.data_path()
coh_stc=mne.read_source_estimate(data_path+'/MEG/sample/coh.stc')
coh_stc.subject='sample' # これは何故か保存されてないみたい
subjects_dir=data_path+'/subjects'
brain=coh_stc.plot('sample','inflated','both',time_label='Coherence %0.1f Hz',subjects_dir=subjects_dir,
                   clim=dict(kind='value',lims=(0.25, 0.4, 0.65)))
brain.show_view('lateral')

信号源空間の結合性の円形表示

In [27]:
% reset
Once deleted, variables cannot be recovered. Proceed (y/[n])? y
In [28]:
# Authors: Martin Luessi <mluessi@nmr.mgh.harvard.edu>
#          Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
#          Nicolas P. Rougier (graph code borrowed from his matplotlib gallery)
#
# License: BSD (3-clause)
In [29]:
import mne
from mne.datasets import sample
from mne.minimum_norm import apply_inverse_epochs,read_inverse_operator
from mne.connectivity import spectral_connectivity
from mne.viz import circular_layout,plot_connectivity_circle
print(__doc__)

data_path=sample.data_path()
subjects_dir=data_path+'/subjects'
fname_inv=data_path+'/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
fname_raw=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw.fif'
fname_event=data_path+'/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'

inverse_operator=read_inverse_operator(fname_inv) # 計算済み
print(type(inverse_operator)) # クラスです
print(inverse_operator.keys()) # いろいろなクラスの集まりです
raw=mne.io.read_raw_fif(fname_raw)
events= mne.read_events(fname_event) # 処理済み
raw.info['bads']+=['MEG 2443']
picks=mne.pick_types(raw.info,meg=True,eeg=False,stim=False,eog=True,exclude='bads')
# クラスepochの作成
event_id,tmin,tmax=1,-0.2,0.5
epochs=mne.Epochs(
    raw,events,event_id,tmin,tmax,picks=picks,baseline=(None,0),
    reject=dict(mag=4e-12,grad=4000e-13,eog=150e-6))
# 信号源空間dSPMの計算 MNE, sLORETAも可
snr=1.0
lambda2=1.0/snr**2
stcs=apply_inverse_epochs(
    epochs,inverse_operator,lambda2,'dSPM',pick_ori="normal",return_generator=True)
print(stcs) # generatorだそうです
# FreeSurferのLabel取得
labels=mne.read_labels_from_annot('sample',parc='aparc',subjects_dir=subjects_dir)
label_colors=[label.color for label in labels]
print(len(labels)) # 68個のラベルです。
# return_generator=Trueが定番です。
src=inverse_operator['src'] #SourceSpacesクラス 要するにメッシュ
label_ts=mne.extract_label_time_course(
    stcs,labels,src,mode='mean_flip',return_generator=True)
print(label_ts) # generatorだそうです。
Built-in functions, exceptions, and other objects.

Noteworthy: None is the `nil' object; Ellipsis represents `...' in slices.
Reading inverse operator decomposition from C:\Users\akira\mne_data\MNE-sample-data/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif...
    Reading inverse operator info...
    [done]
    Reading inverse operator decomposition...
    [done]
    305 x 305 full covariance (kind = 1) found.
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Noise covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 2) found.
    Source covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 6) found.
    Orientation priors read.
    22494 x 22494 diagonal covariance (kind = 5) found.
    Depth priors read.
    Did not find the desired covariance matrix (kind = 3)
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    2 source spaces read
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Source spaces transformed to the inverse solution coordinate frame
<class 'mne.minimum_norm.inverse.InverseOperator'>
['nave', 'methods', 'noisenorm', 'eigen_fields', 'units', 'fmri_prior', 'noise_cov', 'orient_prior', 'mri_head_t', 'sing', 'nsource', 'info', 'src', 'nchan', 'source_cov', 'projs', 'source_nn', 'depth_prior', 'whitener', 'source_ori', 'coord_frame', 'proj', 'eigen_leads_weighted', 'eigen_leads', 'reginv']
Opening raw data file C:\Users\akira\mne_data\MNE-sample-data\MEG\sample\sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Current compensation grade : 0
72 matching events found
Created an SSP operator (subspace dimension = 3)
4 projection items activated
<generator object _apply_inverse_epochs_gen at 0x000000000D6F91B0>
Reading labels from parcellation..
   read 34 labels from C:\Users\akira\mne_data\MNE-sample-data/subjects\sample\label\lh.aparc.annot
   read 34 labels from C:\Users\akira\mne_data\MNE-sample-data/subjects\sample\label\rh.aparc.annot
[done]
68
<generator object _gen_extract_label_time_course at 0x000000000BB41BD0>

α帯域の結合性の計算

1) rawからepochを1つ読み込む
2) SSP適応と基線補正
3) 信号源波形の計算
4) 各部位(label)毎に信号源波形の加算平
5) 結合性の計算し、1)次のepoch
epoch毎に計算することでメモリ消費を抑える

In [30]:
# 結合性の計算
import numpy as np
import matplotlib.pyplot as plt
# % matplotlib inline # いらないみたい
fmin,fmax=8.,13.
sfreq=raw.info['sfreq'] 
con_methods=['pli','wpli2_debiased']
con,freqs,times,n_epochs,n_tapers=spectral_connectivity( # ここで改行すると4文字タブとなる
    label_ts,method=con_methods,mode='multitaper',sfreq=sfreq,fmin=fmin,
    fmax=fmax,faverage=True,mt_adaptive=True,n_jobs=1)
# conは三次元配列のリスト
print(con[0].shape)
print(con[1].shape)
con_res= dict()
for method, c in zip(con_methods,con):
    con_res[method]=c[:,:,0]
# 描画
label_names=[label.name for label in labels]
lh_labels=[name for name in label_names if name.endswith('lh')] # 左半球
label_ypos=list()
for name in lh_labels:
    idx=label_names.index(name)
    ypos=np.mean(labels[idx].pos[:,1])
    label_ypos.append(ypos)
lh_labels=[label for (yp,label) in sorted(zip(label_ypos, lh_labels))]
rh_labels=[label[:-2]+'rh' for label in lh_labels] # 右半球

node_order=list()
node_order.extend(lh_labels[::-1])  # reverse the order
node_order.extend(rh_labels)
node_angles=circular_layout(
    label_names,node_order,start_pos=90,group_boundaries=[0,len(label_names)/2])
# 上位300のみ表示
plot_connectivity_circle(
    con_res['pli'],label_names,n_lines=300,node_angles=node_angles,node_colors=label_colors,
    title='All-to-All Connectivity left-Auditory ' 'Condition (PLI)') # 'A' 'B' = 'AB'
plt.savefig('circle.png',facecolor='black')
# なんかうまく表示されません
fig=plt.figure(num=None,figsize=(8,4),facecolor='black')
no_names=['']*len(label_names)
for ii, method in enumerate(con_methods):
    plot_connectivity_circle(
        con_res[method],no_names,n_lines=300,node_angles=node_angles,node_colors=label_colors,
        title=method,padding=0,fontsize_colorbar=6,fig=fig,subplot=(1,2,ii+1))
plt.show()
Connectivity computation...
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 1
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 3)
    Created the whitener using a full noise covariance matrix (3 small eigenvalues omitted)
    Computing noise-normalization factors (dSPM)...
[done]
Picked 305 channels from the data
Computing inverse...
(eigenleads need to be weighted)...
Processing epoch : 1
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for 2278 connections
    using t=0.000s..0.699s for estimation (106 points)
    frequencies: 8.5Hz..12.7Hz (4 points)
    connectivity scores will be averaged for each band
    using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: PLI, Debiased WPLI Square
    computing connectivity for epoch 1
Processing epoch : 2
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 2
Processing epoch : 3
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 3
Processing epoch : 4
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 4
Processing epoch : 5
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 5
Processing epoch : 6
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 6
Processing epoch : 7
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 7
Processing epoch : 8
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 8
Processing epoch : 9
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 9
Processing epoch : 10
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 10
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 11
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 11
Processing epoch : 12
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 12
Processing epoch : 13
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 13
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 14
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 14
Processing epoch : 15
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 15
Processing epoch : 16
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 16
Processing epoch : 17
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 17
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 18
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 18
Processing epoch : 19
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 19
Processing epoch : 20
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 20
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 21
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 21
Processing epoch : 22
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 22
Processing epoch : 23
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 23
    Rejecting  epoch based on MAG : [u'MEG 1711']
Processing epoch : 24
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 24
Processing epoch : 25
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 25
Processing epoch : 26
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 26
Processing epoch : 27
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 27
Processing epoch : 28
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 28
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 29
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 29
Processing epoch : 30
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 30
Processing epoch : 31
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 31
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 32
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 32
Processing epoch : 33
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 33
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 34
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 34
Processing epoch : 35
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 35
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 36
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 36
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 37
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 37
Processing epoch : 38
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 38
Processing epoch : 39
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 39
Processing epoch : 40
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 40
Processing epoch : 41
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 41
Processing epoch : 42
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 42
Processing epoch : 43
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 43
Processing epoch : 44
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 44
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 45
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 45
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 46
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 46
Processing epoch : 47
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 47
Processing epoch : 48
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 48
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
Processing epoch : 49
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 49
Processing epoch : 50
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 50
Processing epoch : 51
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 51
Processing epoch : 52
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 52
Processing epoch : 53
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 53
Processing epoch : 54
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 54
Processing epoch : 55
Extracting time courses for 68 labels (mode: mean_flip)
    computing connectivity for epoch 55
[done]
    assembling connectivity matrix
[Connectivity computation done]
(68L, 68L, 1L)
(68L, 68L, 1L)