Unable to display the result of spatio-temporal clustering

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

  1. 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?
  2. 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.

Hi all!

I am sorry to bother you, but I would be grateful if you could let me have your answer concerning this matter.
If I should share the variable “clu”, the result of “spatio_temporal_cluster_test()”, with you, I will do so immediately.
If there is anything else you need, please tell me.

Once again, I apologize for any inconvenience this may have caused you, but I am looking forward to your reply.

Best Regard,
Masataka

Hello @masataka-wada,

sorry your question went unanswered for a while. I believe one reason is that the code snippet you presented is rather complex and it’s difficult for readers to understand what exactly you’re doing and when things start to go wrong. A proper minimal working example would probably simplify things…

That said, I just realized that you pass lims=[0, 1, 40]. I’m assuming these limits don’t fit well with your data. Could you simply remove this parameter altogether, such that MNE would automatically try to determine suitable limits? I have the impression that the selection of the midpoint and maximum (1 and 40, respectively) is not appropriate here.

Best wishes,
Richard

Hi, Alex

Thanks for the response!
I apologize that I wrote code that was difficult to read.
I will try to provide minimal working examples from now on.

Thank you very much for your advice!
I have tried to follow your advice, but there was no change.
What I did is following:

adjacency = spatial_src_adjacency(mne.read_source_spaces('fsaverage-ico-5-src.fif'))
T_obs, clusters, cluster_p_values, H0 = clu =\
    spatio_temporal_cluster_test(X, adjacency=adjacency, n_permutations=512, n_jobs=1,
                                 threshold=f_threshold, buffer_size=None)

f = open('clu.pickle','wb')
pickle.dump(clu,f)
f.close

f = open("clu.pickle","rb")
clu = pickle.load(f)

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='auto')

In addition, I have uploaded the variable used in this analysis which I hope you find useful.

Best,
masataka

Hi, @richard

I rechecked my code and variables over and over again.
And then, I noticed that the result of stc_all_cluster_vis.data.shape might be wrong.
That was (20484, 3), but the second dimension should be larger than this, shouldn’t it?

In addition, I think T_obs, which is the result of spatio_temporal_cluster_test, looks also a little bit strange.
At any time point, most are between 10 and 20, but those are obviously too high.

Are these cause the problem that I wrote above?

Best Regard,
Masataka

Quoting from the Returns section of summarize_clusters_stc():

The first time point in this SourceEstimate object is the summation of all the clusters. Subsequent time points contain each individual cluster.

So if your shape is (20484, 3) that means there were 2 significant clusters.

Hi @drammock

I see, Thanks a lot!
The number of significant clusters that I found was 2 as you said.

I will go over the entire process again carefully.
Thanks a lot!

Best,
Masataka