diff --git a/pyproject.toml b/pyproject.toml
index 6d2ab73..ada8a8f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -24,7 +24,7 @@ dependencies = [
   "nitransforms == 23.0.1",
   "niworkflows @ git+https://github.com/nipreps/niworkflows.git@master",
   "pybids >= 0.15.6",
-  "rapidtide == 2.9.5.2",
+  "rapidtide == 2.9.8",
   "sdcflows @ git+https://github.com/nipreps/sdcflows.git@master",
   "smriprep @ git+https://github.com/nipreps/smriprep.git@master",
   "typer",
diff --git a/src/fmripost_rapidtide/config.py b/src/fmripost_rapidtide/config.py
index e28d949..0526b54 100644
--- a/src/fmripost_rapidtide/config.py
+++ b/src/fmripost_rapidtide/config.py
@@ -613,6 +613,8 @@ class workflow(_Config):
     """Generate HCP Grayordinates, accepts either ``'91k'`` (default) or ``'170k'``."""
     dummy_scans = None
     """Set a number of initial scans to be considered nonsteady states."""
+    average_over_runs = True
+    """Whether to average lag maps over runs or not."""
 
 
 class loggers:
diff --git a/src/fmripost_rapidtide/data/io_spec.json b/src/fmripost_rapidtide/data/io_spec.json
index 46dfd29..b272f97 100644
--- a/src/fmripost_rapidtide/data/io_spec.json
+++ b/src/fmripost_rapidtide/data/io_spec.json
@@ -16,15 +16,15 @@
             }
         },
         "derivatives": {
-            "bold_mni152nlin6asym": {
+            "bold_boldref": {
                 "datatype": "func",
                 "echo": null,
                 "part": [
                     "mag",
                     null
                 ],
-                "res": "2",
-                "space": "MNI152NLin6Asym",
+                "res": null,
+                "space": null,
                 "desc": "preproc",
                 "suffix": "bold",
                 "extension": [
@@ -32,15 +32,15 @@
                     ".nii"
                 ]
             },
-            "bold_mask_mni152nlin6asym": {
+            "bold_mask_boldref": {
                 "datatype": "func",
                 "echo": null,
                 "part": [
                     "mag",
                     null
                 ],
-                "res": "2",
-                "space": "MNI152NLin6Asym",
+                "res": null,
+                "space": null,
                 "desc": "brain",
                 "suffix": "mask",
                 "extension": [
@@ -48,16 +48,17 @@
                     ".nii"
                 ]
             },
-            "bold_mask_native": {
+            "boldref": {
                 "datatype": "func",
                 "echo": null,
                 "part": [
                     "mag",
                     null
                 ],
+                "res": null,
                 "space": null,
-                "desc": "brain",
-                "suffix": "mask",
+                "desc": "coreg",
+                "suffix": "boldref",
                 "extension": [
                     ".nii.gz",
                     ".nii"
@@ -79,19 +80,6 @@
                     ".tsv"
                 ]
             },
-            "dseg_mni152nlin6asym": {
-                "datatype": "anat",
-                "task": null,
-                "run": null,
-                "res": "2",
-                "space": "MNI152NLin6Asym",
-                "desc": null,
-                "suffix": "dseg",
-                "extension": [
-                    ".nii.gz",
-                    ".nii"
-                ]
-            },
             "anat_dseg": {
                 "datatype": "anat",
                 "task": null,
@@ -113,6 +101,7 @@
                 "from": "orig",
                 "to": "boldref",
                 "mode": "image",
+                "desc": ["hmc", null],
                 "suffix": "xfm",
                 "extension": ".txt"
             },
@@ -121,6 +110,7 @@
                 "from": "boldref",
                 "to": ["anat", "T1w", "T2w"],
                 "mode": "image",
+                "desc": ["coreg", null],
                 "suffix": "xfm",
                 "extension": ".txt"
             },
@@ -169,9 +159,10 @@
             "name": "threshold",
             "pattern": "(?:^|_)thresh-([a-zA-Z0-9]+)"
         }
-      ],
+    ],
     "patterns": [
         "sub-{subject}[/ses-{session}]/{datatype<func>|func}/sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_echo-{echo}][_part-{part}][_space-{space}][_res-{res}][_desc-{desc}]_{suffix<bold|boldref|dseg|mask>}.{extension<nii|nii.gz|json>|nii.gz}",
+        "sub-{subject}[/ses-{session}]/{datatype<func>|func}/sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_cohort-{cohort}][_seg-{segmentation}][_res-{res}][_stat-{statistic}][_desc-{desc}]_{suffix<boldmap>}{extension<.nii|.nii.gz|.json>|.nii.gz}",
         "sub-{subject}[/ses-{session}]/{datatype<func>|func}/sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_echo-{echo}][_part-{part}][_space-{space}][_res-{res}][_stat-{statistic}][_desc-{desc}]_{suffix<components|mixing>}.{extension<nii|nii.gz|json>|nii.gz}",
         "sub-{subject}[/ses-{session}]/{datatype<func>|func}/sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_echo-{echo}][_part-{part}][_space-{space}][_res-{res}][_stat-{statistic}][_desc-{desc}]_{suffix<components|mixing>}.{extension<tsv|json>|tsv}",
         "sub-{subject}[/ses-{session}]/{datatype<func>|func}/sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_part-{part}][_desc-{desc}]_{suffix<metrics>}.{extension<tsv|json>|tsv}",
diff --git a/src/fmripost_rapidtide/data/rapidtide_spec.yml b/src/fmripost_rapidtide/data/rapidtide_spec.yml
new file mode 100644
index 0000000..4bb987b
--- /dev/null
+++ b/src/fmripost_rapidtide/data/rapidtide_spec.yml
@@ -0,0 +1,118 @@
+autocorr:
+  filename: desc-autocorr_timeseries.nii.gz
+autocorr_json:
+  filename: desc-autocorr_timeseries.json
+cleansimdistdata:
+  filename: desc-cleansimdistdata_info.tsv.gz
+cleansimdistdata_json:
+  filename: desc-cleansimdistdata_info.json
+confoundfiltercleaned:
+  filename: desc-confoundfilterCleaned_bold.nii.gz
+confoundfiltercleaned_json:
+  filename: desc-confoundfilterCleaned_bold.json
+confoundfilterr2hist:
+  filename: desc-confoundfilterR2_hist.tsv.gz
+confoundfilterr2hist_json:
+  filename: desc-confoundfilterR2_hist.json
+confoundfilterr2map:
+  filename: desc-confoundfilterR2_map.nii.gz
+confoundfilterr2map_json:
+  filename: desc-confoundfilterR2_map.json
+corrfitmask:
+  filename: desc-corrfitmask_bold.nii.gz
+corrfitmask_json:
+  filename: desc-corrfitmask_bold.json
+expandedconfounds:
+  filename: desc-expandedconfounds_timeseries.tsv.gz
+expandedconfounds_json:
+  filename: desc-expandedconfounds_timeseries.json
+formattedruntimings:
+  filename: desc-formattedruntimings_info.tsv
+globallaghist:
+  filename: desc-globallag_hist.tsv.gz
+globallaghist_json:
+  filename: desc-globallag_hist.json
+globalmeanmask:
+  filename: desc-globalmean_mask.nii.gz
+globalmeanmask_json:
+  filename: desc-globalmean_mask.json
+initialmovingregressor:
+  filename: desc-initialmovingregressor_timeseries.tsv.gz
+initialmovingregressor_json:
+  filename: desc-initialmovingregressor_timeseries.json
+lagtcgenerator:
+  filename: desc-lagtcgenerator_timeseries.tsv.gz
+lagtcgenerator_json:
+  filename: desc-lagtcgenerator_timeseries.json
+maxcorrhist:
+  filename: desc-maxcorr_hist.tsv.gz
+maxcorrhist_json:
+  filename: desc-maxcorr_hist.json
+maxcorrmap:
+  filename: desc-maxcorr_map.nii.gz
+maxcorrmap_json:
+  filename: desc-maxcorr_map.json
+maxcorrsqmap:
+  filename: desc-maxcorrsq_map.nii.gz
+maxcorrsqmap_json:
+  filename: desc-maxcorrsq_map.json
+maxtimehist:
+  filename: desc-maxtime_hist.tsv.gz
+maxtimehist_json:
+  filename: desc-maxtime_hist.json
+maxtimemap:
+  filename: desc-maxtime_map.nii.gz
+maxtimemap_json:
+  filename: desc-maxtime_map.json
+maxwidthhist:
+  filename: desc-maxwidth_hist.tsv.gz
+maxwidthhist_json:
+  filename: desc-maxwidth_hist.json
+maxwidthmap:
+  filename: desc-maxwidth_map.nii.gz
+maxwidthmap_json:
+  filename: desc-maxwidth_map.json
+meanmap:
+  filename: desc-mean_map.nii.gz
+meanmap_json:
+  filename: desc-mean_map.json
+movingregressor:
+  filename: desc-movingregressor_timeseries.tsv.gz
+movingregressor_json:
+  filename: desc-movingregressor_timeseries.json
+mtthist:
+  filename: desc-MTT_hist.tsv.gz
+mtthist_json:
+  filename: desc-MTT_hist.json
+mttmap:
+  filename: desc-MTT_map.nii.gz
+mttmap_json:
+  filename: desc-MTT_map.json
+nullsimfunchist:
+  filename: desc-nullsimfunc_hist.tsv.gz
+nullsimfunchist_json:
+  filename: desc-nullsimfunc_hist.json
+oversampledmovingregressor:
+  filename: desc-oversampledmovingregressor_timeseries.tsv.gz
+oversampledmovingregressor_json:
+  filename: desc-oversampledmovingregressor_timeseries.json
+preprocessedconfounds:
+  filename: desc-preprocessedconfounds_timeseries.tsv.gz
+preprocessedconfounds_json:
+  filename: desc-preprocessedconfounds_timeseries.json
+processedmask:
+  filename: desc-processed_mask.nii.gz
+processedmask_json:
+  filename: desc-processed_mask.json
+refinemask:
+  filename: desc-refine_mask.nii.gz
+refinemask_json:
+  filename: desc-refine_mask.json
+refinedmovingregressor:
+  filename: desc-refinedmovingregressor_timeseries.tsv.gz
+refinedmovingregressor_json:
+  filename: desc-refinedmovingregressor_timeseries.json
+timepercentilemap:
+  filename: desc-timepercentile_map.nii.gz
+timepercentilemap_json:
+  filename: desc-timepercentile_map.json
diff --git a/src/fmripost_rapidtide/interfaces/rapidtide.py b/src/fmripost_rapidtide/interfaces/rapidtide.py
index 3cd0ffc..71c2dfd 100644
--- a/src/fmripost_rapidtide/interfaces/rapidtide.py
+++ b/src/fmripost_rapidtide/interfaces/rapidtide.py
@@ -5,6 +5,7 @@
 from nipype.interfaces.base import (
     CommandLine,
     CommandLineInputSpec,
+    DynamicTraitedSpec,
     File,
     TraitedSpec,
     traits,
@@ -19,10 +20,11 @@ class _RapidtideInputSpec(CommandLineInputSpec):
         mandatory=True,
         desc='File to denoise',
     )
-    outputname = traits.Str(
+    prefix = traits.Str(
         argstr='%s',
         position=-1,
-        mandatory=True,
+        mandatory=False,
+        genfile=True,
         desc='Output name',
     )
     # Set by the workflow
@@ -253,7 +255,7 @@ class _RapidtideInputSpec(CommandLineInputSpec):
         mandatory=False,
     )
     nprocs = traits.Int(
-        default=1,
+        1,
         usedefault=True,
         argstr='--nprocs %d',
         mandatory=False,
@@ -261,29 +263,214 @@ class _RapidtideInputSpec(CommandLineInputSpec):
 
 
 class _RapidtideOutputSpec(TraitedSpec):
-    delay_map = File(exists=True, desc='3D map of optimal delay times')
-    regressor_file = File(exists=True, desc='Time series of refined regressor')
-    denoised = File(exists=True, desc='Denoised time series')
+    prefix = traits.Str(desc='Directory containing the results, with prefix.')
+    lagtimesfile = File(exists=True, desc='3D map of optimal delay times')
+    lagtcgeneratorfile = File(exists=True, desc='Time series of refined regressor')
+    maskfile = File(exists=True, desc='Mask file (usually called XXX_desc-corrfit_mask.nii.gz)')
 
 
 class Rapidtide(CommandLine):
     """Run the rapidtide command-line interface."""
 
-    _cmd = 'rapidtide --noprogressbar'
+    _cmd = 'rapidtide --noprogressbar --spcalculation --noglm'
     input_spec = _RapidtideInputSpec
     output_spec = _RapidtideOutputSpec
 
+    def _gen_filename(self, name):
+        if name == 'prefix':
+            return os.path.join(os.getcwd(), 'rapidtide')
+
+        return None
+
     def _list_outputs(self):
         outputs = self._outputs().get()
-        out_dir = os.getcwd()
-        outputname = self.inputs.outputname
-        outputs['delay_map'] = os.path.join(out_dir, f'{outputname}_desc-maxtime_map.nii.gz')
-        outputs['regressor_file'] = os.path.join(
-            out_dir,
-            f'{outputname}_desc-refinedmovingregressor_timeseries.tsv.gz',
-        )
+        prefix = self.inputs.prefix
+        outputs['prefix'] = prefix
+        outputs['lagtimesfile'] = f'{prefix}_desc-maxtime_map.nii.gz'
+        outputs['lagtcgeneratorfile'] = f'{prefix}_desc-lagtcgenerator_timeseries.tsv.gz'
+        outputs['maskfile'] = f'{prefix}_desc-corrfit_mask.nii.gz'
+
+        return outputs
+
+
+class _RetroLagTCSInputSpec(CommandLineInputSpec):
+    in_file = File(
+        exists=True,
+        argstr='%s',
+        position=0,
+        mandatory=True,
+        desc='The name of 4D nifti fmri target file.',
+    )
+    maskfile = File(
+        exists=True,
+        argstr='%s',
+        position=1,
+        mandatory=True,
+        desc='The mask file to use (usually called XXX_desc-corrfit_mask.nii.gz)',
+    )
+    lagtimesfile = File(
+        exists=True,
+        argstr='%s',
+        position=2,
+        mandatory=True,
+        desc='The name of the lag times file (usually called XXX_desc-maxtime_map.nii.gz)',
+    )
+    lagtcgeneratorfile = File(
+        exists=True,
+        argstr='%s',
+        position=3,
+        mandatory=True,
+        desc=(
+            'The root name of the lagtc generator file '
+            '(usually called XXX_desc-lagtcgenerator_timeseries)'
+        ),
+    )
+    prefix = traits.Str(
+        argstr='%s',
+        position=4,
+        mandatory=False,
+        genfile=True,
+        desc='Output root.',
+    )
+    glmderivs = traits.Int(
+        0,
+        argstr='--glmderivs %d',
+        mandatory=False,
+        desc='When doing final GLM, include derivatives up to NDERIVS order. Default is 0.',
+        usedefault=True,
+    )
+    nprocs = traits.Int(
+        1,
+        usedefault=True,
+        argstr='--nprocs %d',
+        mandatory=False,
+        desc=(
+            'Use NPROCS worker processes for multiprocessing. '
+            'Setting NPROCS to less than 1 sets the number of worker processes to n_cpus.'
+        ),
+    )
+    numskip = traits.Int(
+        0,
+        argstr='--numskip %d',
+        usedefault=True,
+        mandatory=False,
+        desc='Skip NUMSKIP points at the beginning of the fmri file.',
+    )
+    noprogressbar = traits.Bool(
+        True,
+        argstr='--noprogressbar',
+        mandatory=False,
+        usedefault=True,
+        desc='Will disable showing progress bars (helpful if stdout is going to a file).',
+    )
+    debug = traits.Bool(
+        False,
+        argstr='--debug',
+        mandatory=False,
+        usedefault=True,
+        desc='Output lots of helpful information.',
+    )
+
+
+class _RetroLagTCSOutputSpec(DynamicTraitedSpec):
+    filter_file = File(exists=True, desc='Filter file')
+
+
+class RetroLagTCS(CommandLine):
+    """Run the retrolagtcs command-line interface."""
+
+    _cmd = 'retrolagtcs'
+    input_spec = _RetroLagTCSInputSpec
+    output_spec = _RetroLagTCSOutputSpec
+
+    def _gen_filename(self, name):
+        if name == 'prefix':
+            return os.path.join(os.getcwd(), 'retrolagtcs')
+
+        return None
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        prefix = self.inputs.prefix
+        outputs['filter_file'] = f'{prefix}_desc-lfofilterEV_bold.nii.gz'
+        if self.inputs.glmderivs > 0:
+            for i_deriv in range(self.inputs.glmderivs):
+                outputs[f'filter_file_deriv{i_deriv + 1}'] = (
+                    f'{prefix}_desc-lfofilterEVDeriv{i_deriv + 1}_bold.nii.gz'
+                )
+
+        return outputs
+
+
+class _RetroGLMInputSpec(CommandLineInputSpec):
+    in_file = File(
+        exists=True,
+        argstr='%s',
+        position=0,
+        mandatory=True,
+        desc='File to denoise',
+    )
+    datafileroot = traits.Str(
+        argstr='%s',
+        position=1,
+        mandatory=True,
+        desc=(
+            'The root name of the previously run rapidtide dataset '
+            '(everything up to but not including the underscore.)'
+        ),
+    )
+    prefix = traits.Str(
+        argstr='--alternateoutput %s',
+        mandatory=False,
+        genfile=True,
+        desc='Output name',
+    )
+    glmderivs = traits.Int(
+        0,
+        argstr='--glmderivs %d',
+        usedefault=True,
+        mandatory=False,
+        desc='When doing final GLM, include derivatives up to NDERIVS order.',
+    )
+    nprocs = traits.Int(
+        1,
+        usedefault=True,
+        argstr='--nprocs %d',
+        mandatory=False,
+    )
+    numskip = traits.Int(
+        argstr='--numskip %d',
+        mandatory=False,
+    )
+
+
+class _RetroGLMOutputSpec(TraitedSpec):
+    denoised = File(exists=True, desc='Denoised time series')
+    denoised_json = File(exists=True, desc='Denoised time series metadata')
+
+
+class RetroGLM(CommandLine):
+    """Run the retroglm command-line interface to denoise BOLD with existing rapidtide outputs."""
+
+    _cmd = 'retroglm --noprogressbar'
+    input_spec = _RetroGLMInputSpec
+    output_spec = _RetroGLMOutputSpec
+
+    def _gen_filename(self, name):
+        if name == 'prefix':
+            return os.path.join(os.getcwd(), 'retroglm')
+
+        return None
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        prefix_dir = self.inputs.prefix
+        datafileroot = self.inputs.datafileroot
+        file_prefix = os.path.basename(datafileroot)
         outputs['denoised'] = os.path.join(
-            out_dir,
-            f'{outputname}_desc-lfofilterCleaned_bold.nii.gz',
+            prefix_dir, f'{file_prefix}_desc-lfofilterCleaned_bold.nii.gz'
+        )
+        outputs['denoised_json'] = os.path.join(
+            prefix_dir, f'{file_prefix}_desc-lfofilterCleaned_bold.json'
         )
         return outputs
diff --git a/src/fmripost_rapidtide/interfaces/reportlets.py b/src/fmripost_rapidtide/interfaces/reportlets.py
index 119c5cb..fcf3540 100644
--- a/src/fmripost_rapidtide/interfaces/reportlets.py
+++ b/src/fmripost_rapidtide/interfaces/reportlets.py
@@ -29,7 +29,6 @@
 
 from nipype.interfaces.base import (
     BaseInterfaceInputSpec,
-    Directory,
     File,
     InputMultiObject,
     SimpleInterface,
@@ -163,114 +162,57 @@ def _generate_segment(self):
         )
 
 
-class _ICARapidtideInputSpecRPT(BaseInterfaceInputSpec):
-    in_file = File(
+class _FCInflationPlotInputSpecRPT(BaseInterfaceInputSpec):
+    fcinflation_file = File(
         exists=True,
         mandatory=True,
-        desc='BOLD series input to Rapidtide',
-    )
-    melodic_dir = Directory(
-        exists=True,
-        mandatory=True,
-        desc='MELODIC directory containing the ICA outputs',
-    )
-    rapidtide_noise_ics = File(
-        exists=True,
-        desc='Noise components estimated by Rapidtide, in a comma-separated values file.',
+        desc='FC inflation time series',
     )
     out_report = File(
-        'rapidtide_reportlet.svg',
+        'fcinflation_reportlet.svg',
         usedefault=True,
         desc='Filename for the visual report generated by Nipype.',
     )
-    report_mask = File(
-        exists=True,
-        mandatory=True,
-        desc=(
-            'Mask used to draw the outline on the reportlet. '
-            'If not set the mask will be derived from the data.'
-        ),
-    )
-    compress_report = traits.Bool(
-        True,
-        usedefault=True,
-        desc='Whether to compress the reportlet with SVGO.',
-    )
 
 
-class _ICARapidtideOutputSpecRPT(TraitedSpec):
+class _FCInflationPlotOutputSpecRPT(TraitedSpec):
     out_report = File(
         exists=True,
         desc='Filename for the visual report generated by Nipype.',
     )
 
 
-class ICARapidtideRPT(SimpleInterface):
+class FCInflationPlotRPT(SimpleInterface):
     """Create a reportlet for Rapidtide outputs."""
 
-    input_spec = _ICARapidtideInputSpecRPT
-    output_spec = _ICARapidtideOutputSpecRPT
-
-    def _run_interface(self, runtime):
-        from niworkflows.viz.utils import plot_melodic_components
-
-        out_file = os.path.abspath(self.inputs.out_report)
-
-        plot_melodic_components(
-            melodic_dir=self.inputs.melodic_dir,
-            in_file=self.inputs.in_file,
-            out_file=out_file,
-            compress=self.inputs.compress_report,
-            report_mask=self.inputs.report_mask,
-            noise_components_file=self.inputs.rapidtide_noise_ics,
-        )
-        self._results['out_report'] = out_file
-        return runtime
-
-
-class _ICARapidtideMetricsInputSpecRPT(BaseInterfaceInputSpec):
-    rapidtide_features = File(
-        exists=True,
-        mandatory=True,
-        desc='Rapidtide metrics',
-    )
-    out_report = File(
-        'metrics_reportlet.svg',
-        usedefault=True,
-        desc='Filename for the visual report generated by Nipype.',
-    )
-
-
-class _ICARapidtideMetricsOutputSpecRPT(TraitedSpec):
-    out_report = File(
-        exists=True,
-        desc='Filename for the visual report generated by Nipype.',
-    )
-
-
-class ICARapidtideMetricsRPT(SimpleInterface):
-    """Create a reportlet for Rapidtide outputs."""
-
-    input_spec = _ICARapidtideMetricsInputSpecRPT
-    output_spec = _ICARapidtideMetricsOutputSpecRPT
+    input_spec = _FCInflationPlotInputSpecRPT
+    output_spec = _FCInflationPlotOutputSpecRPT
 
     def _run_interface(self, runtime):
+        import matplotlib.pyplot as plt
         import pandas as pd
         import seaborn as sns
 
-        out_file = os.path.abspath(self.inputs.out_report)
-
-        df = pd.read_table(self.inputs.rapidtide_features)
+        sns.set_theme(style='whitegrid')
 
-        sns.set_theme(style='ticks')
+        out_file = os.path.abspath(self.inputs.out_report)
 
-        g = sns.pairplot(
-            df,
-            hue='classification',
-            vars=['edge_fract', 'csf_fract', 'max_RP_corr', 'HFC'],
-            palette={'rejected': 'red', 'accepted': 'blue'},
-            corner=True,
-        )
-        g.savefig(out_file)
+        df = pd.read_table(self.inputs.fcinflation_file)
+
+        fig, ax = plt.subplots(figsize=(16, 8))
+        palette = ['red', 'lightblue', 'blue']
+        for i_col, col in enumerate(['preprocessed', 'denoised', 'rapidtide']):
+            df[f'{col}_mean_minus_std'] = df[f'{col}_mean'] - df[f'{col}_std']
+            df[f'{col}_mean_plus_std'] = df[f'{col}_mean'] + df[f'{col}_std']
+            sns.lineplot(x='timepoint', y=f'{col}_mean', data=df, ax=ax, color=palette[i_col])
+            ax.fill_between(
+                x=df['timepoint'],
+                y1=df[f'{col}_mean_minus_std'],
+                y2=df[f'{col}_mean_plus_std'],
+                alpha=0.5,
+                color=palette[i_col],
+            )
+        ax.set_xlim(0, df['timepoint'].max())
+        fig.savefig(out_file)
         self._results['out_report'] = out_file
         return runtime
diff --git a/src/fmripost_rapidtide/tests/test_base.py b/src/fmripost_rapidtide/tests/test_base.py
index a36c2e5..d3a8d12 100644
--- a/src/fmripost_rapidtide/tests/test_base.py
+++ b/src/fmripost_rapidtide/tests/test_base.py
@@ -5,19 +5,43 @@
 from fmripost_rapidtide import config
 
 
-def test_init_rapidtide_wf(tmp_path_factory):
-    from fmripost_rapidtide.workflows.rapidtide import init_rapidtide_wf
+def test_init_rapidtide_fit_wf(tmp_path_factory):
+    from fmripost_rapidtide.workflows.rapidtide import init_rapidtide_fit_wf
 
-    tempdir = tmp_path_factory.mktemp('test_init_rapidtide_wf')
+    tempdir = tmp_path_factory.mktemp('test_init_rapidtide_fit_wf')
 
     with mock_config():
         config.execution.output_dir = tempdir / 'out'
         config.execution.work_dir = tempdir / 'work'
-        config.workflow.denoise_method = ['nonaggr', 'orthaggr']
-        config.workflow.melodic_dim = -200
         config.workflow.err_on_warn = False
+        config.workflow.ampthresh = -1.0
+        config.workflow.autorespdelete = False
+        config.workflow.autosync = False
+        config.workflow.bipolar = False
+        config.workflow.confoundderiv = True
+        config.workflow.confoundpowers = 1
+        config.workflow.corrweighting = 'regressor'
+        config.workflow.detrendorder = 1
+        config.workflow.filterband = 'lfo'
+        config.workflow.glmderivs = 1
+        config.workflow.globalpcacomponents = 0.8
+        config.workflow.globalsignalmethod = 'sum'
+        config.workflow.lagmaxthresh = 5.0
+        config.workflow.lagminthresh = 0.5
+        config.workflow.maxpasses = 1
+        config.workflow.numnull = 100
+        config.workflow.numtozero = 0
+        config.workflow.outputlevel = 'min'
+        config.workflow.pcacomponents = 0.8
+        config.workflow.sigmalimit = 1000.0
+        config.workflow.searchrange = [-30, 30]
+        config.workflow.sigmathresh = 100.0
+        config.workflow.simcalcrange = [-1, -1]
+        config.workflow.spatialfilt = 4.0
+        config.workflow.timerange = [-1, -1]
+        config.nipype.omp_nthreads = 1
 
-        wf = init_rapidtide_wf(
+        wf = init_rapidtide_fit_wf(
             bold_file='sub-01_task-rest_bold.nii.gz',
             metadata={'RepetitionTime': 2.0},
             mem_gb={
diff --git a/src/fmripost_rapidtide/tests/test_utils_bids.py b/src/fmripost_rapidtide/tests/test_utils_bids.py
index 5fb4040..3de7dfe 100644
--- a/src/fmripost_rapidtide/tests/test_utils_bids.py
+++ b/src/fmripost_rapidtide/tests/test_utils_bids.py
@@ -80,8 +80,8 @@ def test_collect_derivatives_minimal(minimal_ignore_list):
         # TODO: Add bold_mask_native to the dataset
         # 'bold_mask_native': 'sub-01_task-mixedgamblestask_run-01_desc-brain_mask.nii.gz',
         'bold_mask_native': None,
-        'confounds': None,
-        'hmc': (
+        'bold_confounds': None,
+        'bold_hmc': (
             'sub-01_task-mixedgamblestask_run-01_from-orig_to-boldref_mode-image_desc-hmc_xfm.txt'
         ),
         'boldref2anat': (
@@ -123,8 +123,8 @@ def test_collect_derivatives_full(full_ignore_list):
             'mask.nii.gz'
         ),
         'bold_mask_native': None,
-        'confounds': 'sub-01_task-mixedgamblestask_run-01_desc-confounds_timeseries.tsv',
-        'hmc': (
+        'bold_confounds': 'sub-01_task-mixedgamblestask_run-01_desc-confounds_timeseries.tsv',
+        'bold_hmc': (
             'sub-01_task-mixedgamblestask_run-01_from-orig_to-boldref_mode-image_desc-hmc_xfm.txt'
         ),
         'boldref2anat': (
@@ -148,4 +148,4 @@ def check_expected(subject_data, expected):
             for item, expected_item in zip(subject_data[key], value):
                 assert os.path.basename(item) == expected_item
         else:
-            assert subject_data[key] is value
+            assert subject_data[key] is value, f'Key {key} is not {value}.'
diff --git a/src/fmripost_rapidtide/utils/bids.py b/src/fmripost_rapidtide/utils/bids.py
index c2b0302..ef95aac 100644
--- a/src/fmripost_rapidtide/utils/bids.py
+++ b/src/fmripost_rapidtide/utils/bids.py
@@ -137,6 +137,7 @@ def collect_derivatives(
                 query['to'] = fieldmap_id
 
             item = layout.get(return_type='filename', **query)
+
             if not item:
                 derivs_cache[k] = None
             elif not allow_multiple and len(item) > 1:
@@ -151,13 +152,13 @@ def collect_derivatives(
         spaces_found, bold_outputspaces, bold_mask_outputspaces = [], [], []
         for space in spaces.references:
             # First try to find processed BOLD+mask files in the requested space
-            bold_query = {**entities, **spec['derivatives']['bold_mni152nlin6asym']}
+            bold_query = {**entities, **spec['derivatives']['bold_boldref']}
             bold_query['space'] = space.space
             bold_query = {**bold_query, **space.spec}
             bold_item = layout.get(return_type='filename', **bold_query)
             bold_outputspaces.append(bold_item[0] if bold_item else None)
 
-            mask_query = {**entities, **spec['derivatives']['bold_mask_mni152nlin6asym']}
+            mask_query = {**entities, **spec['derivatives']['bold_boldref']}
             mask_query['space'] = space.space
             mask_query = {**mask_query, **space.spec}
             mask_item = layout.get(return_type='filename', **mask_query)
diff --git a/src/fmripost_rapidtide/workflows/base.py b/src/fmripost_rapidtide/workflows/base.py
index 3aecebc..6ccce51 100644
--- a/src/fmripost_rapidtide/workflows/base.py
+++ b/src/fmripost_rapidtide/workflows/base.py
@@ -134,11 +134,13 @@ def init_single_subject_wf(subject_id: str):
 
     """
     from bids.utils import listify
+    from nipype.interfaces import utility as niu
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
     from niworkflows.interfaces.bids import BIDSInfo
     from niworkflows.interfaces.nilearn import NILEARN_VERSION
 
     from fmripost_rapidtide.interfaces.bids import DerivativesDataSink
+    from fmripost_rapidtide.interfaces.nilearn import MeanImage
     from fmripost_rapidtide.interfaces.reportlets import AboutSummary, SubjectSummary
     from fmripost_rapidtide.utils.bids import collect_derivatives
 
@@ -198,14 +200,14 @@ def init_single_subject_wf(subject_id: str):
             spaces=None,
         )
         # Patch standard-space BOLD files into 'bold' key
-        subject_data['bold'] = listify(subject_data['bold_mni152nlin6asym'])
+        subject_data['bold'] = listify(subject_data['bold_boldref'])
 
-    if not subject_data['bold_mni152nlin6asym']:
+    if not subject_data['bold_boldref']:
         task_id = config.execution.task_id
         raise RuntimeError(
-            f"No MNI152NLin6Asym:res-2 BOLD images found for participant {subject_id} and "
+            f"No boldref:res-native BOLD images found for participant {subject_id} and "
             f"task {task_id if task_id else '<all>'}. "
-            "All workflows require MNI152NLin6Asym:res-2 BOLD images. "
+            "All workflows require boldref:res-native BOLD images. "
             f"Please check your BIDS filters: {config.execution.bids_filters}."
         )
 
@@ -279,23 +281,55 @@ def init_single_subject_wf(subject_id: str):
 """
     workflow.__desc__ += func_pre_desc
 
-    for bold_file in subject_data['bold']:
-        single_run_wf = init_single_run_wf(bold_file)
-        workflow.add_nodes([single_run_wf])
+    denoise_within_run = (len(subject_data['bold']) > 1) and not config.workflow.average_over_runs
+    if not denoise_within_run:
+        # Average the lag map across runs before denoising
+        merge_lag_maps = pe.Node(
+            niu.Merge(len(subject_data['bold'])),
+            name='merge_lag_maps',
+        )
+        average_lag_map = pe.Node(
+            MeanImage(),
+            name='average_lag_map',
+        )
 
-    return clean_datasinks(workflow)
+    for i_run, bold_file in enumerate(subject_data['bold']):
+        fit_single_run_wf = init_fit_single_run_wf(bold_file=bold_file)
+        denoise_single_run_wf = init_denoise_single_run_wf(bold_file=bold_file)
+        workflow.connect([
+            (fit_single_run_wf, denoise_single_run_wf, [
+                ('outputnode.regressor', 'inputnode.regressor'),
+            ]),
+        ])  # fmt:skip
 
+        if denoise_within_run:
+            # Denoise the BOLD data using the run-wise lag map
+            workflow.connect([
+                (fit_single_run_wf, denoise_single_run_wf, [
+                    ('outputnode.delay_map', 'inputnode.delay_map'),
+                ]),
+            ])  # fmt:skip
+        else:
+            # Denoise the BOLD data using the mean lag map
+            workflow.connect([
+                (fit_single_run_wf, merge_lag_maps, [('outputnode.delay_map', f'in{i_run + 1}')]),
+                (merge_lag_maps, average_lag_map, [('out', 'in_file')]),
+                (average_lag_map, denoise_single_run_wf, [('out_file', 'inputnode.delay_map')]),
+            ])  # fmt:skip
+
+    return workflow
 
-def init_single_run_wf(bold_file):
+
+def init_fit_single_run_wf(*, bold_file):
     """Set up a single-run workflow for fMRIPost-rapidtide."""
     from fmriprep.utils.misc import estimate_bold_mem_usage
     from nipype.interfaces import utility as niu
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
+    from fmripost_rapidtide.interfaces.misc import ApplyTransforms
     from fmripost_rapidtide.utils.bids import collect_derivatives, extract_entities
-    from fmripost_rapidtide.workflows.confounds import init_confounds_wf
     from fmripost_rapidtide.workflows.outputs import init_func_fit_reports_wf
-    from fmripost_rapidtide.workflows.rapidtide import init_rapidtide_wf
+    from fmripost_rapidtide.workflows.rapidtide import init_rapidtide_fit_wf
 
     spaces = config.workflow.spaces
     omp_nthreads = config.nipype.omp_nthreads
@@ -344,8 +378,8 @@ def init_single_run_wf(bold_file):
             raise ValueError('Motion parameters cannot be extracted from transforms yet.')
 
     else:
-        # Collect MNI152NLin6Asym:res-2 derivatives
-        # Only derivatives dataset was passed in, so we expected standard-space derivatives
+        # Collect boldref:res-native derivatives
+        # Only derivatives dataset was passed in, so we expected boldref-space derivatives
         functional_cache.update(
             collect_derivatives(
                 raw_dataset=None,
@@ -376,23 +410,39 @@ def init_single_run_wf(bold_file):
         skip_vols = get_nss(functional_cache['confounds'])
 
     # Run rapidtide
-    rapidtide_wf = init_rapidtide_wf(bold_file=bold_file, metadata=bold_metadata, mem_gb=mem_gb)
+    rapidtide_wf = init_rapidtide_fit_wf(
+        bold_file=bold_file,
+        metadata=bold_metadata,
+        mem_gb=mem_gb,
+    )
     rapidtide_wf.inputs.inputnode.confounds = functional_cache['confounds']
     rapidtide_wf.inputs.inputnode.skip_vols = skip_vols
 
-    mni6_buffer = pe.Node(niu.IdentityInterface(fields=['bold', 'bold_mask']), name='mni6_buffer')
+    boldref_buffer = pe.Node(
+        niu.IdentityInterface(fields=['bold', 'bold_mask']),
+        name='boldref_buffer',
+    )
+
+    # Warp the dseg from anatomical space to boldref space
+    dseg_to_boldref = pe.Node(
+        ApplyTransforms(
+            interpolation='GenericLabel',
+            input_image=functional_cache['anat_dseg'],
+            reference_image=functional_cache['boldref'],
+            transforms=[functional_cache['boldref2anat']],
+            invert_transform_flags=[True],
+        ),
+        name='dseg_to_boldref',
+    )
 
-    if ('bold_mni152nlin6asym' not in functional_cache) and ('bold_raw' in functional_cache):
+    if ('bold_boldref' not in functional_cache) and ('bold_raw' in functional_cache):
         # Resample to MNI152NLin6Asym:res-2, for rapidtide denoising
         from fmriprep.workflows.bold.apply import init_bold_volumetric_resample_wf
         from fmriprep.workflows.bold.stc import init_bold_stc_wf
         from niworkflows.interfaces.header import ValidateImage
-        from templateflow.api import get as get_template
-
-        from fmripost_rapidtide.interfaces.misc import ApplyTransforms
 
         workflow.__desc__ += """\
-Raw BOLD series were resampled to MNI152NLin6Asym:res-2, for rapidtide denoising.
+Raw BOLD series were resampled to boldref:res-native, for rapidtide denoising.
 """
 
         validate_bold = pe.Node(
@@ -419,17 +469,7 @@ def init_single_run_wf(bold_file):
         else:
             workflow.connect([(validate_bold, stc_buffer, [('out_file', 'bold_file')])])
 
-        mni6_mask = str(
-            get_template(
-                'MNI152NLin6Asym',
-                resolution=2,
-                desc='brain',
-                suffix='mask',
-                extension=['.nii', '.nii.gz'],
-            )
-        )
-
-        bold_MNI6_wf = init_bold_volumetric_resample_wf(
+        bold_boldref_wf = init_bold_volumetric_resample_wf(
             metadata=bold_metadata,
             fieldmap_id=None,  # XXX: Ignoring the field map for now
             omp_nthreads=omp_nthreads,
@@ -437,71 +477,41 @@ def init_single_run_wf(bold_file):
             jacobian='fmap-jacobian' not in config.workflow.ignore,
             name='bold_MNI6_wf',
         )
-        bold_MNI6_wf.inputs.inputnode.motion_xfm = functional_cache['hmc']
-        bold_MNI6_wf.inputs.inputnode.boldref2fmap_xfm = functional_cache['boldref2fmap']
-        bold_MNI6_wf.inputs.inputnode.boldref2anat_xfm = functional_cache['boldref2anat']
-        bold_MNI6_wf.inputs.inputnode.anat2std_xfm = functional_cache['anat2mni152nlin6asym']
-        bold_MNI6_wf.inputs.inputnode.resolution = '02'
+        bold_boldref_wf.inputs.inputnode.motion_xfm = functional_cache['hmc']
+        bold_boldref_wf.inputs.inputnode.boldref2fmap_xfm = functional_cache['boldref2fmap']
+        bold_boldref_wf.inputs.inputnode.resolution = 'native'
         # use mask as boldref?
-        bold_MNI6_wf.inputs.inputnode.bold_ref_file = functional_cache['bold_mask_native']
-        bold_MNI6_wf.inputs.inputnode.target_mask = mni6_mask
-        bold_MNI6_wf.inputs.inputnode.target_ref_file = mni6_mask
+        bold_boldref_wf.inputs.inputnode.bold_ref_file = functional_cache['boldref']
+        bold_boldref_wf.inputs.inputnode.target_mask = functional_cache['bold_mask_boldref']
+        bold_boldref_wf.inputs.inputnode.target_ref_file = functional_cache['boldref']
 
         workflow.connect([
             # Resample BOLD to MNI152NLin6Asym, may duplicate bold_std_wf above
             # XXX: Ignoring the field map for now
-            # (inputnode, bold_MNI6_wf, [
+            # (inputnode, bold_boldref_wf, [
             #     ('fmap_ref', 'inputnode.fmap_ref'),
             #     ('fmap_coeff', 'inputnode.fmap_coeff'),
             #     ('fmap_id', 'inputnode.fmap_id'),
             # ]),
-            (stc_buffer, bold_MNI6_wf, [('bold_file', 'inputnode.bold_file')]),
-            (bold_MNI6_wf, mni6_buffer, [('outputnode.bold_file', 'bold')]),
+            (stc_buffer, bold_boldref_wf, [('bold_file', 'inputnode.bold_file')]),
+            (bold_boldref_wf, boldref_buffer, [('outputnode.bold_file', 'bold')]),
         ])  # fmt:skip
 
-        # Warp the mask as well
-        mask_to_mni6 = pe.Node(
-            ApplyTransforms(
-                interpolation='GenericLabel',
-                input_image=functional_cache['bold_mask_native'],
-                reference_image=mni6_mask,
-                transforms=[
-                    functional_cache['anat2mni152nlin6asym'],
-                    functional_cache['boldref2anat'],
-                ],
-            ),
-            name='mask_to_mni6',
-        )
-        workflow.connect([(mask_to_mni6, mni6_buffer, [('output_image', 'bold_mask')])])
-
-        # And the dseg
-        dseg_to_mni6 = pe.Node(
-            ApplyTransforms(
-                interpolation='GenericLabel',
-                input_image=functional_cache['anat_dseg'],
-                reference_image=mni6_mask,
-                transforms=[functional_cache['anat2mni152nlin6asym']],
-            ),
-            name='mask_to_mni6',
-        )
-        workflow.connect([(dseg_to_mni6, mni6_buffer, [('output_image', 'dseg')])])
-
-    elif 'bold_mni152nlin6asym' in functional_cache:
+    elif 'bold_boldref' in functional_cache:
         workflow.__desc__ += """\
-Preprocessed BOLD series in MNI152NLin6Asym:res-2 space were collected for rapidtide denoising.
+Preprocessed BOLD series in boldref:res-native space were collected for rapidtide denoising.
 """
-        mni6_buffer.inputs.bold = functional_cache['bold_mni152nlin6asym']
-        mni6_buffer.inputs.bold_mask = functional_cache['bold_mask_mni152nlin6asym']
-        mni6_buffer.inputs.dseg = functional_cache['dseg_mni152nlin6asym']
+        boldref_buffer.inputs.bold = functional_cache['bold_boldref']
+        boldref_buffer.inputs.bold_mask = functional_cache['bold_mask_boldref']
 
     else:
         raise ValueError('No valid BOLD series found for rapidtide denoising.')
 
     workflow.connect([
-        (mni6_buffer, rapidtide_wf, [
-            ('bold', 'inputnode.bold_std'),
-            ('bold_mask', 'inputnode.bold_mask_std'),
-            ('dseg', 'inputnode.dseg_std'),
+        (dseg_to_boldref, rapidtide_wf, [('output_image', 'inputnode.dseg')]),
+        (boldref_buffer, rapidtide_wf, [
+            ('bold', 'inputnode.bold'),
+            ('bold_mask', 'inputnode.bold_mask'),
         ]),
     ])  # fmt:skip
 
@@ -510,36 +520,97 @@ def init_single_run_wf(bold_file):
     func_fit_reports_wf.inputs.inputnode.source_file = bold_file
     func_fit_reports_wf.inputs.inputnode.anat2std_xfm = functional_cache['anat2mni152nlin6asym']
     func_fit_reports_wf.inputs.inputnode.anat_dseg = functional_cache['anat_dseg']
-    workflow.connect([(mni6_buffer, func_fit_reports_wf, [('bold', 'inputnode.bold_mni6')])])
+    workflow.connect([(boldref_buffer, func_fit_reports_wf, [('bold', 'inputnode.bold_mni6')])])
 
-    # Generate confounds
-    confounds_wf = init_confounds_wf(
-        bold_file=bold_file,
-        mem_gb=mem_gb['filesize'],
+    return clean_datasinks(workflow, bold_file=bold_file)
+
+
+def init_denoise_single_run_wf(*, bold_file: str):
+    """Denoise a single run using rapidtide."""
+
+    from nipype.interfaces import utility as niu
+    from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+
+    from fmripost_rapidtide.interfaces.bids import DerivativesDataSink
+    from fmripost_rapidtide.interfaces.rapidtide import RetroGLM
+    from fmripost_rapidtide.workflows.confounds import init_denoising_confounds_wf
+
+    workflow = Workflow(name=_get_wf_name(bold_file, 'rapidtide_denoise'))
+    workflow.__postdesc__ = """\
+Identification and removal of traveling wave artifacts was performed using rapidtide.
+"""
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'bold',
+                'bold_mask',
+                'dseg',
+                'regressor',
+                'lag_map',
+                'skip_vols',
+            ],
+        ),
+        name='inputnode',
+    )
+    denoise_bold = pe.Node(
+        RetroGLM(),
+        name='denoise_bold',
+    )
+    workflow.connect([
+        (inputnode, denoise_bold, [
+            ('bold', 'in_file'),
+            ('bold_mask', 'brainmask'),
+            ('regressor', 'regressor'),
+            ('lag_map', 'lag_map'),
+            ('skip_vols', 'numskip'),
+        ]),
+    ])  # fmt:skip
+
+    ds_denoised_bold = pe.Node(
+        DerivativesDataSink(
+            compress=True,
+            desc='denoised',
+            suffix='bold',
+        ),
+        name='ds_denoised_bold',
+        run_without_submitting=True,
     )
     workflow.connect([
-        (mni6_buffer, confounds_wf, [
+        (denoise_bold, ds_denoised_bold, [
+            ('denoised', 'in_file'),
+            ('denoised_json', 'meta_dict'),
+        ]),
+    ])  # fmt:skip
+
+    # Generate confounds
+    denoising_confounds_wf = init_denoising_confounds_wf()
+    workflow.connect([
+        (inputnode, denoising_confounds_wf, [
             ('bold', 'inputnode.preprocessed_bold'),
             ('bold_mask', 'inputnode.mask'),
         ]),
-        (rapidtide_wf, confounds_wf, [
-            ('outputnode.confound_regressed', 'inputnode.denoised_bold'),
-            ('outputnode.denoised', 'inputnode.rapidtide_bold'),
+        (denoise_bold, denoising_confounds_wf, [
+            ('denoised', 'inputnode.denoised_bold'),
+            ('denoised', 'inputnode.rapidtide_bold'),
         ]),
     ])  # fmt:skip
 
-    return workflow
-
+    return clean_datasinks(workflow, bold_file=bold_file)
 
-def _prefix(subid):
-    return subid if subid.startswith('sub-') else f'sub-{subid}'
 
-
-def clean_datasinks(workflow: pe.Workflow) -> pe.Workflow:
-    """Overwrite ``out_path_base`` of smriprep's DataSinks."""
+def clean_datasinks(workflow: pe.Workflow, bold_file: str) -> pe.Workflow:
+    """Overwrite attributes of DataSinks."""
     for node in workflow.list_node_names():
-        if node.split('.')[-1].startswith('ds_'):
+        node_name = node.split('.')[-1]
+        if node_name.startswith('ds_'):
             workflow.get_node(node).interface.out_path_base = ''
+            workflow.get_node(node).inputs.base_directory = str(config.execution.output_dir)
+            workflow.get_node(node).inputs.source_file = bold_file
+
+        if node_name.startswith('ds_report_'):
+            workflow.get_node(node).inputs.datatype = 'figures'
+
     return workflow
 
 
diff --git a/src/fmripost_rapidtide/workflows/confounds.py b/src/fmripost_rapidtide/workflows/confounds.py
index 9c0201b..d186fee 100644
--- a/src/fmripost_rapidtide/workflows/confounds.py
+++ b/src/fmripost_rapidtide/workflows/confounds.py
@@ -29,10 +29,10 @@
 """
 
 
-def init_confounds_wf(
+def init_denoising_confounds_wf(
     bold_file: str,
     mem_gb: float,
-    name: str = 'confounds_wf',
+    name: str = 'denoising_confounds_wf',
 ):
     from nipype.interfaces import utility as niu
     from nipype.pipeline import engine as pe
@@ -41,6 +41,7 @@ def init_confounds_wf(
     from fmripost_rapidtide.config import DEFAULT_MEMORY_MIN_GB
     from fmripost_rapidtide.interfaces.bids import DerivativesDataSink
     from fmripost_rapidtide.interfaces.confounds import FCInflation
+    from fmripost_rapidtide.interfaces.reportlets import FCInflationPlotRPT
 
     workflow = Workflow(name=name)
 
@@ -130,6 +131,24 @@ def init_confounds_wf(
     )
     workflow.connect([(merge_fci, ds_metrics, [('confounds_metrics', 'in_file')])])
 
+    # Generate reportlets
+    plot_fcinflation = pe.Node(
+        FCInflationPlotRPT(),
+        name='plot_fcinflation',
+    )
+    workflow.connect([(merge_fci, plot_fcinflation, [('confounds_file', 'fcinflation_file')])])
+
+    ds_report_fcinflation = pe.Node(
+        DerivativesDataSink(
+            desc='fcinflation',
+            suffix='bold',
+            extension='.svg',
+        ),
+        name='ds_report_fcinflation',
+        run_without_submitting=True,
+    )
+    workflow.connect([(plot_fcinflation, ds_report_fcinflation, [('out_report', 'in_file')])])
+
     return workflow
 
 
diff --git a/src/fmripost_rapidtide/workflows/rapidtide.py b/src/fmripost_rapidtide/workflows/rapidtide.py
index 79d5473..ea98046 100644
--- a/src/fmripost_rapidtide/workflows/rapidtide.py
+++ b/src/fmripost_rapidtide/workflows/rapidtide.py
@@ -26,11 +26,10 @@
 from nipype.pipeline import engine as pe
 
 from fmripost_rapidtide import config
-from fmripost_rapidtide.interfaces.rapidtide import Rapidtide
 from fmripost_rapidtide.utils.utils import _get_wf_name
 
 
-def init_rapidtide_wf(
+def init_rapidtide_fit_wf(
     *,
     bold_file: str,
     metadata: dict,
@@ -54,9 +53,9 @@ def init_rapidtide_wf(
             :graph2use: orig
             :simple_form: yes
 
-            from fmripost_rapidtide.workflows.bold.confounds import init_rapidtide_wf
+            from fmripost_rapidtide.workflows.rapidtide import init_rapidtide_fit_wf
 
-            wf = init_rapidtide_wf(
+            wf = init_rapidtide_fit_wf(
                 bold_file="fake.nii.gz",
                 metadata={"RepetitionTime": 1.0},
             )
@@ -70,11 +69,11 @@ def init_rapidtide_wf(
 
     Inputs
     ------
-    bold_std
+    bold
         BOLD series in template space
-    bold_mask_std
+    bold_mask
         BOLD series mask in template space
-    dseg_std
+    dseg
         Tissue segmentation in template space
     confounds
         fMRIPrep-formatted confounds file, which must include the following columns:
@@ -91,20 +90,21 @@ def init_rapidtide_wf(
     from nipype.interfaces.base import Undefined
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
+    from fmripost_rapidtide.interfaces.bids import DerivativesDataSink
     from fmripost_rapidtide.interfaces.nilearn import SplitDseg
+    from fmripost_rapidtide.interfaces.rapidtide import Rapidtide
 
     workflow = Workflow(name=_get_wf_name(bold_file, 'rapidtide'))
     workflow.__postdesc__ = """\
-Automatic removal of motion artifacts using independent component analysis
-[Rapidtide, @rapidtide] was performed on the *preprocessed BOLD on MNI152NLin6Asym space*.
+Identification and removal of traveling wave artifacts was performed using rapidtide.
 """
 
     inputnode = pe.Node(
         niu.IdentityInterface(
             fields=[
-                'bold_std',
-                'bold_mask_std',
-                'dseg_std',
+                'bold',
+                'bold_mask',
+                'dseg',
                 'confounds',
                 'skip_vols',
             ],
@@ -117,8 +117,6 @@ def init_rapidtide_wf(
             fields=[
                 'delay_map',
                 'regressor_file',
-                'confound_regressed',
-                'denoised',
             ],
         ),
         name='outputnode',
@@ -129,47 +127,46 @@ def init_rapidtide_wf(
         SplitDseg(),
         name='split_tissues',
     )
-    workflow.connect([(inputnode, split_tissues, [('dseg_std', 'dseg')])])
+    workflow.connect([(inputnode, split_tissues, [('dseg', 'dseg')])])
 
     # Run the Rapidtide classifier
     # XXX: simcalcrange is converted to list of strings
     rapidtide = pe.Node(
         Rapidtide(
-            outputname='rapidtide',
-            datatstep=metadata['RepetitionTime'],
+            ampthresh=config.workflow.ampthresh,
+            autorespdelete=config.workflow.autorespdelete,
             autosync=config.workflow.autosync,
+            bipolar=config.workflow.bipolar,
+            confoundderiv=config.workflow.confoundderiv,
+            confoundfile=config.workflow.confoundfile or Undefined,
+            confoundpowers=config.workflow.confoundpowers,
+            convergencethresh=config.workflow.convergencethresh or Undefined,
+            corrweighting=config.workflow.corrweighting,
+            datatstep=metadata['RepetitionTime'],
+            detrendorder=config.workflow.detrendorder,
             filterband=config.workflow.filterband,
             filterfreqs=config.workflow.filterfreqs or Undefined,
             filterstopfreqs=config.workflow.filterstopfreqs or Undefined,
-            numnull=config.workflow.numnull,
-            detrendorder=config.workflow.detrendorder,
-            spatialfilt=config.workflow.spatialfilt,
-            confoundfile=config.workflow.confoundfile or Undefined,
-            confoundpowers=config.workflow.confoundpowers,
-            confoundderiv=config.workflow.confoundderiv,
-            globalsignalmethod=config.workflow.globalsignalmethod,
+            fixdelay=config.workflow.fixdelay or Undefined,
+            glmderivs=config.workflow.glmderivs,
+            glmsourcefile=config.workflow.glmsourcefile or Undefined,
             globalpcacomponents=config.workflow.globalpcacomponents,
+            globalsignalmethod=config.workflow.globalsignalmethod,
+            lagmaxthresh=config.workflow.lagmaxthresh,
+            lagminthresh=config.workflow.lagminthresh,
+            maxpasses=config.workflow.maxpasses,
+            nprocs=config.nipype.omp_nthreads,
+            numnull=config.workflow.numnull,
             numtozero=config.workflow.numtozero,
-            timerange=[int(float(i)) for i in config.workflow.timerange],
-            corrweighting=config.workflow.corrweighting,
-            simcalcrange=[int(float(i)) for i in config.workflow.simcalcrange],
-            fixdelay=config.workflow.fixdelay or Undefined,
+            outputlevel=config.workflow.outputlevel,
+            pcacomponents=config.workflow.pcacomponents,
             searchrange=[int(float(i)) for i in config.workflow.searchrange],
             sigmalimit=config.workflow.sigmalimit,
-            bipolar=config.workflow.bipolar,
-            lagminthresh=config.workflow.lagminthresh,
-            lagmaxthresh=config.workflow.lagmaxthresh,
-            ampthresh=config.workflow.ampthresh,
             sigmathresh=config.workflow.sigmathresh,
-            pcacomponents=config.workflow.pcacomponents,
-            convergencethresh=config.workflow.convergencethresh or Undefined,
-            maxpasses=config.workflow.maxpasses,
-            glmsourcefile=config.workflow.glmsourcefile or Undefined,
-            glmderivs=config.workflow.glmderivs,
-            outputlevel=config.workflow.outputlevel,
+            simcalcrange=[int(float(i)) for i in config.workflow.simcalcrange],
+            spatialfilt=config.workflow.spatialfilt,
             territorymap=config.workflow.territorymap or Undefined,
-            autorespdelete=config.workflow.autorespdelete,
-            nprocs=config.nipype.omp_nthreads,
+            timerange=[int(float(i)) for i in config.workflow.timerange],
         ),
         name='rapidtide',
         mem_gb=mem_gb['filesize'] * 6,
@@ -177,8 +174,8 @@ def init_rapidtide_wf(
     )
     workflow.connect([
         (inputnode, rapidtide, [
-            ('bold_std', 'in_file'),
-            ('bold_mask_std', 'brainmask'),
+            ('bold', 'in_file'),
+            ('bold_mask', 'brainmask'),
             ('confounds', 'motionfile'),
             ('skip_vols', 'numskip'),
         ]),
@@ -189,13 +186,278 @@ def init_rapidtide_wf(
             ('gm', 'refineinclude'),  # GM mask for refinement
             ('gm', 'offsetinclude'),  # GM mask for offset calculation
         ]),
-        (rapidtide, outputnode, [
-            ('delay_map', 'delay_map'),
-            ('regressor_file', 'regressor_file'),
-            ('denoised', 'denoised'),
-        ])
     ])  # fmt:skip
 
-    # Generate figures for report
+    ds_delay_map = pe.Node(
+        DerivativesDataSink(
+            compress=True,
+            desc='delay',
+            suffix='boldmap',
+        ),
+        name='ds_delay_map',
+        run_without_submitting=True,
+    )
+    workflow.connect([
+        (rapidtide, ds_delay_map, [
+            ('maxtimemap', 'in_file'),
+            ('maxtimemap_json', 'meta_dict'),
+        ]),
+        (ds_delay_map, outputnode, [('out_file', 'delay_map')]),
+    ])  # fmt:skip
+
+    ds_regressor = pe.Node(
+        DerivativesDataSink(
+            desc='regressor',
+            suffix='timeseries',
+            extension='.tsv',
+        ),
+        name='ds_regressor',
+        run_without_submitting=True,
+    )
+    workflow.connect([
+        (rapidtide, ds_regressor, [
+            ('lagtcgenerator', 'in_file'),
+            ('lagtcgenerator_json', 'meta_dict'),
+        ]),
+        (ds_regressor, outputnode, [('out_file', 'regressor_file')]),
+    ])  # fmt:skip
+
+    return workflow
+
+
+def init_rapidtide_denoise_wf(
+    *,
+    bold_file: str,
+    metadata: dict,
+    mem_gb: dict,
+):
+    """Build a workflow that runs `Rapidtide`_.
+
+    This workflow wraps `Rapidtide`_ to characterize and remove the traveling wave artifact.
+
+    The following steps are performed:
+
+    #. Remove non-steady state volumes from the bold series.
+    #. Run rapidtide
+    #. Collect rapidtide outputs
+    #. Generate a confounds file with the rapidtide outputs
+
+    .. _Rapidtide: https://rapidtide.readthedocs.io/
+
+    Workflow Graph
+        .. workflow::
+            :graph2use: orig
+            :simple_form: yes
+
+            from fmripost_rapidtide.workflows.rapidtide import init_rapidtide_fit_wf
+
+            wf = init_rapidtide_fit_wf(
+                bold_file="fake.nii.gz",
+                metadata={"RepetitionTime": 1.0},
+            )
+
+    Parameters
+    ----------
+    bold_file
+        BOLD series used as name source for derivatives
+    metadata : :obj:`dict`
+        BIDS metadata for BOLD file
+
+    Inputs
+    ------
+    bold
+        BOLD series in template space
+    bold_mask
+        BOLD series mask in template space
+    dseg
+        Tissue segmentation in template space
+    confounds
+        fMRIPrep-formatted confounds file, which must include the following columns:
+        "trans_x", "trans_y", "trans_z", "rot_x", "rot_y", "rot_z".
+    skip_vols
+        number of non steady state volumes
+
+    Outputs
+    -------
+    denoised_bold
+    confounds_file
+    """
+
+    from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+
+    from fmripost_rapidtide.interfaces.bids import DerivativesDataSink
+    from fmripost_rapidtide.interfaces.rapidtide import RetroGLM
+
+    workflow = Workflow(name=_get_wf_name(bold_file, 'denoise'))
+    workflow.__postdesc__ = """\
+Identification and removal of traveling wave artifacts was performed using rapidtide.
+"""
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'bold',
+                'bold_mask',
+                'rapidtide_dir',
+                'skip_vols',
+            ],
+        ),
+        name='inputnode',
+    )
+
+    # Remove the traveling wave artifact
+    retroglm = pe.Node(
+        RetroGLM(
+            nprocs=config.nipype.omp_nthreads,
+            glmderivs=config.workflow.glmderivs,
+        ),
+        name='retroglm',
+        mem_gb=mem_gb['filesize'] * 6,
+        n_procs=config.nipype.omp_nthreads,
+    )
+    workflow.connect([
+        (inputnode, retroglm, [
+            ('bold', 'in_file'),
+            ('bold_mask', 'brainmask'),
+            ('rapidtide_dir', 'datafileroot'),
+            ('skip_vols', 'numskip'),
+        ]),
+    ])  # fmt:skip
+
+    ds_denoised_bold = pe.Node(
+        DerivativesDataSink(
+            compress=True,
+            desc='denoised',
+            suffix='bold',
+        ),
+        name='ds_denoised_bold',
+        run_without_submitting=True,
+    )
+    workflow.connect([
+        (retroglm, ds_denoised_bold, [
+            ('denoised', 'in_file'),
+            ('denoised_json', 'meta_dict'),
+        ]),
+    ])  # fmt:skip
+
+    return workflow
+
+
+def init_rapidtide_confounds_wf(
+    *,
+    bold_file: str,
+    metadata: dict,
+    mem_gb: dict,
+):
+    """Build a workflow that runs `Rapidtide`_.
+
+    This workflow wraps `Rapidtide`_ to characterize and remove the traveling wave artifact.
+
+    The following steps are performed:
+
+    #. Remove non-steady state volumes from the bold series.
+    #. Run rapidtide
+    #. Collect rapidtide outputs
+    #. Generate a confounds file with the rapidtide outputs
+
+    .. _Rapidtide: https://rapidtide.readthedocs.io/
+
+    Workflow Graph
+        .. workflow::
+            :graph2use: orig
+            :simple_form: yes
+
+            from fmripost_rapidtide.workflows.rapidtide import init_rapidtide_fit_wf
+
+            wf = init_rapidtide_fit_wf(
+                bold_file="fake.nii.gz",
+                metadata={"RepetitionTime": 1.0},
+            )
+
+    Parameters
+    ----------
+    bold_file
+        BOLD series used as name source for derivatives
+    metadata : :obj:`dict`
+        BIDS metadata for BOLD file
+
+    Inputs
+    ------
+    bold
+        BOLD series in template space
+    bold_mask
+        BOLD series mask in template space
+    dseg
+        Tissue segmentation in template space
+    confounds
+        fMRIPrep-formatted confounds file, which must include the following columns:
+        "trans_x", "trans_y", "trans_z", "rot_x", "rot_y", "rot_z".
+    skip_vols
+        number of non steady state volumes
+
+    Outputs
+    -------
+    denoised_bold
+    confounds_file
+    """
+
+    from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+
+    from fmripost_rapidtide.interfaces.bids import DerivativesDataSink
+    from fmripost_rapidtide.interfaces.rapidtide import RetroLagTCS
+
+    workflow = Workflow(name=_get_wf_name(bold_file, 'denoise'))
+    workflow.__postdesc__ = """\
+Identification and removal of traveling wave artifacts was performed using rapidtide.
+"""
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'bold',
+                'bold_mask',
+                'delay_map',
+                'regressor',
+                'skip_vols',
+            ],
+        ),
+        name='inputnode',
+    )
+
+    # Generate the traveling wave artifact voxel-wise regressor
+    retrolagtcs = pe.Node(
+        RetroLagTCS(
+            nprocs=config.nipype.omp_nthreads,
+            glmderivs=config.workflow.glmderivs,
+        ),
+        name='retrolagtcs',
+        mem_gb=mem_gb['filesize'] * 6,
+        n_procs=config.nipype.omp_nthreads,
+    )
+    workflow.connect([
+        (inputnode, retrolagtcs, [
+            ('bold', 'in_file'),
+            ('bold_mask', 'maskfile'),
+            ('delay_map', 'lagtimesfile'),
+            ('regressor', 'lagtcgeneratorfile'),
+            ('skip_vols', 'numskip'),
+        ]),
+    ])  # fmt:skip
+
+    ds_voxelwise_regressor = pe.Node(
+        DerivativesDataSink(
+            compress=True,
+            desc='sLFO',
+            suffix='timeseries',
+        ),
+        name='ds_voxelwise_regressor',
+        run_without_submitting=True,
+    )
+    workflow.connect([
+        (retrolagtcs, ds_voxelwise_regressor, [
+            ('denoised', 'in_file'),
+            ('denoised_json', 'meta_dict'),
+        ]),
+    ])  # fmt:skip
 
     return workflow
diff --git a/tests/data/ds000005-fmriprep b/tests/data/ds000005-fmriprep
new file mode 160000
index 0000000..caeb5ac
--- /dev/null
+++ b/tests/data/ds000005-fmriprep
@@ -0,0 +1 @@
+Subproject commit caeb5acc567ddc94846bfa9e801088a7be3085c0