Merge pull request #39 from nipreps/add_extract_tacs_petprep

This PR implements support for extracting time activity curves from the PET data, given a specific segmentation. If a --ref-mask is specified, it will also extract a time activity curve for this region. If PVC correction is chosen, it will extract the time activity curve from the PVC corrected data.

Other relevant changes:

Expanded the outputs documentation with explanations of the TACs workflows. It now shows how both segmentation-based TACs and reference-region TACs are generated and named
Implemented a new ExtractRefTAC interface to compute TACs from reference masks and included it in a dedicated init_pet_ref_tacs_wf workflow

Author: mnoergaard

docs/outputs.rst

@@ -330,14 +330,30 @@ corresponding PET time series::
   662.82715 0.0 487.37302 1.01591 15.281561 0.0333937 0.3328022893  -0.2220965269 -0.0912891436 0.2326688125  0.279138129 -0.111878887  0.16901660629999998 0.0550480212  0.1798747037  -0.25383302620000003  0.1646403629  0.3953613889  0.0 0.010164  -0.0103568  0.0424513 0.0 -0.0  0.00019174
 
 **Time activity curves**.
-For each PET file, mean uptake values are extracted from the anatomical segmentation.
-The resulting table has ``FrameTimesStart`` and ``FrameTimesEnd`` columns followed
-by one column per region::
+The workflow :func:`petprep.workflows.pet.tacs.init_pet_tacs_wf` extracts mean uptake
+from an anatomical segmentation. The resulting table has ``FrameTimesStart`` and
+``FrameTimesEnd`` columns followed by one column per region::
 
   sub-<subject_label>/
     pet/
       sub-<subject_label>_[specifiers]_desc-preproc_timeseries.tsv
 
+If partial volume correction is applied, the filenames also include the
+``_pvc-<method>`` entity, indicating the algorithm used.
+
+If a reference mask is specified, the workflow
+:func:`petprep.workflows.pet.ref_tacs.init_pet_ref_tacs_wf` extracts a separate
+table containing the mean uptake within that region::
+
+  sub-<subject_label>/
+    pet/
+      sub-<subject_label>_[specifiers]_desc-<seg>_ref-<ref>_timeseries.tsv
+
+The ``ref`` entity captures the reference region identifier provided via the
+:ref:`CLI options <cli_refmask>` ``--ref-mask-name`` and ``--ref-mask-index``.
+When partial volume correction is performed, the ``_pvc-<method>`` entity is
+also included.
+
 Confounds
 ---------
 The :abbr:`BOLD (blood-oxygen level dependent)` signal measured with fMRI is a mixture of fluctuations

docs/usage.rst

@@ -203,6 +203,8 @@ method with a 5 mm PSF::
     $ petprep /data/bids_root /out participant \
         --pvc-tool petpvc --pvc-method GTM --pvc-psf 5
 
+.. _cli_refmask:
+
 Reference region masks
 ----------------------
 *PETPrep* can build masks for reference regions used in quantification.

petprep/cli/workflow.py

@@ -151,15 +151,15 @@ def build_workflow(config_file, retval):
     missing = check_deps(retval['workflow'])
     if missing:
         build_log.critical(
-            'Cannot run fMRIPrep. Missing dependencies:%s',
+            'Cannot run PETPrep. Missing dependencies:%s',
             '\n\t* '.join([''] + [f'{cmd} (Interface: {iface})' for iface, cmd in missing]),
         )
         retval['return_code'] = 127  # 127 == command not found.
         return retval
 
     config.to_filename(config_file)
     build_log.info(
-        'fMRIPrep workflow graph with %d nodes built successfully.',
+        'PETPrep workflow graph with %d nodes built successfully.',
         len(retval['workflow']._get_all_nodes()),
     )
     retval['return_code'] = 0
@@ -194,7 +194,7 @@ def build_boilerplate(config_file, workflow):
 
         bib_text = data.load.readable('boilerplate.bib').read_text()
         citation_files['bib'].write_text(
-            bib_text.replace('fMRIPrep <version>', f'fMRIPrep {config.environment.version}')
+            bib_text.replace('PETPrep <version>', f'PETPrep {config.environment.version}')
         )
 
         # Generate HTML file resolving citations
