#### libraries ####
import os
import mne
import numpy as np
import matplotlib.pyplot as plt
import pickle
import pandas as pd
from scipy import signal
from scipy.io import wavfile
# from pybv import write_brainvision
from pyprep.prep_pipeline import PrepPipeline
from mne_icalabel import label_componentsLaryngeal Preprocessing Pipeline
1 Libraries
2 Trigger latency correction
Fixing trigger timing & renaming events
2.1 parameters
# directory
input_dir = os.getcwd() + '/../raw data/eeg_brainvision/'
output_dir = os.getcwd() + '/../preprocessed/1_trigger_lag_corrected/'
# create a folder if the folder doesn't exist
os.makedirs(output_dir, exist_ok=True)
# create a dictionary for blocks and block markers
block_dict = {
'dorsal_highStan_lowDevi': 1000, # standard: dorsal_high; deviant: dorsal_low
'dorsal_lowStan_highDevi': 2000, # standard: dorsal_low; deviant: dorsal_high
'glottal_highStan_lowDevi': 3000, # standard: glottal_high; deviant: glottal_low
'glottal_lowStan_highDevi': 4000, # standard: glottal_low; deviant: glottal_high
}
# exclude participants with missing data
exclude_ppts = [
'0027', # vot, incomplete eeg recording due to wrong sound library setting
'0028', # f0, missing block 3 and 4
'0030', # missing block 3 and 4
'0034', # missing block 3
'0040', # f0, missing block 3 and 4
'0046', # f0, missing block 2,3,4
]
# trigger searching window (actual trigger time based on audio - trigger time in the data)
t_left = -0.01
t_right = 0.4 # long latency due to using pygame sound library#######################################################
#### create a dictionary for trigger codes and descriptions ####
df = pd.read_csv("vot_mapping.txt", delimiter='\t')
VOT_mapping = dict(zip(df['code'], df['description'])) # key = 'id', value = 'name'
df = pd.read_csv("f0_mapping.txt", delimiter='\t')
F0_mapping = dict(zip(df['code'], df['description'])) # key = 'id', value = 'name'
#######################################################
# initialize a dictionary to save bad stims
all_bad_stim_dict = {}
# read in eeg data files
all_files = os.listdir(input_dir)
#### for each file, create a dictionary to store each block and the indices of trials of that block ####
for file in all_files:
if file.endswith('.vhdr') and (file.split('_')[1] not in exclude_ppts) and (file.split('.')[0]+ '_corr.fif' not in os.listdir(output_dir)):
# read in vhdr files
raw = mne.io.read_raw_brainvision(input_dir + file, preload = True)
# If the aux channel is not named 'StimTrak', change the channel name to 'StimTrak'
if raw.info['ch_names'][31] != 'StimTrak':
# change the channel name
raw.rename_channels({raw.info['ch_names'][31]: 'StimTrak'})
# specify the audio channel
raw.set_channel_types({'StimTrak':'misc'})
# extract the info of experiment version
exp_ver = file.split('.')[0].split('_')[2]
# extract sampling rate
eeg_sfreq = raw.info['sfreq']
#############################################################################################
#### create a dictionary for each stim's standard version and deviant version of markers ####
# intialize a dictionary for triggers
trigger_dict = {}
# read in trigger_codes file, and create a standard trigger code and a deviant triggre code for each stimulus
with open('../experiment programs/'+exp_ver+'_eeg/trigger_codes.txt','r') as tf:
for line in tf:
# read in the current line
line = line.replace('\n','')
# separate fileNames and triggerMarker
label, marker = line.split('\t')
# convert trigger markers to integer
marker = int(marker)
# create label for stims used as standards
trigger_dict[marker] = label + '-s'
# add 100 for deviant marker
marker_deviant = marker + 100
# create label for stims used as deviants
trigger_dict[marker_deviant] = label + '-d'
#############################################################################################
########################
#### get audio data ####
# initialize dictionaries
audio = {}
lengths = {}
# del trigger_dict[99999]
# for each trigger code and label in the trigger dictionary
for marker,label in trigger_dict.items():
# extract file name, up to the -2 character of the label
file_name = label[:-2]
# if not already in audio dictionary, get the info of the audio file
if file_name not in audio:
# get sample rate and data of the audio file
sampleRate, data = wavfile.read('../experiment programs/'+exp_ver+'_eeg/stimuli/{}.wav'.format(file_name))
# calculate sound file length
lengths[file_name] = len(data)/sampleRate
# reduce the sampling rate of the audio file by the factor of int(sampleRate/eeg_sfreq)
data_downsampled = signal.decimate(data, int(sampleRate/eeg_sfreq), ftype='fir')
# add info the audio dictionary
audio[file_name] = data_downsampled
########################
##################################################################
#### recode trigger marker to reflect condition and stim info ####
# for each stimulus, mark the block info
events_from_annot, event_dict = mne.events_from_annotations(raw, verbose='WARNING')
# remove the New Segment marker (which has a trigger value of 99999) and pause marker (222) and continue marker (223)
events_from_annot = events_from_annot[events_from_annot[:, 2] < 200]
# initialize the train for each standards+deviant sequence
train = np.array([]).astype(int)
all_block = {}
isStandard = True # whether the standard in a train has been noted
# loop over each trigger
for i in range(events_from_annot.shape[0]):
# add the current token to the train
train = np.append(train,i)
# if the trigger value is smaller than 100 and the standard toggle is true
if (events_from_annot[i,2]<100) & isStandard:
# split the file name
trigger_splitted = trigger_dict[events_from_annot[i,2]].split('_')
# get the standard category name
block = trigger_splitted[1] + '_' + trigger_splitted[2] + 'Stan'
# toggle standard
isStandard = False
# if the trigger value is over 100
elif events_from_annot[i,2]>100:
# split the file name
trigger_splitted = trigger_dict[events_from_annot[i,2]].split('_')
# append the deviant stim category
block = block + '_' + trigger_splitted[2] + 'Devi'
# if the block is not present in all block
if block not in all_block:
# add the new block and the token idx
all_block[block] = train[2:] # [2:] to exclude the first 2 standards
else:
# add the token idx to the existing block
all_block[block] = np.concatenate([all_block[block],train[2:]],axis = None)
# reset train
train = np.array([]).astype(int)
# toggle standard note
isStandard = True
# loop over each block and its trigger index
for k,v in all_block.items():
# recode the trigger value to reflect block and stim category
events_from_annot[tuple(v),2] += block_dict[k]
# exclude the first 2 standards
events_from_annot = events_from_annot[events_from_annot[:,2]>1000]
##################################################################
######################################################
#### cross correction to detect the trigger delay ####
# initialize delay list
delays = np.array([])
# initialize bad stim list
bad_stim = []
corr_results = []
# loop over each event
for i in range(events_from_annot.shape[0]):
# get current event info [time, duration, annotation]
event = events_from_annot[i]
# get the onset latency (s) of the event
time = event[0]/eeg_sfreq
# get the file name of the event
name = trigger_dict[event[2]%100].split('-')[0]
# get the data from the sound channel
audio_eeg = raw.get_data(
picks = ['StimTrak'],
tmin = time + t_left,
tmax = time + lengths[name] + t_right,
)[0]
# actual stimulus
audio_stim = audio[name]
# Z-score normalization (subtract mean, divide by std)
audio_eeg = (audio_eeg - np.mean(audio_eeg)) / np.std(audio_eeg)
audio_stim = (audio_stim - np.mean(audio_stim)) / np.std(audio_stim)
# cross-correlation
corr = signal.correlate(audio_eeg, audio_stim, mode='full')
# Normalize cross-correlation
corr = corr / (np.linalg.norm(audio_eeg) * np.linalg.norm(audio_stim))
# Find peak correlation value
max_corr = np.max(corr)
# if the maximum correction (sum of products) is less than a threshold (empirical threshold, 0.5 is generally good)
if max_corr < 0.5:
# mark the stim bad
bad_stim.append(i)
# append the maximum correlation
corr_results.append(max_corr)
# the lags for cross-correlation
lags = signal.correlation_lags(
audio_eeg.size,
audio_stim.size,
mode="full")
# get the lag of the maximum cross-correlation
lag = lags[np.argmax(corr)] + t_left*eeg_sfreq
# save the lag for non-starting events
delays = np.append(delays,lag)
######################################################
##################################################################################################
#### plot the wave from the stim track and the eeg channel of the token with the minimum corr ####
min_corr = np.argmin(corr_results)
# get current event info [time, duration, annotation]
event = events_from_annot[min_corr]
# get the onset latency (s) of the event
time = event[0]/eeg_sfreq
# get the file name of the event
name = trigger_dict[event[2]%100].split('-')[0]
# get the data from the sound channel
audio_eeg = raw.get_data(
picks = ['StimTrak'],
tmin = time + t_left,
tmax = time + lengths[name] + t_right,
)[0]
# actual stimulus
audio_stim = audio[name]
# Z-score normalization (subtract mean, divide by std)
audio_eeg = (audio_eeg - np.mean(audio_eeg)) / np.std(audio_eeg)
audio_stim = (audio_stim - np.mean(audio_stim)) / np.std(audio_stim)
# plot
fig, ax = plt.subplots()
ax.plot(audio_eeg, label = 'StimTrak', alpha = 0.6)
ax.plot(audio_stim, label = 'wave', alpha = 0.6)
ax.set_title(file)
ax.legend()
fig.savefig(output_dir + file.split('.')[0] + "_minCor.png", dpi=300, bbox_inches='tight')
##################################################################################################
################################################
#### correct for trigger lag and save files ####
# add number of bad stim info
all_bad_stim_dict[file] = len(bad_stim)
# remove items associated with bad stims from the event list
events_from_annot = np.delete(events_from_annot, bad_stim, 0)
# remove items associated with bad stims from the delay list
delays = np.delete(delays, bad_stim, 0)
# add delay back to the onset latency of each event
events_from_annot[:,0] = events_from_annot[:,0] + delays
# convert individual event marker to conditions
# events_from_annot[:,2] = events_from_annot[:,2] - events_from_annot[:,2]%100
# create annotations
annot_from_events = mne.annotations_from_events(
events = events_from_annot,
event_desc = eval(exp_ver + '_mapping'),
sfreq = eeg_sfreq
)
# set annotations
raw.set_annotations(annot_from_events)
# drop the audio channel in data
raw.drop_channels(['StimTrak'])
# save as a file-into-file data
raw.save(output_dir + file.split('.')[0]+ '_corr.fif')
# save lag data
np.savetxt(output_dir + file.replace('.vhdr', '_delays.txt'), delays, fmt='%i')
################################################
# save the number of bad stims of all participant
with open(output_dir + 'bad_stim.txt', 'a') as f:
for key, value in all_bad_stim_dict.items():
if value > 0:
f.write(key + '\t' + str(value) + '\n')3 Bad channel correction
- filtering
- resampling
- remove line noise
- bad channel detection & repairing
- add back reference channel TP9
3.1 parameters
# set directory
input_dir = os.getcwd() + '/../preprocessed/1_trigger_lag_corrected/'
output_dir = os.getcwd() + '/../preprocessed/2_bad_channel_corrected/'
# create a folder if the folder doesn't exist
os.makedirs(output_dir, exist_ok=True)
# filter cutoff frequencies (low/high)
f_low = 1
f_high = 100
# resampling frequency
f_res = 250
# line frequency
line_freq = 60
# preprocessing parameters
prep_params = {
"ref_chs": 'eeg',
"reref_chs": 'eeg', # average re-reference
"line_freqs": np.arange(line_freq, f_res/2, line_freq),
}
# create a montage file for the pipeline
montage = mne.channels.make_standard_montage("standard_1020")
# interpolation method
# method=dict(eeg="spline")#####################################################
#### Preprocessing (filtering, resampling, bad channel detection/interpoloation, re-reference) ####
#####################################################
# get all file namesin the folder
all_input = os.listdir(input_dir)
all_output = os.listdir(output_dir)
for file in all_input:
if file.endswith("corr.fif") and (file.split('.')[0]+ '_prep.fif' not in all_output):
# read in file
raw = mne.io.read_raw_fif(input_dir + file, preload=True)
# set channel type
raw.set_channel_types({'Fp1':'eog', 'Fp2':'eog'})
# filter
raw.filter(l_freq = f_low, h_freq = f_high)
#### cut off the beginning and ending part ####
# get the onset of the first and the last event ####
events_from_annot, event_dict = mne.events_from_annotations(raw, verbose='WARNING')
# define the beginning time (in seconds)
crop_start = events_from_annot[0][0]/raw.info['sfreq'] - 10
# define the ending time (in seconds)
crop_end = events_from_annot[-1][0]/raw.info['sfreq'] + 10
# crop the data
raw.crop(
tmin=max(crop_start, raw.times[0]),
tmax=min(crop_end, raw.times[-1])
)
# resample
raw.resample(sfreq = f_res)
# read in channel location info
raw.set_montage(montage)
#### Use PrePipeline to preprocess ####
'''
1. detect and interpolate bad channels
2. remove line noise
3. re-reference
'''
# apply pyprep
prep = PrepPipeline(raw, prep_params, montage, random_state=42)
prep.fit()
# export a txt file for the interpolated channel info
with open(output_dir + 'bad_channel.txt', 'a+') as f:
_ =f.write(
file + ':\n' +
"- Bad channels original: {}".format(prep.noisy_channels_original["bad_all"]) + '\n' +
"- Bad channels after robust average reference: {}".format(prep.interpolated_channels) + '\n' +
"- Bad channels after interpolation: {}".format(prep.still_noisy_channels) + '\n'
)
# save the pyprep preprocessed data
raw = prep.raw
# add back the reference channel
raw = mne.add_reference_channels(raw,'TP9')
# add the channel loc info (for the newly added reference channel)
raw.set_montage(montage)
# save
raw.save(output_dir + file.split('.')[0]+ '_prep.fif')4 ICA bad trial correction
4.1 parameters
# directory
input_dir = os.getcwd() + '/../preprocessed/2_bad_channel_corrected/'
output_dir = os.getcwd() + '/../preprocessed/3_ica/'
# create a folder if the folder doesn't exist
os.makedirs(output_dir, exist_ok=True)
# up to which IC you want to consider
ic_upto = 15
# ic_upto = 99# get all file names in the folder
all_input = os.listdir(input_dir)
all_output = os.listdir(output_dir)
# initialize a dictionary for files
for file in all_input:
if file.endswith("prep.fif") and (file.split('.')[0]+ '_ica.fif' not in all_output):
# read in file
raw = mne.io.read_raw_fif(input_dir + file, preload=True)
# make a filtered file copy ICA. It works better on signals with 1 Hz high-pass filtered and 100 Hz low-pass filtered
raw_filt = raw.copy().filter(l_freq = 1, h_freq = 100)
# apply a common average referencing, to comply with the ICLabel requirements
raw_filt.set_eeg_reference("average")
# initialize ica parameters
ica = mne.preprocessing.ICA(
# n_components=0.999999,
max_iter='auto', # n-1
# use ‘extended infomax’ method for fitting the ICA, to comply with the ICLabel requirements
method = 'infomax',
fit_params = dict(extended=True),
random_state = 42,
)
# get ica solution
ica.fit(raw_filt, picks = ['eeg'])
# save ica solutions
ica.save(output_dir + file.split('.')[0]+ '_icaSolution.fif')
# use ICLabel for automatic IC labeling
ic_labels = label_components(raw_filt, ica, method="iclabel")
# save
with open(output_dir + file.split('.')[0]+ '_icLabels.pickle', 'wb') as f:
pickle.dump(ic_labels, f)
# exclude bad IC
labels = ic_labels["labels"]
exclude_idx = [
idx for idx, label in enumerate(labels) if idx<ic_upto and label not in ["brain", "other"]
]
# ica.apply() changes the Raw object in-place
ica.apply(raw, exclude=exclude_idx)
# record the bad ICs in bad_ICs.txt
with open(output_dir + '/bad_ICs.txt', 'a+') as f:
_ = f.write(file + '\t' + str(exclude_idx) + '\n')
# save data
raw.save(output_dir + file.split('.')[0]+ '_ica.fif')5 ERP Segmentation
segmenting continuous EEG into epochs - re-reference - segmentation
5.1 parameters
# directory
input_dir = os.getcwd() + '/../preprocessed/3_ica/'
output_dir = os.getcwd() + '/../preprocessed/4_erp_epochs/' # for ERP
# create a folder if the folder doesn't exist
os.makedirs(output_dir, exist_ok=True)
# epoch window:
erp_t_start = -0.2; erp_t_end = 0.8
baseline = (-0.2, 0)
# criteria to reject epoch
# reject_criteria = dict(eeg = 100e-6) # 100 µV
# reject_criteria = dict(eeg = 150e-6) # 150 µV
reject_criteria = dict(eeg=200e-6) # 200 µV# initialize a list for participants with too many bad trials
too_many_bad_trial_participants = []
# get file names
all_input = os.listdir(input_dir)
all_output = os.listdir(output_dir)
# re-reference, then epoch
for file in all_input:
if file.endswith("ica.fif") and (file.split('_')[2] + '_' + file.split('_')[1] + '_epo.fif' not in all_output):
# read in data
raw = mne.io.read_raw_fif(input_dir + file, preload = True)
# average-mastoids re-reference
raw.set_eeg_reference(ref_channels = ['TP9', 'TP10'])
#### this is for source calculation ####
# filter the data, optional
# raw = raw.filter(l_freq=None, h_freq=30)
# sphere = mne.make_sphere_model('auto', 'auto', raw.info)
# src = mne.setup_volume_source_space(sphere=sphere, exclude=30., pos=15.)
# forward = mne.make_forward_solution(raw.info, trans=None, src=src, bem=sphere)
# raw = raw.set_eeg_reference('REST', forward=forward)
########################################
# pick EEG channels
# picks = mne.pick_types(raw.info, eeg = True)
# get event info for segmentation
events_from_annot, event_dict = mne.events_from_annotations(raw, verbose='WARNING')
# segmentation for ERP
epochs = mne.Epochs(
raw,
events = events_from_annot, event_id = event_dict,
tmin = erp_t_start, tmax = erp_t_end,
# apply baseline correction
baseline = baseline,
# remove epochs that meet the rejection criteria
reject = reject_criteria,
preload = True,
)
##########################################################
#### remove 0-trial events, and log segmentation info ####
ppt = file.split('_')[2] + '_' + file.split('_')[1]
for k, v in event_dict.items():
# good trial count
trial_count = len(epochs[k])
# remove 0 trial event
if trial_count==0:
del epochs.event_id[k]
# good trial rate
goodTrial_rate = round( trial_count/sum(events_from_annot[:,2]==v), 2 )
# record epoch summary
with open(output_dir + 'epoch_summary.txt', 'a+') as f:
_ =f.write(ppt + '\t' + k + '\t' + str(trial_count) + '\t' + str(goodTrial_rate) + '\n')
# mark a participant bad if any condition has fewer than 1/2 trials
if ( goodTrial_rate < 0.5 ):
# mark the participant file as bad
if ppt not in too_many_bad_trial_participants:
too_many_bad_trial_participants.append(ppt)
##########################################################
# save single participant file
epochs.save(output_dir + ppt + '_epo.fif', overwrite=True)
# export the record of bad participants
with open(output_dir + 'too_many_bad_trial_participants.txt', 'w') as f:
for item in too_many_bad_trial_participants:
f.write(item + '\n')6 ERP average
6.1 parameters
# directory
input_dir = os.getcwd() + '/../preprocessed/4_erp_epochs/'
output_dir = os.getcwd() + '/../preprocessed/5_averaged/'
# create a folder if the folder doesn't exist
os.makedirs(output_dir, exist_ok=True)#### get ERP ####
# get file names
all_input = os.listdir(input_dir)
all_output = os.listdir(output_dir)
# initialize a dictionary to store data
all_evokeds = {}
# for each file
for file in all_input:
if file.endswith("_epo.fif"):
# extract participant number
ppt = file.split('_')[0] + '_' + file.split('_')[1]
# read in data
epochs = mne.read_epochs(input_dir + file, preload = True)
# average | get ERP for each condition
evoked = epochs.average(by_event_type=True)
# initialize dictionary for single-participant ERP
all_evokeds[ppt] = {}
# add key for each condition for analysis
for cond in evoked:
# append the evoked data to the dictioncary of evoked data
all_evokeds[ppt][cond.comment] = cond
# Saving the ERP data:
with open(output_dir + '/all_evokeds.pkl', 'wb') as f:
pickle.dump(all_evokeds, f)
del all_evokeds