run_envcorr_adapted.py

import mne
import os.path as op

import numpy as np

from mne.connectivity import envelope_correlation
from mne.preprocessing import compute_proj_ecg, compute_proj_eog

# Adapted from:
# https://github.com/SherazKhan/mne-camcan/blob/master/examples/run_envcorr.py
# Changes:
# - switch to oct6, remove overlap from epochs, fix mag coil types,
#   use morph_labels, switch to faster hilbert, get rid of get_stc
#   get rid of bct, switch to fsaverage, remove tSSS, change flip to mean,
#   modify method, modify SNR

subject = 'CC110033'
subjects_dir = 'recons'
inv_fname = op.join('rest', 'sub-CC110033', 'meg', 'sub-CC11033-inv.fif')
raw_fname = op.join('rest', 'sub-' + subject, 'meg', 'rest_raw.fif')

##############################################################################
# Load and preprocess the data

raw = mne.io.read_raw_fif(raw_fname)
raw.load_data()
raw.info['bads'] += [u'MEG2113', u'MEG1941', u'MEG1412', u'MEG2331']
raw.filter(None, 40)
projs_ecg, ecg_events = compute_proj_ecg(raw, n_grad=1, n_mag=2)
projs_eog, eog_events = compute_proj_eog(raw, n_grad=1, n_mag=2)
raw.info['projs'] = projs_ecg + projs_eog
raw.apply_proj()
cov = mne.compute_raw_covariance(raw, tmin=0, tmax=None)
raw.filter(14, 30)
reject = dict(grad=1500e-13, mag=2e-12)
duration = 30.  # seconds
events = mne.make_fixed_length_events(raw, duration=duration)
epochs = mne.Epochs(raw, events, tmin=0, tmax=duration,
                    baseline=None, preload=True, reject=reject, decim=10)
print(len(epochs), len(epochs.drop_log) - len(epochs))

##############################################################################
# Load and morph labels

mne.datasets.fetch_aparc_sub_parcellation(subjects_dir=subjects_dir)
labels_orig = mne.read_labels_from_annot('fsaverage', 'aparc_sub', 'both',
                                         subjects_dir=subjects_dir)
labels_orig = [label for label in labels_orig
               if not label.name.startswith('unknown')]
labels = mne.morph_labels(labels_orig, subject, subjects_dir=subjects_dir)

##############################################################################
# Apply the inverse

inv = mne.minimum_norm.read_inverse_operator(inv_fname)
stcs = mne.minimum_norm.apply_inverse_epochs(
    epochs, inv, 1. / 9., method='MNE', pick_ori="normal",
    return_generator=True)
labels_data = mne.extract_label_time_course(
    stcs, labels, inv['src'], return_generator=True, verbose=True)

##############################################################################
# Compute the envelope correlation

corr = envelope_correlation(labels_data)
deg = mne.connectivity.degree(corr, 0.15)
stc = mne.labels_to_stc(labels_orig, deg, subject='fsaverage')
stc = stc.in_label(mne.Label(np.arange(10242), hemi='lh') +
                   mne.Label(np.arange(10242), hemi='rh'))
lims = len(labels_orig) * np.array([0.2, 0.4, 0.6])
lims = [75, 85, 95]
brain = stc.plot(subject='fsaverage', hemi='split', colormap='gnuplot',
                 views=['lat', 'med'],
                 clim=dict(kind='percent', lims=lims),
                 surface='inflated', subjects_dir=subjects_dir,
                 smoothing_steps=25)
brain.save_image('run_envcorr_adapted.png')