ENH: store confound timeseries data

.circleci/circle_bold.txt

@@ -30,8 +30,14 @@ sub-ds205s03/figures/sub-ds205s03_task-view_run-02_desc-stdev_bold.svg
 sub-ds205s03/figures/sub-ds205s03_task-view_run-02_desc-zoomed_bold.svg
 sub-ds205s03/func
 sub-ds205s03/func/sub-ds205s03_task-functionallocalizer_run-01_bold.json
+sub-ds205s03/func/sub-ds205s03_task-functionallocalizer_run-01_timeseries.json
+sub-ds205s03/func/sub-ds205s03_task-functionallocalizer_run-01_timeseries.tsv
 sub-ds205s03/func/sub-ds205s03_task-view_run-01_bold.json
+sub-ds205s03/func/sub-ds205s03_task-view_run-01_timeseries.json
+sub-ds205s03/func/sub-ds205s03_task-view_run-01_timeseries.tsv
 sub-ds205s03/func/sub-ds205s03_task-view_run-02_bold.json
+sub-ds205s03/func/sub-ds205s03_task-view_run-02_timeseries.json
+sub-ds205s03/func/sub-ds205s03_task-view_run-02_timeseries.tsv
 sub-ds205s03_task-functionallocalizer_run-01_bold.html
 sub-ds205s03_task-view_run-01_bold.html
 sub-ds205s03_task-view_run-02_bold.html
@@ -60,8 +66,14 @@ sub-ds205s07/figures/sub-ds205s07_task-view_run-02_desc-stdev_bold.svg
 sub-ds205s07/figures/sub-ds205s07_task-view_run-02_desc-zoomed_bold.svg
 sub-ds205s07/func
 sub-ds205s07/func/sub-ds205s07_task-functionallocalizer_run-01_bold.json
+sub-ds205s07/func/sub-ds205s07_task-functionallocalizer_run-01_timeseries.json
+sub-ds205s07/func/sub-ds205s07_task-functionallocalizer_run-01_timeseries.tsv
 sub-ds205s07/func/sub-ds205s07_task-view_run-01_bold.json
+sub-ds205s07/func/sub-ds205s07_task-view_run-01_timeseries.json
+sub-ds205s07/func/sub-ds205s07_task-view_run-01_timeseries.tsv
 sub-ds205s07/func/sub-ds205s07_task-view_run-02_bold.json
+sub-ds205s07/func/sub-ds205s07_task-view_run-02_timeseries.json
+sub-ds205s07/func/sub-ds205s07_task-view_run-02_timeseries.tsv
 sub-ds205s07_task-functionallocalizer_run-01_bold.html
 sub-ds205s07_task-view_run-01_bold.html
 sub-ds205s07_task-view_run-02_bold.html
@@ -90,8 +102,14 @@ sub-ds205s09/figures/sub-ds205s09_task-view_acq-RL_run-01_desc-stdev_bold.svg
 sub-ds205s09/figures/sub-ds205s09_task-view_acq-RL_run-01_desc-zoomed_bold.svg
 sub-ds205s09/func
 sub-ds205s09/func/sub-ds205s09_task-view_acq-LR_run-01_bold.json
+sub-ds205s09/func/sub-ds205s09_task-view_acq-LR_run-01_timeseries.json
+sub-ds205s09/func/sub-ds205s09_task-view_acq-LR_run-01_timeseries.tsv
 sub-ds205s09/func/sub-ds205s09_task-view_acq-LR_run-02_bold.json
+sub-ds205s09/func/sub-ds205s09_task-view_acq-LR_run-02_timeseries.json
+sub-ds205s09/func/sub-ds205s09_task-view_acq-LR_run-02_timeseries.tsv
 sub-ds205s09/func/sub-ds205s09_task-view_acq-RL_run-01_bold.json
+sub-ds205s09/func/sub-ds205s09_task-view_acq-RL_run-01_timeseries.json
+sub-ds205s09/func/sub-ds205s09_task-view_acq-RL_run-01_timeseries.tsv
 sub-ds205s09_task-view_acq-LR_run-01_bold.html
 sub-ds205s09_task-view_acq-LR_run-02_bold.html
 sub-ds205s09_task-view_acq-RL_run-01_bold.html