@@ -205,7 +205,7 @@ def build_boilerplate(config_file, workflow):
             str(citation_files['bib']),
             '--citeproc',
             '--metadata',
-            'pagetitle="fMRIPrep citation boilerplate"',
+            'pagetitle="PETPrep citation boilerplate"',
             str(citation_files['md']),
             '-o',
             str(citation_files['html']),

petprep/interfaces/__init__.py

@@ -3,7 +3,7 @@
 from niworkflows.interfaces.bids import DerivativesDataSink as _DDSink
 
 from .cifti import GeneratePetCifti
-from .tacs import ExtractTACs
+from .tacs import ExtractTACs, ExtractRefTAC
 
 
 class DerivativesDataSink(_DDSink):
@@ -14,4 +14,5 @@ class DerivativesDataSink(_DDSink):
     'DerivativesDataSink',
     'GeneratePetCifti',
     'ExtractTACs',
+    'ExtractRefTAC',
 )

petprep/interfaces/tacs.py

@@ -3,9 +3,14 @@
 import nibabel as nb
 import numpy as np
 import pandas as pd
-from nipype.interfaces.base import BaseInterfaceInputSpec, File, SimpleInterface, TraitedSpec
+from nipype.interfaces.base import (
+    BaseInterfaceInputSpec,
+    File,
+    SimpleInterface,
+    TraitedSpec,
+)
 from nipype.utils.filemanip import fname_presuffix
-from niworkflows.utils.timeseries import _nifti_timeseries
+from nipype.interfaces.base import traits
 
 
 class _ExtractTACsInputSpec(BaseInterfaceInputSpec):
@@ -81,4 +86,58 @@ def _run_interface(self, runtime):
         return runtime
 
 
-__all__ = ('ExtractTACs',)
+class _ExtractRefTACInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc="PET file in anatomical space")
+    mask_file = File(
+        exists=True, mandatory=True, desc="Reference mask in anatomical space"
+    )
+    ref_mask_name = traits.Str(mandatory=True, desc="Name of reference region")
+    metadata = File(exists=True, mandatory=True, desc="PET JSON metadata file")
+
+
+class _ExtractRefTACOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc="Reference region time activity curve")
+
+
+class ExtractRefTAC(SimpleInterface):
+    """Extract a time activity curve from a reference mask."""
+
+    input_spec = _ExtractRefTACInputSpec
+    output_spec = _ExtractRefTACOutputSpec
+
+    def _run_interface(self, runtime):
+        pet_img = nb.load(self.inputs.in_file)
+        pet_data = pet_img.get_fdata()
+        if pet_img.ndim == 3:
+            pet_data = pet_data[..., np.newaxis]
+
+        mask = nb.load(self.inputs.mask_file).get_fdata() > 0
+
+        with open(self.inputs.metadata) as f:
+            metadata = json.load(f)
+
+        frame_times = metadata.get("FrameTimesStart", [])
+        frame_durations = metadata.get("FrameDuration", [])
+
+        if len(frame_times) != len(frame_durations):
+            raise ValueError("FrameTimesStart and FrameDuration must have equal length")
+
+        timeseries = pet_data[mask, :].mean(axis=0)
+        frame_times_end = np.add(frame_times, frame_durations).tolist()
+        df = pd.DataFrame({self.inputs.ref_mask_name: timeseries})
+        df.insert(0, "FrameTimesEnd", frame_times_end)
+        df.insert(0, "FrameTimesStart", list(frame_times))
+
+        out_file = fname_presuffix(
+            self.inputs.in_file,
+            suffix="_timeseries.tsv",
+            newpath=runtime.cwd,
+            use_ext=False,
+        )
+        df.to_csv(out_file, sep="\t", index=False, na_rep="n/a")
+
+        self._results["out_file"] = out_file
+        return runtime
+
+
+__all__ = ("ExtractTACs", "ExtractRefTAC")

petprep/interfaces/tests/test_tacs.py

@@ -1,4 +1,7 @@
+import gzip
 import json
+import pickle
+from pathlib import Path
 
 import nibabel as nb
 import numpy as np
@@ -7,14 +10,19 @@
 from nipype.pipeline import engine as pe
 from nipype.pipeline.engine.nodes import NodeExecutionError
 
