Add carpet plots to output (#49)

* Try to get carpetplots working.

* Update aroma.py

* Keep working.

* Update aroma.py

* Update env.yml

* Update fetch_templates.py

* Simplify things.

* Fix things up a bit.

* Update base.py

* Update base.py

* Work on carpet plots.

env.yml

@@ -29,6 +29,9 @@ dependencies:
   - fsl-melodic=2111.3
   - fsl-miscmaths=2203.2
   - fsl-susan=2111.0
+  # Workflow dependencies: ANTs
+  - ants=2.5
+  # Other dependencies
   - pip
   - pip:
     - templateflow==24.0

scripts/fetch_templates.py

@@ -25,11 +25,11 @@ def fetch_MNI2009():
     """
     template = 'MNI152NLin2009cAsym'
 
-    tf.get(template, resolution=(1, 2), desc=None, suffix='T1w')
-    tf.get(template, resolution=(1, 2), desc='brain', suffix='mask')
+    # tf.get(template, resolution=(1, 2), desc=None, suffix='T1w')
+    # tf.get(template, resolution=(1, 2), desc='brain', suffix='mask')
     tf.get(template, resolution=1, atlas=None, desc='carpet', suffix='dseg')
-    tf.get(template, resolution=2, desc='fMRIPrep', suffix='boldref')
-    tf.get(template, resolution=1, label='brain', suffix='probseg')
+    # tf.get(template, resolution=2, desc='fMRIPrep', suffix='boldref')
+    # tf.get(template, resolution=1, label='brain', suffix='probseg')
 
 
 def fetch_MNI6():
@@ -48,6 +48,10 @@ def fetch_MNI6():
     tf.get(template, resolution=(1, 2), desc='brain', suffix='mask')
     # CIFTI
     tf.get(template, resolution=2, atlas='HCP', suffix='dseg')
+    # Transform from MNI152NLin2009cAsym to MNI152NLin6Asym
+    tf.get(
+        template, mode='image', suffix='xfm', extension='.h5', **{'from': 'MNI152NLin2009cAsym'}
+    )
 
 
 def fetch_OASIS():
@@ -108,7 +112,7 @@ def fetch_fsLR():
 
 
 def fetch_all():
-    # fetch_MNI2009()
+    fetch_MNI2009()
     fetch_MNI6()
     # fetch_OASIS()
     # fetch_fsaverage()

src/fmripost_aroma/cli/workflow.py

@@ -90,11 +90,12 @@ def build_workflow(config_file, retval):
     # Called with reports only
     if config.execution.reports_only:
         build_log.log(25, 'Running --reports-only on participants %s', ', '.join(subject_list))
-        retval['return_code'] = generate_reports(
+        failed_reports = generate_reports(
             subject_list=config.execution.participant_label,
             output_dir=config.execution.output_dir,
             run_uuid=config.execution.run_uuid,
         )
+        retval['return_code'] = len(failed_reports)
         return retval
 
     # Build main workflow

src/fmripost_aroma/data/reports-spec-func.yml

@@ -21,11 +21,14 @@ sections:
       anatomical white matter mask, which appears as a red contour.
     static: false
     subtitle: Alignment of functional and anatomical MRI data (coregistration)
-  - bids:
-      datatype: figures
-      desc: '[aggr,nonaggr,orthaggr]carpetplot'
-      extension: [.html]
-      suffix: bold
+  - bids: {datatype: figures, desc: preprocCarpetplot, suffix: bold}
+    title: Preprocessed BOLD
+  - bids: {datatype: figures, desc: nonaggrCarpetplot, suffix: bold}
+    title: Non-Aggressively Denoised BOLD
+  - bids: {datatype: figures, desc: aggrCarpetplot, suffix: bold}
+    title: Aggressively Denoised BOLD
+  - bids: {datatype: figures, desc: orthaggrCarpetplot, suffix: bold}
+    title: Othogonalized Aggressively Denoised BOLD
 - name: About
   nested: true
   reportlets:

src/fmripost_aroma/data/reports-spec.yml

@@ -21,11 +21,14 @@ sections:
       anatomical white matter mask, which appears as a red contour.
     static: false
     subtitle: Alignment of functional and anatomical MRI data (coregistration)
-  - bids:
-      datatype: figures
-      desc: '[aggr,nonaggr,orthaggr]carpetplot'
-      extension: [.html]
-      suffix: bold
+  - bids: {datatype: figures, desc: preprocCarpetplot, suffix: bold}
+    title: Preprocessed BOLD
+  - bids: {datatype: figures, desc: nonaggrCarpetplot, suffix: bold}
+    title: Non-Aggressively Denoised BOLD
+  - bids: {datatype: figures, desc: aggrCarpetplot, suffix: bold}
+    title: Aggressively Denoised BOLD
+  - bids: {datatype: figures, desc: orthaggrCarpetplot, suffix: bold}
+    title: Othogonalized Aggressively Denoised BOLD
 - name: About
   nested: true
   reportlets:

src/fmripost_aroma/interfaces/reportlets.py

@@ -67,7 +67,7 @@
 ABOUT_TEMPLATE = """\t<ul>
 \t\t<li>fMRIPost-AROMA version: {version}</li>
 \t\t<li>fMRIPost-AROMA command: <code>{command}</code></li>
-\t\t<li>Date preprocessed: {date}</li>
+\t\t<li>Date postprocessed: {date}</li>
 \t</ul>
 </div>
 """

src/fmripost_aroma/tests/test_base.py

@@ -33,5 +33,8 @@ def test_init_denoise_wf(tmp_path_factory):
         config.execution.output_dir = tempdir / 'out'
         config.execution.work_dir = tempdir / 'work'
 
-        wf = init_denoise_wf(bold_file='sub-01_task-rest_bold.nii.gz')
+        wf = init_denoise_wf(
+            bold_file='sub-01_task-rest_bold.nii.gz',
+            metadata={'RepetitionTime': 2.0},
+        )
         assert wf.name == 'denoise_task_rest_wf'

src/fmripost_aroma/workflows/aroma.py

@@ -388,7 +388,7 @@ def init_ica_aroma_wf(
     return workflow
 
 
-def init_denoise_wf(bold_file):
+def init_denoise_wf(bold_file, metadata):
     """Build a workflow that denoises a BOLD series using AROMA confounds.
 
     This workflow performs the denoising in the requested output space(s).
@@ -397,6 +397,7 @@ def init_denoise_wf(bold_file):
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
     from fmripost_aroma.interfaces.confounds import ICADenoise
+    from fmripost_aroma.workflows.confounds import init_carpetplot_wf
 
     workflow = Workflow(name=_get_wf_name(bold_file, 'denoise'))
 
@@ -405,6 +406,7 @@ def init_denoise_wf(bold_file):
             fields=[
                 'bold_file',
                 'bold_mask',
+                'confounds_file',
                 'mixing',
                 'classifications',
                 'skip_vols',
@@ -427,6 +429,22 @@ def init_denoise_wf(bold_file):
         ]),
     ])  # fmt:skip
 
