diff --git a/aslprep/cli/aggregate_qc.py b/aslprep/cli/aggregate_qc.py
index e2e730460..15679349f 100644
--- a/aslprep/cli/aggregate_qc.py
+++ b/aslprep/cli/aggregate_qc.py
@@ -2,46 +2,307 @@
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """Aggregate QC measures across all subjects in dataset."""
 
-import os
 from pathlib import Path
 
-import pandas as pd
+from aslprep._warnings import logging
 
 
-def get_parser():
+def _build_parser():
     """Build parser object."""
     from argparse import ArgumentParser, RawTextHelpFormatter
+    from functools import partial
+
+    def _path_exists(path, parser):
+        """Ensure a given path exists."""
+        if path is None or not Path(path).exists():
+            raise parser.error(f'Path does not exist: <{path}>.')
+        return Path(path).absolute()
 
     parser = ArgumentParser(description=__doc__, formatter_class=RawTextHelpFormatter)
 
-    parser.add_argument('aslprep_dir', action='store', type=Path, help='aslprep output dir')
+    PathExists = partial(_path_exists, parser=parser)
+
+    parser.add_argument(
+        'aslprep_dir',
+        action='store',
+        type=PathExists,
+        help='Path to ASLPrep derivatives dataset.',
+    )
+    parser.add_argument(
+        'output_dir',
+        action='store',
+        type=Path,
+        help='Path to output directory where QC measures will be saved.',
+    )
+    parser.add_argument(
+        'analysis_level',
+        choices=['group'],
+        help='processing stage to be run, only "group" in the case of this CLI.',
+    )
+
+    g_other = parser.add_argument_group('Other options')
+    g_other.add_argument(
+        '-w',
+        '--work-dir',
+        action='store',
+        type=Path,
+        default=Path('work').absolute(),
+        help='path where intermediate results should be stored',
+    )
 
-    parser.add_argument('output_prefix', action='store', type=str, help='output prefix for group')
     return parser
 
 
