diff --git a/xcp_d/interfaces/plotting.py b/xcp_d/interfaces/plotting.py index 7abf8b8ce..780e5407e 100644 --- a/xcp_d/interfaces/plotting.py +++ b/xcp_d/interfaces/plotting.py @@ -1,7 +1,6 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Plotting interfaces.""" -import json import os import matplotlib.pyplot as plt @@ -33,7 +32,7 @@ from xcp_d.utils.filemanip import fname_presuffix from xcp_d.utils.modified_data import compute_fd from xcp_d.utils.plotting import FMRIPlot, plot_fmri_es, surf_data_from_cifti -from xcp_d.utils.qcmetrics import compute_dvars, compute_registration_qc +from xcp_d.utils.qcmetrics import compute_dvars from xcp_d.utils.write_save import read_ndata LOGGER = logging.getLogger("nipype.interface") @@ -205,22 +204,11 @@ def _run_interface(self, runtime): class _QCPlotsInputSpec(BaseInterfaceInputSpec): - name_source = File( - exists=False, - mandatory=True, - desc=( - "Preprocessed BOLD file. Used to find files. " - "In the case of the concatenation workflow, " - "this may be a nonexistent file " - "(i.e., the preprocessed BOLD file, with the run entity removed)." - ), - ) bold_file = File( exists=True, mandatory=True, desc="Preprocessed BOLD file, after dummy scan removal. Used in carpet plot.", ) - dummy_scans = traits.Int(mandatory=True, desc="Dummy time to drop") temporal_mask = traits.Either( File(exists=True), Undefined, @@ -247,19 +235,9 @@ class _QCPlotsInputSpec(BaseInterfaceInputSpec): # Inputs used only for nifti data seg_file = File(exists=True, mandatory=False, desc="Seg file for nifti") - anat_brainmask = File( - exists=True, - mandatory=False, - desc="Anatomically-derived brain mask in anatomical space.", - ) - template_mask = File(exists=True, mandatory=False, desc="Template mask") - bold2T1w_mask = File(exists=True, mandatory=False, desc="BOLD mask in MNI space") - bold2temp_mask = File(exists=True, mandatory=False, desc="BOLD mask in anatomical space") class _QCPlotsOutputSpec(TraitedSpec): - qc_file = File(exists=True, desc="QC TSV file.") - qc_metadata = File(exists=True, desc="Sidecar JSON for QC TSV file.") raw_qcplot = File(exists=True, desc="qc plot before regression") clean_qcplot = File(exists=True, desc="qc plot after regression") @@ -276,12 +254,10 @@ class QCPlots(SimpleInterface): .. doctest:: qcplots = QCPlots() qcplots.inputs.cleaned_file = datafile - qcplots.inputs.name_source = rawbold qcplots.inputs.bold_file = rawbold qcplots.inputs.TR = TR qcplots.inputs.temporal_mask = temporalmask qcplots.inputs.mask_file = mask - qcplots.inputs.dummy_scans = dummy_scans qcplots.run() .. testcleanup:: >>> tmpdir.cleanup() @@ -303,22 +279,14 @@ def _run_interface(self, runtime): head_radius=self.inputs.head_radius, ) - # Get rmsd - rmsd = confounds_df["rmsd"].to_numpy() - # Determine number of dummy volumes and load temporal mask - dummy_scans = self.inputs.dummy_scans if isdefined(self.inputs.temporal_mask): censoring_df = pd.read_table(self.inputs.temporal_mask) tmask_arr = censoring_df["framewise_displacement"].values else: tmask_arr = np.zeros(preproc_fd_timeseries.size, dtype=int) - num_censored_volumes = int(tmask_arr.sum()) - num_retained_volumes = int((tmask_arr == 0).sum()) - # Apply temporal mask to interpolated/full data - rmsd_censored = rmsd[tmask_arr == 0] postproc_fd_timeseries = preproc_fd_timeseries[tmask_arr == 0] # get QC plot names @@ -391,173 +359,6 @@ def _run_interface(self, runtime): ) plt.close(postproc_fig) - # Get the different components in the bold file name - # eg: ['sub-colornest001', 'ses-1'], etc. - _, bold_file_name = os.path.split(self.inputs.name_source) - bold_file_name_components = bold_file_name.split("_") - - # Fill out dictionary with entities from filename - qc_values_dict = {} - for entity in bold_file_name_components[:-1]: - qc_values_dict[entity.split("-")[0]] = entity.split("-")[1] - - # Calculate QC measures - mean_fd = np.mean(preproc_fd_timeseries) - mean_fd_post_censoring = np.mean(postproc_fd_timeseries) - mean_relative_rms = np.nanmean(rmsd_censored) # first value can be NaN if no dummy scans - mean_dvars_before_processing = np.mean(dvars_before_processing) - mean_dvars_after_processing = np.mean(dvars_after_processing) - fd_dvars_correlation_initial = np.corrcoef(preproc_fd_timeseries, dvars_before_processing)[ - 0, 1 - ] - fd_dvars_correlation_final = np.corrcoef(postproc_fd_timeseries, dvars_after_processing)[ - 0, 1 - ] - rmsd_max_value = np.nanmax(rmsd_censored) - - # A summary of all the values - qc_values_dict.update( - { - "mean_fd": [mean_fd], - "mean_fd_post_censoring": [mean_fd_post_censoring], - "mean_relative_rms": [mean_relative_rms], - "max_relative_rms": [rmsd_max_value], - "mean_dvars_initial": [mean_dvars_before_processing], - "mean_dvars_final": [mean_dvars_after_processing], - "num_dummy_volumes": [dummy_scans], - "num_censored_volumes": [num_censored_volumes], - "num_retained_volumes": [num_retained_volumes], - "fd_dvars_correlation_initial": [fd_dvars_correlation_initial], - "fd_dvars_correlation_final": [fd_dvars_correlation_final], - } - ) - - qc_metadata = { - "mean_fd": { - "LongName": "Mean Framewise Displacement", - "Description": ( - "Average framewise displacement without any motion parameter filtering. " - "This value includes high-motion outliers, but not dummy volumes. " - "FD is calculated according to the Power definition." - ), - "Units": "mm / volume", - "Term URL": "https://doi.org/10.1016/j.neuroimage.2011.10.018", - }, - "mean_fd_post_censoring": { - "LongName": "Mean Framewise Displacement After Censoring", - "Description": ( - "Average framewise displacement without any motion parameter filtering. " - "This value does not include high-motion outliers or dummy volumes. " - "FD is calculated according to the Power definition." - ), - "Units": "mm / volume", - "Term URL": "https://doi.org/10.1016/j.neuroimage.2011.10.018", - }, - "mean_relative_rms": { - "LongName": "Mean Relative Root Mean Squared", - "Description": ( - "Average relative root mean squared calculated from motion parameters, " - "after removal of dummy volumes and high-motion outliers. " - "Relative in this case means 'relative to the previous scan'." - ), - "Units": "arbitrary", - }, - "max_relative_rms": { - "LongName": "Maximum Relative Root Mean Squared", - "Description": ( - "Maximum relative root mean squared calculated from motion parameters, " - "after removal of dummy volumes and high-motion outliers. " - "Relative in this case means 'relative to the previous scan'." - ), - "Units": "arbitrary", - }, - "mean_dvars_initial": { - "LongName": "Mean DVARS Before Postprocessing", - "Description": ( - "Average DVARS (temporal derivative of root mean squared variance over " - "voxels) calculated from the preprocessed BOLD file, after dummy scan removal." - ), - "TermURL": "https://doi.org/10.1016/j.neuroimage.2011.02.073", - }, - "mean_dvars_final": { - "LongName": "Mean DVARS After Postprocessing", - "Description": ( - "Average DVARS (temporal derivative of root mean squared variance over " - "voxels) calculated from the denoised BOLD file." - ), - "TermURL": "https://doi.org/10.1016/j.neuroimage.2011.02.073", - }, - "num_dummy_volumes": { - "LongName": "Number of Dummy Volumes", - "Description": ( - "The number of non-steady state volumes removed from the time series by XCP-D." - ), - }, - "num_censored_volumes": { - "LongName": "Number of Censored Volumes", - "Description": ( - "The number of high-motion outlier volumes censored by XCP-D. " - "This does not include dummy volumes." - ), - }, - "num_retained_volumes": { - "LongName": "Number of Retained Volumes", - "Description": ( - "The number of volumes retained in the denoised dataset. " - "This does not include dummy volumes or high-motion outliers." - ), - }, - "fd_dvars_correlation_initial": { - "LongName": "FD-DVARS Correlation Before Postprocessing", - "Description": ( - "The Pearson correlation coefficient between framewise displacement and DVARS " - "(temporal derivative of root mean squared variance over voxels), " - "after removal of dummy volumes, but before removal of high-motion outliers." - ), - }, - "fd_dvars_correlation_final": { - "LongName": "FD-DVARS Correlation After Postprocessing", - "Description": ( - "The Pearson correlation coefficient between framewise displacement and DVARS " - "(temporal derivative of root mean squared variance over voxels), " - "after postprocessing. " - "The FD time series is unfiltered, but censored. " - "The DVARS time series is calculated from the denoised BOLD data." - ), - }, - } - - if self.inputs.bold2T1w_mask: # If a bold mask in T1w is provided - # Compute quality of registration - registration_qc, registration_metadata = compute_registration_qc( - bold2t1w_mask=self.inputs.bold2T1w_mask, - anat_brainmask=self.inputs.anat_brainmask, - bold2template_mask=self.inputs.bold2temp_mask, - template_mask=self.inputs.template_mask, - ) - qc_values_dict.update(registration_qc) # Add values to dictionary - qc_metadata.update(registration_metadata) - - # Convert dictionary to df and write out the qc file - df = pd.DataFrame(qc_values_dict) - self._results["qc_file"] = fname_presuffix( - self.inputs.cleaned_file, - suffix="qc_bold.tsv", - newpath=runtime.cwd, - use_ext=False, - ) - df.to_csv(self._results["qc_file"], index=False, header=True, sep="\t") - - # Write out the metadata file - self._results["qc_metadata"] = fname_presuffix( - self.inputs.cleaned_file, - suffix="qc_bold.json", - newpath=runtime.cwd, - use_ext=False, - ) - with open(self._results["qc_metadata"], "w") as fo: - json.dump(qc_metadata, fo, indent=4, sort_keys=True) - return runtime diff --git a/xcp_d/interfaces/utils.py b/xcp_d/interfaces/utils.py index 86c7821f8..796fdbe4c 100644 --- a/xcp_d/interfaces/utils.py +++ b/xcp_d/interfaces/utils.py @@ -1,5 +1,11 @@ """Miscellaneous utility interfaces.""" +import json +import os + +import h5py +import numpy as np +import pandas as pd from nipype import logging from nipype.interfaces.base import ( BaseInterfaceInputSpec, @@ -8,11 +14,16 @@ SimpleInterface, TraitedSpec, Undefined, + isdefined, traits, traits_extension, ) -from xcp_d.utils.modified_data import downcast_to_32 +from xcp_d.utils.confounds import load_motion +from xcp_d.utils.filemanip import fname_presuffix +from xcp_d.utils.modified_data import compute_fd, downcast_to_32 +from xcp_d.utils.qcmetrics import compute_dvars, compute_registration_qc +from xcp_d.utils.write_save import read_ndata LOGGER = logging.getLogger("nipype.interface") @@ -151,3 +162,399 @@ def _run_interface(self, runtime): outlist.append(item) self._results["outlist"] = outlist return runtime + + +class _LINCQCInputSpec(BaseInterfaceInputSpec): + name_source = File( + exists=False, + mandatory=True, + desc=( + "Preprocessed BOLD file. Used to find files. " + "In the case of the concatenation workflow, " + "this may be a nonexistent file " + "(i.e., the preprocessed BOLD file, with the run entity removed)." + ), + ) + bold_file = File( + exists=True, + mandatory=True, + desc="Preprocessed BOLD file, after dummy scan removal. Used in carpet plot.", + ) + dummy_scans = traits.Int(mandatory=True, desc="Dummy time to drop") + temporal_mask = traits.Either( + File(exists=True), + Undefined, + desc="Temporal mask", + ) + fmriprep_confounds_file = File( + exists=True, + mandatory=True, + desc="fMRIPrep confounds file, after dummy scans removal", + ) + cleaned_file = File( + exists=True, + mandatory=True, + desc="Processed file, after denoising and censoring.", + ) + TR = traits.Float(mandatory=True, desc="Repetition time, in seconds.") + head_radius = traits.Float(mandatory=True, desc="Head radius for FD calculation, in mm.") + bold_mask_inputspace = traits.Either( + None, + File(exists=True), + mandatory=True, + desc=( + "Mask file from NIfTI. May be None, for CIFTI processing. " + "The mask is in the same space as the BOLD data, which may not be the same as the " + "bold_mask_stdspace file. " + "Used to load the masked BOLD data. Not used for QC metrics." + ), + ) + + # Inputs used only for nifti data + anat_mask_anatspace = File( + exists=True, + mandatory=False, + desc=( + "Anatomically-derived brain mask in anatomical space. " + "Used to calculate coregistration QC metrics." + ), + ) + template_mask = File( + exists=True, + mandatory=False, + desc=( + "Template's official brain mask. " + "This matches the space of bold_mask_stdspace, " + "but does not necessarily match the space of bold_mask_inputspace. " + "Used to calculate normalization QC metrics." + ), + ) + bold_mask_anatspace = File( + exists=True, + mandatory=False, + desc="BOLD mask in anatomical space. Used to calculate coregistration QC metrics.", + ) + bold_mask_stdspace = File( + exists=True, + mandatory=False, + desc=( + "BOLD mask in template space. " + "This matches the space of template_mask, " + "but does not necessarily match the space of bold_mask_inputspace. " + "Used to calculate normalization QC metrics." + ), + ) + + +class _LINCQCOutputSpec(TraitedSpec): + qc_file = File(exists=True, desc="QC TSV file.") + qc_metadata = File(exists=True, desc="Sidecar JSON for QC TSV file.") + + +class LINCQC(SimpleInterface): + """Calculate QC metrics used by the LINC lab.""" + + input_spec = _LINCQCInputSpec + output_spec = _LINCQCOutputSpec + + def _run_interface(self, runtime): + # Load confound matrix and load motion without motion filtering + confounds_df = pd.read_table(self.inputs.fmriprep_confounds_file) + preproc_motion_df = load_motion( + confounds_df.copy(), + TR=self.inputs.TR, + motion_filter_type=None, + ) + preproc_fd = compute_fd(confound=preproc_motion_df, head_radius=self.inputs.head_radius) + rmsd = confounds_df["rmsd"].to_numpy() + + # Determine number of dummy volumes and load temporal mask + dummy_scans = self.inputs.dummy_scans + if isdefined(self.inputs.temporal_mask): + censoring_df = pd.read_table(self.inputs.temporal_mask) + tmask_arr = censoring_df["framewise_displacement"].values + else: + tmask_arr = np.zeros(preproc_fd.size, dtype=int) + + num_censored_volumes = int(tmask_arr.sum()) + num_retained_volumes = int((tmask_arr == 0).sum()) + + # Apply temporal mask to interpolated/full data + rmsd_censored = rmsd[tmask_arr == 0] + postproc_fd = preproc_fd[tmask_arr == 0] + + dvars_before_processing = compute_dvars( + datat=read_ndata( + datafile=self.inputs.bold_file, + maskfile=self.inputs.bold_mask_inputspace, + ), + )[1] + dvars_after_processing = compute_dvars( + datat=read_ndata( + datafile=self.inputs.cleaned_file, + maskfile=self.inputs.bold_mask_inputspace, + ), + )[1] + if preproc_fd.size != dvars_before_processing.size: + raise ValueError(f"FD {preproc_fd.size} != DVARS {dvars_before_processing.size}\n") + + # Get the different components in the bold file name + # eg: ['sub-colornest001', 'ses-1'], etc. + _, bold_file_name = os.path.split(self.inputs.name_source) + bold_file_name_components = bold_file_name.split("_") + + # Fill out dictionary with entities from filename + qc_values_dict = {} + for entity in bold_file_name_components[:-1]: + qc_values_dict[entity.split("-")[0]] = entity.split("-")[1] + + # Calculate QC measures + mean_fd = np.mean(preproc_fd) + mean_fd_post_censoring = np.mean(postproc_fd) + mean_relative_rms = np.nanmean(rmsd_censored) # first value can be NaN if no dummy scans + mean_dvars_before_processing = np.mean(dvars_before_processing) + mean_dvars_after_processing = np.mean(dvars_after_processing) + fd_dvars_correlation_initial = np.corrcoef(preproc_fd, dvars_before_processing)[0, 1] + fd_dvars_correlation_final = np.corrcoef(postproc_fd, dvars_after_processing)[0, 1] + rmsd_max_value = np.nanmax(rmsd_censored) + + # A summary of all the values + qc_values_dict.update( + { + "mean_fd": [mean_fd], + "mean_fd_post_censoring": [mean_fd_post_censoring], + "mean_relative_rms": [mean_relative_rms], + "max_relative_rms": [rmsd_max_value], + "mean_dvars_initial": [mean_dvars_before_processing], + "mean_dvars_final": [mean_dvars_after_processing], + "num_dummy_volumes": [dummy_scans], + "num_censored_volumes": [num_censored_volumes], + "num_retained_volumes": [num_retained_volumes], + "fd_dvars_correlation_initial": [fd_dvars_correlation_initial], + "fd_dvars_correlation_final": [fd_dvars_correlation_final], + } + ) + + qc_metadata = { + "mean_fd": { + "LongName": "Mean Framewise Displacement", + "Description": ( + "Average framewise displacement without any motion parameter filtering. " + "This value includes high-motion outliers, but not dummy volumes. " + "FD is calculated according to the Power definition." + ), + "Units": "mm / volume", + "Term URL": "https://doi.org/10.1016/j.neuroimage.2011.10.018", + }, + "mean_fd_post_censoring": { + "LongName": "Mean Framewise Displacement After Censoring", + "Description": ( + "Average framewise displacement without any motion parameter filtering. " + "This value does not include high-motion outliers or dummy volumes. " + "FD is calculated according to the Power definition." + ), + "Units": "mm / volume", + "Term URL": "https://doi.org/10.1016/j.neuroimage.2011.10.018", + }, + "mean_relative_rms": { + "LongName": "Mean Relative Root Mean Squared", + "Description": ( + "Average relative root mean squared calculated from motion parameters, " + "after removal of dummy volumes and high-motion outliers. " + "Relative in this case means 'relative to the previous scan'." + ), + "Units": "arbitrary", + }, + "max_relative_rms": { + "LongName": "Maximum Relative Root Mean Squared", + "Description": ( + "Maximum relative root mean squared calculated from motion parameters, " + "after removal of dummy volumes and high-motion outliers. " + "Relative in this case means 'relative to the previous scan'." + ), + "Units": "arbitrary", + }, + "mean_dvars_initial": { + "LongName": "Mean DVARS Before Postprocessing", + "Description": ( + "Average DVARS (temporal derivative of root mean squared variance over " + "voxels) calculated from the preprocessed BOLD file, after dummy scan removal." + ), + "TermURL": "https://doi.org/10.1016/j.neuroimage.2011.02.073", + }, + "mean_dvars_final": { + "LongName": "Mean DVARS After Postprocessing", + "Description": ( + "Average DVARS (temporal derivative of root mean squared variance over " + "voxels) calculated from the denoised BOLD file." + ), + "TermURL": "https://doi.org/10.1016/j.neuroimage.2011.02.073", + }, + "num_dummy_volumes": { + "LongName": "Number of Dummy Volumes", + "Description": ( + "The number of non-steady state volumes removed from the time series by XCP-D." + ), + }, + "num_censored_volumes": { + "LongName": "Number of Censored Volumes", + "Description": ( + "The number of high-motion outlier volumes censored by XCP-D. " + "This does not include dummy volumes." + ), + }, + "num_retained_volumes": { + "LongName": "Number of Retained Volumes", + "Description": ( + "The number of volumes retained in the denoised dataset. " + "This does not include dummy volumes or high-motion outliers." + ), + }, + "fd_dvars_correlation_initial": { + "LongName": "FD-DVARS Correlation Before Postprocessing", + "Description": ( + "The Pearson correlation coefficient between framewise displacement and DVARS " + "(temporal derivative of root mean squared variance over voxels), " + "after removal of dummy volumes, but before removal of high-motion outliers." + ), + }, + "fd_dvars_correlation_final": { + "LongName": "FD-DVARS Correlation After Postprocessing", + "Description": ( + "The Pearson correlation coefficient between framewise displacement and DVARS " + "(temporal derivative of root mean squared variance over voxels), " + "after postprocessing. " + "The FD time series is unfiltered, but censored. " + "The DVARS time series is calculated from the denoised BOLD data." + ), + }, + } + + if self.inputs.bold_mask_anatspace: # If a bold mask in T1w is provided + # Compute quality of registration + registration_qc, registration_metadata = compute_registration_qc( + bold_mask_anatspace=self.inputs.bold_mask_anatspace, + anat_mask_anatspace=self.inputs.anat_mask_anatspace, + bold_mask_stdspace=self.inputs.bold_mask_stdspace, + template_mask=self.inputs.template_mask, + ) + qc_values_dict.update(registration_qc) # Add values to dictionary + qc_metadata.update(registration_metadata) + + # Convert dictionary to df and write out the qc file + df = pd.DataFrame(qc_values_dict) + self._results["qc_file"] = fname_presuffix( + self.inputs.cleaned_file, + suffix="qc_bold.tsv", + newpath=runtime.cwd, + use_ext=False, + ) + df.to_csv(self._results["qc_file"], index=False, header=True, sep="\t") + + # Write out the metadata file + self._results["qc_metadata"] = fname_presuffix( + self.inputs.cleaned_file, + suffix="qc_bold.json", + newpath=runtime.cwd, + use_ext=False, + ) + with open(self._results["qc_metadata"], "w") as fo: + json.dump(qc_metadata, fo, indent=4, sort_keys=True) + + return runtime + + +class _ABCCQCInputSpec(BaseInterfaceInputSpec): + filtered_motion = File( + exists=True, + mandatory=True, + desc="", + ) + TR = traits.Float(mandatory=True, desc="Repetition Time") + + +class _ABCCQCOutputSpec(TraitedSpec): + qc_file = File(exists=True, desc="ABCC QC HDF5 file.") + + +class ABCCQC(SimpleInterface): + """Create an HDF5-format file containing a DCAN-format dataset. + + Notes + ----- + The metrics in the file are: + + - ``FD_threshold``: a number >= 0 that represents the FD threshold used to calculate + the metrics in this list. + - ``frame_removal``: a binary vector/array the same length as the number of frames + in the concatenated time series, indicates whether a frame is removed (1) or not (0) + - ``format_string`` (legacy): a string that denotes how the frames were excluded. + This uses a notation devised by Avi Snyder. + - ``total_frame_count``: a whole number that represents the total number of frames + in the concatenated series + - ``remaining_frame_count``: a whole number that represents the number of remaining + frames in the concatenated series + - ``remaining_seconds``: a whole number that represents the amount of time remaining + after thresholding + - ``remaining_frame_mean_FD``: a number >= 0 that represents the mean FD of the + remaining frames + """ + + input_spec = _ABCCQCInputSpec + output_spec = _ABCCQCOutputSpec + + def _run_interface(self, runtime): + TR = self.inputs.TR + + self._results["qc_file"] = fname_presuffix( + self.inputs.filtered_motion, + suffix="qc_bold.hdf5", + newpath=runtime.cwd, + use_ext=False, + ) + + # Load filtered framewise_displacement values from file + filtered_motion_df = pd.read_table(self.inputs.filtered_motion) + fd = filtered_motion_df["framewise_displacement"].values + + with h5py.File(self._results["qc_file"], "w") as dcan: + for thresh in np.linspace(0, 1, 101): + thresh = np.around(thresh, 2) + + dcan.create_dataset( + f"/dcan_motion/fd_{thresh}/skip", + data=0, + dtype="float", + ) + dcan.create_dataset( + f"/dcan_motion/fd_{thresh}/binary_mask", + data=(fd > thresh).astype(int), + dtype="float", + ) + dcan.create_dataset( + f"/dcan_motion/fd_{thresh}/threshold", + data=thresh, + dtype="float", + ) + dcan.create_dataset( + f"/dcan_motion/fd_{thresh}/total_frame_count", + data=len(fd), + dtype="float", + ) + dcan.create_dataset( + f"/dcan_motion/fd_{thresh}/remaining_total_frame_count", + data=len(fd[fd <= thresh]), + dtype="float", + ) + dcan.create_dataset( + f"/dcan_motion/fd_{thresh}/remaining_seconds", + data=len(fd[fd <= thresh]) * TR, + dtype="float", + ) + dcan.create_dataset( + f"/dcan_motion/fd_{thresh}/remaining_frame_mean_FD", + data=(fd[fd <= thresh]).mean(), + dtype="float", + ) + + return runtime diff --git a/xcp_d/interfaces/workbench.py b/xcp_d/interfaces/workbench.py index fc99b3dac..df879037e 100644 --- a/xcp_d/interfaces/workbench.py +++ b/xcp_d/interfaces/workbench.py @@ -1224,74 +1224,6 @@ def _list_outputs(self): return outputs -class _CiftiChangeMappingInputSpec(CommandLineInputSpec): - """Input specification for the CiftiChangeMapping command.""" - - data_cifti = File( - exists=True, - mandatory=True, - argstr="%s", - position=0, - desc="The cifti file to use the data from.", - ) - direction = traits.Enum( - "ROW", - "COLUMN", - mandatory=True, - argstr="%s", - position=1, - desc="Which mapping to parcellate (integer, ROW, or COLUMN)", - ) - cifti_out = File( - name_source=["label"], - name_template="converted_%s.dscalar.nii", - keep_extension=False, - argstr="%s", - position=2, - desc="The output cifti file.", - ) - scalar = traits.Bool( - False, - usedefault=True, - argstr="-scalar", - position=3, - desc="Set the mapping to scalar", - ) - - -class _CiftiChangeMappingOutputSpec(TraitedSpec): - """Output specification for the CiftiChangeMapping command.""" - - cifti_out = File(exists=True, desc="output CIFTI file") - - -class CiftiChangeMapping(WBCommand): - """Convert to scalar, copy mapping, etc. - - Take an existing cifti file and change one of the mappings. - Exactly one of -series, -scalar, or -from-cifti must be specified. - The direction can be either an integer starting from 1, or the strings 'ROW' or 'COLUMN'. - - Examples - -------- - >>> ccdft = CiftiChangeMapping() - >>> ccdft.inputs.data_cifti = "tpl-fsLR_atlas-Gordon_den-32k_dseg.dlabel.nii" - >>> ccdft.inputs.direction = "ROW" - >>> ccdft.inputs.cifti_out = "out.dscalar.nii" - >>> ccdft.inputs.scalar = True - >>> ccdft.cmdline - wb_command -cifti-change-mapping \ - tpl-fsLR_atlas-Gordon_den-32k_dseg.dlabel.nii \ - ROW \ - out.dscalar.nii \ - -scalar - """ - - input_spec = _CiftiChangeMappingInputSpec - output_spec = _CiftiChangeMappingOutputSpec - _cmd = "wb_command -cifti-change-mapping" - - class _CiftiMathInputSpec(CommandLineInputSpec): """Input specification for the CiftiMath command.""" diff --git a/xcp_d/utils/qcmetrics.py b/xcp_d/utils/qcmetrics.py index 340e521c8..32efe925a 100644 --- a/xcp_d/utils/qcmetrics.py +++ b/xcp_d/utils/qcmetrics.py @@ -1,17 +1,18 @@ """Quality control metrics.""" -import h5py import nibabel as nb import numpy as np -import pandas as pd from nipype import logging -from xcp_d.utils.doc import fill_doc - LOGGER = logging.getLogger("nipype.utils") -def compute_registration_qc(bold2t1w_mask, anat_brainmask, bold2template_mask, template_mask): +def compute_registration_qc( + bold_mask_anatspace, + anat_mask_anatspace, + bold_mask_stdspace, + template_mask, +): """Compute quality of registration metrics. This function will calculate a series of metrics, including: @@ -25,14 +26,14 @@ def compute_registration_qc(bold2t1w_mask, anat_brainmask, bold2template_mask, t Parameters ---------- - bold2t1w_mask : :obj:`str` - Path to the BOLD mask in anatomical (T1w or T2w) space. - anat_brainmask : :obj:`str` + bold_mask_anatspace : :obj:`str` + Path to the BOLD brain mask in anatomical (T1w or T2w) space. + anat_mask_anatspace : :obj:`str` Path to the anatomically-derived brain mask in anatomical space. - bold2template_mask : :obj:`str` - Path to the BOLD mask in template space. + bold_mask_stdspace : :obj:`str` + Path to the BOLD brain mask in template space. template_mask : :obj:`str` - Path to the template mask. + Path to the template's official brain mask. Returns ------- @@ -41,18 +42,18 @@ def compute_registration_qc(bold2t1w_mask, anat_brainmask, bold2template_mask, t qc_metadata : dict Metadata describing the QC measures. """ - bold2t1w_mask_arr = nb.load(bold2t1w_mask).get_fdata() - t1w_mask_arr = nb.load(anat_brainmask).get_fdata() - bold2template_mask_arr = nb.load(bold2template_mask).get_fdata() + bold_mask_anatspace_arr = nb.load(bold_mask_anatspace).get_fdata() + anat_mask_anatspace_arr = nb.load(anat_mask_anatspace).get_fdata() + bold_mask_stdspace_arr = nb.load(bold_mask_stdspace).get_fdata() template_mask_arr = nb.load(template_mask).get_fdata() reg_qc = { - "coreg_dice": [dice(bold2t1w_mask_arr, t1w_mask_arr)], - "coreg_correlation": [pearson(bold2t1w_mask_arr, t1w_mask_arr)], - "coreg_overlap": [overlap(bold2t1w_mask_arr, t1w_mask_arr)], - "norm_dice": [dice(bold2template_mask_arr, template_mask_arr)], - "norm_correlation": [pearson(bold2template_mask_arr, template_mask_arr)], - "norm_overlap": [overlap(bold2template_mask_arr, template_mask_arr)], + "coreg_dice": [dice(bold_mask_anatspace_arr, anat_mask_anatspace_arr)], + "coreg_correlation": [pearson(bold_mask_anatspace_arr, anat_mask_anatspace_arr)], + "coreg_overlap": [overlap(bold_mask_anatspace_arr, anat_mask_anatspace_arr)], + "norm_dice": [dice(bold_mask_stdspace_arr, template_mask_arr)], + "norm_correlation": [pearson(bold_mask_stdspace_arr, template_mask_arr)], + "norm_overlap": [overlap(bold_mask_stdspace_arr, template_mask_arr)], } qc_metadata = { "coreg_dice": { @@ -290,107 +291,3 @@ def compute_dvars( dvars_stdz = np.insert(dvars_stdz, 0, 0) return dvars_nstd, dvars_stdz - - -def make_abcc_qc_file(filtered_motion, TR): - """Make DCAN HDF5 file from single motion file. - - NOTE: This is a Node function. - - Parameters - ---------- - filtered_motion_file : :obj:`str` - File from which to extract information. - TR : :obj:`float` - Repetition time. - - Returns - ------- - dcan_df_file : :obj:`str` - Name of the HDF5-format file that is created. - """ - import os - - from xcp_d.utils.qcmetrics import make_dcan_df - - dcan_df_file = os.path.abspath("desc-abcc_qc.hdf5") - - make_dcan_df(filtered_motion, dcan_df_file, TR) - return dcan_df_file - - -@fill_doc -def make_dcan_df(filtered_motion, name, TR): - """Create an HDF5-format file containing a DCAN-format dataset. - - Parameters - ---------- - %(filtered_motion)s - name : :obj:`str` - Name of the HDF5-format file to be created. - %(TR)s - - Notes - ----- - The metrics in the file are: - - - ``FD_threshold``: a number >= 0 that represents the FD threshold used to calculate - the metrics in this list. - - ``frame_removal``: a binary vector/array the same length as the number of frames - in the concatenated time series, indicates whether a frame is removed (1) or not (0) - - ``format_string`` (legacy): a string that denotes how the frames were excluded. - This uses a notation devised by Avi Snyder. - - ``total_frame_count``: a whole number that represents the total number of frames - in the concatenated series - - ``remaining_frame_count``: a whole number that represents the number of remaining - frames in the concatenated series - - ``remaining_seconds``: a whole number that represents the amount of time remaining - after thresholding - - ``remaining_frame_mean_FD``: a number >= 0 that represents the mean FD of the - remaining frames - """ - LOGGER.debug(f"Generating DCAN file: {name}") - - # Load filtered framewise_displacement values from file - filtered_motion_df = pd.read_table(filtered_motion) - fd = filtered_motion_df["framewise_displacement"].values - - with h5py.File(name, "w") as dcan: - for thresh in np.linspace(0, 1, 101): - thresh = np.around(thresh, 2) - - dcan.create_dataset( - f"/dcan_motion/fd_{thresh}/skip", - data=0, - dtype="float", - ) - dcan.create_dataset( - f"/dcan_motion/fd_{thresh}/binary_mask", - data=(fd > thresh).astype(int), - dtype="float", - ) - dcan.create_dataset( - f"/dcan_motion/fd_{thresh}/threshold", - data=thresh, - dtype="float", - ) - dcan.create_dataset( - f"/dcan_motion/fd_{thresh}/total_frame_count", - data=len(fd), - dtype="float", - ) - dcan.create_dataset( - f"/dcan_motion/fd_{thresh}/remaining_total_frame_count", - data=len(fd[fd <= thresh]), - dtype="float", - ) - dcan.create_dataset( - f"/dcan_motion/fd_{thresh}/remaining_seconds", - data=len(fd[fd <= thresh]) * TR, - dtype="float", - ) - dcan.create_dataset( - f"/dcan_motion/fd_{thresh}/remaining_frame_mean_FD", - data=(fd[fd <= thresh]).mean(), - dtype="float", - ) diff --git a/xcp_d/workflows/plotting.py b/xcp_d/workflows/plotting.py index 1acbc0cc4..b4e554e95 100644 --- a/xcp_d/workflows/plotting.py +++ b/xcp_d/workflows/plotting.py @@ -11,8 +11,8 @@ from xcp_d.interfaces.bids import DerivativesDataSink from xcp_d.interfaces.plotting import QCPlots, QCPlotsES from xcp_d.interfaces.report import FunctionalSummary +from xcp_d.interfaces.utils import ABCCQC, LINCQC from xcp_d.utils.doc import fill_doc -from xcp_d.utils.qcmetrics import make_abcc_qc_file from xcp_d.utils.utils import get_bold2std_and_t1w_xfms, get_std2bold_xfms @@ -157,7 +157,6 @@ def init_qc_report_wf( n_procs=omp_nthreads, mem_gb=1, ) - workflow.connect([ (inputnode, warp_boldmask_to_t1w, [ ("bold_mask", "input_image"), @@ -179,7 +178,6 @@ def init_qc_report_wf( n_procs=omp_nthreads, mem_gb=1, ) - workflow.connect([ (inputnode, warp_boldmask_to_mni, [("bold_mask", "input_image")]), (get_native2space_transforms, warp_boldmask_to_mni, [ @@ -220,7 +218,6 @@ def init_qc_report_wf( ), name="get_std2native_transform", ) - workflow.connect([(inputnode, get_mni_to_bold_xfms, [("name_source", "bold_file")])]) # Use MNI152NLin2009cAsym tissue-type segmentation file for carpet plots. @@ -268,26 +265,24 @@ def init_qc_report_wf( n_procs=omp_nthreads, mem_gb=3, ) - workflow.connect([ (inputnode, warp_dseg_to_bold, [("boldref", "reference_image")]), (add_xfm_to_nlin6asym, warp_dseg_to_bold, [("out", "transforms")]), ]) # fmt:skip if config.workflow.linc_qc: - qc_report = pe.Node( - QCPlots( + make_linc_qc = pe.Node( + LINCQC( TR=TR, head_radius=head_radius, template_mask=nlin2009casym_brain_mask, ), - name="qc_report", + name="make_linc_qc", mem_gb=2, n_procs=omp_nthreads, ) - workflow.connect([ - (inputnode, qc_report, [ + (inputnode, make_linc_qc, [ ("name_source", "name_source"), ("preprocessed_bold", "bold_file"), ("censored_denoised_bold", "cleaned_file"), @@ -295,19 +290,18 @@ def init_qc_report_wf( ("temporal_mask", "temporal_mask"), ("dummy_scans", "dummy_scans"), ]), - (qc_report, outputnode, [("qc_file", "qc_file")]), + (make_linc_qc, outputnode, [("qc_file", "qc_file")]), ]) # fmt:skip if config.workflow.file_format == "nifti": workflow.connect([ - (inputnode, qc_report, [("bold_mask", "mask_file")]), - (warp_dseg_to_bold, qc_report, [("output_image", "seg_file")]), - (warp_boldmask_to_t1w, qc_report, [("output_image", "bold2T1w_mask")]), - (warp_boldmask_to_mni, qc_report, [("output_image", "bold2temp_mask")]), - (warp_anatmask_to_t1w, qc_report, [("output_image", "anat_brainmask")]), + (inputnode, make_linc_qc, [("bold_mask", "bold_mask_inputspace")]), + (warp_boldmask_to_t1w, make_linc_qc, [("output_image", "bold_mask_anatspace")]), + (warp_boldmask_to_mni, make_linc_qc, [("output_image", "bold_mask_stdspace")]), + (warp_anatmask_to_t1w, make_linc_qc, [("output_image", "anat_mask_anatspace")]), ]) # fmt:skip else: - qc_report.inputs.mask_file = None + make_linc_qc.inputs.bold_mask_inputspace = None ds_qc_metadata = pe.Node( DerivativesDataSink( @@ -321,40 +315,60 @@ def init_qc_report_wf( name="ds_qc_metadata", run_without_submitting=True, ) - workflow.connect([ (inputnode, ds_qc_metadata, [("name_source", "source_file")]), - (qc_report, ds_qc_metadata, [("qc_metadata", "in_file")]), + (make_linc_qc, ds_qc_metadata, [("qc_metadata", "in_file")]), + ]) # fmt:skip + + make_qc_plots_nipreps = pe.Node( + QCPlots(TR=TR, head_radius=head_radius), + name="make_qc_plots_nipreps", + mem_gb=2, + n_procs=omp_nthreads, + ) + workflow.connect([ + (inputnode, make_qc_plots_nipreps, [ + ("preprocessed_bold", "bold_file"), + ("censored_denoised_bold", "cleaned_file"), + ("fmriprep_confounds_file", "fmriprep_confounds_file"), + ("temporal_mask", "temporal_mask"), + ]), ]) # fmt:skip - ds_report_preprocessing = pe.Node( + if config.workflow.file_format == "nifti": + workflow.connect([ + (inputnode, make_qc_plots_nipreps, [("bold_mask", "mask_file")]), + (warp_dseg_to_bold, make_qc_plots_nipreps, [("output_image", "seg_file")]), + ]) # fmt:skip + else: + make_qc_plots_nipreps.inputs.mask_file = None + + ds_preproc_qc_plot_nipreps = pe.Node( DerivativesDataSink( base_directory=output_dir, desc="preprocessing", datatype="figures", ), - name="ds_report_preprocessing", + name="ds_preproc_qc_plot_nipreps", run_without_submitting=False, ) - workflow.connect([ - (inputnode, ds_report_preprocessing, [("name_source", "source_file")]), - (qc_report, ds_report_preprocessing, [("raw_qcplot", "in_file")]), + (inputnode, ds_preproc_qc_plot_nipreps, [("name_source", "source_file")]), + (make_qc_plots_nipreps, ds_preproc_qc_plot_nipreps, [("raw_qcplot", "in_file")]), ]) # fmt:skip - ds_report_postprocessing = pe.Node( + ds_postproc_qc_plot_nipreps = pe.Node( DerivativesDataSink( base_directory=output_dir, desc="postprocessing", datatype="figures", ), - name="ds_report_postprocessing", + name="ds_postproc_qc_plot_nipreps", run_without_submitting=False, ) - workflow.connect([ - (inputnode, ds_report_postprocessing, [("name_source", "source_file")]), - (qc_report, ds_report_postprocessing, [("clean_qcplot", "in_file")]), + (inputnode, ds_postproc_qc_plot_nipreps, [("name_source", "source_file")]), + (make_qc_plots_nipreps, ds_postproc_qc_plot_nipreps, [("clean_qcplot", "in_file")]), ]) # fmt:skip functional_qc = pe.Node( @@ -363,10 +377,9 @@ def init_qc_report_wf( run_without_submitting=False, mem_gb=2, ) - workflow.connect([ (inputnode, functional_qc, [("name_source", "bold_file")]), - (qc_report, functional_qc, [("qc_file", "qc_file")]), + (make_linc_qc, functional_qc, [("qc_file", "qc_file")]), ]) # fmt:skip ds_report_qualitycontrol = pe.Node( @@ -378,7 +391,6 @@ def init_qc_report_wf( name="ds_report_qualitycontrol", run_without_submitting=False, ) - workflow.connect([ (inputnode, ds_report_qualitycontrol, [("name_source", "source_file")]), (functional_qc, ds_report_qualitycontrol, [("out_report", "in_file")]), @@ -388,19 +400,13 @@ def init_qc_report_wf( workflow.add_nodes([outputnode]) if config.workflow.abcc_qc: - make_abcc_qc_file_node = pe.Node( - Function( - input_names=["filtered_motion", "TR"], - output_names=["dcan_df_file"], - function=make_abcc_qc_file, - ), - name="make_abcc_qc_file_node", + make_abcc_qc = pe.Node( + ABCCQC(TR=TR), + name="make_abcc_qc", + mem_gb=2, + n_procs=omp_nthreads, ) - make_abcc_qc_file_node.inputs.TR = TR - - workflow.connect([ - (inputnode, make_abcc_qc_file_node, [("filtered_motion", "filtered_motion")]), - ]) # fmt:skip + workflow.connect([(inputnode, make_abcc_qc, [("filtered_motion", "filtered_motion")])]) ds_abcc_qc = pe.Node( DerivativesDataSink( @@ -413,22 +419,20 @@ def init_qc_report_wf( name="ds_abcc_qc", run_without_submitting=True, ) - workflow.connect([ (inputnode, ds_abcc_qc, [("name_source", "source_file")]), - (make_abcc_qc_file_node, ds_abcc_qc, [("dcan_df_file", "in_file")]), + (make_abcc_qc, ds_abcc_qc, [("qc_file", "in_file")]), ]) # fmt:skip # Generate preprocessing and postprocessing carpet plots. - plot_execsummary_carpets_dcan = pe.Node( + make_qc_plots_es = pe.Node( QCPlotsES(TR=TR, standardize=config.workflow.params == "none"), - name="plot_execsummary_carpets_dcan", + name="make_qc_plots_es", mem_gb=2, n_procs=omp_nthreads, ) - workflow.connect([ - (inputnode, plot_execsummary_carpets_dcan, [ + (inputnode, make_qc_plots_es, [ ("preprocessed_bold", "preprocessed_bold"), ("denoised_interpolated_bold", "denoised_interpolated_bold"), ("filtered_motion", "filtered_motion"), @@ -439,46 +443,38 @@ def init_qc_report_wf( if config.workflow.file_format == "nifti": workflow.connect([ - (inputnode, plot_execsummary_carpets_dcan, [("bold_mask", "mask")]), - (warp_dseg_to_bold, plot_execsummary_carpets_dcan, [ - ("output_image", "seg_data"), - ]), + (inputnode, make_qc_plots_es, [("bold_mask", "mask")]), + (warp_dseg_to_bold, make_qc_plots_es, [("output_image", "seg_data")]), ]) # fmt:skip - ds_preproc_execsummary_carpet_dcan = pe.Node( + ds_preproc_qc_plot_es = pe.Node( DerivativesDataSink( base_directory=output_dir, dismiss_entities=["den"], datatype="figures", desc="preprocESQC", ), - name="ds_preproc_execsummary_carpet_dcan", + name="ds_preproc_qc_plot_es", run_without_submitting=True, ) - workflow.connect([ - (inputnode, ds_preproc_execsummary_carpet_dcan, [("name_source", "source_file")]), - (plot_execsummary_carpets_dcan, ds_preproc_execsummary_carpet_dcan, [ - ("before_process", "in_file"), - ]), + (inputnode, ds_preproc_qc_plot_es, [("name_source", "source_file")]), + (make_qc_plots_es, ds_preproc_qc_plot_es, [("before_process", "in_file")]), ]) # fmt:skip - ds_postproc_execsummary_carpet_dcan = pe.Node( + ds_postproc_qc_plot_es = pe.Node( DerivativesDataSink( base_directory=output_dir, dismiss_entities=["den"], datatype="figures", desc="postprocESQC", ), - name="ds_postproc_execsummary_carpet_dcan", + name="ds_postproc_qc_plot_es", run_without_submitting=True, ) - workflow.connect([ - (inputnode, ds_postproc_execsummary_carpet_dcan, [("name_source", "source_file")]), - (plot_execsummary_carpets_dcan, ds_postproc_execsummary_carpet_dcan, [ - ("after_process", "in_file"), - ]), + (inputnode, ds_postproc_qc_plot_es, [("name_source", "source_file")]), + (make_qc_plots_es, ds_postproc_qc_plot_es, [("after_process", "in_file")]), ]) # fmt:skip return workflow