Build and adjust a pipeline

A pipeline is a list of named steps. Each step is a small scikit-learn style estimator: its parameters are set when it is created, fit learns what it needs from a recording, transform applies it, and afterwards it reports what it did (qc_) and draws it (plot).

This tutorial builds the recommended MEG pipeline for a Neuromag system step by step on MNE’s sample recording, looks at every step, and then changes it.

import tempfile
from pathlib import Path

import mne

import meeg_utils as meu
from meeg_utils import steps as S

mne.set_log_level("ERROR")

data_path = mne.datasets.sample.data_path()
raw = mne.io.read_raw_fif(data_path / "MEG" / "sample" / "sample_audvis_raw.fif")
raw.pick(["meg", "eog", "stim"]).load_data()
raw.info["line_freq"] = 60.0
meu.io.detect_system(raw)

'neuromag'

The steps

Neuromag systems come with a cross-talk and a fine-calibration file, which SSS needs. In a BIDS dataset they are found automatically; here we pass them.

cross_talk = data_path / "SSS" / "ct_sparse_mgh.fif"
calibration = data_path / "SSS" / "sss_cal_mgh.dat"

pipe = meu.Pipeline(
    [
        ("bads", S.BadChannels("maxwell", cross_talk=cross_talk, calibration=calibration)),
        ("sss", S.Maxwell(cross_talk=cross_talk, calibration=calibration)),
        ("highpass", S.Filter(0.1, None)),
        ("resample", S.Resample(250.0)),
        ("line_noise", S.LineNoise("zapline-plus")),
        ("ica", S.ICA(20, picks="meg", method="infomax", labeler="megnet")),
    ]
)
pipe

#	Step	Configuration (changed parameters)
0	bads	BadChannels(calibration=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/sss_cal_mgh.dat'), cross_talk=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/ct_sparse_mgh.fif'), method='maxwell')
1	sss	Maxwell(calibration=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/sss_cal_mgh.dat'), cross_talk=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/ct_sparse_mgh.fif'))
2	highpass	Filter(h_freq=None, l_freq=0.1)
3	resample	Resample(sfreq=250.0)
4	line_noise	LineNoise()
5	ica	ICA(labeler='megnet', n_components=20, picks='meg')

This is the meg-erp preset; meu.Pipeline.preset("meg-erp", system="neuromag") builds the same pipeline and finds the two files in the BIDS dataset of the recording.

clean = pipe.fit_transform(raw)

Bad channels

mne.preprocessing.find_bad_channels_maxwell() scores every channel against its SSS reconstruction. MEG 2443 was already marked bad in the file:

print(pipe.qc_["bads"]["meg"])
figs = pipe["bads"].plot("scores")

{'bads': ['MEG 2313'], 'n_bads': 1, 'fraction': 0.0032679738562091504, 'by_criterion': {'noisy': ['MEG 2313'], 'flat': []}}

SSS

Signal-space separation keeps the fields from inside the head and rebuilds the bad channels. The spectra before and after:

print({k: pipe.qc_["sss"][k] for k in ("reconstructed", "n_basis", "power_change_db")})
figs = pipe["sss"].plot("psd")

{'reconstructed': ['MEG 2443', 'MEG 2313'], 'n_basis': 72, 'power_change_db': {'mag': np.float64(-9.28), 'grad': np.float64(-1.97)}}

Line noise

ZapLine-plus removes the 60 Hz components and their harmonics, chunk by chunk, without notching the spectrum:

figs = pipe["line_noise"].plot("psd")

ICA with MEGnet

MEGnet labels the 20 components (it was trained on 20 Infomax components of 250 Hz data, which is why the pipeline resamples first). Components whose artifact probability is at least 0.8 are removed:

print(pipe.qc_["ica"]["labels"])
figs = pipe["ica"].plot("labels")

['brain/other', 'eye blink', 'brain/other', 'brain/other', 'heart beat', 'brain/other', 'brain/other', 'brain/other', 'brain/other', 'brain/other', 'brain/other', 'brain/other', 'brain/other', 'eye movement', 'brain/other', 'brain/other', 'brain/other', 'brain/other', 'heart beat', 'brain/other']

figs = pipe["ica"].plot("components")

The excluded components, with their spectrum and time course (these figures need the data):

figs = pipe["ica"].plot("properties", inst=raw)

• 
• 

If you disagree with a label, correct it and transform again; the decomposition is kept:

pipe["ica"].relabel({4: "brain"})
clean = pipe.transform(raw)

Change the pipeline

Parameters are addressed as <step>__<parameter>, as in scikit-learn:

pipe.set_params(highpass__l_freq=1.0, ica__threshold=0.9)
pipe["highpass"]

