WIP: ICA clean raw data

config.py

@@ -569,7 +569,7 @@
 # AUTOMATIC REJECTION OF ARTIFACTS
 # --------------------------------
 
-reject: Optional[dict] = None
+reject: Optional[Union[dict, Literal['auto']]] = None
 """
 The rejection limits to mark epochs as bads.
 This allows to remove strong transient artifacts.
@@ -579,17 +579,17 @@
 
 Pass ``None`` to avoid automated epoch rejection based on amplitude.
 
-???+ note "Note"
-    These numbers tend to vary between subjects.. You might want to consider
-    using the autoreject method by Jas et al. 2018.
-    See https://autoreject.github.io
-
+Pass ``'auto'`` if you want to automate the estimation of the reject
+parameter using AutoReject [Jas et al. 2017] (See https://autoreject.github.io).
+AutoReject is useful as the optimal rejection thresholds tend to vary between
+subjects.
 
 ???+ example "Example"
     ```python
     reject = {'grad': 4000e-13, 'mag': 4e-12, 'eog': 150e-6}
     reject = {'grad': 4000e-13, 'mag': 4e-12, 'eeg': 200e-6}
     reject = None
+    reject = 'auto'
     ```
 """
 
@@ -760,6 +760,17 @@
     ```
 """
 
+fixed_length_epochs_duration: Optional[float] = None
+"""
+Duration of epochs in seconds.
+"""
+
+fixed_length_epochs_overlap: Optional[float] = None
+"""
+Overlap between epochs in seconds. This is used if the task is ``rest``
+and when the annotations do not contain any stimulation or behavior events.
+"""
+
 baseline: Optional[Tuple[Optional[float], Optional[float]]] = (None, 0)
 """
 Specifies which time interval to use for baseline correction of epochs;
@@ -797,6 +808,11 @@
     ```
 """
 
+no_epoching: bool = False
+"""
+Do not perform epoching on the data.
+"""
+
 ###############################################################################
 # ARTIFACT REMOVAL
 # ----------------
@@ -1365,6 +1381,14 @@ def get_t1_from_meeg(bids_path):
            f'are: "FLASH", "T1", and "auto".')
     raise ValueError(msg)
 
+if no_epoching:
+    epochs_tmin = 0
+    fixed_length_epochs_duration = fixed_length_epochs_duration or 1
+    fixed_length_epochs_overlap = fixed_length_epochs_overlap or 0.0
+
+if task == 'rest':
+    no_epoching = True
+
 
 ###############################################################################
 # Helper functions
@@ -1488,6 +1512,9 @@ def get_reject() -> dict:
     if reject is None:
         return dict()
 
+    if reject == 'auto':
+        return dict()
+
     reject_ = reject.copy()  # Avoid clash with global variable.
 
     if ch_types == ['eeg']:

scripts/preprocessing/03-make_epochs.py

@@ -16,7 +16,6 @@
 import mne
 from mne.parallel import parallel_func
 from mne_bids import BIDSPath
-
 import config
 from config import gen_log_message, on_error, failsafe_run
 
@@ -61,7 +60,22 @@ def run_epochs(subject, session=None):
     else:
         raw = mne.concatenate_raws(raw_list)
 
-    events, event_id = mne.events_from_annotations(raw)
+    # Compute events for rest tasks
+    if config.no_epoching:
+        stop = raw.times[-1] - config.fixed_length_epochs_duration
+        # duration = config.epochs_tmax - config.epochs_tmin
+        assert config.epochs_tmin == 0., "epochs_tmin must be 0 for rest"
+        assert config.fixed_length_epochs_overlap is not None, \
+            "epochs_overlap cannot be None for rest"
+        events = mne.make_fixed_length_events(
+            raw, id=3000, start=0,
+            duration=config.fixed_length_epochs_duration,
+            overlap=config.fixed_length_epochs_overlap,
+            stop=stop)
+        event_id = dict(rest=3000)
+    else:  # Events for task runs
+        events, event_id = mne.events_from_annotations(raw)
+
     if "eeg" in config.ch_types:
         projection = True if config.eeg_reference == 'average' else False
         raw.set_eeg_reference(config.eeg_reference, projection=projection)
@@ -88,14 +102,17 @@ def run_epochs(subject, session=None):
 
     # Epoch the data
     msg = (f'Creating epochs with duration: '
-           f'[{config.epochs_tmin}, {config.epochs_tmin}] sec')
+           f'[{config.epochs_tmin}, {config.epochs_tmax}] sec')
     logger.info(gen_log_message(message=msg, step=3, subject=subject,
                                 session=session))
+
+    reject = config.get_reject()
+
     epochs = mne.Epochs(raw, events=events, event_id=event_id,
                         tmin=config.epochs_tmin, tmax=config.epochs_tmax,
                         proj=True, baseline=None,
                         preload=False, decim=config.decim,
-                        reject=config.get_reject(),
+                        reject=reject,
                         reject_tmin=config.reject_tmin,
                         reject_tmax=config.reject_tmax,
                         metadata=metadata,

scripts/preprocessing/04a-run_ica.py

@@ -105,15 +105,20 @@ def make_epochs_for_ica(raw, subject, session):
                             root=config.deriv_root,
                             check=False)
     epochs = mne.read_epochs(epochs_fname)
-    selection = epochs.selection
+
+    # Not sure why it is necessary to recreate epochs when we already have them
+    if config.no_epoching:
+        return epochs
 
     # Now, create new epochs, and only keep the ones we kept in step 3.
     # Because some events present in event_id may disappear entirely from the
     # data, we pass `on_missing='ignore'` to mne.Epochs. Also note that we do
     # not pass the `reject` parameter here.
 
+    selection = epochs.selection
     events, event_id = mne.events_from_annotations(raw)
     events = events[selection]
+
     epochs_ica = mne.Epochs(raw, events=events, event_id=event_id,
                             tmin=epochs.tmin, tmax=epochs.tmax,
                             baseline=None,

scripts/preprocessing/05a-apply_ica.py

@@ -45,17 +45,34 @@ def apply_ica(subject, session):
                              root=config.deriv_root,
                              check=False)
 
-    fname_epo_in = bids_basename.copy().update(suffix='epo', extension='.fif')
-    fname_epo_out = bids_basename.copy().update(
-        suffix='epo', processing='clean', extension='.fif')
+    if config.no_epoching:
+
+        fname_in = {
+            run:
+            bids_basename.copy().update(run=run, processing='filt',
+                                        suffix='raw', check=False,
+                                        extension='.fif')
+            for run in config.get_runs()
+        }
+
+        fname_out = {
+            run:
+            bids_basename.copy().update(run=run, processing='clean',
+                                        suffix='raw', check=False,
+                                        extension='.fif')
+            for run in config.get_runs()
+        }
+
+    else:
+        fname_in = bids_basename.copy().update(suffix='epo', extension='.fif')
+        fname_out = bids_basename.copy().update(
+            suffix='epo', processing='clean', extension='.fif')
+
     fname_ica = bids_basename.copy().update(suffix='ica', extension='.fif')
     fname_ica_components = bids_basename.copy().update(
         processing='ica', suffix='components', extension='.tsv')
 
-    # Load epochs to reject ICA components.
-    epochs = mne.read_epochs(fname_epo_in, preload=True)
-
-    msg = f'Input: {fname_epo_in}, Output: {fname_epo_out}'
+    msg = f'Input: {fname_in}, Output: {fname_out}'
     logger.info(gen_log_message(message=msg, step=5, subject=subject,
                                 session=session))
 
@@ -76,41 +93,69 @@ def apply_ica(subject, session):
                    .loc[tsv_data['status'] == 'bad', 'component']
                    .to_list())
 
-    # Compare ERP/ERF before and after ICA artifact rejection. The evoked
-    # response is calculated across ALL epochs, just like ICA was run on
-    # all epochs, regardless of their respective experimental condition.
-    #
-    # Note that up until now, we haven't actually rejected any ICs from the
-    # epochs.
-    evoked = epochs.average()
-
-    # Plot source time course
-    fig = ica.plot_sources(evoked, show=config.interactive)
-    report.add_figs_to_section(figs=fig,
-                               captions='All ICs - Source time course')
-
-    # Plot original & corrected data
-    fig = ica.plot_overlay(evoked, show=config.interactive)
-    report.add_figs_to_section(figs=fig,
-                               captions=f'Evoked response (across all epochs) '
-                                        f'before and after cleaning via ICA '
-                                        f'({len(ica.exclude)} ICs removed)')
-    report.save(report_fname, overwrite=True, open_browser=False)
-
-    # Now actually reject the components.
-    msg = f'Rejecting ICs: {ica.exclude}'
-    logger.info(gen_log_message(message=msg, step=5, subject=subject,
-                                session=session))
-    epochs_cleaned = ica.apply(epochs.copy())  # Copy b/c works in-place!
-    epochs_cleaned.apply_baseline(config.baseline)
+    if config.no_epoching:
 
-    msg = 'Saving cleaned epochs.'
-    logger.info(gen_log_message(message=msg, step=5, subject=subject,
-                                session=session))
-    epochs_cleaned.save(fname_epo_out, overwrite=True)
+        msg = f'Rejecting ICs: {ica.exclude}'
+        logger.info(gen_log_message(message=msg, step=5, subject=subject,
+                                    session=session))
+
+        # Load and apply ICA to each run separately
+        for run in fname_in:
+
+            raw = mne.io.read_raw_fif(fname_in[run], preload=True)
+            raw_cleaned = ica.apply(raw.copy())  # Not sure if copy needed
+
+            msg = 'Saving cleaned raw.'
+            logger.info(gen_log_message(message=msg, step=5, subject=subject,
+                                        session=session))
+            raw_cleaned.save(fname_out[run], overwrite=True)
+
+            fig = ica.plot_sources(raw, show=config.interactive)
+            report.add_figs_to_section(figs=fig, captions='All ICs - run '
+                                                          f'{run}')
+
+            # Does it make sense to plot_overlay for raw data? I don't think so
+
+        report.save(report_fname, overwrite=True, open_browser=False)
+
+    else:
+
+        # Load epochs to reject ICA components.
+        epochs = mne.read_epochs(fname_in, preload=True)
+
+        # Compare ERP/ERF before and after ICA artifact rejection. The evoked
+        # response is calculated across ALL epochs, just like ICA was run on
+        # all epochs, regardless of their respective experimental condition.
+        #
+        # Note that up until now, we haven't actually rejected any ICs from the
+        # epochs.
+        evoked = epochs.average()
+
+        # Plot source time course
+        fig = ica.plot_sources(evoked, show=config.interactive)
+        report.add_figs_to_section(figs=fig,
+                                   captions='All ICs - Source time course')
+
+        # Plot original & corrected data
+        fig = ica.plot_overlay(evoked, show=config.interactive)
+        report.add_figs_to_section(figs=fig,
+                                   captions=f'Evoked response (across all '
+                                            f'epochs) before and after '
+                                            f'cleaning via ICA '
+                                            f'({len(ica.exclude)} ICs removed)'
+                                   )
+        report.save(report_fname, overwrite=True, open_browser=False)
+
+        epochs_cleaned = ica.apply(epochs.copy())  # Copy b/c works in-place!
+        epochs_cleaned.apply_baseline(config.baseline)
+
+        msg = 'Saving cleaned epochs.'
+        logger.info(gen_log_message(message=msg, step=5, subject=subject,
+                                    session=session))
+        epochs_cleaned.save(fname_out, overwrite=True)
 
-    if config.interactive:
-        epochs_cleaned.plot_image(combine='gfp', sigma=2., cmap="YlGnBu_r")
+        if config.interactive:
+            epochs_cleaned.plot_image(combine='gfp', sigma=2., cmap="YlGnBu_r")
 
 
 def main():