Errors when morphing source estimates for group analysis

[masataka-wada]

Hi, all!

• MNE-Python version: 0.23.4
• operating system: Win10

I have a problem with morphing source estimates for group analysis.
The dataset consists of 80 subjects out of which 7 subjects have this problem.
I have checked the Freesurfer reconstruction of these subjects and they look OK.

What I am doing right now is as follows:

fsave_vertices = [s['vertno'] for s in src]
morph = mne.compute_source_morph(stc, subject, 'fsaverage',spacing=fsave_vertices, subjects_dir=subjects_dir)

I get the following error:

  File "C:\Users\masaw\Documents\Python Scripts\run\M_Source_group_analysis.py", line 95, in <module>
    morph = mne.compute_source_morph(stc, subject, 'fsaverage',

  File "<decorator-gen-327>", line 22, in compute_source_morph

  File "C:\Users\masaw\anaconda3\envs\mne\lib\site-packages\mne\morph.py", line 248, in compute_source_morph
    morph_mat = _compute_morph_matrix(

  File "C:\Users\masaw\anaconda3\envs\mne\lib\site-packages\mne\morph.py", line 1112, in _compute_morph_matrix
    morpher.append(_hemi_morph(

  File "C:\Users\masaw\anaconda3\envs\mne\lib\site-packages\mne\morph.py", line 1152, in _hemi_morph
    mm = maps[vertices_to] * mm

  File "C:\Users\masaw\anaconda3\envs\mne\lib\site-packages\scipy\sparse\_index.py", line 33, in __getitem__
    row, col = self._validate_indices(key)

  File "C:\Users\masaw\anaconda3\envs\mne\lib\site-packages\scipy\sparse\_index.py", line 138, in _validate_indices
    row = self._asindices(row, M)

  File "C:\Users\masaw\anaconda3\envs\mne\lib\site-packages\scipy\sparse\_index.py", line 170, in _asindices
    raise IndexError('index (%d) out of range' % max_indx)

IndexError: index (165866) out of range

Any help would be appreciated.

Thanks.

[agramfort (Alexandre Gramfort)]

that’s weird. Are you sure the stc, src used are all from the subject processed here?

if so can you share with me a folder and a script to replicate the pb?

thx
Alex

[masataka-wada]

Hi, Alex!

Thanks for your comment!
I pre-processed the data through EEGLAB, but it is not the cause of this problem because other data pre-processed through EEGLAB don’t have this problem.
The other process was conducted through MNE using the following code.

import os.path as op
import os
import glob
import mne
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 = 8
n_subjects2 = 5
datas_path = "E:\\MNE-forum\\"
EEG_path = os.path.join(datas_path,'EEG-datas')
EEG_file = glob.glob(os.path.join(EEG_path, '*.set'))
subjects_dir = os.path.join(datas_path, 'subjects_dir')

for sub_num in range(len(EEG_file)):
    #set montage
    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().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
    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)
        
    # compute_source_morph
    fsave_vertices = [s['vertno'] for s in src]
    morph = mne.compute_source_morph(stc, subject, 'fsaverage',
                                     spacing=fsave_vertices, subjects_dir=subjects_dir)

You can download the data which brings about the problem.
https://drive.google.com/drive/folders/1ti8D4b5rxyEbVZ5xVx6t4yY6MGu5ocgw?usp=sharing

Thanks!!

[agramfort (Alexandre Gramfort)]

what do you get with:

mne.viz.plot_alignment(epochs.info, trans=trans, subject=subject, subjects_dir=subjects_dir, src=src, eeg=True)

just to make sure your head model looks good.

Alex

[masataka-wada]

Thanks for the important advice.

This is the results of the following code,
mne.viz.plot_alignment(evoked.info, trans=trans, subject=subject, subjects_dir=subjects_dir, src=src, eeg=['projected', 'original'])

And this is the results of the following code,

mne.viz.plot_bem(subject=subject, subjects_dir=subjects_dir, 
                 brain_surfaces=[],
                 orientation='coronal')

Finally, I have tried to check the 3D bem surfaces using freesufer with the following code,
mne freeview_bem_surfaces -s subject
And the results are these.

But, they look good.

Masataka

[agramfort (Alexandre Gramfort)]

Ok I think i got it.

You wrote:

fsave_vertices = [s[“vertno”] for s in src]

and you passed this to spacing parameter in compute_source_morph.

the spacing param in compute_source_morph is relative to the target brain ie fsaverage.
So you should pass vertices defined on fsaverage or just set spacing=5 (ie the default)
to end up with an ico5 spacing in fsaverage space

hope this helps
Alex

[masataka-wada]

Hi Alex,

Your excellent advice helped me to solve my problem.
Thank you very much for your kindness.

Best regards,
Masataka