-from petprep.interfaces.tacs import ExtractTACs
+from petprep.interfaces.tacs import ExtractRefTAC, ExtractTACs
+from petprep.workflows.pet.ref_tacs import init_pet_ref_tacs_wf
+from petprep.workflows.pet.tacs import init_pet_tacs_wf
 
 
 def test_ExtractTACs(tmp_path):
-    pet_data = np.stack([
-        np.ones((2, 2, 2)),
-        np.ones((2, 2, 2)) * 2,
-    ], axis=-1)
+    pet_data = np.stack(
+        [
+            np.ones((2, 2, 2)),
+            np.ones((2, 2, 2)) * 2,
+        ],
+        axis=-1,
+    )
     pet_file = tmp_path / 'pet.nii.gz'
     nb.Nifti1Image(pet_data, np.eye(4)).to_filename(pet_file)
 
@@ -47,10 +55,13 @@ def test_ExtractTACs(tmp_path):
 
 
 def test_ExtractTACs_mismatched_meta(tmp_path):
-    pet_data = np.stack([
-        np.ones((2, 2, 2)),
-        np.ones((2, 2, 2)) * 2,
-    ], axis=-1)
+    pet_data = np.stack(
+        [
+            np.ones((2, 2, 2)),
+            np.ones((2, 2, 2)) * 2,
+        ],
+        axis=-1,
+    )
     pet_file = tmp_path / 'pet.nii.gz'
     nb.Nifti1Image(pet_data, np.eye(4)).to_filename(pet_file)
 
@@ -76,4 +87,149 @@ def test_ExtractTACs_mismatched_meta(tmp_path):
     )
 
     with pytest.raises(NodeExecutionError):