+    preproc_carpetplot_wf = init_carpetplot_wf(
+        mem_gb=config.DEFAULT_MEMORY_MIN_GB,
+        metadata=metadata,
+        cifti_output=False,
+        name='preproc_carpet_wf',
+    )
+    preproc_carpetplot_wf.inputs.inputnode.desc = 'preprocCarpetplot'
+    workflow.connect([
+        (inputnode, preproc_carpetplot_wf, [
+            ('bold_mask', 'inputnode.bold_mask'),
+            ('confounds_file', 'inputnode.confounds_file'),
+            ('skip_vols', 'inputnode.dummy_scans'),
+        ]),
+        (inputnode, preproc_carpetplot_wf, [('bold_file', 'inputnode.bold')]),
+    ])  # fmt:skip
+
     for denoise_method in config.workflow.denoise_method:
         denoise = pe.Node(
             ICADenoise(method=denoise_method),
@@ -476,6 +494,22 @@ def init_denoise_wf(bold_file):
             (add_non_steady_state, ds_denoised, [('bold_add', 'in_file')]),
         ])  # fmt:skip
 
+        carpetplot_wf = init_carpetplot_wf(
+            mem_gb=config.DEFAULT_MEMORY_MIN_GB,
+            metadata=metadata,
+            cifti_output=False,
+            name=f'{denoise_method}_carpet_wf',
+        )
+        carpetplot_wf.inputs.inputnode.desc = f'{denoise_method}Carpetplot'
+        workflow.connect([
+            (inputnode, carpetplot_wf, [
+                ('bold_mask', 'inputnode.bold_mask'),
+                ('confounds_file', 'inputnode.confounds_file'),
+                ('skip_vols', 'inputnode.dummy_scans'),
+            ]),
+            (ds_denoised, carpetplot_wf, [('out_file', 'inputnode.bold')]),
+        ])  # fmt:skip
+
     return workflow
 
 

