Skip to content

Commit

Permalink
Merge pull request #1166 from psadil/enh/store-qc-timeseries
Browse files Browse the repository at this point in the history 
ENH: store confound timeseries data
  • Loading branch information
@effigies
effigies committed Dec 6, 2023
2 parents db3d3cb + 92cb229 commit 41fdc70
Show file tree
Hide file tree
Showing 4 changed files with 241 additions and 5 deletions.
18 changes: 18 additions & 0 deletions .circleci/circle_bold.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion mriqc/interfaces/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand All @@ -47,6 +47,7 @@ class DerivativesDataSink(_DDSink):
"DerivativesDataSink",
"EnsureSize",
"FunctionalQC",
"GatherTimeseries",
"Harmonize",
"IQMFileSink",
"RotationMask",
Expand Down
181 changes: 181 additions & 0 deletions mriqc/interfaces/functional.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,8 @@
traits,
Undefined,
)
from nipype.utils.misc import normalize_mc_params
import pandas as pd


class FunctionalQCInputSpec(BaseInterfaceInputSpec):
Expand Down Expand Up @@ -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
Expand Down
44 changes: 40 additions & 4 deletions mriqc/workflows/functional/base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -302,6 +308,7 @@ def compute_iqms(name="ComputeIQMs"):
"fd_thres",
"in_tsnr",
"metadata",
"mpars",
"exclude_index",
"subject",
"session",
Expand Down Expand Up @@ -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"),
Expand All @@ -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

Expand All @@ -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")]),
Expand All @@ -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

Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 41fdc70

Please sign in to comment.