Today I revisited the whole brain volume analysis to find the spatial regions of evoked potentials.
The code is below. The strategy was a bit different this time. Rather than comparing difference
waves, I applied the inverse operator to each epoch, then performed a t-test on each voxel at each
time point, with false discovery rate correction. I save the t-value map for viewing in freeview.
Finally, I show the peak activation of the right amygdala for comparison with previous results.
The threshold for significance after correction is 3.35. Consistent with yesterday's results, there
is a brief window of significant activity around 1.9 s after stimulus presentation.
On Monday, I'll extract voxels from the right amygdala and average them together, to look at their
waveform and compare to yesterday's analysis.
%matplotlib inline
%load ~/PythonEEG/spm_volume.py
import os.path as op
import numpy as np
from scipy import stats
import nibabel as nib
import mne
from mne.datasets import spm_face
from mne.minimum_norm import make_inverse_operator, apply_inverse_epochs
data_path = spm_face.data_path()
subject = 'spm'
subjects_dir = data_path + '/subjects'
subject_dir = subjects_dir + '/' + subject
meg_dir = data_path + '/MEG/' + subject
epo_fname = meg_dir + '/SPM_CTF_MEG_example_faces1_3D_epo.fif'
mgz_fname = subject_dir + '/mri/brain.mgz'
bem_fname = subject_dir + '/bem/spm-5120-5120-5120-bem-sol.fif'
mri_fname = meg_dir + '/SPM_CTF_MEG_example_faces1_3D_raw-trans.fif'
fwd_fname = meg_dir + '/SPM_CTF_MEG_example_faces1_3D_vol-7-fwd.fif'
# load the epoched data
epochs = mne.read_epochs(epo_fname)
######################################
# read or compute the forward solution
if op.exists(fwd_fname):
# read the forward solution
fwd = mne.read_forward_solution(fwd_fname)
else:
# make volume source space
src = mne.setup_volume_source_space(subject, pos=7.0, mri=mgz_fname,
bem=bem_fname)
# make forward solution
fwd = mne.make_forward_solution(epochs.info, mri=mri_fname, src=src,
bem=bem_fname)
# write the source space to file
mne.write_forward_solution(fwd_fname, fwd, overwrite=True)
##########################
# compute inverse solution
snr = 5.0
lambda2 = 1.0 / snr ** 2
method = 'dSPM'
# estimate noise covariance
noise_cov = mne.compute_covariance(epochs.crop(None, 0, copy=True))
# make the inverse operator
inv = make_inverse_operator(epochs.info, fwd, noise_cov, loose=None, depth=None)
# apply the operator to the epoched data
stcs = apply_inverse_epochs(epochs.crop(0, None, copy=False), inv, lambda2, method)
# create an empty array to store data
X = np.zeros((len(epochs), stcs[0].shape[0], stcs[0].shape[1]))
# iterate through the epochs
for i, s in enumerate(stcs):
X[i] = s.data
# compare the two conditions
t, p = stats.ttest_ind(X[epochs.events[:,2]==1], X[epochs.events[:,2]==2])
# correct for multiple comparisons
sig, pcorr = mne.stats.fdr_correction(p)
# save t-values as source estimate volume
stc = s.copy()
stc._data = t
img = stc.save_as_volume(fwd['src'], mri_resolution=False)
nib.save(img, 'spm_volume_tmap.nii.gz')
from IPython.display import Image
Image('/Users/Alan/Desktop/Screen Shot 2014-05-30 at 5.08.56 PM.png')
