Add support for autoreject (local) before running ICA

docs/source/v1.5.md.inc

@@ -16,9 +16,12 @@ All users are encouraged to update.
 - Reduced logging when reports are created and saved (#799 by @hoechenberger)
 - Added [`"picard-extended_infomax"`][mne_bids_pipeline._config.ica_algorithm] ICA algorithm to perform "extended Infomax"-like ICA decomposition using Picard (#801 by @hoechenberger)
 - Added support for using "local" [`autoreject`](https://autoreject.github.io) to find (and repair) bad channels on a
-  per-epochs basis; this can be enabled by setting [`reject`][mne_bids_pipeline._config.reject] to `"autoreject_local"`.
-  The behavior can further be controlled via the new setting
+  per-epochs basis as the last preprocessing step; this can be enabled by setting [`reject`][mne_bids_pipeline._config.reject]
+  to `"autoreject_local"`. The behavior can further be controlled via the new setting
   [`autoreject_n_interpolate`][mne_bids_pipeline._config.autoreject_n_interpolate]. (#807 by @hoechenberger)
+- Added support for "local" and "global"  [`autoreject`](https://autoreject.github.io) to find (and repair) bad channels on a per-epochs
+  and global basis before submitting them to ICA fitting. This can be enabled by setting [`ica_reject`][mne_bids_pipeline._config.ica_reject]
+  to `"autoreject_local"` or `"autoreject_global"`, respectively. (#810 by @hoechenberger)
 - Website documentation tables can now be sorted (e.g., to find examples that use a specific feature) (#808 by @larsoner)
 
 [//]: # (### :warning: Behavior changes)

mne_bids_pipeline/_config.py

@@ -1233,13 +1233,25 @@
 
 # Rejection based on ICA
 # ~~~~~~~~~~~~~~~~~~~~~~
-ica_reject: Optional[Dict[str, float]] = None
+ica_reject: Optional[
+    Union[Dict[str, float], Literal["autoreject_global", "autoreject_local"]]
+] = None
 """
-Peak-to-peak amplitude limits to exclude epochs from ICA fitting. Epochs exceeding these
-limits will be excluded from ICA fitting. This allows you to remove strong transient
-artifacts, which could negatively affect ICA performance.
+Peak-to-peak amplitude limits to exclude epochs from ICA fitting. This allows you to
+remove strong transient artifacts from the epochs used for fitting ICA, which could
+negatively affect ICA performance. 
+
+The parameter values are the same as for [`reject`][mne_bids_pipeline._config.reject].
 
-This will also be applied to ECG and EOG epochs created during preprocessing.
+If not using `"autoreject_local"` or `"autoreject_global"`, the rejection limits
+will also be applied to the ECG and EOG epochs created to find heart beats and ocular
+artifacts.
+
+???+ info
+    MNE-BIDS-Pipeline will automatically try to detect EOG and ECG artifacts in
+    your data, and remove them. For this to work properly, it is recommended
+    to **not** specify rejection thresholds for EOG and ECG channels here –
+    otherwise, ICA won't be able to "see" these artifacts.
 
 ???+ info 
     This setting is applied only to the epochs that are used for **fitting** ICA. The
@@ -1249,19 +1261,13 @@
     contain large-amplitude artifacts. Those epochs can then be rejected by using
     the [`reject`][mne_bids_pipeline._config.reject] parameter.
 
-The BIDS Pipeline will automatically try to detect EOG and ECG artifacts in
-your data, and remove them. For this to work properly, it is recommended
-to **not** specify rejection thresholds for EOG and ECG channels here –
-otherwise, ICA won't be able to "see" these artifacts.
-
-If `None` (default), do not apply artifact rejection. If a dictionary,
-manually specify peak-to-peak rejection thresholds (see examples).
-
 ???+ example "Example"
     ```python
     ica_reject = {'grad': 10e-10, 'mag': 20e-12, 'eeg': 400e-6}
     ica_reject = {'grad': 15e-10}
-    ica_reject = None  # no rejection
+    ica_reject = None  # no rejection before fitting ICA
+    ica_reject = "autoreject_global"  # find global (per channel type) PTP thresholds before fitting ICA
+    ica_reject = "autoreject_local"  # find local (per channel) thresholds and repair epochs before fitting ICA
     ```
 """
 

mne_bids_pipeline/steps/preprocessing/_06a_run_ica.py

@@ -16,6 +16,7 @@
 
 import pandas as pd
 import numpy as np
+import autoreject
 
 import mne
 from mne.report import Report
@@ -191,7 +192,7 @@ def detect_bad_components(
     ica: mne.preprocessing.ICA,
     ch_names: Optional[List[str]],
     subject: str,
-    session: str,
+    session: Optional[str],
 ) -> Tuple[List[int], np.ndarray]:
     artifact = which.upper()
     msg = f"Performing automated {artifact} artifact detection …"
@@ -390,33 +391,51 @@ def run_ica(
     epochs_ecg = ecg_epochs_all_runs
     epochs_eog = eog_epochs_all_runs
 
-    del epochs_all_runs, eog_epochs_all_runs, ecg_epochs_all_runs
+    del epochs_all_runs, eog_epochs_all_runs, ecg_epochs_all_runs, run
 
     # Set an EEG reference
     if "eeg" in cfg.ch_types:
         projection = True if cfg.eeg_reference == "average" else False
         epochs.set_eeg_reference(cfg.eeg_reference, projection=projection)
 
-    # Reject epochs based on peak-to-peak rejection thresholds
-    ica_reject = _get_reject(
-        subject=subject,
-        session=session,
-        reject=cfg.ica_reject,
-        ch_types=cfg.ch_types,
-        param="ica_reject",
-    )
+    if cfg.ica_reject == "autoreject_local":
+        msg = "Using autoreject to find and repair bad epochs before fitting ICA"
+        logger.info(**gen_log_kwargs(message=msg))
 
-    msg = f"Using PTP rejection thresholds: {ica_reject}"
-    logger.info(**gen_log_kwargs(message=msg))
+        ar = autoreject.AutoReject(
+            n_interpolate=cfg.autoreject_n_interpolate,
+            random_state=cfg.random_state,
+            n_jobs=exec_params.n_jobs,
+            verbose=False,
+        )
+        ar.fit(epochs)
+        epochs, reject_log = ar.transform(epochs, return_log=True)
+        msg = (
+            f"autoreject marked {reject_log.bad_epochs.sum()} epochs as bad "
+            f"(cross-validated n_interpolate limit: {ar.n_interpolate_})"
+        )
+        logger.info(**gen_log_kwargs(message=msg))
+    else:
+        # Reject epochs based on peak-to-peak rejection thresholds
+        ica_reject = _get_reject(
+            subject=subject,
+            session=session,
+            reject=cfg.ica_reject,
+            ch_types=cfg.ch_types,
+            param="ica_reject",
+        )
+
+        msg = f"Using PTP rejection thresholds: {ica_reject}"
+        logger.info(**gen_log_kwargs(message=msg))
 
-    epochs.drop_bad(reject=ica_reject)
-    if epochs_eog is not None:
-        epochs_eog.drop_bad(reject=ica_reject)
-    if epochs_ecg is not None:
-        epochs_ecg.drop_bad(reject=ica_reject)
+        epochs.drop_bad(reject=ica_reject)
+        if epochs_eog is not None:
+            epochs_eog.drop_bad(reject=ica_reject)
+        if epochs_ecg is not None:
+            epochs_ecg.drop_bad(reject=ica_reject)
 
     # Now actually perform ICA.
-    msg = "Calculating ICA solution."
+    msg = f"Calculating ICA solution using method: {cfg.ica_algorithm}."
     logger.info(**gen_log_kwargs(message=msg))
     ica = fit_ica(cfg=cfg, epochs=epochs, subject=subject, session=session)
 
@@ -497,6 +516,24 @@ def run_ica(
     eog_scores = None if len(eog_scores) == 0 else eog_scores
 
     with _agg_backend():
+        if cfg.ica_reject == "autoreject_local":
+            caption = (
+                f"Autoreject was run to produce cleaner epochs before fitting ICA. "
+                f"{reject_log.bad_epochs.sum()} epochs were rejected because more than "
+                f"{ar.n_interpolate_} channels were bad (cross-validated n_interpolate "
+                f"limit; excluding globally bad and non-data channels, shown in white)."
+            )
+            report.add_figure(
+                fig=reject_log.plot(
+                    orientation="horizontal", aspect="auto", show=False
+                ),
+                title="Epochs: Autoreject cleaning",
+                caption=caption,
+                tags=("ica", "epochs", "autoreject"),
+                replace=True,
+            )
+            del caption
+
         report.add_epochs(
             epochs=epochs,
             title="Epochs used for ICA fitting",
@@ -536,7 +573,7 @@ def run_ica(
 def get_config(
     *,
     config: SimpleNamespace,
-    subject: Optional[str] = None,
+    subject: str,
     session: Optional[str] = None,
 ) -> SimpleNamespace:
     cfg = SimpleNamespace(
@@ -551,6 +588,7 @@ def get_config(
         ica_reject=config.ica_reject,
         ica_eog_threshold=config.ica_eog_threshold,
         ica_ctps_ecg_threshold=config.ica_ctps_ecg_threshold,
+        autoreject_n_interpolate=config.autoreject_n_interpolate,
         random_state=config.random_state,
         ch_types=config.ch_types,
         l_freq=config.l_freq,

mne_bids_pipeline/steps/preprocessing/_08_ptp_reject.py

@@ -116,7 +116,7 @@ def drop_ptp(
             epochs=epochs,
         )
 
-        if cfg.spatial_filter == "ica":
+        if cfg.spatial_filter == "ica" and cfg.ica_reject != "autoreject_local":
             ica_reject = _get_reject(
                 subject=subject,
                 session=session,
@@ -156,9 +156,13 @@ def drop_ptp(
         )
         logger.warning(**gen_log_kwargs(message=msg))
     elif n_epochs_after_reject == 0:
+        rejection_type = (
+            cfg.reject
+            if cfg.reject in ["autoreject_global", "autoreject_local"]
+            else "PTP-based"
+        )
         raise RuntimeError(
-            "No epochs remaining after peak-to-peak-based "
-            "rejection. Cannot continue."
+            f"No epochs remaining after {rejection_type} rejection. Cannot continue."
         )
 
     msg = "Saving cleaned, baseline-corrected epochs …"

mne_bids_pipeline/tests/configs/config_ERP_CORE.py

@@ -75,17 +75,22 @@
     spatial_filter = None
     reject = "autoreject_local"
     autoreject_n_interpolate = [2, 4]
+elif task == "P3":  # test autoreject local with ICA (before and after ICA cleaning)
+    spatial_filter = "ica"
+    ica_reject = "autoreject_local"
+    reject = "autoreject_local"
+    autoreject_n_interpolate = [4]  # Only for testing!
 else:
+    spatial_filter = "ica"
     ica_reject = dict(eeg=350e-6, eog=500e-6)
     reject = "autoreject_global"
 
-    spatial_filter = "ica"
-    ica_max_iterations = 1000
-    ica_eog_threshold = 2
-    ica_decim = 2  # speed up ICA fitting
+# These settings are only used for the cases where spatial_filter is not None
+ica_max_iterations = 1000
+ica_eog_threshold = 2
+ica_decim = 2  # speed up ICA fitting
 
 run_source_estimation = False
-
 on_rename_missing_events = "ignore"
 
 parallel_backend = "dask"