+def parse_args(args=None, namespace=None):
+    """Parse args and run further checks on the command line."""
+    parser = _build_parser()
+    opts = parser.parse_args(args, namespace)
+
+    output_dir = opts.output_dir
+    input_dir = opts.aslprep_dir
+    work_dir = opts.work_dir
+
+    output_dir.mkdir(exist_ok=True, parents=True)
+    work_dir.mkdir(exist_ok=True, parents=True)
+
+    log_dir = output_dir / 'logs'
+    log_dir.mkdir(exist_ok=True, parents=True)
+
+    build_log = logging.getLogger()
+
+    # Wipe out existing work_dir
+    if opts.clean_workdir and work_dir.exists():
+        from niworkflows.utils.misc import clean_directory
+
+        build_log.info(f'Clearing previous aslprep working directory: {work_dir}')
+        if not clean_directory(work_dir):
+            build_log.warning(f'Could not clear all contents of working directory: {work_dir}')
+
+    return opts
+
+
+def build_workflow(opts, retval):
+    """Build the workflow."""
+    from bids.layout import BIDSLayout
+    from nipype.interfaces import utility as niu
+    from nipype.pipeline import engine as pe
+    from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+
+    from aslprep.data import load as load_data
+    from aslprep.interfaces.bids import DerivativesDataSink
+    from aslprep.interfaces.group import (
+        AggregateCBFQC,
+        MeanCBF,
+        PlotAggregatedCBFQC,
+        StatisticalMapRPT,
+    )
+
+    aslprep_config = load_data('aslprep_bids_config.json')
+    layout = BIDSLayout(
+        opts.aslprep_dir,
+        validate=False,
+        config=['bids', 'derivatives', aslprep_config],
+    )
+    qc_data = collect_aslprep_qc_derivatives(layout)
+
+    workflow = Workflow(name='aslprep_group_wf')
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'cbf_qc',
+                'cbf_img',
+            ],
+        ),
+        name='inputnode',
+    )
+    inputnode.inputs.cbf_qc = qc_data['cbf_qc']
+    inputnode.inputs.cbf_img = qc_data['cbf_img']
+
+    # Aggregate QC metrics
+    aggregate_cbf_qc = pe.Node(
+        AggregateCBFQC(),
+        name='aggregate_cbf_qc',
+    )
+    workflow.connect([(inputnode, aggregate_cbf_qc, [('cbf_qc', 'in_files')])])
+
+    # Plot aggregated QC metrics
+    plot_cbf_qc = pe.Node(
+        PlotAggregatedCBFQC(),
+        name='plot_cbf_qc',
+    )
+    workflow.connect([(aggregate_cbf_qc, plot_cbf_qc, [('out_file', 'in_file')])])
+
+    # Save aggregated QC metrics
+    ds_report_qei = pe.MapNode(
+        DerivativesDataSink(
+            base_directory=opts.output_dir,
+            datatype='figures',
+            suffix='cbf',
+        ),
+        iterfield=['in_file', 'desc'],
+        name='ds_report_qei',
+    )
+    workflow.connect([
+        (plot_cbf_qc, ds_report_qei, [
+            ('out_files', 'in_file'),
+            ('measures', 'desc'),
+        ]),
+    ])  # fmt:skip
+
+    # Average CBF maps across subjects
+    average_cbf_maps = pe.Node(
+        MeanCBF(),
+        name='average_cbf_maps',
+    )
+    workflow.connect([(inputnode, average_cbf_maps, [('cbf_img', 'in_files')])])
+
+    # Plot average and SD CBF maps
+    plot_mean_cbf = pe.Node(
+        StatisticalMapRPT(cmap='viridis'),
+        name='plot_mean_cbf',
+    )
+    workflow.connect([
+        (inputnode, plot_mean_cbf, [('template', 'underlay')]),
+        (average_cbf_maps, plot_mean_cbf, [('mean_map', 'overlay')]),
+    ])  # fmt:skip
+
+    ds_report_mean_cbf = pe.Node(
+        DerivativesDataSink(
+            base_directory=opts.output_dir,
+            datatype='figures',
+            desc='mean',
+            suffix='cbf',
+        ),
+        name='ds_report_mean_cbf',
+    )
+    workflow.connect([(plot_mean_cbf, ds_report_mean_cbf, [('out_report', 'in_file')])])
+
+    plot_sd_cbf = pe.Node(
+        StatisticalMapRPT(cmap='viridis'),
+        name='plot_sd_cbf',
+    )
+    workflow.connect([
+        (inputnode, plot_sd_cbf, [('template', 'underlay')]),
+        (average_cbf_maps, plot_sd_cbf, [('sd_map', 'overlay')]),
+    ])  # fmt:skip
+
+    ds_report_sd_cbf = pe.Node(
+        DerivativesDataSink(
+            base_directory=opts.output_dir,
+            datatype='figures',
+            desc='sd',
+            suffix='cbf',
+        ),
+        name='ds_report_sd_cbf',
+    )
+    workflow.connect([(plot_sd_cbf, ds_report_sd_cbf, [('out_report', 'in_file')])])
+
+    return workflow
+
+
+def collect_aslprep_qc_derivatives(layout):
+    """Collect ASLPrep QC derivatives."""
+    qc_data = {}
+    qc_files = layout.get(
+        desc='qualitycontrol',
+        suffix='cbf',
+        extension='.tsv',
+        return_type='file',
+    )
+    qc_data['cbf_qc'] = qc_files
+    for qc_file in qc_files:
+        # Collect associated files
+        # Standard-space CBF map
+        # We also need the space used so we can grab the right template
+        ...
+
+    return qc_data
+
+
+def build_boilerplate(opts, aslprep_wf):
+    """Build the boilerplate."""
+    pass
+
+
 def main():
