Preprocessing

Preprocessing.py

+# -*- coding: utf-8 -*-
+"""
+Updated Oct 18 2022
+
+@author: Qianliang Li (glia@dtu.dk)
+
+The following preprocessing steps were employed:
+    1. Load the data
+    2. 1-100Hz Bandpass filtering
+    3. 50Hz Notch filtering
+    4. Retrieving event information about eyes open/closed states
+    5. Epoching to 4s non-overlapping segments
+    6. Visual inspection of the data
+        Bad channels and non-ocular artifacts are removed
+    7. Robust re-reference to common average (without bad channels)
+    8. Interpolation of bad channels
+    9. ICA artifact removal of ocular and ECG artifacts
+    10. AutoReject guided final visual inspection
+
+Due to privacy issues of clinical data, our data is not publically available.
+For demonstration purposes, I am using a publically available EEG dataset
+and treating it as resting-state eyes open/closed
+Link to the demonstration data: www.bci2000.org
+"""
+
+# Set working directory
+import os
+wkdir = "/home/glia/EEG"
+os.chdir(wkdir)
+
+# Load all libraries from the Preamble
+from Preamble import *
+
+# To demonstrate the script, a publically available EEG dataset are used
+# EEG recordings from 2 subjects are used as an example
+# The EEGBCI200 is task based, but we will treat it as "resting-state"
+# And the 2 runs as Eyes Closed and Eyes Open
+n_subjects = 2
+Subject_id = [1,2]
+# Download and get filenames
+from mne.datasets import eegbci
+files = []
+for i in range(n_subjects):
+    raw_fnames = eegbci.load_data(Subject_id[i], [1,2]) # The first 2 runs
+    files.append(raw_fnames)
+
+# # Original code to get filenames in a folder
+# data_path = "EEG_folder"
+# # Get filenames
+# files = []
+# for r, d, f in os.walk(data_path):
+#     for file in f:
+#         if ".bdf" in file:
+#             files.append(os.path.join(r, file))
+
+# Eye status
+anno_to_event = {'Eyes Closed': 1, 'Eyes Open': 2} # manually defined event id
+eye_status = list(anno_to_event.keys())
+n_eye_status = len(eye_status)
+
+# Epoch settings
+epoch_len = 4 # length of epochs in seconds
+n_epochs_trial = int(60/epoch_len) # number of epochs in a trial
+n_trials = 2 # 1 eyes closed followed by 1 eyes open
+
+# Montage settings
+montage = mne.channels.make_standard_montage('standard_1005')
+#montage = mne.channels.read_custom_montage(filename) # custom montage file
+
+# %% Load, filter and epoch (Steps 1 to 5)
+# Pre-allocate memory
+epochs = [0]*n_subjects
+for i in range(n_subjects):
+    # MNE python supports many EEG formats. Make sure to use the proper one
+    raw = mne.io.concatenate_raws([mne.io.read_raw_edf(f, preload=True) for f in files[i]])
+    # Fix EEGBCI channel names
+    eegbci.standardize(raw)
+    # Set montage
+    raw.set_montage(montage)
+    # Only use the EEG channels
+    raw.pick_types(meg=False, eeg=True, stim=False)
+    # Bandpass filter (will also detrend, median = 0)
+    raw.filter(1, 70, fir_design="firwin", verbose=0) # Due to sfreq = 160 I cannot lowpass at 100Hz
+    # raw.filter(1, 100, fir_design="firwin", verbose=0) # original line
+    # Notch filter to remove power-line noise (50Hz is common in Europe)
+    raw.notch_filter(50, fir_design="firwin", verbose=0)
+    # Epoch to 4 sec
+    # The first 60s are treated as Eyes Closed and following 60s as Eyes Open
+    event = np.zeros((int(2*n_epochs_trial),3), dtype=int) # manually make event array (n_events,3)
+    event[:,0] = np.array(np.linspace(0,2*60-epoch_len,int(2*n_epochs_trial))*int(raw.info["sfreq"]), dtype=int) # first column, the time for events
+    # Hardcoded based on data format
+    event[:n_epochs_trial,2] = 1 # Eyes closed
+    event[n_epochs_trial:,2] = 2 # Eyes open
+    # Make the epochs. (The -1/int(raw[n].info["sfreq"]) is needed because python start with 0, so 30000 points is 0 to 29999
+    epochs[i] = mne.Epochs(raw, event, event_id=anno_to_event,
+                 tmin=0, tmax=epoch_len-1/int(raw.info["sfreq"]),baseline=None, verbose=0).load_data()
+    
+    print("Subject:{} finished epoching ({}/{})".format(Subject_id[i],i+1,n_subjects))
+
+# Compute the number of epochs in each trial
+n_epochs = pd.DataFrame(np.zeros((n_subjects,n_trials)))
+for i in range(n_subjects):
+    for t in range(n_trials):
+        try:
+            n_epochs_ec = len(epochs[i][t:int((t+1)*n_epochs_trial)][eye_status[0]])
+        except:
+            n_epochs_ec = 0
+        try: 
+            n_epochs_eo = len(epochs[i][t:int((t+1)*n_epochs_trial)][eye_status[1]])
+        except:
+            n_epochs_eo = 0
+        n_epochs.iloc[i,t] = np.max([n_epochs_ec,n_epochs_eo])
+# Calculate cumulative sum which is used for later indexing
+cumsum_n_epochs = np.cumsum(n_epochs, axis=1)
+# Insert 0 as first column
+cumsum_n_epochs.insert(0,"Start",0)
+
+# %% First visual inspection, re-referencing and ch interpolation (Steps 6 to 8)
+# When in doubt, I did not exclude the data
+# Ocular artifacts are not removed, as we will use ICA to correct for those
+# Bad epochs are dropped
+# Bad channels are removed and interpolated on a trial basis
+# i.e. for each 1min eyes open or eyes closed trial
+
+# Visualize each
+epochs[0].plot(scalings=200e-6, n_epochs=15, n_channels=24)
+epochs[1].plot(scalings=200e-6, n_epochs=15, n_channels=24)
+
+# Pre-allocate memory
+bad_channels = [0]*len(epochs)
+bad_epochs = [0]*len(epochs)
+
+# Manually found bad data segments
+#bad_channels[0] = [0] # If all trials are fine I do not make a line with bad channels
+bad_epochs[0] = [6] # Epoch 6
+bad_channels[1] = ["T7",0] # T7 is bad in first eyes closed trial
+bad_epochs[1] = [12,13] # 2 bad epochs
+
+# Cautionary note: bad_epochs are using index. Not the epoch number in MNE-viewer
+# They are often the same, but not always, like in this demonstration
+# where there are no Epoch 15 because it was not full length data due to
+# the way I defined the event times, so when you drop Epoch idx 15, it is the
+# epoch that MNE labeled as 16 in this example!
+
+# Pre-allocate memory
+cleaned_epochs = epochs.copy()
+
+# Interpolate bad channels and re-reference to robust common average
+bad_ch_counter = 0 # count how many bad channels for summary
+bad_ch_idx = np.array(bad_channels,dtype=object).nonzero()[0]
+for i in range(len(bad_ch_idx)):
+    n = bad_ch_idx[i]
+    n_epochs_trial = int(60/epoch_len)
+    # Get epoch number for each subject
+    epoch_number = np.arange(0,len(epochs[n]))
+    # Pre-allocate memory
+    temp_epochs = [0]*n_trials
+    for trial in range(n_trials):
+        # Retrieve trial
+        trial_idx = (epoch_number >= cumsum_n_epochs.iloc[n,trial]) & (epoch_number < cumsum_n_epochs.iloc[n,trial+1])
+        temp_epochs[trial] = epochs[n].copy().drop(np.invert(trial_idx))
+        # Set bad channel
+        trial_bad_ch = bad_channels[bad_ch_idx[i]][trial]
+        if trial_bad_ch == 0:
+            # Do not perform interpolation, only re-referencing to average
+            temp_epochs[trial].set_eeg_reference(ref_channels="average", verbose=0)
+        else:
+            if type(trial_bad_ch) == str: # fix if only one ch is provided
+                trial_bad_ch = [trial_bad_ch] # make it to list
+            temp_epochs[trial].info["bads"] = trial_bad_ch
+            # Re-reference to Cz
+            temp_epochs[trial].set_eeg_reference(ref_channels="average", verbose=0)
+            # Interpolate the bad channels
+            temp_epochs[trial].interpolate_bads(reset_bads=True)
+            # Increase counter
+            bad_ch_counter += len(trial_bad_ch)
+    # Concatenate temporary epoch and save
+    cleaned_epochs[n] = mne.concatenate_epochs(temp_epochs, add_offset=False)
+    # Notice that an offset is still added when using to_data_frame!
+
+# Re-reference all other data (that did not have bad channels) to common avg
+good_ch_idx = np.where(np.array(bad_channels,dtype=object) == 0)[0]
+for i in range(len(good_ch_idx)):
+    n = good_ch_idx[i]
+    cleaned_epochs[n].set_eeg_reference(ref_channels="average", verbose=0)
+
+# Drop bad epochs
+bad_epoch_counter = 0
+bad_epochs_idx = np.array(bad_epochs,dtype=object).nonzero()[0]
+for i in range(len(bad_epochs_idx)):
+    n = bad_epochs_idx[i]
+    subject_bad_epochs = bad_epochs[n]
+    cleaned_epochs[n].drop(subject_bad_epochs)
+    bad_epoch_counter += len(subject_bad_epochs)
+
+# Summarize how many bad channels and epochs there were manually defined
+(bad_ch_counter/(n_trials*epochs[0].info["nchan"]*n_subjects))*100 # 0.4% bad channels
+(bad_epoch_counter/(141+150*(n_subjects-1)))*100 # 1% bad epochs rejected
+
+
+# %% ICA is performed to remove eye blinks, ECG and EOG artifacts
+# Make list to contain all ICAs
+ica = [0]*len(cleaned_epochs)
+# Make ICA objects
+for n in range(len(cleaned_epochs)):
+    # Matrix rank is -1 because I used common average reference
+    # If any channels were interpolated the rank is further reduced
+    if any([i in [n] for i in bad_ch_idx]):
+        inter_pol_ch = [i for i in bad_channels[n] if i != 0] # remove all 0
+        # Make a flat list to use np.unique
+        flat_list = []
+        for sublist in inter_pol_ch:
+            if type(sublist) == str: # fix if only one ch is provided
+                sublist = [sublist]
+            for item in sublist:
+                flat_list.append(item)
+        n_inter_pol_ch = len(np.unique(flat_list))
+        matrix_rank = cleaned_epochs[n].info["nchan"]-1-n_inter_pol_ch
+    else:
+        matrix_rank = cleaned_epochs[n].info["nchan"]-1
+
+    ica[n] = mne.preprocessing.ICA(method="fastica", random_state=42, verbose=0,
+                                       max_iter=500, n_components=matrix_rank)
+    ica[n].fit(cleaned_epochs[n])
+    print("{} out of {} ICAs processed".format(n+1,len(cleaned_epochs)))
+
+# Plot the components for visual inspection
+def ica_analysis(n):
+    plt.close("all")
+    # Plot original
+    cleaned_epochs[n].plot(scalings=200e-6, n_epochs = 10)
+    # Plot ICA - compontents
+    ica[n].plot_sources(cleaned_epochs[n], start = 0, stop = 10)
+    ica[n].plot_sources(cleaned_epochs[n], start = 5, stop = 8) # zoomed in helps for ECG artifact recognition
+    ica[n].plot_components(picks=np.arange(0,20))
+
+# Manual ICA decomposition visualization
+n = 1; ica_analysis(n) # manually change n
+
+# # Plot specific component for further inspection of the marked artifacts
+# ica[n].plot_properties(cleaned_epochs[n], picks = artifacts[n])
+
+# Pre-allocate memeory
+artifacts = [0]*len(cleaned_epochs)
+# Manually determine artifacts
+artifacts[0] = [0, 5, 12] # eye blinks, eye blinks, eye movement
+artifacts[1] = [0, 2] # eye blinks, eye movement
+
+# Remove the artifact components from the signal
+corrected_epochs = cleaned_epochs.copy()
+
+for n in range(len(cleaned_epochs)):
+    # Define the components with artifacts    
+    ica[n].exclude = artifacts[n]
+    # Remove on corrected data
+    ica[n].apply(corrected_epochs[n].load_data())
+
+# Inspect how the ICA worked
+n=0; corrected_epochs[n].plot(scalings=200e-6, n_epochs = 10)
+
+# %% Detect bad epochs automatically using AutoReject
+"""
+AutoReject will use cross-validation to determine optimal Peak-to-Peak threshold
+This threshold will be determined for each channel for each subject
+The threshold is used to mark whether each ch and epoch are bad
+Interpolation is performed on bad epoch/ch segments
+If many neighboring ch are bad, then the algorithm will score the worst ch
+(based on peak-to-peak amplitude) and only interpolate the worst
+If too many ch in one epoch are bad, the epoch will be rejected
+
+I am not using Autoreject directly to interpolate, but as a guide.
+The marked bad ch and epochs are then manually inspected and determined
+whether they should be dropped or not.
+The algorithm is run for each eye status separately
+"""
+
+# Suppress plots
+# import matplotlib
+# matplotlib.use("Agg")
+
+save_path = "./Autoreject_overview" # for autoreject overview
+# Pre-allocate memory
+reject_log = [0]*len(corrected_epochs)
+dropped_epochs = [0]*len(corrected_epochs)
+ar = [0]*len(corrected_epochs)
+mean_threshold = [0]*len(corrected_epochs)
+for i in range(len(corrected_epochs)):
+    reject_log0 = [0]*n_eye_status
+    ar0 = [0]*n_eye_status
+    mean_threshold0 = [0]*n_eye_status
+    drop_epochs0 = [0]*n_eye_status
+    for e in range(n_eye_status):
+        ee = eye_status[e]
+        # Initialize class
+        ar0[e] = AutoReject(consensus=np.linspace(0,1,11), cv=10, n_jobs=8,
+                            verbose=False, random_state=42)
+        # Fit to data - but do not transform
+        ar0[e].fit(corrected_epochs[i][ee])
+        # Get rejection log
+        reject_log0[e] = ar0[e].get_reject_log(corrected_epochs[i][ee])
+        # Plot and save Autorejected epochs
+        fig = reject_log0[e].plot(orientation="horizontal", show=False)
+        fig.savefig(os.path.join(save_path,"AR_" + str(Subject_id_concat[i]) + "_" + str(ee) + ".png"))
+        # Close figure window
+        plt.close(fig)
+        # Save mean peak-to-peak voltage threshold used
+        mean_threshold0[e] = np.mean(list(ar0[e].threshes_.values()))
+        # Save suggested dropped epochs
+        drop_epochs0[e] = reject_log0[e].bad_epochs
+    # Concatenate dropped epochs
+    drop_epochs1 = np.concatenate(drop_epochs0)
+    # Save
+    dropped_epochs[i] = drop_epochs1.nonzero()[0]
+    reject_log[i] = reject_log0
+    ar[i] = ar0
+    mean_threshold[i] = mean_threshold0
+    print("{} out of {} subjects finished autorejecting".format(i+1,len(corrected_epochs)))
+
+# Overview of dropped epochs
+Drop_epochs_df = pd.DataFrame.from_records(dropped_epochs) # convert to dataframe
+Drop_epochs_df.insert(0, "Subject_id", Subject_id)
+
+# Re-enable plots
+# import matplotlib.pyplot as plt
+# %matplotlib qt
+
+print(Drop_epochs_df) # no further bad epochs after ICA
+# But for demonstration purposes we will add one more bad epochs in subject 2
+
+### Dropped epochs are used as guide to manually inspect potential thresholded bad epochs
+# Visualize
+n=1 ; corrected_epochs[n].plot(scalings=100e-6, n_channels=31, n_epochs=15)
+
+bad_epochs_ar = dropped_epochs # make copy before modifying
+# Manually modify when appropriate
+# n = 0 - agreed to all
+bad_epochs_ar[1] = np.array([8], dtype="int64") # added 1 more than suggested
+# n = 1: T7 is an outlier and bad
+
+# Combine dropped epochs from first visual inspection and Autoreject
+bad_epochs_comb = [0]*len(corrected_epochs)
+for i in range(len(corrected_epochs)):
+    # Retrieve bad epochs from first manual inspection
+    bad1 = np.array(bad_epochs[i])
+    # If there are no dropped epochs from first inspection, just use second
+    if type(bad_epochs[i]) == int:
+        bad_epochs_comb[i] = bad_epochs_ar[i]
+        continue
+    else:
+        # Retrieve bad epochs from manual of AR
+        bad2 = bad_epochs_ar[i].copy()
+        # Fix index due to iterative dropped epoch
+        # E.g. if 45 is dropped, then all epochs above 45 will have 1 index lower and should be fixed
+        for drops in range(len(bad1)):
+            bad2[bad2 >= bad1[drops]] += 1
+        # Concatenate from both
+        bad_epochs_comb[i] = np.sort(np.concatenate([bad1,bad2]))
+
+# Convert to dataframe
+bad_epochs_comb_df = pd.DataFrame.from_records(bad_epochs_comb)
+bad_epochs_comb_df.insert(0, "Subject_ID", Subject_id)
+# Save
+bad_epochs_comb_df.to_pickle("./Preprocessing/dropped_epochs.pkl")
+
+# Drop bad epochs
+final_epochs = corrected_epochs.copy()
+bad_epochs_df = pd.DataFrame.from_records(bad_epochs_ar)
+bad_epoch_counter2 = 0
+bad_epochs_idx2 = np.where(bad_epochs_df.iloc[:,0].notnull())[0]
+for i in range(len(bad_epochs_idx2)):
+    n = bad_epochs_idx2[i]
+    subject_bad_epochs = bad_epochs_ar[n]
+    final_epochs[n].drop(subject_bad_epochs)
+    bad_epoch_counter2 += len(subject_bad_epochs)
+
+# Summarize how many bad channels and epochs there were manually defined
+(bad_epoch_counter2+bad_epoch_counter)/(141+150*(n_subjects-1))*100 # in total 1.37% epochs rejected
+
+# Re-reference and interpolate
+def re_reference_interpolate(input_epochs, n, bads):
+    # Re-reference to Cz from average
+    input_epochs[n].set_eeg_reference(ref_channels=["Cz"], verbose=0)
+    # Set bads
+    input_epochs[n].info["bads"] = bads
+    # Re-reference to average without bads
+    input_epochs[n].set_eeg_reference(ref_channels="average", verbose=0)
+    # Interpolate
+    input_epochs[n].interpolate_bads(reset_bads=True)
+
+re_reference_interpolate(final_epochs, 1, ["T7"]) # n = 1
+
+# Save all bad channels from first visual inspection and autorejection guided
+bad_channels[1] = ["T7"]*n_trials
+# To keep track of which channels were interpolated
+with open("./Preprocessing/bad_ch.pkl", "wb") as filehandle:
+    pickle.dump(bad_channels, filehandle)
+
+
+# %% Save the preprocessed epochs
+save_path = "./PreprocessedData"
+for n in range(len(corrected_epochs)):
+    final_epochs[n].save(fname = os.path.join(save_path,str("{}_preprocessed".format(Subject_id[n]) + "-epo.fif")),
+                    overwrite=True, verbose=0)
+
+# %% Evaluate dropped epochs to determine threshold for exclusion
+dropped_epochs_df = np.load("./Preprocessing/dropped_epochs.pkl", allow_pickle=True)
+
+# Convert the df to number of dropped epochs
+# In the original data I had multiple datasets and combined the dropped df here
+All_drop_epoch_df = dropped_epochs_df # originally using pd.concat
+
+def drop_epoch_counter(row,df_row):
+    res = np.sum(df_row.iloc[row,1:].notnull())
+    return res
+
+Number_dropped_epochs = [0]*All_drop_epoch_df.shape[0]
+for i in range(All_drop_epoch_df.shape[0]):
+    Number_dropped_epochs[i] = drop_epoch_counter(i,All_drop_epoch_df)
+
+All_drop_epoch_df2 = pd.DataFrame({"Subject_ID":All_drop_epoch_df["Subject_ID"],
+                                   "Number_drop_epoch":Number_dropped_epochs})
+
+# Plot histogram
+bin_seq = range(0,np.max(All_drop_epoch_df2["Number_drop_epoch"])+1)
+All_drop_epoch_df2.hist(column="Number_drop_epoch",figsize=(12,8),bins=bin_seq)
+
+# View all subjects with more than 15 dropped epochs
+All_drop_epoch_df2.loc[All_drop_epoch_df2["Number_drop_epoch"]>=15,:]
+# Number of dropped subjects using above 30 as cutoff
+len(list(All_drop_epoch_df2.loc[All_drop_epoch_df2["Number_drop_epoch"]>30,"Subject_ID"]))
+# Using cutoff for exclusion: 20% or more dropped epochs