- MNE version: 1.0.3
- operating system: Windows 10
Dear MNE users,
Thanks a lot for all the documentation and tutorials offered on the MNE website.
I need to analyse some EEG data from patients with disorders of consciousness. The experiment had a very simple design with two events: patients were asked to try to move their thumb or to relax (in alternating order), during a time course of 5 minutes (36 events move / 36 events relax).
I would like to perform an ERD analysis (similar to this one), to see whether there are changes in power in the mu (8-13Hz) and beta (13-30Hz) bands over the sensorimotor cortex, when patients were asked to move.
Even though there are many helpful examples on the MNE website on how to preprocess EEG data, I’m always a bit confused when it comes to joining the different steps.
I’m trying to build my own pipeline using autoreject and ICA to clean the data, before moving to the ERD analysis. Could someone tell me if I’m following all the right steps, and in the right order?
Steps:
- definying montage
- applying low and high pass filters (1-40Hz)
- finding events / epoching: do I need to use a baseline of -0.5 when epoching around the events? So that when I analyse the ERD of the epochs “move” I can directly see what changes from before time 0?
- bad electrodes detection using RANSAC (which are excluded from ICA): I wanted to automate bad channels detection using RANSAC from autoreject and then passing the list of bad channels to the epoched file so that ICA won’t consider them. Is it correct?
- bad epochs detection (before ICA): as suggested in this autoreject example, I
first need to run autoreject, then ICA then autoreject again. Have I understood correctly?
- ICA without bad electrodes (as marked by RANSAC) / bad epochs (defined with the first round of autoreject)
- I need to visually inspect the ICA components and mark eye blinks
- apply ICA on the epochs
- re run autoreject on the epochs, including all channels (even those that were marked as bad by ransac)
- run the baseline correction
- save the epochs.
I made some sort of preliminary code with all those steps. Could someone tell me if what I’m doing is correct? Do I need to add a rejection threshold somewhere in the code? I have also tried to identify eye blinks using create_eog_epochs after having defined 2 electrodes as EOG electrodes. However it does not seem to be finding any EOG events: why could that be?
I appreciate any tip / answer. Thanks a lot!
Here’s the code.
def preprocess(session_name, output_folder, raw, montage, stim_channel, events_id = dict(),
tmin = -.5, tmax = 4, filter = (1, 40), reject=dict(eeg=300e-6), baseline = None, save = False):
print('Preprocessing: ' + session_name)
print('Output folder: ' + output_folder)
import os
import autoreject
import time
import mne
import matplotlib.pyplot as plt
import numpy as np
# Specifying montage
montage = mne.channels.make_standard_montage(montage)
# Filtering Data
raw_filtered = raw.copy().filter(filter[0], filter[1])
# Reference on average
raw_filtered.set_eeg_reference(ref_channels='average')
picks_eeg = mne.pick_types(raw.info, meg=False, eeg=True, eog=True, stim=False,
exclude='bads')
# Finding events
events = mne.find_events(raw, stim_channel='STI 014', verbose=True) # or STIM
# Epoching data
epochs = mne.Epochs(raw_filtered, events, events_id, tmin, tmax, baseline=None,
reject=None, verbose=False, detrend=0, preload=True)
# Bad channel detection (ransac)
from autoreject import Ransac # this uses ransac from PREP
print('Detecting bad channels with Ransac')
print(time.strftime("%H:%M:%S", time.localtime()))
ransac = Ransac(verbose=True, picks=picks_eeg, n_jobs=-1)
epochs_clean = ransac.fit_transform(epochs)
bad_chs = ransac.bad_chs_ # list with bad channels according to RANSAC
print('Bad channels detected: ' + str(len(bad_chs)))
print(time.strftime("%H:%M:%S", time.localtime()))
# Bad epochs rejection (autoreject before ICA)
print('Bad epochs detection (autorject)')
ar = autoreject.AutoReject(n_interpolate=[1, 2, 3, 4], random_state=11,
n_jobs=-1, verbose=True)
ar.fit(epochs[:20]) # fit on a few epochs to save time
print(time.strftime("%H:%M:%S", time.localtime()))
epochs_ar, reject_log = ar.transform(epochs, return_log=True)
reject_log.plot('horizontal', show = False)
plt.savefig(output_folder + session_name + '_autoreject_before_ica.jpg', dpi = 160)
print(f'Rejected Epochs: {(reject_log.bad_epochs == True).sum()}')
print(time.strftime("%H:%M:%S", time.localtime()))
# ICA without bad epochs / bad channelss
epochs.info['bads'] = bad_chs
print('Fitting ICA')
ica = mne.preprocessing.ICA(random_state=99)
ica.fit(epochs[~reject_log.bad_epochs])
print('Completed.')
print(time.strftime("%H:%M:%S", time.localtime()))
# ICA visual inspection
print('Visually inspect results and select eye-blinks components...')
ica.plot_components(picks = np.arange(0,30,1))
print(ica.exclude) # list of components excluded
ica.plot_components(ica.exclude, show = False)
plt.savefig(output_folder + session_name + '_ica_bad_components.jpg', dpi = 160)
ica.plot_overlay(epochs.average(), exclude=ica.exclude, show = False)
plt.savefig(output_folder + session_name + '_ica_overlay.jpg', dpi = 160)
ica.apply(epochs, exclude=ica.exclude)
# Autoreject after high pass filter and ICA
epochs.info['bads'] = []
ar = autoreject.AutoReject(n_interpolate=[1, 2, 3, 4], random_state=11,
n_jobs=-1, verbose=True)
ar.fit(epochs[:20]) # fit on the first 20 epochs to save time
epochs_ar, reject_log = ar.transform(epochs, return_log=True)
epochs[reject_log.bad_epochs].plot(scalings=dict(eeg=100e-6))
reject_log.plot('horizontal', show = False)
plt.savefig(output_folder + session_name + '_autoreject_after_ica.jpg', dpi = 160)
# We will visualize the cleaned average data and compare it against the bad segments.
evoked_bad = epochs[reject_log.bad_epochs].average()
plt.figure()
plt.plot(evoked_bad.times, evoked_bad.data.T * 1e6, 'r', zorder=-1, show = False)
epochs_ar.average().plot(axes=plt.gca(), show = False)
plt.savefig(output_folder + session_name + '_epochs_ar_average_badseg.jpg', dpi = 160)
# Re-Referencing after ICA
mne.set_eeg_reference(epochs_ar, ref_channels='average', copy = False)
epochs_ar.apply_baseline()
# Saving cleaned epochs
if save == True:
epochs_ar.save(output_folder + session_name + '_cleaned_epochs.mff' ,overwrite=True)
return epochs_ar