-    """Run the workflow."""
-    opts = get_parser().parse_args()
+    """Entry point."""
+    import gc
+    import sys
+    from multiprocessing import Manager, Process
+    from os import EX_SOFTWARE
+
+    from aslprep.utils.bids import write_bidsignore, write_derivative_description
+
+    opts = parse_args()
+    input_dir = opts.aslprep_dir
+    output_dir = opts.output_dir
+
+    logger = logging.getLogger()
+
+    # CRITICAL Call build_workflow(config_file, retval) in a subprocess.
+    # Because Python on Linux does not ever free virtual memory (VM), running the
+    # workflow construction jailed within a process preempts excessive VM buildup.
+    with Manager() as mgr:
+        retval = mgr.dict()
+        p = Process(target=build_workflow, args=(opts, retval))
+        p.start()
+        p.join()
+        retval = dict(retval.items())  # Convert to base dictionary
+
+        if p.exitcode:
+            retval['return_code'] = p.exitcode
+
+    exitcode = retval.get('return_code', 0)
+    aslprep_wf = retval.get('workflow', None)
+
+    exitcode = exitcode or (aslprep_wf is None) * EX_SOFTWARE
+    if exitcode != 0:
+        sys.exit(exitcode)
+
+    # Generate boilerplate
+    with Manager() as mgr:
+        p = Process(target=build_boilerplate, args=(opts, aslprep_wf))
+        p.start()
+        p.join()
+
+    # Clean up master process before running workflow, which may create forks
+    gc.collect()
 
-    allsubj_dir = os.path.abspath(opts.aslprep_dir)
-    outputfile = os.getcwd() + '/' + str(opts.output_prefix) + '_allsubjects_qc.tsv'
+    logger.log(25, 'ASLPrep-Group started!')
+    errno = 1  # Default is error exit unless otherwise set
+    try:
+        aslprep_wf.run({'plugin': 'Linear'})
+    except Exception as e:
+        logger.critical('ASLPrep-Group failed: %s', e)
+        raise
+    else:
+        logger.log(25, 'ASLPrep-Group finished successfully!')
+        errno = 0
+    finally:
+        from aslprep import data
+        from aslprep.reports.core import generate_group_report
 
-    qclist = []
-    for r, _, f in os.walk(allsubj_dir):
-        for filex in f:
-            if filex.endswith('desc-qualitycontrol_cbf.tsv'):
-                qclist.append(r + '/' + filex)
+        # Generate reports phase
+        failed_reports = generate_group_report(
+            output_dir,
+            bootstrap_file=data.load('reports-spec-group.yml'),
+        )
+        write_derivative_description(input_dir, output_dir)
+        write_bidsignore(output_dir)
 
-    datax = pd.read_table(qclist[0])
-    for i in range(1, len(qclist)):
-        dy = pd.read_table(qclist[i])
-        datax = pd.concat([datax, dy])
+        if failed_reports:
+            msg = (
+                'Report generation was not successful for the following participants '
+                f': {", ".join(failed_reports)}.'
+            )
+            logger.error(msg)
 
-    datax.to_csv(outputfile, index=None, sep='\t')
+        sys.exit(int((errno + len(failed_reports)) > 0))
 
 
 if __name__ == '__main__':
     raise RuntimeError(
-        'this should be run after running aslprep;\nit required installation of  aslprep'
+        'aslprep/cli/aggregate_qc.py should not be run directly;\n'
+        'Please `pip install` aslprep and use the `aslprep-group` command'
     )