mriqc/interfaces/__init__.py

@@ -32,7 +32,7 @@
 )
 from mriqc.interfaces.bids import IQMFileSink
 from mriqc.interfaces.common import ConformImage, EnsureSize
-from mriqc.interfaces.functional import FunctionalQC, Spikes
+from mriqc.interfaces.functional import FunctionalQC, GatherTimeseries, Spikes
 from mriqc.interfaces.webapi import UploadIQMs
 
 
@@ -47,6 +47,7 @@ class DerivativesDataSink(_DDSink):
     "DerivativesDataSink",
     "EnsureSize",
     "FunctionalQC",
+    "GatherTimeseries",
     "Harmonize",
     "IQMFileSink",
     "RotationMask",

mriqc/interfaces/functional.py

@@ -38,6 +38,8 @@
     traits,
     Undefined,
 )
+from nipype.utils.misc import normalize_mc_params
+import pandas as pd
 
 
 class FunctionalQCInputSpec(BaseInterfaceInputSpec):
@@ -309,6 +311,185 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class GatherTimeseriesInputSpec(TraitedSpec):
+    dvars = File(exists=True, mandatory=True, desc='file containing DVARS')
+    fd = File(exists=True, mandatory=True, desc='input framewise displacement')
+    mpars = File(exists=True, mandatory=True, desc='input motion parameters')
+    mpars_source = traits.Enum(
+        "FSL",
+        "AFNI",
+        "SPM",
+        "FSFAST",
+        "NIPY",
+        desc="Source of movement parameters",
+        mandatory=True,
+    )
+    outliers = File(
+        exists=True,
+        mandatory=True,
+        desc="input file containing timeseries of AFNI's outlier count")
+    quality = File(
+        exists=True,
+        mandatory=True,
+        desc="input file containing AFNI's Quality Index")
+
+
+class GatherTimeseriesOutputSpec(TraitedSpec):
+    timeseries_file = File(desc='output confounds file')
+    timeseries_metadata = traits.Dict(desc='Metadata dictionary describing columns')
+
+
+class GatherTimeseries(SimpleInterface):
+    """
+    Gather quality metrics that are timeseries into one TSV file
+
+    """
+    input_spec = GatherTimeseriesInputSpec
+    output_spec = GatherTimeseriesOutputSpec
+
+    def _run_interface(self, runtime):
+
+        # motion parameters
+        mpars = np.apply_along_axis(
+            func1d=normalize_mc_params,
+            axis=1,
+            arr=np.loadtxt(self.inputs.mpars),  # mpars is N_t x 6
+            source=self.inputs.mpars_source,
+        )
+        timeseries = pd.DataFrame(
+            mpars,
+            columns=[
+                "trans_x",
+                "trans_y",
+                "trans_z",
+                "rot_x",
+                "rot_y",
+                "rot_z"
+            ])
+
+        # DVARS
+        dvars = pd.read_csv(
+            self.inputs.dvars,
+            delim_whitespace=True,
+            skiprows=1,  # column names have spaces
+            header=None,
+            names=["dvars_std", "dvars_nstd", "dvars_vstd"])
+        dvars.index = pd.RangeIndex(1, timeseries.index.max() + 1)
+
+        # FD
+        fd = pd.read_csv(
+            self.inputs.fd,
+            delim_whitespace=True,
+            header=0,
+            names=["framewise_displacement"])
+        fd.index = pd.RangeIndex(1, timeseries.index.max() + 1)
+
+        # AQI
+        aqi = pd.read_csv(
+            self.inputs.quality,
+            delim_whitespace=True,
+            header=None,
+            names=["aqi"])
+
+        # Outliers
+        aor = pd.read_csv(
+            self.inputs.outliers,
+            delim_whitespace=True,
+            header=None,
+            names=["aor"])
+
+        timeseries = pd.concat((timeseries, dvars, fd, aqi, aor), axis=1)
+
+        timeseries_file = op.join(runtime.cwd, "timeseries.tsv")
+
+        timeseries.to_csv(timeseries_file, sep='\t', index=False, na_rep='n/a')
+
+        self._results['timeseries_file'] = timeseries_file
+        self._results['timeseries_metadata'] = _build_timeseries_metadata()
+        return runtime
+
+
+def _build_timeseries_metadata():
+    return {
+        "trans_x": {
+            "LongName": "Translation Along X Axis",
+            "Description": "Estimated Motion Parameter",
+            "Units": "mm"
+        },
+        "trans_y": {
+            "LongName": "Translation Along Y Axis",
+            "Description": "Estimated Motion Parameter",
+            "Units": "mm"
+        },
+        "trans_z": {
+            "LongName": "Translation Along Z Axis",
+            "Description": "Estimated Motion Parameter",
+            "Units": "mm",
+        },
+        "rot_x": {
+            "LongName": "Rotation Around X Axis",
+            "Description": "Estimated Motion Parameter",
+            "Units": "rad"
+        },
+        "rot_y": {
+            "LongName": "Rotation Around X Axis",
+            "Description": "Estimated Motion Parameter",
+            "Units": "rad"
+        },
+        "rot_z": {
+            "LongName": "Rotation Around X Axis",
+            "Description": "Estimated Motion Parameter",
+            "Units": "rad"
+        },
+        "dvars_std": {
+            "LongName": "Derivative of RMS Variance over Voxels, Standardized",
+            "Description": (
+                "Indexes the rate of change of BOLD signal across"
+                "the entire brain at each frame of data, normalized with the"
+                "standard deviation of the temporal difference time series"
+            )
+        },
+        "dvars_nstd": {
+            "LongName": (
+                "Derivative of RMS Variance over Voxels,"
+                "Non-Standardized"
+            ),
+            "Description": (
+                "Indexes the rate of change of BOLD signal across"
+                "the entire brain at each frame of data, not normalized."
+            )
+        },
+        "dvars_vstd": {
+            "LongName": "Derivative of RMS Variance over Voxels, Standardized",
+            "Description": (
+                "Indexes the rate of change of BOLD signal across"
+                "the entire brain at each frame of data, normalized across"
+                "time by that voxel standard deviation across time,"
+                "before computing the RMS of the temporal difference"
+            )
+        },
+        "framewise_displacement": {
+            "LongName": "Framewise Displacement",
+            "Description": (
+                "A quantification of the estimated bulk-head"
+                "motion calculated using formula proposed by Power (2012)"
+            ),
+            "Units": "mm"
+        },
+        "aqi": {
+            "LongName": "AFNI's Quality Index",
+            "Description": "Mean quality index as computed by AFNI's 3dTqual"
+        },
+        "aor": {
+            "LongName": "AFNI's Fraction of Outliers per Volume",
+            "Description": (
+                "Mean fraction of outliers per fMRI volume as"
+                "given by AFNI's 3dToutcount"
+            )
+        }
+    }
+
+
 def find_peaks(data):
     t_z = [data[:, :, i, :].mean(axis=0).mean(axis=0) for i in range(data.shape[2])]
     return t_z

