Plot the spatial extent of a cluster (as those returned from the cluster-based permutation stats) on a brain.

cluster_plotting_example.py

"""
Plot the spatial extent of a cluster (as those returned from the cluster-based
permutation stats) on a brain.

Author: Marijn van Vliet <[email protected]>
"""

import mne
import numpy as np

# Load raw data (mark a bad channel)
path = mne.datasets.sample.data_path() / "MEG" / "sample"
subjects_dir = mne.datasets.sample.data_path() / "subjects"
raw = mne.io.read_raw_fif(path / "sample_audvis_filt-0-40_raw.fif", preload=True)
raw.info["bads"] = ["MEG 2443"]

# Clean blinks a little
regression_model = mne.preprocessing.EOGRegression(picks="meg").fit(raw)
raw_clean = regression_model.apply(raw)

# Create epochs and evoked
epochs = mne.Epochs(
    raw_clean,
    events=mne.find_events(raw_clean, stim_channel="STI 014"),
    event_id=dict(vis_l=3, vis_r=4),  # trials showing visual stimuli
    tmin=-0.2,
    tmax=0.5,
    reject=dict(grad=4000e-13, mag=4e-12),
    preload=True,
)
evoked = epochs.average()

# Load pre-computed inverse operator for source estimation
inverse_operator = mne.minimum_norm.read_inverse_operator(
    path / "sample_audvis-meg-oct-6-meg-inv.fif"
)
stc = mne.minimum_norm.apply_inverse(evoked, inverse_operator, lambda2=1 / 9)

# morph to fsaverage as that's what you are using to plot things.
src_fsaverage = mne.read_source_spaces(
    subjects_dir / "fsaverage" / "bem" / "fsaverage-ico-5-src.fif"
)
morph = mne.compute_source_morph(
    stc,
    subject_from="sample",
    subject_to="fsaverage",
    src_to=src_fsaverage,
    subjects_dir=subjects_dir,
)
stc = morph.apply(stc)

# Create a cluster like those you would get out of the permutation test function.
# We're not going to do statistics in this example, but instead just threshold the
# values to get a cluster.
X = stc.data.T  # this is how the data is fed to the statistics function
cluster = np.where(X > 8)

# A cluster is a tuple of (time_idx, vertices_idx)
time_idx, vertices_idx = cluster

##
# Ok, now for the good part. Start plotting various things. Note that I'm plotting
# everything on the "sample" brain while you are working on "fsaverage".
subjects_dir = mne.datasets.sample.data_path() / "subjects"

# A cluster is defined both in space and time. If we want to plot the boundaries of the
# cluster in space, we must choose a specific time for which to show the boundaries (as
# they change over time). Let's pick 0.1 seconds here, as that's when we have the most
# activity in the sample data we are plotting.
chosen_time = 0.1  # seconds
chosen_time_idx = np.searchsorted(stc.times, chosen_time)

# Select only the vertex indices at the chosen time
vertices_idx = [v for v, t in zip(vertices_idx, time_idx) if t == chosen_time_idx]

# Let's create an anatomical label containing these vertex indices.
# Problem 1): a label must be defined for either the left or right hemisphere. It cannot
# span both hemispheres. So we must filter the vertices based on their hemisphere.
# Problem 2): we have vertex *indices* that need to be transformed into proper vertex
# numbers. Not every vertex in the original high-resolution brain mesh is a source point
# in the source estimate. Do draw nice smooth curves, we need to interpolate the vertex
# indices.

# Both problems can be solved by accessing the vertices defined in the source space
# object. Since we're plotting on the fsaverage brain, we must use the fsaverage source
# space.
src_lh, src_rh = src_fsaverage

# filter the vertices_idx based on their hemisphere
lh_verts, rh_verts = src_lh["vertno"], src_rh["vertno"]
n_lh_verts = len(lh_verts)
cluster_lh_verts = [lh_verts[v] for v in vertices_idx if v < n_lh_verts]
cluster_rh_verts = [rh_verts[v - n_lh_verts] for v in vertices_idx if v >= n_lh_verts]

# vertices must be unique and in increasing order
cluster_lh_verts = np.unique(cluster_lh_verts)
cluster_rh_verts = np.unique(cluster_rh_verts)

# We are now ready to create the anatomical label objects
lh_label = mne.Label(cluster_lh_verts, hemi="lh")
rh_label = mne.Label(cluster_rh_verts, hemi="rh")

# Interpolate the vertices in each label to the full resolution mesh
lh_label = lh_label.fill(src_fsaverage)
rh_label = rh_label.fill(src_fsaverage)

# Plot the evoked response and the two anatomical labels
brain = stc.plot(
    subject="fsaverage",
    subjects_dir=subjects_dir,
    hemi="both",
    initial_time=chosen_time,
)
# borders=1 means show them as a 1px outline instead of a filled shape
brain.add_label(lh_label, borders=1, color="magenta")
brain.add_label(rh_label, borders=1, color="magenta")