Add the code from PR "Add support for MNE-ICALabel #812"

.circleci/config.yml

@@ -1081,7 +1081,7 @@ jobs:
       - run:
           name: Install dependencies
           command: |
-            pip install -ve .[docs]
+            pip install -ve .[docs] onnxruntime
       - run:
           name: Build documentation
           command: |

docs/source/dev.md.inc

@@ -2,6 +2,7 @@
 
 ### :new: New features & enhancements
 
+- Support for using `mne-icalabel` to automatically label ICA components. This requires the `mne-icalabel` package to be installed. (#1018 and #812 by @jschepers, @behinger, @hoechenberger, and @larsoner)
 - It is now possible to use separate MRIs for each session within a subject, as in longitudinal studies. This is achieved by creating separate "subject" folders for each subject-session combination, with the naming convention `sub-XXX_ses-YYY`, in the freesurfer `SUBJECTS_DIR`. (#987 by @drammock)
 - New config option [`allow_missing_sessions`][mne_bids_pipeline._config.allow_missing_sessions] allows to continue when not all sessions are present for all subjects. (#1000 by @drammock)
 - New config option [`mf_extra_kws`][mne_bids_pipeline._config.mf_extra_kws] passes additional keyword arguments to `mne.preprocessing.maxwell_filter`. (#1038 by @drammock)

mne_bids_pipeline/_config.py

@@ -13,6 +13,7 @@
     DigMontageType,
     FloatArrayLike,
     PathLike,
+    UniqueSequence,
 )
 
 # %%
@@ -1430,6 +1431,11 @@
     us so we can discuss.
 """
 
+ica_h_freq: float | None = None
+"""
+The cutoff frequency of the low-pass filter to apply before running ICA.
+"""
+
 ica_max_iterations: int = 500
 """
 Maximum number of iterations to decompose the data into independent
@@ -1476,19 +1482,70 @@
 `1` or `None` to not perform any decimation.
 """
 
+ica_use_ecg_detection: bool = True
+"""
+Whether to use the MNE ECG detection on the ICA components.
+"""
+
 ica_ecg_threshold: float = 0.1
 """
 The cross-trial phase statistics (CTPS) threshold parameter used for detecting
 ECG-related ICs.
 """
 
+ica_use_eog_detection: bool = True
+"""
+Whether to use the MNE EOG detection on the ICA components.
+"""
+
 ica_eog_threshold: float = 3.0
 """
 The threshold to use during automated EOG classification. Lower values mean
 that more ICs will be identified as EOG-related. If too low, the
 false-alarm rate increases dramatically.
 """
 
+ica_use_icalabel: bool = False
+"""
+Whether to use MNE-ICALabel to automatically label ICA components. Only available for
+EEG data.
+
+!!! info
+    Using MNE-ICALabel mandates that you also set:
+    ```python
+    eeg_reference = "average"
+    ica_l_freq = 1
+    ica_h_freq = 100
+    ```
+
+!!! info
+    Using this requires `mne-icalabel` package, which in turn will require you to
+    install a suitable runtime (`onnxruntime` or `pytorch`).
+"""
+
+ica_icalabel_include: Annotated[
+    UniqueSequence[
+        Literal[
+            "brain",
+            "muscle artifact",
+            "eye blink",
+            "heart beat",
+            "line noise",
+            "channel noise",
+            "other",
+        ]
+    ],
+    Len(1, 7),
+] = ("brain", "other")
+"""
+Which independent components (ICs) to keep based on the labels given by ICLabel.
+Possible labels are:
+
+```
+["brain", "muscle artifact", "eye blink", "heart beat", "line noise", "channel noise", "other"]
+```
+"""  # noqa: E501
+
 # ### Amplitude-based artifact rejection
 #
 # ???+ info "Good Practice / Advice"

mne_bids_pipeline/_config_import.py

@@ -426,6 +426,27 @@ def _check_config(config: SimpleNamespace, config_path: PathLike | None) -> None
                 f"but got shape {destination.shape}"
             )
 
+    # From: https://github.com/mne-tools/mne-bids-pipeline/pull/812
+    # MNE-ICALabel
+    if config.ica_use_icalabel:
+        pre = "When using MNE-ICALabel, you must set"
+        if config.ica_l_freq != 1.0 or config.ica_h_freq != 100.0:
+            raise ValueError(
+                f"{pre} ica_l_freq=1 and h_freq=100, "
+                f"but got: ica_l_freq={config.ica_l_freq} and "
+                f"ica_h_freq={config.ica_h_freq}"
+            )
+        if config.eeg_reference != "average":
+            raise ValueError(
+                f'{pre} eeg_reference="average", but got: '
+                f"eeg_reference={config.eeg_reference}"
+            )
+        if config.ica_algorithm not in ("picard-extended_infomax", "extended_infomax"):
+            raise ValueError(
+                f'{pre} ica_algorithm="picard-extended_infomax" or "extended_infomax", '
+                f"but got: ica_algorithm={repr(config.ica_algorithm)}"
+            )
+
 
 def _default_factory(key: str, val: Any) -> Any:
     # convert a default to a default factory if needed, having an explicit
@@ -435,6 +456,7 @@ def _default_factory(key: str, val: Any) -> Any:
         {"custom": (8, 24.0, 40)},  # decoding_csp_freqs
         ["evoked"],  # inverse_targets
         [4, 8, 16],  # autoreject_n_interpolate
+        ("brain", "other"),  # ica_icalabel_include
     ]
 
     def default_factory() -> Any:

mne_bids_pipeline/steps/preprocessing/_06a1_fit_ica.py

@@ -79,6 +79,14 @@ def run_ica(
     """Run ICA."""
     import matplotlib.pyplot as plt
 
+    if cfg.ica_use_icalabel:
+        # The ICALabel network was trained on extended-Infomax ICA decompositions fit
+        # on data flltered between 1 and 100 Hz.
+        assert cfg.ica_algorithm in ["picard-extended_infomax", "extended_infomax"]
+        assert cfg.ica_l_freq == 1.0
+        assert cfg.ica_h_freq == 100.0
+        assert cfg.eeg_reference == "average"
+
     raw_fnames = [in_files.pop(f"raw_run-{run}") for run in cfg.runs]
     out_files = dict()
     bids_basename = raw_fnames[0].copy().update(processing=None, split=None, run=None)
@@ -105,20 +113,38 @@ def run_ica(
         # Sanity check – make sure we're using the correct data!
         if cfg.raw_resample_sfreq is not None:
             assert np.allclose(raw.info["sfreq"], cfg.raw_resample_sfreq)
-        if cfg.l_freq is not None:
-            assert np.allclose(raw.info["highpass"], cfg.l_freq)
 
         if idx == 0:
-            if cfg.ica_l_freq is None:
+            # We have to do some gymnastics here to permit for example 128 Hz-sampled
+            # data to be used with mne-icalabel, which wants data low-pass filtered
+            # at 100 Hz
+            h_freq = cfg.ica_h_freq
+            nyq = raw.info["sfreq"] / 2.0
+            if h_freq is not None and h_freq >= nyq:
                 msg = (
-                    f"Not applying high-pass filter (data is already filtered, "
-                    f"cutoff: {raw.info['highpass']} Hz)."
+                    f"Low-pass filter cutoff {h_freq} Hz is higher "
+                    f"than Nyquist {nyq} Hz"
                 )
+                if cfg.ica_use_icalabel:
+                    msg += ", setting to None for compatibility with MNE-ICALabel."
+                    logger.warning(**gen_log_kwargs(message=msg))
+                    h_freq = None
+                else:
+                    raise ValueError(msg)
+            msg = ""
+            if cfg.ica_l_freq is not None and h_freq is not None:
+                msg = (
+                    f"Applying band-pass filter with {cfg.ica_l_freq}-{h_freq} "
+                    "Hz cutoffs"
+                )
+            elif cfg.ica_l_freq is not None:
+                msg = f"Applying high-pass filter with {cfg.ica_l_freq} Hz cutoff"
+            elif h_freq is not None:
+                msg = f"Applying low-pass filter with {h_freq} Hz cutoff"
+            if cfg.ica_l_freq is not None or h_freq is not None:
                 logger.info(**gen_log_kwargs(message=msg))
-            else:
-                msg = f"Applying high-pass filter with {cfg.ica_l_freq} Hz cutoff …"
-                logger.info(**gen_log_kwargs(message=msg))
-                raw.filter(l_freq=cfg.ica_l_freq, h_freq=None, n_jobs=1)
+                raw.filter(l_freq=cfg.ica_l_freq, h_freq=h_freq, n_jobs=1)
+            del nyq, h_freq
 
         # Only keep the subset of the mapping that applies to the current run
         event_id = event_name_to_code_map.copy()
@@ -166,7 +192,13 @@ def run_ica(
     # Set an EEG reference
     if "eeg" in cfg.ch_types:
         projection = True if cfg.eeg_reference == "average" else False
+        if not projection:
+            msg = "Applying average reference to EEG epochs used for ICA fitting."
+            logger.info(**gen_log_kwargs(message=msg))
+
         epochs.set_eeg_reference(cfg.eeg_reference, projection=projection)
+        if cfg.ica_use_icalabel:
+            epochs.apply_proj()  # Apply the reference projection
 
     ar_reject_log = ar_n_interpolate_ = None
     if cfg.ica_reject == "autoreject_local":
@@ -333,16 +365,17 @@ def get_config(
         conditions=config.conditions,
         runs=get_runs(config=config, subject=subject),
         task_is_rest=config.task_is_rest,
+        ica_h_freq=config.ica_h_freq,
         ica_l_freq=config.ica_l_freq,
         ica_algorithm=config.ica_algorithm,
         ica_n_components=config.ica_n_components,
         ica_max_iterations=config.ica_max_iterations,
         ica_decim=config.ica_decim,
         ica_reject=config.ica_reject,
+        ica_use_icalabel=config.ica_use_icalabel,
         autoreject_n_interpolate=config.autoreject_n_interpolate,
         random_state=config.random_state,
         ch_types=config.ch_types,
-        l_freq=config.l_freq,
         epochs_decim=config.epochs_decim,
         raw_resample_sfreq=config.raw_resample_sfreq,
         event_repeated=config.event_repeated,