mriqc/workflows/functional/base.py

@@ -177,7 +177,8 @@ def fmri_qc_workflow(name="funcMRIQC"):
         (sanitize, iqmswf, [("out_file", "inputnode.in_ras")]),
         (mean, iqmswf, [("out_file", "inputnode.epi_mean")]),
         (hmcwf, iqmswf, [("outputnode.out_file", "inputnode.hmc_epi"),
-                         ("outputnode.out_fd", "inputnode.hmc_fd")]),
+                         ("outputnode.out_fd", "inputnode.hmc_fd"),
+                         ("outputnode.mpars", "inputnode.mpars")]),
         (tsnr, iqmswf, [("tsnr_file", "inputnode.in_tsnr")]),
         (non_steady_state_detector, iqmswf, [("n_volumes_to_discard", "inputnode.exclude_index")]),
         # Feed reportlet generation
@@ -282,7 +283,12 @@ def compute_iqms(name="ComputeIQMs"):
     from nipype.algorithms.confounds import ComputeDVARS
     from nipype.interfaces.afni import OutlierCount, QualityIndex
 
-    from mriqc.interfaces import FunctionalQC, IQMFileSink
+    from mriqc.interfaces import (
+        DerivativesDataSink,
+        FunctionalQC,
+        IQMFileSink,
+        GatherTimeseries
+    )
     from mriqc.interfaces.reports import AddProvenance
     from mriqc.interfaces.transitional import GCOR
     from mriqc.workflows.utils import _tofloat, get_fwhmx
