diff --git a/.circleci/circle_bold.txt b/.circleci/circle_bold.txt index 6c28affc..586f0af1 100644 --- a/.circleci/circle_bold.txt +++ b/.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 diff --git a/mriqc/interfaces/__init__.py b/mriqc/interfaces/__init__.py index 645ce7b2..de5cac73 100644 --- a/mriqc/interfaces/__init__.py +++ b/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", diff --git a/mriqc/interfaces/functional.py b/mriqc/interfaces/functional.py index dd72af0f..0f644e24 100644 --- a/mriqc/interfaces/functional.py +++ b/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 diff --git a/mriqc/workflows/functional/base.py b/mriqc/workflows/functional/base.py index 3db1cb22..abe4ed1b 100644 --- a/mriqc/workflows/functional/base.py +++ b/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