Pipelines

A pipeline is a Julia function that applies a fixed sequence of processing steps to a recording. Pipelines make the analysis reproducible: every file of a study is processed with exactly the same steps and parameters, and the pipeline script documents the analysis.

Example pipelines

NeuroAnalyzer comes with example pipelines in the pipelines/test_pipeline.jl file of the package. Copy the file to your pipelines folder (~/NeuroAnalyzer/pipelines) and use it as a starting point:

using NeuroAnalyzer
src = joinpath(pkgdir(NeuroAnalyzer), "pipelines", "test_pipeline.jl")
dst = joinpath(homedir(), "NeuroAnalyzer", "pipelines", "my_pipeline.jl")
cp(src, dst)

Loading the pipelines:

include(joinpath(homedir(), "NeuroAnalyzer", "pipelines", "my_pipeline.jl"))

The file contains four pipelines:

Pipeline	Input	Steps	Output
`test_pipeline()`	continuous EEG	power line noise removal, band-pass filtering, common average reference	processed `NEURO` object
`preprocess_pipeline()`	continuous EEG	channel locations, DC removal, power line noise removal, band-pass filtering, bad channel detection and interpolation, common average reference, optional ICA	preprocessed `NEURO` object
`spectral_pipeline()`	preprocessed continuous EEG	epoching, epoch rejection, absolute and relative band powers, alpha peak frequency, aperiodic (1/f) offset and exponent	`DataFrame` (one row per channel), rejected epochs
`erp_pipeline()`	preprocessed continuous EEG with markers	epoching around a marker, epoch rejection, baseline correction, averaging, positive and negative peaks	ERP object, rejected epochs, `DataFrame` of peaks

All pipelines work on a copy: the original object is not modified.

Preprocessing

eeg = import_recording("files/eeg.edf")
eeg = preprocess_pipeline(eeg; locs_file = "files/standard-10-20-cap19-elmiko.ced")
save_neuro(eeg; file_name = "files/eeg_clean.hdf")

The most important keyword arguments:

locs_file: channel locations file; required for the interpolation of bad channels (without locations, bad channels are only marked in eeg.header.recording[:bad_channel] and excluded from the average reference)
pl_frq: power line frequency (default: taken from the recording header)
hp, lp: band-pass filter cutoffs in Hz (default: 0.5 and 45 Hz)
bad_channels, amp_t: bad channel detection and its amplitude threshold in μV
ica_n, ica_remove: number of ICA components and the components to remove

ICA components to remove (e.g. eye blinks) must be chosen after inspecting them, so ICA is a two-step procedure:

eeg = preprocess_pipeline(eeg_raw; locs_file = "files/standard-10-20-cap19-elmiko.ced")
ica = ica_decompose(eeg; ch = "eeg", n = 15)
plot_icatopo(eeg; ch = "eeg", ic = ica.ic, ic_mw = ica.ic_mw, tpos = 10)
# components 1 and 3 are eye blinks
eeg = preprocess_pipeline(eeg_raw; locs_file = "files/standard-10-20-cap19-elmiko.ced", ica_n = 15, ica_remove = [1, 3])

Resting-state spectral analysis

r = spectral_pipeline(eeg; ep_len = 4, amp_t = 150)
r.df          # band powers, relative powers, alpha peak frequency, aperiodic offset and exponent
r.rejected    # rejected epochs

Frequency bands are set with bands (names as in band_frq(), default: delta, theta, alpha, beta and gamma 1) and the frequency range of the total power and of the aperiodic fit with flim (default: 1–40 Hz).

Event-related potentials

r = erp_pipeline(eeg; marker = "S  1", offset = 0.2, ep_len = 1.0, amp_t = 100)
plot_erp(r.erp; ch = "Cz")
r.peaks       # latency (s) and amplitude of the positive and negative peak in every channel

The marker value must be one of the values in eeg.markers[!, :value]. The epochs start offset seconds before the marker; the pre-stimulus interval is used as the baseline. The peaks are searched for within twin (default: from the marker to the end of the epoch).

Processing many files

Processing all EDF files in a folder and saving the results:

using DataFrames, CSV

files = filter(endswith(".edf"), readdir("data"; join = true))
results = DataFrame()
for f in files
    eeg = import_recording(f)
    eeg = preprocess_pipeline(eeg; locs_file = "files/standard-10-20-cap19-elmiko.ced")
    save_neuro(eeg; file_name = replace(f, ".edf" => "_clean.hdf"))
    r = spectral_pipeline(eeg)
    r.df[!, :file] .= basename(f)
    append!(results, r.df)
end
CSV.write("data/spectral_results.csv", results)

Tip: Save the preprocessed objects: the analysis steps can then be repeated or changed without preprocessing the recordings again.

Tip: A file that cannot be processed (e.g. all epochs rejected) stops the loop. To skip such files, wrap the loop body in try … catch e … end and record the error.

Writing a pipeline

A pipeline is an ordinary Julia function that takes a NEURO object and keyword arguments for the parameters, and returns the processed object or the results. A minimal pipeline:

function my_pipeline(eeg::NeuroAnalyzer.NEURO; hp::Real = 1, lp::Real = 40)::NeuroAnalyzer.NEURO
    obj = deepcopy(eeg)
    remove_dc!(obj; ch = "eeg")
    NeuroAnalyzer.filter!(obj; ch = "eeg", fprototype = :butterworth, ftype = :hp, cutoff = hp, order = 4)
    NeuroAnalyzer.filter!(obj; ch = "eeg", fprototype = :fir, ftype = :lp, cutoff = lp, bw = 2)
    reference_avg!(obj)
    return obj
end

Recommendations:

work on a copy (deepcopy()) or use the non-mutating functions, so that the original recording is not modified
pass all parameters as keyword arguments with sensible defaults instead of hard-coding them
check the input (e.g. continuous vs. epoched data) and throw ArgumentError with a clear message
keep one pipeline per analysis step (preprocessing, spectral analysis, ERP, …), so that the steps can be combined
every processing step is recorded in obj.history, so the processed object documents how it was created