@@ -302,6 +308,7 @@ def compute_iqms(name="ComputeIQMs"):
                 "fd_thres",
                 "in_tsnr",
                 "metadata",
+                "mpars",
                 "exclude_index",
                 "subject",
                 "session",
@@ -365,6 +372,13 @@ def compute_iqms(name="ComputeIQMs"):
         iterfield=["in_epi", "in_hmc", "in_tsnr", "in_dvars", "in_fwhm"],
     )
 
+    timeseries = pe.MapNode(
+        GatherTimeseries(mpars_source="AFNI"),
+        name="timeseries",
+        mem_gb=mem_gb * 3,
+        iterfield=["dvars", "outliers", "quality", "fd"]
+    )
+
     # fmt: off
     workflow.connect([
         (inputnode, dvnode, [("hmc_epi", "in_file"),
@@ -385,7 +399,11 @@ def compute_iqms(name="ComputeIQMs"):
         (dvnode, measures, [("out_all", "in_dvars")]),
         (fwhm, measures, [(("fwhm", _tofloat), "in_fwhm")]),
         (dvnode, outputnode, [("out_all", "out_dvars")]),
-        (outliers, outputnode, [("out_file", "outliers")])
+        (outliers, outputnode, [("out_file", "outliers")]),
+        (outliers, timeseries, [("out_file", "outliers")]),
+        (quality, timeseries, [("out_file", "quality")]),
+        (dvnode, timeseries, [("out_all", "dvars")]),
+        (inputnode, timeseries, [("hmc_fd", "fd"), ("mpars", "mpars")]),
     ])
     # fmt: on
 
@@ -408,6 +426,17 @@ def compute_iqms(name="ComputeIQMs"):
         iterfield=["in_file", "root", "metadata", "provenance"],
     )
 
+    # Save timeseries TSV file
+    ds_timeseries = pe.MapNode(
+        DerivativesDataSink(
+            base_directory=str(config.execution.output_dir),
+            suffix="timeseries"
+        ),
+        name="ds_timeseries",
+        run_without_submitting=True,
+        iterfield=["in_file", "source_file", "meta_dict"],
+    )
+
     # fmt: off
     workflow.connect([
         (inputnode, addprov, [("in_file", "in_file")]),
@@ -426,6 +455,9 @@ def compute_iqms(name="ComputeIQMs"):
         (quality, datasink, [(("out_file", _parse_tqual), "aqi")]),
         (measures, datasink, [("out_qc", "root")]),
         (datasink, outputnode, [("out_file", "out_file")]),
+        (inputnode, ds_timeseries, [("in_file", "source_file")]),
+        (timeseries, ds_timeseries, [("timeseries_file", "in_file"),
+                                     ("timeseries_metadata", "meta_dict")]),
     ])
     # fmt: on
 
@@ -509,7 +541,10 @@ def hmc(name="fMRI_HMC", omp_nthreads=None):
         name="inputnode",
     )
 
-    outputnode = pe.Node(niu.IdentityInterface(fields=["out_file", "out_fd"]), name="outputnode")
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=["out_file", "out_fd", "mpars"]),
+        name="outputnode",
+    )
 
     # calculate hmc parameters
     estimate_hm = pe.Node(
@@ -539,6 +574,7 @@ def hmc(name="fMRI_HMC", omp_nthreads=None):
         (estimate_hm, apply_hmc, [("oned_matrix_save", "in_xfm")]),
         (apply_hmc, outputnode, [("out", "out_file")]),
         (estimate_hm, fdnode, [("oned_file", "in_file")]),
+        (estimate_hm, outputnode, [("oned_file", "mpars")]),
         (fdnode, outputnode, [("out_file", "out_fd")]),
     ])
     # fmt: on