Stage 3-b: ICA artefact decomposition

Author: Hadi Zaatiti <hadi.zaatiti@nyu.edu>

What & why

Stage 3-b of the KIT general pipeline decomposes every stage-2 recording with independent component analysis (ICA) so that stereotyped artefacts (heartbeat, eye blinks and movements, muscle) can be removed before the data is filtered, epoched and averaged. It follows the recipe of the COGNESTIC course (08_EEG_MEG/A_1_Artifact_Correction.ipynb, itself based on the MNE-Python ICA tutorial):

  1. fit the ICA on a copy of the data high-pass filtered at 1 Hz (slow drifts are detrimental to the decomposition), with ICA(n_components=0.90, method="fastica", max_iter="auto", random_state=97);

  2. look at the solution: plot_sources, plot_components, plot_properties, plot_overlay, the explained variance;

  3. let the data propose artefact components: find_bads_ecg on a synthetic ECG (correlation, threshold auto), find_bads_eog when an EOG or proxy channel exists, with plot_scores and the sources of the heartbeat-locked average as evidence;

  4. set ica.exclude and apply the solution to the unfiltered data.

As in stage 3-a, no component is removed without a human: the automatic matches are proposals, pre-marked in MNE’s source browser next to the figures above; the reviewer confirms, rejects or extends them and closes the window, and that is the decision. Stage 4 then applies the reviewed solution to the unfiltered recording right before filtering, so no second copy of the raw data is written.

The decomposition uses the 207 KIT head gradiometers (MNE type mag) minus the bad channels decided by stage 3-a. The reference magnetometers (ref_meg, including MEG 092) never take part.

What you need to provide

  • The stage-2 outputs (kit2fiff_summary.csv and the FIFs).

  • The stage 3-a records: a recording without one is decomposed with the FIF’s own info['bads']; a recording whose candidates are still pending is decomposed with only its human-provided bad channels (the stage warns in both cases).

  • A desktop session with a Qt binding and mne-qt-browser for the review; scikit-learn for FastICA (both in the neurowaves recipe). Fitting needs memory: the 207 channels of a 42-minute KidLang recording take about 8 GB while they are loaded and filtered.

Configuration file

The ica block of the YAML config (defaults shown; the KidLang settings are the same):

ica:
  enabled: true
  fit:
    l_freq: 1.0            # high-pass of the copy the ICA is fitted on (COGNESTIC)
    h_freq: 40             # low-pass of that copy (null = none); see step 1
    n_components: 0.9      # fraction of variance to keep (< 1) or a number of components
    method: fastica
    max_iter: auto
    random_state: 97       # fixed seed: reruns give identical components
    decim: 5               # use every 5th sample when fitting (speed)
    max_duration_s: 0      # fit on the first N seconds only (0 = whole recording)
    reject_t: null         # peak-to-peak rejection while fitting, in T (e.g. 5.0e-12)
  ecg:
    enabled: true
    method: correlation    # as in COGNESTIC ('ctps' is MNE's other option)
    measure: zscore        # 'zscore' = outliers among the components; 'correlation' = absolute value
    threshold: auto        # 3 standard deviations with zscore; e.g. 0.9 with correlation
    ch_name: null          # a real ECG channel name, if one was recorded
  eog:
    enabled: false         # KIT records no EOG; set the channels to use (e.g. MISC lines)
    channels: []
    threshold: 3.0
  muscle:
    enabled: false         # find_bads_muscle (MNE >= 1.1), opt-in
    threshold: 0.5
  manual_exclude: {}       # components always removed, per subject: {Y0312: [3, 7]}
  review:
    n_cuts: 8              # cuts of the high-passed data shown in the browser
    cut_duration_s: 30
    skip_edges_s: 5
    resample_hz: 150
    window_s: 10           # seconds per browser page
    antialiasing: false
    use_opengl: false
    whole_recording: false # true = browse the entire high-passed recording (slow to open)
    max_properties: 8      # property figures opened for the proposed components
  save_plots: true
  save_cleaned: false      # also write <entities>_desc-icaclean_raw.fif (large)

manual_exclude is keyed by subject label with or without sub- and applies to every recording of that subject. reports.scripts.3-b-kit_ica enables or disables the report refresh after the stage.

How to run

From the repository root with the neurowaves environment activated:

python -m pipeline.mne_pipelines.kit_general_pipelines.3-b-kit_ica \
    -c pipeline/mne_pipelines/kit_general_pipelines/pipeline_config_files/kidlang_config.yml

Options:

Flag

Effect

-c, --config PATH

YAML config (always pass the dataset’s file).

(no flag)

Interactive run: fit (or reuse the saved solution), propose, open the figures and the source browser for every recording without a saved review. While you review one recording, the next one is fitted in a background process, so after the first there is no wait.

--no-gui

