-
Notifications
You must be signed in to change notification settings - Fork 3
/
participant.py
246 lines (201 loc) · 8.81 KB
/
participant.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
import numpy as np
import pandas as pd
from mne import Epochs
from mne.time_frequency import tfr_morlet
from .averaging import compute_evokeds
from .epoching import (compute_single_trials, get_bad_channels, get_bad_epochs,
get_events, match_log_to_epochs, read_log,
update_skip_log_rows)
from .io import (read_eeg, save_clean, save_df, save_epochs, save_evokeds,
save_montage, save_report)
from .preprocessing import (add_heog_veog, apply_montage, correct_besa,
correct_ica, interpolate_bad_channels)
from .report import create_report
from .tfr import compute_single_trials_tfr, subtract_evoked
def participant_pipeline(
raw_file,
log_file,
besa_file=None,
bad_channels=None,
auto_bad_channels=None,
skip_log_rows=None,
skip_log_conditions=None,
downsample_sfreq=None,
veog_channels='auto',
heog_channels='auto',
montage='easycap-M1',
ref_channels='average',
ica_method=None,
ica_n_components=None,
highpass_freq=0.1,
lowpass_freq=40.0,
triggers=None,
triggers_column=None,
epochs_tmin=-0.5,
epochs_tmax=1.5,
baseline=(-0.2, 0.0),
reject_peak_to_peak=200.0,
components={'name': [], 'tmin': [], 'tmax': [], 'roi': []},
average_by=None,
perform_tfr=False,
tfr_subtract_evoked=False,
tfr_freqs=np.linspace(4.0, 40.0, num=37),
tfr_cycles=np.linspace(2.0, 20.0, num=37),
tfr_mode='percent',
tfr_baseline=(-0.45, -0.05),
tfr_components={
'name': [], 'tmin': [], 'tmax': [], 'fmin': [], 'fmax': [], 'roi': []},
clean_dir=None,
epochs_dir=None,
trials_dir=None,
evokeds_dir=None,
chanlocs_dir=None,
tfr_dir=None,
report_dir=None,
to_df=True,
):
"""Process EEG data for a single participant.
Performs preprocessing and computes single trial mean amplitudes for ERP
components of interest as well as averaged waveforms.
Parameters & returns
--------------------
See `Usage <../usage.html>`_ for the pipeline input arguments and outputs.
"""
# Backup input arguments for re-use
config = locals()
# Read raw data
raw, participant_id = read_eeg(raw_file)
# Create backup of the raw data for the HTML report
if report_dir is not None:
dirty = raw.copy()
# Downsample
if downsample_sfreq is not None:
sfreq = raw.info['sfreq']
downsample_sfreq = float(downsample_sfreq)
print(f'Downsampling from {sfreq} Hz to {downsample_sfreq} Hz')
raw.resample(downsample_sfreq)
# Add EOG channels
raw = add_heog_veog(raw, veog_channels, heog_channels)
# Apply custom or standard montage
apply_montage(raw, montage)
# Handle any bad channels
raw, interpolated_channels = interpolate_bad_channels(
raw, bad_channels, auto_bad_channels)
# Re-reference to a set of channels or the average
_ = raw.set_eeg_reference(ref_channels)
# Do ocular correction with BESA and/or ICA
if besa_file is not None:
raw = correct_besa(raw, besa_file)
if ica_method is not None:
raw, ica = correct_ica(raw, ica_method, ica_n_components)
else:
ica = None
# Filtering
filt = raw.copy().filter(highpass_freq, lowpass_freq, n_jobs=1, picks='eeg')
# Determine events and the corresponding (selection of) triggers
events, event_id = get_events(filt, triggers)
# Epoching including baseline correction
if baseline is not None:
baseline = tuple(baseline)
epochs = Epochs(filt, events, event_id, epochs_tmin, epochs_tmax, baseline,
preload=True, on_missing='warn')
# Automatically detect bad channels and interpolate if necessary
if bad_channels == 'auto' and auto_bad_channels is None:
auto_bad_channels = get_bad_channels(epochs)
config['auto_bad_channels'] = auto_bad_channels
if auto_bad_channels != []:
print('Restarting with interpolation of bad channels')
return participant_pipeline(**config)
# Add bad ICA components to config
if ica is not None:
if ica_n_components is None or ica_n_components < 1.0:
config['auto_ica_n_components'] = int(ica.n_components_)
config['auto_ica_bad_components'] = [int(x) for x in ica.exclude]
# Drop the last sample to produce a nice even number
_ = epochs.crop(tmin=None, tmax=epochs_tmax, include_tmax=False)
print(epochs.__str__().replace(u"\u2013", "-"))
# Read behavioral log file and match to the epochs
skip_log_rows = update_skip_log_rows(skip_log_rows, epochs)
log = read_log(log_file, skip_log_rows, skip_log_conditions)
if triggers_column is not None:
log, missing_ixs = match_log_to_epochs(epochs, log, triggers_column)
config['auto_missing_epochs'] = missing_ixs
epochs.metadata = log
epochs.metadata.insert(0, column='participant_id', value=participant_id)
# If log file was provided as a DataFrame, convert to dict for config
if isinstance(log_file, pd.DataFrame):
log_file = log_file.astype(object) # Convert NaNs to null for JSON
log_file = log_file.where(pd.notnull(log_file), None)
config['log_file'] = log_file.to_dict(orient='list')
# Get indices of bad epochs
bad_ixs = get_bad_epochs(epochs, reject_peak_to_peak)
config['auto_rejected_epochs'] = bad_ixs
# Compute single trial mean ERP amplitudes and add to metadata
trials = compute_single_trials(epochs, components, bad_ixs)
# Compute evokeds
evokeds, evokeds_df = compute_evokeds(
epochs, average_by, bad_ixs, participant_id)
# Save cleaned continuous data
if clean_dir is not None:
save_clean(filt, clean_dir, participant_id)
# Save channel locations
if chanlocs_dir is not None:
save_montage(epochs, chanlocs_dir)
# Save epochs as data frame and/or MNE object
if epochs_dir is not None:
save_epochs(epochs, epochs_dir, participant_id, to_df)
# Save evokeds as data frame and/or MNE object
if evokeds_dir is not None:
save_evokeds(evokeds, evokeds_df, evokeds_dir, participant_id, to_df)
# Create and save HTML report
if report_dir is not None:
dirty.info['bads'] = interpolated_channels
report = create_report(participant_id, dirty, ica, filt, events,
event_id, epochs, evokeds)
save_report(report, report_dir, participant_id)
# Time-frequency analysis
if perform_tfr:
# Epoching again without filtering
epochs_unfilt = Epochs(raw, events, event_id, epochs_tmin, epochs_tmax,
baseline, preload=True, on_missing='warn',
verbose=False)
# Drop the last sample to produce a nice even number
_ = epochs_unfilt.crop(tmin=None, tmax=epochs_tmax, include_tmax=False)
# Copy original metadata
epochs_unfilt.metadata = epochs.metadata.copy()
# Optionally subtract evoked activity
# See, e.g., https://doi.org/10.1016/j.neuroimage.2006.02.034
if tfr_subtract_evoked:
epochs_unfilt = subtract_evoked(epochs_unfilt, average_by, evokeds)
# Morlet wavelet convolution
print('Doing time-frequency transform with Morlet wavelets')
tfr = tfr_morlet(epochs_unfilt, tfr_freqs, tfr_cycles, use_fft=True,
return_itc=False, n_jobs=1, average=False)
# First, divisive baseline correction using the full epoch
# See https://doi.org/10.3389/fpsyg.2011.00236
if tfr_mode is not None:
tfr_modes = \
['ratio', 'logratio', 'percent', 'zscore', 'zlogratio']
assert tfr_mode in tfr_modes, \
f'`tfr_baseline_mode` must be one of {tfr_modes}'
tfr.apply_baseline(baseline=(None, None), mode=tfr_mode)
# Second, additive baseline correction using the prestimulus interval
if tfr_baseline is not None:
tfr_baseline = tuple(tfr_baseline)
tfr.apply_baseline(baseline=tfr_baseline, mode='mean')
# Reduce numerical precision to reduce object size
tfr.data = np.float32(tfr.data)
# Add single trial mean power to metadata
trials = compute_single_trials_tfr(tfr, tfr_components, bad_ixs)
# Save single trial data (again)
if trials_dir is not None:
save_df(trials, trials_dir, participant_id, suffix='trials')
# Compute evoked power
tfr_evokeds, tfr_evokeds_df = compute_evokeds(
tfr, average_by, bad_ixs, participant_id)
# Save evoked power
if tfr_dir is not None:
save_evokeds(
tfr_evokeds, tfr_evokeds_df, tfr_dir, participant_id, to_df)
return trials, evokeds, evokeds_df, config, tfr_evokeds, tfr_evokeds_df
return trials, evokeds, evokeds_df, config