src/fmripost_aroma/workflows/base.py

@@ -201,15 +201,6 @@ def init_single_subject_wf(subject_id: str):
         # Patch standard-space BOLD files into 'bold' key
         subject_data['bold'] = listify(subject_data['bold_mni152nlin6asym'])
 
-    if not subject_data['bold_mni152nlin6asym']:
-        task_id = config.execution.task_id
-        raise RuntimeError(
-            f"No MNI152NLin6Asym:res-2 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. "
-            f"Please check your BIDS filters: {config.execution.bids_filters}."
-        )
-
     # Make sure we always go through these two checks
     if not subject_data['bold']:
         task_id = config.execution.task_id
@@ -406,10 +397,11 @@ def init_single_run_wf(bold_file):
 
     if config.workflow.denoise_method:
         # Now denoise the output-space BOLD data using ICA-AROMA
-        denoise_wf = init_denoise_wf(bold_file=bold_file)
+        denoise_wf = init_denoise_wf(bold_file=bold_file, metadata=bold_metadata)
         denoise_wf.inputs.inputnode.skip_vols = skip_vols
         denoise_wf.inputs.inputnode.space = 'MNI152NLin6Asym'
         denoise_wf.inputs.inputnode.res = '2'
+        denoise_wf.inputs.inputnode.confounds_file = functional_cache['confounds']
 
         workflow.connect([
             (mni6_buffer, denoise_wf, [
@@ -422,6 +414,12 @@ def init_single_run_wf(bold_file):
             ]),
         ])  # fmt:skip
 
+    # Fill-in datasinks seen so far
+    for node in workflow.list_node_names():
+        if node.split('.')[-1].startswith('ds_'):
+            workflow.get_node(node).inputs.base_directory = config.execution.output_dir
+            workflow.get_node(node).inputs.source_file = bold_file
+
     return workflow
 
 

src/fmripost_aroma/workflows/confounds.py

