I would like to plot only the significant vertices in source space but there is no masking option as far as I can see. I wonder if there is maybe a valid reason not to only plot significantly different vertices in source space?
A second question is related to the source estimate: I would like to work with single trial power to estimate brain-behavior relationships. Can I get the source power on a single trial level?
a_source_power, freqs = apply_dics_csd(csd_a, filter_a)
b_source_power, freqs = apply_dics_csd(csd_b, filter_b)
stc = a_source_power - b_source_power
stc_all.append(stc)
contrast_array = np.array([s.data for s in stc_all])
contrast_array = np.squeeze(contrast_array) # get rid of 0 dimension
# run permutation test over source space
t, p, ho = mne.stats.permutation_t_test(contrast_array)
print(sum(p < 0.05))
# take mean over first dimension (subjects)
contrast_mean = np.mean(contrast_array, axis=0)
# convert array back to source estimate structure
contrast_stc = mne.SourceEstimate(contrast_mean,
stc_all[0].vertices, tmin=tmin,
tstep=0.004)
vertno_max, time_max = contrast_stc.get_peak()
surfer_kwargs = dict(
hemi='both', subjects_dir=subjects_dir, subject=subject,
views='flat', initial_time=tmax, time_unit='s',
size=(800, 800), smoothing_steps=10, surface='flat')
brain = contrast_stc.plot(**surfer_kwargs)
brain.add_foci(vertno_max, coords_as_verts=True, hemi='rh',
color='blue',
scale_factor=0.6, alpha=0.5)
brain.add_text(0.1, 0.9,
'dics (beamforming) (plus location of maximal activation)',
'title',
font_size=14)`
Masking an STC is possible. Make a copy of the STC, set the data value of all non-significant vertices to zero, and then use stc_to_label() to create a custom label (example). Then go back to your original STC, and use stc.in_label() to restrict it to the label before plotting.
Maybe @britta-wstnr can help with the second question?
Thank you for your help Dan. I created now the label based on my mask. I have a question though:
Here is the code I implemented (stc_all is a list that holds the the source estimate from a power difference between conditions estimate using dics beamformer for each participant):
contrast_array = np.array([s.data for s in stc_all])
# dims = 43, 8196
contrast_array = np.squeeze(contrast_array) # get rid of 0 dimension
# run permutation test over source space
t, p, ho = mne.stats.permutation_t_test(contrast_array)
print(sum(p < 0.05))
# 15 significantly different vertices
p_bool = p < 0.05
# get the indices of the significant vertices
ind = list(np.where(p < 0.05))
only_sig = contrast_array[:, ind[0]]
# take mean over first dimension (subjects)
contrast_mean = np.mean(contrast_array, axis=0)
# TODO: convert to list comprehension
# now we mask all the non significant vertices
mask = []
for index, value in enumerate(contrast_mean):
if p_bool[index] == 0:
t = 0
mask.append(t)
else:
t = value
mask.append(t)
mask_array = np.array(mask)
# convert mask back to source space
stc_masked = mne.SourceEstimate(mask_array,
stc_all[0].vertices, tmin=-0.4,
tstep=0.004)
# now create the functional label from the masked stc
beta_label = mne.stc_to_label(stc_masked, src=src,
subjects_dir=subjects_dir,
connected=True, smooth=True)
# beta label is a list of size 2: first list entry is left hemisphere , second entry is right hemisphere
# only vertices in rh: 6 labels with with 34 to 301 verices?
# convert array back to source estimate structure
contrast_stc = mne.SourceEstimate(contrast_mean,
stc_all[0].vertices, tmin=tmin,
tstep=0.004)
# now mask the contrast (first entry in beta label is hemisphere:
# 0 = lh, 1 = rh, second index is the label to show, ideally show all
stc_masked = contrast_stc.in_label(beta_label[1][0])
my question is now:
the permutation test over all vertices (8000 something) tells me there are 15 significant vertices. Now the label includes 6 labels in the right hemisphere with 34 to 300 vertices each. How does this work?
Am I missing something or is there something wrong in my code?
ok, so contrast_array is 43 participants x 8196 vertices (and you’ve already collapsed across the time dimension?), and the values are “estimated within-subject difference in power between experimental conditions”.
I think maybe a simpler approach than the one I mentioned before is to skip the whole process of making a label; depending on what you want to do this might not be enough though:
contrast_stc = stc_all[0].copy()
contrast_stc.data = np.where(p < 0.05, contrast_all.mean(axis=0), 0)
contrast_stc.plot()
# if you need the label object anyway:
# beta_label = mne.stc_to_label(contrast_stc, ...)
I think the reason that the labels have way more vertices than you expect them to have is that the labels are being created from a higher-resolution source space than the one implicitly contained in the STC (because you specified smooth=True). Also, regarding your code comment about “ideally show all labels”: you can combine the 6 RH labels in beta[1] with the + operator, to make one (spatially non-contiguous) label.