Fit, propose and save the figures only; exclude nothing (review_status: pending). Needs no display; use it to precompute the solutions, the GUI run then opens straight into the review.

--review-all

Redo the review of reviewed recordings, starting from the previous decision.

--reset-review

Discard the saved reviews; recordings become pending again.

--overwrite

Refit even when a solution with the same input, settings and bad channels exists (reviews are kept but worth redoing).

--no-plots

Skip the saved figures (the live figures of the review still open).

--save-cleaned

Also write <entities>_desc-icaclean_raw.fif, the unfiltered recording with the reviewed components removed (as large as the input). Stage 4 does not need it.

A solution is reused when the FIF (path, size, modification time), the fit/ecg/eog/muscle settings and the bad-channel list are unchanged; changing any of them refits or re-proposes on the next run. The VS Code entries (Kidlang) ICA: fit + review (GUI), (Kidlang) ICA: fit only (no GUI) and (Kidlang) ICA: redo all reviews (GUI) run the three variants.

Step by step

Step 1: fit on a high-passed copy

The 207 gradiometers minus the bad channels (plus the configured EOG or proxy channels, if any) are loaded, band-pass filtered between fit.l_freq and fit.h_freq (1 to 40 Hz), and the ICA is fitted on them with decim (every n-th sample), optionally on the first max_duration_s seconds and with peak-to-peak rejection reject_t. COGNESTIC only high-passes at 1 Hz, but its sample data was already low-passed at 40 Hz; KIT raw data is broadband at 1000 Hz, and without the low-pass the sensor noise above 40 Hz costs dozens of extra components (measured on a 42-minute KidLang recording: 90 components for 10 minutes of broadband data and a fit that did not finish in 17 minutes on the whole recording, against 23 components and 43 s with the 1 to 40 Hz band). decim subsamples without an anti-alias filter, which is another reason to low-pass first. n_components: 0.9 keeps as many components as needed to explain 90 % of the variance; the fixed random_state makes the solution reproducible. The explained-variance ratio of the solution and of each component is computed on a 60 s crop and stored. Loading and filtering one recording takes about a minute and 8 GB of memory.

Step 2: the figures

Saved under figures/ and, during the review, opened live:

  • plot_components: the topographies of all components (several pages of 20); live, clicking a topography opens its properties;

  • plot_scores: the match score of every component for each proposer, proposed components in red;

  • plot_properties of each proposed component (up to review.max_properties): topography, source segments, spectrum, variance;

  • plot_sources of the ECG-evoked response (the sources averaged around the detected heartbeats): heartbeat components show a clear QRS shape;

  • plot_overlay: the data before (red) and after (black) removing the proposals, and after the review the same with the final exclusions (overlay-final).

Step 3: proposals

  • ECG: find_bads_ecg correlates every component with the ECG signal. KIT records no ECG channel, so MNE builds a synthetic one from the average of the MEG channels (ch_name: null), which contains the cardiac field clearly enough. With measure: zscore (default) the correlations are z-scored and threshold: auto proposes the components more than 3 standard deviations from the others; with measure: correlation the threshold is an absolute correlation (for example 0.9). The heartbeat often appears as two components (the cardiac field is a travelling dipole). On KIT data the synthetic ECG is a weak reference (the cardiac field is a small part of the MEG average), so a heartbeat component may not stand out statistically: the stage therefore also prints and stores the three highest ECG correlations as hints, and the ECG-evoked sources figure shows which components carry the QRS shape. measure: correlation with threshold: 0.5 proposes more liberally.

  • EOG: find_bads_eog against the configured channels, only when eog.enabled and eog.channels are set. Without an EOG channel the blink components are not proposed, but they are easy to spot in the review: a frontal, dipolar topography and a blink-shaped time course.

  • Muscle: find_bads_muscle when enabled (rarely needed for MEG).

  • manual_exclude components are always removed and pre-marked.

The union is the list of candidates; they are proposals and nothing more.

Step 4: human review and decision

The stage prints the proposals with their scores, opens the figures of step 2 and MNE’s source browser (plot_sources, Qt backend from mne-qt-browser) on a review copy of the high-passed data: review.n_cuts cuts of review.cut_duration_s seconds at evenly spaced positions (4 minutes in total by default), resampled to review.resample_hz, with an annotation giving the position of each cut in the recording (whole_recording: true shows everything, slowly). The proposals and the manual components are pre-marked.

In the browser:

  • click a component name to toggle its exclusion (excluded = grey);

  • right-click a name to open its properties; arrow keys scroll, Page Up/Page Down change the component page;

  • close the window to decide. The components excluded at that moment, plus manual_exclude, are the final exclusions (review_status: reviewed); the solution is saved with them (ica.save), the scores table is updated and the final overlay is drawn.

With --no-gui the review does not run: review_status stays pending, exclude holds only the manual components, and stage 4 warns that proposals are waiting. Proposals that appear after a review (changed proposer settings) are put to the reviewer on the next interactive run; a headless run records them as unseen_candidates.