@@ -54,17 +54,11 @@ def init_carpetplot_wf(
     Inputs
     ------
     bold
-        BOLD image, after the prescribed corrections (STC, HMC and SDC)
-        when available.
+        BOLD image, in MNI152NLin6Asym space + 2mm resolution.
     bold_mask
-        BOLD series mask
+        BOLD series mask in same space as ``bold``.
     confounds_file
         TSV of all aggregated confounds
-    boldref2anat_xfm
-        Affine matrix that maps the BOLD reference space into alignment with
-        the anatomical (T1w) space
-    std2anat_xfm
-        ANTs-compatible affine-and-warp transform file
     cifti_bold
         BOLD image in CIFTI format, to be used in place of volumetric BOLD
     crown_mask
@@ -79,29 +73,26 @@ def init_carpetplot_wf(
     out_carpetplot
         Path of the generated SVG file
     """
+    from fmriprep.interfaces.confounds import FMRISummary
     from nipype.interfaces import utility as niu
     from nipype.pipeline import engine as pe
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
     from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms
     from templateflow.api import get as get_template
 
     from fmripost_aroma.config import DEFAULT_MEMORY_MIN_GB
-    from fmripost_aroma.interfaces import DerivativesDataSink
-    from fmripost_aroma.interfaces.confounds import FMRISummary
+    from fmripost_aroma.interfaces.bids import DerivativesDataSink
 
     inputnode = pe.Node(
         niu.IdentityInterface(
             fields=[
                 'bold',
                 'bold_mask',
                 'confounds_file',
-                'boldref2anat_xfm',
-                'std2anat_xfm',
                 'cifti_bold',
-                'crown_mask',
-                'acompcor_mask',
                 'dummy_scans',
-            ]
+                'desc',
+            ],
         ),
         name='inputnode',
     )
@@ -113,10 +104,12 @@ def init_carpetplot_wf(
         FMRISummary(
             tr=metadata['RepetitionTime'],
             confounds_list=[
-                ('global_signal', None, 'GS'),
-                ('csf', None, 'CSF'),
-                ('white_matter', None, 'WM'),
-                ('std_dvars', None, 'DVARS'),
+                ('trans_x', 'mm', 'x'),
+                ('trans_y', 'mm', 'y'),
+                ('trans_z', 'mm', 'z'),
+                ('rot_x', 'deg', 'pitch'),
+                ('rot_y', 'deg', 'roll'),
+                ('rot_z', 'deg', 'yaw'),
                 ('framewise_displacement', 'mm', 'FD'),
             ],
         ),
@@ -125,7 +118,6 @@ def init_carpetplot_wf(
     )
     ds_report_bold_conf = pe.Node(
         DerivativesDataSink(
-            desc='carpetplot',
             datatype='figures',
             extension='svg',
             dismiss_entities=('echo', 'den', 'res'),
@@ -137,10 +129,8 @@ def init_carpetplot_wf(
 
     parcels = pe.Node(niu.Function(function=_carpet_parcellation), name='parcels')
     parcels.inputs.nifti = not cifti_output
-    # List transforms
-    mrg_xfms = pe.Node(niu.Merge(2), name='mrg_xfms')
 
-    # Warp segmentation into EPI space
+    # Warp segmentation into MNI152NLin6Asym space
     resample_parc = pe.Node(
         ApplyTransforms(
             dimension=3,
@@ -153,7 +143,17 @@ def init_carpetplot_wf(
                     extension=['.nii', '.nii.gz'],
                 )
             ),
-            invert_transform_flags=[True, False],
+            transforms=[
+                str(
+                    get_template(
+                        template='MNI152NLin6Asym',
+                        mode='image',
+                        suffix='xfm',
+                        extension='.h5',
+                        **{'from': 'MNI152NLin2009cAsym'},
+                    ),
+                ),
+            ],
             interpolation='MultiLabel',
             args='-u int',
         ),
@@ -165,28 +165,22 @@ def init_carpetplot_wf(
         workflow.connect(inputnode, 'cifti_bold', conf_plot, 'in_cifti')
 
     workflow.connect([
-        (inputnode, mrg_xfms, [
-            ('boldref2anat_xfm', 'in1'),
-            ('std2anat_xfm', 'in2'),
-        ]),
         (inputnode, resample_parc, [('bold_mask', 'reference_image')]),
-        (inputnode, parcels, [('crown_mask', 'crown_mask')]),
-        (inputnode, parcels, [('acompcor_mask', 'acompcor_mask')]),
         (inputnode, conf_plot, [
             ('bold', 'in_nifti'),
             ('confounds_file', 'confounds_file'),
             ('dummy_scans', 'drop_trs'),
         ]),
-        (mrg_xfms, resample_parc, [('out', 'transforms')]),
         (resample_parc, parcels, [('output_image', 'segmentation')]),
         (parcels, conf_plot, [('out', 'in_segm')]),
+        (inputnode, ds_report_bold_conf, [('desc', 'desc')]),
         (conf_plot, ds_report_bold_conf, [('out_file', 'in_file')]),
         (conf_plot, outputnode, [('out_file', 'out_carpetplot')]),
     ])  # fmt:skip
     return workflow
 
 
-def _carpet_parcellation(segmentation, crown_mask, acompcor_mask, nifti=False):
+def _carpet_parcellation(segmentation, nifti=False):
     """Generate the union of two masks."""
     from pathlib import Path
 
@@ -202,9 +196,9 @@ def _carpet_parcellation(segmentation, crown_mask, acompcor_mask, nifti=False):
     lut[255] = 5 if nifti else 0  # Cerebellum
     # Apply lookup table
     seg = lut[np.uint16(img.dataobj)]
-    seg[np.bool_(nb.load(crown_mask).dataobj)] = 6 if nifti else 2
+    # seg[np.bool_(nb.load(crown_mask).dataobj)] = 6 if nifti else 2
     # Separate deep from shallow WM+CSF
-    seg[np.bool_(nb.load(acompcor_mask).dataobj)] = 4 if nifti else 1
+    # seg[np.bool_(nb.load(acompcor_mask).dataobj)] = 4 if nifti else 1
 
     outimg = img.__class__(seg.astype('uint8'), img.affine, img.header)
     outimg.set_data_dtype('uint8')