-        node.run()
+        node.run()
+
+
+def test_tacs_workflow(tmp_path):
+    """Workflow passes resampled PET to ExtractTACs."""
+    pet_data = np.stack(
+        [
+            np.ones((2, 2, 2)),
+            np.ones((2, 2, 2)) * 2,
+        ],
+        axis=-1,
+    )
+    pet_file = tmp_path / 'pet.nii.gz'
+    nb.Nifti1Image(pet_data, np.eye(4)).to_filename(pet_file)
+
+    seg_data = np.tile([[1, 2], [1, 2]], (2, 1, 1)).astype('int16')
+    seg_file = tmp_path / 'seg.nii.gz'
+    nb.Nifti1Image(seg_data, np.eye(4)).to_filename(seg_file)
+
+    dseg_tsv = tmp_path / 'seg.tsv'
+    pd.DataFrame({'index': [1, 2], 'name': ['A', 'B']}).to_csv(dseg_tsv, sep='\t', index=False)
+
+    meta_json = tmp_path / 'pet.json'
+    meta_json.write_text(json.dumps({'FrameTimesStart': [0, 1], 'FrameDuration': [1, 1]}))
+
+    wf = init_pet_tacs_wf()
+    wf.base_dir = str(tmp_path)
+    wf.config['execution']['remove_unnecessary_outputs'] = False
+    wf.inputs.inputnode.pet_anat = str(pet_file)
+    wf.inputs.inputnode.segmentation = str(seg_file)
+    wf.inputs.inputnode.dseg_tsv = str(dseg_tsv)
+    wf.inputs.inputnode.metadata = str(meta_json)
+
+    wf.run()
+
+    resampled_pet = tmp_path / 'pet_tacs_wf' / 'resample_pet' / 'pet_resampled.nii.gz'
+    tac_inputs_file = tmp_path / 'pet_tacs_wf' / 'tac' / '_inputs.pklz'
+    with gzip.open(tac_inputs_file, 'rb') as f:
+        inputs = pickle.load(f)
+
+    assert inputs['in_file'] == str(resampled_pet)
+    assert Path(tmp_path / 'pet_tacs_wf' / 'tac' / 'pet_resampled_timeseries.tsv').exists()
+
+
+def test_ExtractRefTAC(tmp_path):
+    pet_data = np.stack(
+        [
+            np.ones((2, 2, 2)),
+            np.ones((2, 2, 2)) * 2,
+        ],
+        axis=-1,
+    )
+    pet_file = tmp_path / "pet.nii.gz"
+    nb.Nifti1Image(pet_data, np.eye(4)).to_filename(pet_file)
+
+    mask_data = np.zeros((2, 2, 2), dtype="int16")
+    mask_data[0] = 1
+    mask_file = tmp_path / "mask.nii.gz"
+    nb.Nifti1Image(mask_data, np.eye(4)).to_filename(mask_file)
+
+    meta_json = tmp_path / "pet.json"
+    meta_json.write_text(
+        json.dumps({"FrameTimesStart": [0, 1], "FrameDuration": [1, 1]})
+    )
+
+    node = pe.Node(
+        ExtractRefTAC(
+            in_file=str(pet_file),
+            mask_file=str(mask_file),
+            ref_mask_name="ref",
+            metadata=str(meta_json),
+        ),
+        name="tac",
+        base_dir=tmp_path,
+    )
+    res = node.run()
+
+    out = pd.read_csv(res.outputs.out_file, sep="\t")
+    assert list(out.columns) == ["FrameTimesStart", "FrameTimesEnd", "ref"]
+    assert np.allclose(out["ref"], [1, 2])
+
+
+def test_ExtractRefTAC_mismatched_meta(tmp_path):
+    pet_data = np.stack(
+        [
+            np.ones((2, 2, 2)),
+            np.ones((2, 2, 2)) * 2,
+        ],
+        axis=-1,
+    )
+    pet_file = tmp_path / "pet.nii.gz"
+    nb.Nifti1Image(pet_data, np.eye(4)).to_filename(pet_file)
+
+    mask_data = np.zeros((2, 2, 2), dtype="int16")
+    mask_data[0] = 1
+    mask_file = tmp_path / "mask.nii.gz"
+    nb.Nifti1Image(mask_data, np.eye(4)).to_filename(mask_file)
+
+    meta_json = tmp_path / "pet.json"
+    meta_json.write_text(json.dumps({"FrameTimesStart": [0], "FrameDuration": [1, 1]}))
+
+    node = pe.Node(
+        ExtractRefTAC(
+            in_file=str(pet_file),
+            mask_file=str(mask_file),
+            ref_mask_name="ref",
+            metadata=str(meta_json),
+        ),
+        name="tac2",
+        base_dir=tmp_path,
+    )
+
+    with pytest.raises(NodeExecutionError):
+        node.run()
+
+
+def test_ref_tacs_workflow_mismatched_meta(tmp_path):
+    """Workflow should fail with inconsistent metadata."""
+    pet_data = np.stack(
+        [
+            np.ones((2, 2, 2)),
+            np.ones((2, 2, 2)) * 2,
+        ],
+        axis=-1,
+    )
+    pet_file = tmp_path / 'pet.nii.gz'
+    nb.Nifti1Image(pet_data, np.eye(4)).to_filename(pet_file)
+
+    mask_data = np.zeros((2, 2, 2), dtype='int16')
+    mask_data[0] = 1
+    mask_file = tmp_path / 'mask.nii.gz'
+    nb.Nifti1Image(mask_data, np.eye(4)).to_filename(mask_file)
+
+    meta_json = tmp_path / 'pet.json'
+    meta_json.write_text(json.dumps({'FrameTimesStart': [0], 'FrameDuration': [1, 1]}))
+
+    wf = init_pet_ref_tacs_wf()
+    wf.base_dir = str(tmp_path)
+    wf.config['execution']['remove_unnecessary_outputs'] = False
+    wf.inputs.inputnode.pet_anat = str(pet_file)
+    wf.inputs.inputnode.mask_file = str(mask_file)
+    wf.inputs.inputnode.metadata = str(meta_json)
+    wf.inputs.inputnode.ref_mask_name = 'ref'
+
+    with pytest.raises(NodeExecutionError):
+        wf.run()

petprep/workflows/pet/__init__.py

@@ -20,11 +20,13 @@
 from .registration import init_pet_reg_wf
 from .resampling import init_pet_surf_wf
 from .tacs import init_pet_tacs_wf
+from .ref_tacs import init_pet_ref_tacs_wf
 
 __all__ = [
     'init_pet_confs_wf',
     'init_pet_hmc_wf',
     'init_pet_reg_wf',
     'init_pet_surf_wf',
     'init_pet_tacs_wf',
+    'init_pet_ref_tacs_wf',
 ]