Outputs

<bids_root>/derivatives/3-b-kit_ica/
  config_<timestamp>.yml
  ica_summary.csv                       # one row per recording, whole cohort, rewritten after every recording
  sub-<id>/
    ica_log.txt
    [ses-<id>/]
      sub-<id>..._desc-ica_ica.fif        # the ICA solution incl. exclude (mne.preprocessing.read_ica)
      sub-<id>..._desc-ica.json           # proposals, scores, review, settings, provenance
      sub-<id>..._desc-ica_scores.tsv     # one row per component
      sub-<id>..._desc-icareview_raw.fif  # review data, kept only until the review is done
      sub-<id>..._desc-icaclean_raw.fif   # only with --save-cleaned
      figures/
        sub-<id>..._desc-ica_components-<n>.png
        sub-<id>..._desc-ica_scores-ecg.png        (and -eog, -muscle when enabled)
        sub-<id>..._desc-ica_properties-ICA<nnn>.png
        sub-<id>..._desc-ica_ecg_sources.png
        sub-<id>..._desc-ica_overlay.png
        sub-<id>..._desc-ica_overlay-final.png

The JSON record holds exclude and review_status (what stage 4 uses), candidates, proposals (per proposer: indices, scores, method, threshold), manual_exclude, unseen_candidates, review (pre-marked, decision, timestamp, browser backend, data shown, started_from), n_components, explained_variance, variance_per_component, channels_fit, bad_channels (names, status, record used), fit_settings, detection_settings, input (FIF identity), ica_fif, scores_tsv, figures, fit_duration_s and mne_version.

The scores TSV has, per component: variance_ratio, ecg_score, eog_score, muscle_score, the *_candidate flags, candidate, manual and excluded.

ica_summary.csv has one row per recording: the BIDS entities, in_fif, n_channels_fit, bad_channels, bad_channel_status, fit_duration_s, n_components, explained_variance, the *_candidates lists, candidates, manual_exclude, review_status, excluded, n_excluded, ica_fif, record_json, scores_tsv, cleaned_fif, status, error and mne_version.

How later stages use the result

utilities.final_ica_for_recording(bids_root, entities, source_fif) returns (ica_fif, exclude, status, record_path) with status ok, pending_review (only the manual components in exclude), stale_input (record made from another version of the FIF), missing_solution or no_record. Stage 4 calls it after resolving the bad channels, warns on pending_review and stale_input, applies ica.apply with the exclusions before the notch and band-pass filters (COGNESTIC applies the solution fitted on filtered data to the original data; filtering afterwards is a linear operation and commutes with it), stores ica_excluded, ica_source, ica_record and ica_fif in the provenance sidecar and the ica_excluded / ica_source columns of filter_summary.csv, and reports an existing filtered FIF whose exclusions differ as existing_output_unverified (rerun stage 4 with --overwrite). Nothing is applied when the record excludes no component.

Report

The section 3-b. ICA artefact decomposition of the Quarto report explains the four steps with the parameter values of the run, shows the cohort table (reviewed and pending counts, components removed per recording) and, per recording, a collapsed block with the proposals, the decision, the scores of the proposed and excluded components and the figures. The report is refreshed in the background after every recording.

Edge cases & safeguards

  • Short recordings or no channels: a recording shorter than 10 s, or one with fewer than two gradiometers left after the bad channels, is refused with a clear error (status: error); the review cuts of a short recording shrink to what exists.

  • FIF changed since the fit (stage 2 rerun): the solution is refitted on the next run; stage 4 reports stale_input until then.

  • Bad-channel list changed (stage 3-a rerun): the solution is refitted, since the channel set changed; the previous review is discarded with it.

  • FastICA does not converge: MNE warns; raise fit.max_iter or lower n_components.

  • Interrupted run: finished recordings keep their solution, record and summary row; a background fit of the next recording is terminated and redone on the next run.

  • No Qt binding / no display: the interactive run stops at start; use --no-gui.

Troubleshooting

  • Fitting is slow or yields very many components: make sure fit.h_freq is set (40 Hz); broadband data needs several times more components and the fit can run for tens of minutes. fit.max_duration_s (for example 600 s) fits on the first 10 minutes only. The console prints the fit time and the number of components.

  • No ECG proposal although a heartbeat component is visible: with measure: zscore a component is proposed only when it stands out from the others; try measure: correlation with threshold: 0.9, or just mark the component in the review (the ECG-evoked sources figure shows which ones carry the QRS shape).

  • Too many proposals: raise the threshold, or ignore them in the review; only what is marked at close is removed.

  • The browser is slow: same remedies as for stage 3-a (review.antialiasing: false, review.use_opengl: true on a machine with a GPU, mne-qt-browser installed, no debugger).