Hi, all!
- MNE-Python version: 0.23.4
- operating system: Win10
I have a problem with plotting the results of spatio-temporal clustering.
I think I could run the function āmne.stats.spatio_temporal_cluster_testā without wny error.
However, the result of āmne.stats.summarize_clusters_stc.plotā does not seem to be working.
The main problems are as follows
- the time window of my data is 0.1-0.6(s), but the time window of the seek bar in the plot is -0.1-9.9(s). why is there such a difference?
- There is no color signals which represent the significant cluster, instead, the brain image appears completely white.
The detail of them is like this.
What I am doing right now is as follows:
import os.path as op
import os
import glob
import numpy as np
from scipy import stats as stats
import mne
from mne import spatial_src_adjacency
from mne.stats import spatio_temporal_cluster_test, summarize_clusters_stc
from mne.cov import compute_covariance
from mne.minimum_norm import make_inverse_operator, apply_inverse
from mne import make_forward_solution
import collections
import re
print(__doc__)
# setting
snr = 1.0 # use smaller SNR for raw data
method = 'dSPM'
parc = 'aparc' # the parcellation to use, e.g., 'aparc' 'aparc.a2009s'
lambda2 = 1.0 / snr ** 2
n_subjects1 = 30
n_subjects2 = 60
tmin=-0.1
tmax=0.6
# read eeg
datas_path = "/media/owner/HDPH-UT/"
subjects_dir = "/home/owner/fs-output"
#read EEG data
EEG_path1 = os.path.join(datas_path,'final_analysis', 'PAS', '30-HC','single3')
EEG_path2 = os.path.join(datas_path,'final_analysis', 'PAS', '31-MDD1','single3')
for gro_num in range(2):
if gro_num == 0:
EEG_path = EEG_path1
elif gro_num == 1:
EEG_path = EEG_path2
EEG_file = glob.glob(os.path.join(EEG_path1, '*.set'))
for sub_num in range(len(EEG_file)):
epochs = mne.read_epochs_eeglab(EEG_file[sub_num])
epochs2 = epochs.copy().rename_channels({'FP1':'Fp1', 'FPZ':'Fpz', 'FP2':'Fp2', 'FZ':'Fz', 'FCZ':'FCz','CZ':'Cz', 'CPZ':'CPz', 'PZ':'Pz', 'POZ':'POz','CB1':'O1w', 'O1':'Ozw','OZ': 'O2w','O2': 'BAD1','CB2': 'BAD2'})
epochs3 = epochs2.copy().rename_channels({'O1w':'O1', 'Ozw':'Oz','O2w': 'O2'})
epochs4 = epochs3.drop_channels(['BAD1', 'BAD2'])
od = collections.OrderedDict()
for chan_name in epochs4.ch_names:
od[chan_name] = epochs4.get_montage()._get_ch_pos()[chan_name]/1000
posit = mne.channels.make_dig_montage(od)
epochs5 = epochs4.copy().set_montage(posit)
epochs_1020 = epochs5.copy().crop(tmin=tmin, tmax=tmax).set_montage('standard_1020')
# read MRI
subjects_dir_list2 = os.listdir(path=subjects_dir)
subjects_dir_list3 = [s for s in subjects_dir_list2 if re.match('.*T1w_MPR.*nii', s)]
subjects_dir_list4 = [s for s in subjects_dir_list2 if re.match('.*navigation.*nii', s)]
subjects_dir_list = subjects_dir_list3 + subjects_dir_list4
EEG_id = re.search(r'\d\d\d\d\d+', EEG_file[sub_num])
for mri_num in range(len(subjects_dir_list)):
mri_id = re.search(r'\d\d\d\d\d+', subjects_dir_list[mri_num])
if EEG_id[0] == mri_id[0]:
subject = subjects_dir_list[mri_num]
break
trans_fname = op.join(subjects_dir, subject, '%s-trans.fif' % subject)
trans = mne.read_trans(trans_fname)
bem_dir = os.path.join(subjects_dir,subject, 'bem')
src = mne.read_source_spaces(os.path.join(subjects_dir,subject, 'bem','%s-oct6-src.fif' % subject))
fname_model = op.join(bem_dir, '%s-5120-5120-5120-bem.fif' % subject)
fname_bem = op.join(bem_dir, '%s-5120-5120-5120-bem-sol.fif' % subject)
# create noise_covāforward_modelāinverse_operatorāapply_inverse_epochsć§SourceEstimation(stc)ćä½ę
noise_cov = compute_covariance(epochs_1020, tmax=0. ,method='auto')
fwd = make_forward_solution(epochs_1020.info, trans=trans, src=src, bem=fname_bem,
mindist=5.0)
inverse_operator = make_inverse_operator(epochs_1020.info, fwd, noise_cov)
evoked = epochs_1020.copy().set_eeg_reference(projection=True).average()
stc = apply_inverse(evoked, inverse_operator=inverse_operator,lambda2=lambda2, method=method)
stc.resample(200, npad='auto')
# compute_source_morph
#fsave_vertices = [s['vertno'] for s in src]
morph = mne.compute_source_morph(stc, subject, 'fsaverage',
spacing=5, subjects_dir=subjects_dir)
stc = morph.apply(stc)
n_vertices_fsave, n_times = stc.data.shape
tstep = stc.tstep * 1000 # convert to milliseconds
if sub_num == 0:
if gro_num == 0:
np.random.seed(0)
X1 = np.random.randn(n_vertices_fsave, n_times, n_subjects1) * 10
elif gro_num == 1:
np.random.seed(0)
X2 = np.random.randn(n_vertices_fsave, n_times, n_subjects2) * 10
if gro_num == 0:
X1[:, :, sub_num] = stc.data[:, :]
if gro_num == 1:
X2[:, :, sub_num] = stc.data[:, :]
del(src, stc, fwd, inverse_operator, morph, noise_cov)
# transpose
X1 = np.abs(X1) # only magnitude
X2 = np.abs(X2) # only magnitude
X1 = np.transpose(X1, [2, 1, 0])
X2 = np.transpose(X2, [2, 1, 0])
X = [X1, X2]
del(X1, X2)
# statistics
fsave_bem_dir = "/home/owner/mne_data/MNE-sample-data/subjects/fsaverage/bem"
print('Computing adjacency.')
src = mne.read_source_spaces(os.path.join(fsave_bem_dir,'fsaverage-ico-5-src.fif'))
adjacency = spatial_src_adjacency(src)
p_threshold = 0.05
f_threshold = stats.distributions.f.ppf(1. - p_threshold / 2.,
n_subjects1 - 1, n_subjects2 - 1)
print('Clustering.')
T_obs, clusters, cluster_p_values, H0 = clu =\
spatio_temporal_cluster_test(X, adjacency=adjacency, n_jobs=1,
threshold=f_threshold, buffer_size=None)
good_cluster_inds = np.where(cluster_p_values < 0.05)[0]
# plot
fsave_vertices = [np.arange(n_vertices_fsave/2), np.arange(n_vertices_fsave/2)]
stc_all_cluster_vis = summarize_clusters_stc(clu, tstep=tstep, tmin=-0.1,
vertices=fsave_vertices,
subject='fsaverage')
brain = stc_all_cluster_vis.plot('fsaverage', hemi='both',
views='lateral', subjects_dir=subjects_dir, spacing='ico5',
time_label='temporal extent (ms)',
clim=dict(kind='value', lims=[0, 1, 40]))
Any help would be appreciated.
Thanks.
Masataka.