WIP: ICA clean raw data

config.py

@@ -579,12 +579,6 @@
 
 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
-
-
 ???+ example "Example"
     ```python
     reject = {'grad': 4000e-13, 'mag': 4e-12, 'eog': 150e-6}
@@ -760,6 +754,18 @@
     ```
 """
 
+fixed_length_epochs_duration: Optional[float] = 1
+"""
+Duration of epochs in seconds. Used if `no_epoching` was set, in order to
+perform ICA cleaning.
+"""
+
+fixed_length_epochs_overlap: Optional[float] = 0
+"""
+Overlap between epochs in seconds. Used if `no_epoching` was set, in order to
+perform ICA cleaning.
+"""
+
 baseline: Optional[Tuple[Optional[float], Optional[float]]] = (None, 0)
 """
 Specifies which time interval to use for baseline correction of epochs;
@@ -797,6 +803,11 @@
     ```
 """
 
+no_epoching: bool = False
+"""
+Do not perform epoching on the data.
+"""
+
 ###############################################################################
 # ARTIFACT REMOVAL
 # ----------------
@@ -1365,6 +1376,12 @@ def get_t1_from_meeg(bids_path):
            f'are: "FLASH", "T1", and "auto".')
     raise ValueError(msg)
 
+if no_epoching:
+    epochs_tmin = 0
+
+if task == 'rest':
+    no_epoching = True
+
 
 ###############################################################################
 # Helper functions

run.py

@@ -184,7 +184,9 @@ def process(config: PathLike,
         script_paths = [*SCRIPT_PATHS['init'], *script_paths]
 
     for script_path in script_paths:
-        step_name = script_path.name.replace('.py', '')[3:]
+        step_name = script_path.name.replace('.py', '')
+        if '-' in step_name:
+            step_name = '-'.join(step_name.split('-')[1:])
         logger.info(f'Now running: {step_name}')
         _run_script(script_path, config, root_dir, subject, session, task, run)
         logger.info(f'Successfully finished running: {step_name}')

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,28 @@ def make_epochs_for_ica(raw, subject, session):
                             root=config.deriv_root,
                             check=False)
     epochs = mne.read_epochs(epochs_fname)
-    selection = epochs.selection
+
+    # Should be factorized as it is the same as in 03-make_epochs.py:64,
+    # but not obvious how
+    if config.no_epoching:
+        stop = raw.times[-1] - config.fixed_length_epochs_duration
+        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, event_id = mne.events_from_annotations(raw)
 
     # 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.
 
-    events, event_id = mne.events_from_annotations(raw)
+    selection = epochs.selection
     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)
 
-    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)
+    if config.no_epoching:
+
+        # 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():