Filter(h_freq=None, l_freq=1.0)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Steps can be added, replaced and removed by name, and a pipeline can be sliced:

pipe.insert_after("line_noise", "lowpass", S.Filter(None, 40.0))
pipe.remove("resample")
list(pipe.named_steps), pipe[:2]

(['bads', 'sss', 'highpass', 'line_noise', 'lowpass', 'ica'], Pipeline(steps=[('bads',
                 BadChannels(calibration=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/sss_cal_mgh.dat'),
                             cross_talk=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/ct_sparse_mgh.fif'),
                             method='maxwell')),
                ('sss',
                 Maxwell(calibration=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/sss_cal_mgh.dat'),
                         cross_talk=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/ct_sparse_mgh.fif')))]))

A pipeline is saved as YAML, to rerun it or to publish it with a paper:

fname = Path(tempfile.mkdtemp()) / "pipeline.yaml"
pipe.to_yaml(fname)
print(fname.read_text())

meeg_utils_version: 0.2.0
copy: true
steps:
- name: bads
  class: meeg_utils.steps.channels:BadChannels
  params:
    calibration: /home/runner/mne_data/MNE-sample-data/SSS/sss_cal_mgh.dat
    cross_talk: /home/runner/mne_data/MNE-sample-data/SSS/ct_sparse_mgh.fif
    h_freq: 40.0
    limit: 7.0
    method: maxwell
    origin: auto
    psd: false
    random_state: 42
    ransac: true
    reject_by_annotation: null
    robust_reference: true
- name: sss
  class: meeg_utils.steps.maxwell:Maxwell
  params:
    bad_condition: error
    calibration: /home/runner/mne_data/MNE-sample-data/SSS/sss_cal_mgh.dat
    coord_frame: head
    cross_talk: /home/runner/mne_data/MNE-sample-data/SSS/ct_sparse_mgh.fif
    destination: null
    ext_order: 3
    head_pos: null
    int_order: 8
    origin: auto
    regularize: in
    remove_chpi: true
    skip_by_annotation:
    - edge
    - bad_acq_skip
    st_correlation: 0.98
    st_duration: null
- name: highpass
  class: meeg_utils.steps.filtering:Filter
  params:
    filter_length: auto
    fir_design: firwin
    fir_window: hamming
    h_freq: null
    h_trans_bandwidth: auto
    iir_params: null
    l_freq: 1.0
    l_trans_bandwidth: auto
    method: fir
    n_jobs: null
    phase: zero
    picks: null
- name: line_noise
  class: meeg_utils.steps.line_noise:LineNoise
  params:
    adaptive_params: null
    fline: null
    method: zapline-plus
    n_harmonics: null
    n_remove: auto
    nfft: 1024
- name: lowpass
  class: meeg_utils.steps.filtering:Filter
  params:
    filter_length: auto
    fir_design: firwin
    fir_window: hamming
    h_freq: 40.0
    h_trans_bandwidth: auto
    iir_params: null
    l_freq: null
    l_trans_bandwidth: auto
    method: fir
    n_jobs: null
    phase: zero
    picks: null
- name: ica
  class: meeg_utils.steps.ica:ICA
  params:
    decim: null
    exclude_labels: null
    fit_h_freq: 100.0
    fit_l_freq: 1.0
    fit_params: null
    labeler: megnet
    max_iter: auto
    method: infomax
    n_components: 20
    picks: meg
    random_state: 42
    reject_by_annotation: true
    threshold: 0.9

meu.Pipeline.from_yaml(fname) reads it back, and meu run pipeline.yaml --sources bids/ --out derivatives/ runs it from the command line.

MEG and EEG together

ByChannelType runs separate steps on each channel type and writes the results back into one recording, e.g. SSS for the MEG and PREP plus an average reference for the EEG:

combined = meu.Pipeline(
    [
        (
            "by_type",
            S.ByChannelType(
                {
                    "meg": [("sss", S.Maxwell(cross_talk=cross_talk, calibration=calibration))],
                    "eeg": [("bads", S.BadChannels("prep")), ("reference", S.Reference())],
                }
            ),
        ),
        ("highpass", S.Filter(0.1, None)),
    ]
)
combined

#	Step	Configuration (changed parameters)
0	by_type	ByChannelType(branches={'eeg': [('bads', BadChannels(method='prep')), ('reference', Reference())], 'meg': [('sss', Maxwell(calibration=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/sss_cal_mgh.dat'), cross_talk=PosixPath('/home/runner/mne_data/MNE-sample-data/SSS/ct_sparse_mgh.fif')))]})
1	highpass	Filter(h_freq=None, l_freq=0.1)