diff --git a/aslprep/interfaces/group.py b/aslprep/interfaces/group.py
new file mode 100644
index 000000000..3419e7113
--- /dev/null
+++ b/aslprep/interfaces/group.py
@@ -0,0 +1,244 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""Interfaces for calculating and collecting confounds."""
+
+import os
+
+import matplotlib.pyplot as plt
+import pandas as pd
+import seaborn as sns
+from nipype.interfaces.base import (
+    BaseInterfaceInputSpec,
+    File,
+    SimpleInterface,
+    TraitedSpec,
+    isdefined,
+    traits,
+)
+
+
+class _AggregateCBFQCInputSpec(BaseInterfaceInputSpec):
+    in_files = traits.List(File(exists=True), desc='list of QC tsv files')
+
+
+class _AggregateCBFQCOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='aggregated QC file')
+
+
+class AggregateCBFQC(SimpleInterface):
+    """Aggregate CBF QC metrics across runs, sessions, and subjects."""
+
+    input_spec = _AggregateCBFQCInputSpec
+    output_spec = _AggregateCBFQCOutputSpec
+
+    def _run_interface(self, runtime):
+        in_files = self.inputs.in_files
+        dfs = []
+        for in_file in in_files:
+            temp_df = pd.read_table(in_file)
+            temp_df['source_file'] = os.path.basename(in_file)
+            dfs.append(temp_df)
+
+        df = pd.concat(dfs)
+
+        self._results['out_file'] = os.path.join(runtime.cwd, 'aggregated_qc.tsv')
+        df.to_csv(self._results['out_file'], index=None, sep='\t')
+
+        return runtime
+
+
+class _PlotAggregatedCBFQCInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, desc='aggregated QC file')
+
+
+class _PlotAggregatedCBFQCOutputSpec(TraitedSpec):
+    out_files = traits.List(File(exists=True), desc='list of output figures')
+    measures = traits.List(traits.String, desc='list of measures that have been plotted')
+
+
+class PlotAggregatedCBFQC(SimpleInterface):
+    """Plot aggregated CBF QC metrics across runs, sessions, and subjects."""
+
+    input_spec = _PlotAggregatedCBFQCInputSpec
+    output_spec = _PlotAggregatedCBFQCOutputSpec
+
+    def _run_interface(self, runtime):
+        from bids.layout import parse_file_entities
+
+        df = pd.read_table(self.inputs.in_file)
+        all_filenames = df['source_file'].unique()
+        all_entities = [list(parse_file_entities(f).keys()) for f in all_filenames]
+        # flatten list of lists
+        all_entities = [item for sublist in all_entities for item in sublist]
+        all_entities = list(set(all_entities))
+
+        cols_to_ignore = all_entities + ['source_file']
+        measures = [c for c in df.columns if c not in cols_to_ignore]
+
+        entities_to_ignore = ['subject', 'session', 'run']
+        entities_to_group = [e for e in all_entities if e not in entities_to_ignore]
+
+        for measure in measures:
+            figure_file = self.plot_aggregated_qc_metrics(
+                df, measure, entities_to_group, runtime.cwd
+            )
+            self._results['out_files'].append(figure_file)
+            self._results['measures'].append(measure)
+
+        return runtime
+
+    def plot_aggregated_qc_metrics(self, df, measure, entity_columns, cwd):
+        """Plot aggregated QC metrics for a given measure and entity columns."""
+        # Identify entities with varying values in dataset
+        grouping_columns = []
+        for col in entity_columns:
+            if df[col].unique().size > 1:
+                grouping_columns.append(col)
+                df[col] = df[col].fillna('NONE')
+
+        # Create a new column that indexes combinations of grouping columns
+        df['group'] = ''
+        for col in grouping_columns:
+            df['group'] += col + '-' + df[col] + '_'
+
+        # Remove trailing underscore
+        df['group'] = df['group'].str.rstrip('_')
+
+        kwargs = {}
+        if len(grouping_columns):
+            kwargs = {'hue': 'group'}
+
+        fig, ax = plt.subplots()
+        sns.stripplot(data=df, y=measure, ax=ax, figure=fig, **kwargs)
+        figure_file = os.path.join(cwd, f'{measure}.svg')
+        fig.savefig(figure_file)
+        return figure_file
+
+
+class _MeanCBFInputSpec(BaseInterfaceInputSpec):
+    in_files = traits.List(File(exists=True), desc='list of CBF files')
+
+
+class _MeanCBFOutputSpec(TraitedSpec):
+    mean_map = File(exists=True, desc='mean CBF map')
+    sd_map = File(exists=True, desc='SD CBF map')
+
+
+class MeanCBF(SimpleInterface):
+    """Calculate mean and SD CBF maps dataset."""
+
+    input_spec = _MeanCBFInputSpec
+    output_spec = _MeanCBFOutputSpec
+
+    def _run_interface(self, runtime):
+        from nilearn import image
+
+        mean_map = image.mean_img(self.inputs.in_files)
+        sd_map = image.math_img('np.std(img, axis=3)', img=self.inputs.in_files)
+
+        mean_map_file = os.path.join(runtime.cwd, 'mean_cbf.nii.gz')
+        sd_map_file = os.path.join(runtime.cwd, 'sd_cbf.nii.gz')
+
+        mean_map.to_filename(mean_map_file)
+        sd_map.to_filename(sd_map_file)
+
+        self._results['mean_map'] = mean_map_file
+        self._results['sd_map'] = sd_map_file
+
+        return runtime
+
+
+class _StatisticalMapInputSpecRPT(BaseInterfaceInputSpec):
+    overlay = File(
+        exists=True,
+        mandatory=True,
+        desc='FC inflation time series',
+    )
+    underlay = File(
+        exists=True,
+        mandatory=True,
+        desc='Underlay image',
+    )
+    mask = File(
+        exists=True,
+        mandatory=False,
+        desc='Mask image',
+    )
+    cmap = traits.Str(
+        'viridis',
+        desc='Colormap',
+        usedefault=True,
+    )
+    out_report = File(
+        'statistical_map_report.svg',
+        usedefault=True,
+        desc='Filename for the visual report generated by Nipype.',
+    )
+
+
+class _StatisticalMapOutputSpecRPT(TraitedSpec):
+    out_report = File(
+        exists=True,
+        desc='Filename for the visual report generated by Nipype.',
+    )
+
+
+class StatisticalMapRPT(SimpleInterface):
+    """Create a reportlet for Rapidtide outputs."""
+
+    input_spec = _StatisticalMapInputSpecRPT
+    output_spec = _StatisticalMapOutputSpecRPT
+
+    def _run_interface(self, runtime):
+        from uuid import uuid4
+
+        from nilearn import image, masking, plotting
+        from nireports._vendored.svgutils.transform import fromstring
+        from nireports.reportlets.utils import compose_view, cuts_from_bbox, extract_svg
+
+        out_file = os.path.abspath(self.inputs.out_report)
+
+        if isdefined(self.inputs.mask):
+            mask_img = image.load_img(self.inputs.mask)
+            overlay_img = masking.unmask(
+                masking.apply_mask(self.inputs.overlay, self.inputs.mask),
+                self.inputs.mask,
+            )
+            # since the moving image is already in the fixed image space we
+            # should apply the same mask
+            underlay_img = image.load_img(self.inputs.underlay)
+        else:
+            overlay_img = image.load_img(self.inputs.overlay)
+            underlay_img = image.load_img(self.inputs.underlay)
+            mask_img = image.threshold_img(overlay_img, 1e-3)
+
+        n_cuts = 7
+        cuts = cuts_from_bbox(mask_img, cuts=n_cuts)
+        order = ('z', 'x', 'y')
+        out_files = []
+
+        # Plot each cut axis
+        plot_params = {}
+        for mode in list(order):
+            plot_params['display_mode'] = mode
+            plot_params['cut_coords'] = cuts[mode]
+            plot_params['title'] = None
+            plot_params['cmap'] = self.inputs.cmap
+
+            # Generate nilearn figure
+            display = plotting.plot_stat_map(
+                overlay_img,
+                bg_img=underlay_img,
+                **plot_params,
+            )
+
+            svg = extract_svg(display, compress=False)
+            display.close()
+
+            # Find and replace the figure_1 id.
+            svg = svg.replace('figure_1', f'{mode}-{uuid4()}', 1)
+            out_files.append(fromstring(svg))
+
+        compose_view(bg_svgs=out_files, fg_svgs=None, out_file=out_file)
+        self._results['out_report'] = out_file
+        return runtime