ENH: Write out raw data and SSP events

docs/source/v1.6.md.inc

@@ -6,6 +6,8 @@
 
 - Added [`regress_artifact`][mne_bids_pipeline._config.regress_artifact] to allow artifact regression (e.g., of MEG reference sensors in KIT systems) (#837 by @larsoner)
 - Chosen `reject` parameters are now saved in the generated HTML reports (#839 by @larsoner)
+- Added saving of clean raw data in addition to epochs (#840 by @larsoner)
+- Added saving of detected blink and cardiac events used to calculate SSP projectors (#840 by @larsoner)
 
 [//]: # (### :warning: Behavior changes)
 
@@ -21,6 +23,7 @@
 - Fix bug where EEG `reject` params were not used for `ch_types = ["meg", "eeg"]` (#839 by @larsoner)
 - Fix bug where implicit `mf_reference_run` could change across invocations of `mne_bids_pipeline`, breaking caching (#839 by @larsoner)
 - Fix bug where `--no-cache` had no effect (#839 by @larsoner)
+- Fix bug where raw, empty-room, and custom noise covariances were errantly calculated on data without ICA or SSP applied (#840 by @larsoner)
 
 ### :medical_symbol: Code health
 

mne_bids_pipeline/_config.py

@@ -1161,14 +1161,14 @@
 ways using the configuration options you can find below.
 """
 
-min_ecg_epochs: int = 5
+min_ecg_epochs: Annotated[int, Ge(1)] = 5
 """
-Minimal number of ECG epochs needed to compute SSP or ICA rejection.
+Minimal number of ECG epochs needed to compute SSP projectors.
 """
 
-min_eog_epochs: int = 5
+min_eog_epochs: Annotated[int, Ge(1)] = 5
 """
-Minimal number of EOG epochs needed to compute SSP or ICA rejection.
+Minimal number of EOG epochs needed to compute SSP projectors.
 """
 
 

mne_bids_pipeline/_config_utils.py

@@ -484,7 +484,7 @@ def get_noise_cov_bids_path(
         task=cfg.task,
         acquisition=cfg.acq,
         run=None,
-        processing=cfg.proc,
+        processing="clean",
         recording=cfg.rec,
         space=cfg.space,
         suffix="cov",
@@ -638,3 +638,24 @@ def _pl(x, *, non_pl="", pl="s"):
     """Determine if plural should be used."""
     len_x = x if isinstance(x, (int, np.generic)) else len(x)
     return non_pl if len_x == 1 else pl
+
+
+def _proj_path(
+    *,
+    cfg: SimpleNamespace,
+    subject: str,
+    session: Optional[str],
+) -> BIDSPath:
+    return BIDSPath(
+        subject=subject,
+        session=session,
+        task=cfg.task,
+        acquisition=cfg.acq,
+        recording=cfg.rec,
+        space=cfg.space,
+        datatype=cfg.datatype,
+        root=cfg.deriv_root,
+        extension=".fif",
+        suffix="proj",
+        check=False,
+    )

mne_bids_pipeline/steps/preprocessing/_05_regress_artifact.py

@@ -82,7 +82,7 @@ def run_regress_artifact(
     model.apply(raw, copy=False)
     if projs:
         raw.add_proj(projs)
-    raw.save(out_files[in_key], overwrite=True)
+    raw.save(out_files[in_key], overwrite=True, split_size=cfg._raw_split_size)
     _update_for_splits(out_files, in_key)
     model.save(out_files["regress"], overwrite=True)
     assert len(in_files) == 0, in_files.keys()

mne_bids_pipeline/steps/preprocessing/_06b_run_ssp.py

@@ -7,13 +7,15 @@
 from typing import Optional
 
 import mne
+import numpy as np
 from mne import compute_proj_epochs, compute_proj_evoked
-from mne.preprocessing import create_ecg_epochs, create_eog_epochs
+from mne.preprocessing import find_ecg_events, find_eog_events
 from mne_bids import BIDSPath
 
 from ..._config_utils import (
     _bids_kwargs,
     _pl,
+    _proj_path,
     get_runs,
     get_sessions,
     get_subjects,
@@ -25,6 +27,11 @@
 from ..._run import _prep_out_files, _update_for_splits, failsafe_run, save_logs
 
 
+def _find_ecg_events(raw: mne.io.Raw, ch_name: Optional[str]) -> np.ndarray:
+    """Wrap find_ecg_events to use the same defaults as create_ecg_events."""
+    return find_ecg_events(raw, ch_name=ch_name, l_freq=8, h_freq=16)[0]
+
+
 def get_input_fnames_run_ssp(
     *,
     cfg: SimpleNamespace,
@@ -69,14 +76,7 @@ def run_ssp(
     # compute SSP on all runs of raw
     raw_fnames = [in_files.pop(f"raw_run-{run}") for run in cfg.runs]
 
-    # when saving proj, use run=None
-    out_files = dict()
-    out_files["proj"] = (
-        raw_fnames[0]
-        .copy()
-        .update(run=None, suffix="proj", split=None, processing=None, check=False)
-    )
-
+    out_files = dict(proj=_proj_path(cfg=cfg, subject=subject, session=session))
     msg = (
         f"Input{_pl(raw_fnames)} ({len(raw_fnames)}): "
         f'{raw_fnames[0].basename}{_pl(raw_fnames, pl=" ...")}'
@@ -93,7 +93,7 @@ def run_ssp(
     projs = dict()
     proj_kinds = ("ecg", "eog")
     rate_names = dict(ecg="heart", eog="blink")
-    epochs_fun = dict(ecg=create_ecg_epochs, eog=create_eog_epochs)
+    events_fun = dict(ecg=_find_ecg_events, eog=find_eog_events)
     minimums = dict(ecg=cfg.min_ecg_epochs, eog=cfg.min_eog_epochs)
     rejects = dict(ecg=cfg.ssp_reject_ecg, eog=cfg.ssp_reject_eog)
     avg = dict(ecg=cfg.ecg_proj_from_average, eog=cfg.eog_proj_from_average)
@@ -111,17 +111,38 @@ def run_ssp(
         projs[kind] = []
         if not any(n_projs[kind].values()):
             continue
-        proj_epochs = epochs_fun[kind](
-            raw,
-            ch_name=ch_name[kind],
-            decim=cfg.epochs_decim,
-        )
-        n_orig = len(proj_epochs.selection)
+        events = events_fun[kind](raw=raw, ch_name=ch_name[kind])
+        n_orig = len(events)
         rate = n_orig / raw.times[-1] * 60
         bpm_msg = f"{rate:5.1f} bpm"
         msg = f"Detected {rate_names[kind]} rate: {bpm_msg}"
         logger.info(**gen_log_kwargs(message=msg))
-        # Enough to start
+        # Enough to create epochs
+        if len(events) < minimums[kind]:
+            msg = (
+                f"No {kind.upper()} projectors computed: got "
+                f"{len(events)} original events < {minimums[kind]} {bpm_msg}"
+            )
+            logger.warning(**gen_log_kwargs(message=msg))
+            continue
+        out_files[f"events_{kind}"] = (
+            out_files["proj"]
+            .copy()
+            .update(suffix=f"{kind}-eve", split=None, check=False, extension=".txt")
+        )
+        mne.write_events(out_files[f"events_{kind}"], events, overwrite=True)
+        proj_epochs = mne.Epochs(
+            raw,
+            events=events,
+            event_id=events[0, 2],
+            tmin=-0.5,
+            tmax=0.5,
+            proj=False,
+            baseline=(None, None),
+            reject_by_annotation=True,
+            preload=True,
+            decim=cfg.epochs_decim,
+        )
         if len(proj_epochs) >= minimums[kind]:
             reject_ = _get_reject(
                 subject=subject,
@@ -134,7 +155,6 @@ def run_ssp(
             proj_epochs.drop_bad(reject=reject_)
         # Still enough after rejection
         if len(proj_epochs) >= minimums[kind]:
-            proj_epochs.apply_baseline((None, None))
             use = proj_epochs.average() if avg[kind] else proj_epochs
             fun = compute_proj_evoked if avg[kind] else compute_proj_epochs
             desc_prefix = (

mne_bids_pipeline/steps/preprocessing/_07_make_epochs.py

@@ -54,7 +54,7 @@ def get_input_fnames_epochs(
         extension=".fif",
         datatype=cfg.datatype,
         root=cfg.deriv_root,
-        processing="filt",
+        processing=cfg.processing,
     ).update(suffix="raw", check=False)
 
     # Generate a list of raw data paths (i.e., paths of individual runs)
@@ -276,7 +276,7 @@ def _get_events(cfg, subject, session):
         acquisition=cfg.acq,
         recording=cfg.rec,
         space=cfg.space,
-        processing="filt",
+        processing=cfg.processing,
         suffix="raw",
         extension=".fif",
         datatype=cfg.datatype,
@@ -322,6 +322,7 @@ def get_config(
         rest_epochs_overlap=config.rest_epochs_overlap,
         _epochs_split_size=config._epochs_split_size,
         runs=get_runs(config=config, subject=subject),
+        processing="filt" if config.regress_artifact is None else "regress",
         **_bids_kwargs(config=config),
     )
     return cfg