"""
Generate XCP-stype quality control files.
.. seealso::
`User Guide: Generate eXtensible Connectivity Pipeline-style quality control files <https://fcp-indi.github.io/docs/user/xcpqc>`_
Columns
-------
sub : str
subject label :footcite:`BIDS21`
ses : str
session label :footcite:`BIDS21`
task : str
task label :footcite:`BIDS21`
run : int
run index :footcite:`BIDS21`
desc : str
description :footcite:`BIDS21`
regressors : str
'Name' of regressors in the current fork
space : str
space label :footcite:`BIDS21`
meanFD : float
mean Jenkinson framewise displacement :footcite:`xcp_22,Jenk02` :func:`CPAC.generate_motion_statistics.calculate_FD_J` after preprocessing
relMeansRMSMotion : float
"mean value of RMS motion" :footcite:`xcp_22,Ciri19`
relMaxRMSMotion : float
"maximum vaue of RMS motion" :footcite:`xcp_22,Ciri19`
meanDVInit : float
"mean DVARS" :footcite:`xcp_22,Ciri19`
meanDVFinal : float
"mean DVARS" :footcite:`xcp_22,Ciri19`
nVolCensored : int
"total number of volume(s) censored :footcite:`Ciri19`
nVolsRemoved : int
number of volumes in derivative minus number of volumes in original
functional scan
motionDVCorrInit : float
"correlation of RMS and DVARS before regresion" :footcite:`Ciri19`
motionDVCorrFinal : float
"correlation of RMS and DVARS after regresion" :footcite:`Ciri19`
coregDice : float
"Coregsitration of Functional and T1w:[…] Dice index" :footcite:`xcp_22,Ciri19`
coregJaccard : float
"Coregsitration of Functional and T1w:[…] Jaccard index" :footcite:`xcp_22,Ciri19`
coregCrossCorr : float
"Coregsitration of Functional and T1w:[…] cross correlation" :footcite:`xcp_22,Ciri19`
coregCoverag : float
"Coregsitration of Functional and T1w:[…] Coverage index" :footcite:`xcp_22,Ciri19`
normDice : float
"Normalization of T1w/Functional to Template:[…] Dice index" :footcite:`xcp_22,Ciri19`
normJaccard : float
"Normalization of T1w/Functional to Template:[…] Jaccard index" :footcite:`xcp_22,Ciri19`
normCrossCorr : float
"Normalization of T1w/Functional to Template:[…] cross correlation" :footcite:`xcp_22,Ciri19`
normCoverage : float
"Normalization of T1w/Functional to Template:[…] Coverage index" :footcite:`xcp_22,Ciri19`
""" # pylint: disable=line-too-long
from io import BufferedReader
import os
import re
from bids.layout import parse_file_entities
import numpy as np
import pandas as pd
import nibabel as nib
from nipype.interfaces import afni, fsl
from CPAC.generate_motion_statistics.generate_motion_statistics import (
DVARS_strip_t0,
ImageTo1D,
)
from CPAC.pipeline import nipype_pipeline_engine as pe
from CPAC.pipeline.nodeblock import nodeblock
from CPAC.qc.qcmetrics import regisQ
from CPAC.utils.interfaces.function import Function
motion_params = [
"dvars",
"framewise-displacement-jenkinson",
"desc-movementParametersUnfiltered_motion",
"desc-movementParameters_motion",
]
def _connect_motion(wf, nodes, strat_pool, qc_file, pipe_num):
"""
Connect the motion metrics to the workflow.
Parameters
----------
wf : nipype.pipeline.engine.Workflow
The workflow to connect the motion metrics to.
nodes : dict
Dictionary of nodes already collected from the strategy pool.
strat_pool : CPAC.pipeline.engine.ResourcePool
The current strategy pool.
qc_file : nipype.pipeline.engine.Node
A function node with the function ``generate_xcp_qc``.
pipe_num : int
Returns
-------
wf : nipype.pipeline.engine.Workflow
"""
# pylint: disable=invalid-name, too-many-arguments
try:
nodes = {**nodes, "censor-indices": strat_pool.node_data("censor-indices")}
wf.connect(
nodes["censor-indices"].node,
nodes["censor-indices"].out,
qc_file,
"censor_indices",
)
except LookupError:
qc_file.inputs.censor_indices = []
cal_DVARS = pe.Node(
ImageTo1D(method="dvars"),
name=f"cal_DVARS_{pipe_num}",
mem_gb=0.4,
mem_x=(739971956005215 / 151115727451828646838272, "in_file"),
throttle=True,
)
cal_DVARS_strip = pe.Node(
Function(
input_names=["file_1D"],
output_names=["out_file"],
function=DVARS_strip_t0,
as_module=True,
),
name=f"cal_DVARS_strip_{pipe_num}",
)
wf.connect(
[
(
nodes["desc-preproc_bold"].node,
cal_DVARS,
[(nodes["desc-preproc_bold"].out, "in_file")],
),
(
nodes["space-bold_desc-brain_mask"].node,
cal_DVARS,
[(nodes["space-bold_desc-brain_mask"].out, "mask")],
),
(cal_DVARS, cal_DVARS_strip, [("out_file", "file_1D")]),
(cal_DVARS_strip, qc_file, [("out_file", "dvars_after")]),
*[
(nodes[node].node, qc_file, [(nodes[node].out, node.replace("-", "_"))])
for node in motion_params
if node in nodes
],
]
)
return wf
[docs]
def dvcorr(dvars, fdj):
"""Correlate DVARS and FD-J."""
dvars = np.loadtxt(dvars)
fdj = np.loadtxt(fdj)
if len(dvars) != len(fdj) - 1:
msg = (
"len(DVARS) should be 1 less than len(FDJ), but their respective "
f"lengths are {len(dvars)} and {len(fdj)}."
)
raise ValueError(msg)
return np.corrcoef(dvars, fdj[1:])[0, 1]
# This function is for a function node for which
# Nipype will connect many other nodes as inputs
[docs]
def generate_xcp_qc( # noqa: PLR0913
sub,
ses,
task,
run,
desc,
regressors,
bold2t1w_mask,
t1w_mask,
bold2template_mask,
template_mask,
original_func,
final_func,
movement_parameters,
dvars,
censor_indices,
framewise_displacement_jenkinson,
dvars_after,
template,
):
"""
Generate an RBC-style QC CSV.
Parameters
----------
sub : str
subject ID
ses : str
session ID
task : str
task ID
run : str or int
run ID
desc : str
description string
regressors : str
'Name' of regressors in fork
original_func : str
path to original 'bold' image
final_bold : str
path to 'space-template_desc-preproc_bold' image
bold2t1w_mask : str
path to bold-to-T1w transform applied to space-bold_desc-brain_mask
with space-T1w_desc-brain_mask reference
t1w_mask : str
path to space-T1w_desc-brain_mask
bold2template_mask : str
path to space-template_desc-bold_mask
template_mask : str
path to T1w-brain-template-mask or EPI-template-mask
movement_parameters: str
path to movement parameters
dvars : str
path to DVARS before motion correction
censor_indices : list
list of indices of censored volumes
framewise_displacement_jenkinson : str
path to framewise displacement (Jenkinson) before motion correction
dvars_after : str
path to DVARS on final 'bold' image
template : str
path to registration template
Returns
-------
str
path to space-template_desc-xcp_quality TSV
"""
columns = (
"sub,ses,task,run,desc,regressors,space,meanFD,relMeansRMSMotion,"
"relMaxRMSMotion,meanDVInit,meanDVFinal,nVolCensored,nVolsRemoved,"
"motionDVCorrInit,motionDVCorrFinal,coregDice,coregJaccard,"
"coregCrossCorr,coregCoverage,normDice,normJaccard,normCrossCorr,"
"normCoverage".split(",")
)
images = {
"original_func": nib.load(original_func),
"final_func": nib.load(final_func),
}
# `sub` through `space`
from_bids = {
"sub": sub,
"ses": ses,
"task": task,
"run": run,
"desc": desc,
"regressors": regressors,
"space": os.path.basename(template).split(".", 1)[0].split("_", 1)[0],
}
if from_bids["space"].startswith("tpl-"):
from_bids["space"] = from_bids["space"][4:]
# `nVolCensored` & `nVolsRemoved`
n_vols_censored = len(censor_indices) if censor_indices is not None else "unknown"
shape_params = {
"nVolCensored": n_vols_censored,
"nVolsRemoved": images["original_func"].shape[3]
- images["final_func"].shape[3],
}
if isinstance(final_func, BufferedReader):
final_func = final_func.name
qc_filepath = os.path.join(os.getcwd(), "xcpqc.tsv")
desc_span = re.search(r"_desc-.*_", final_func)
if desc_span:
desc_span = desc_span.span()
final_func = "_".join([final_func[: desc_span[0]], final_func[desc_span[1] :]])
del desc_span
# `meanFD (Jenkinson)`
power_params = {"meanFD": np.mean(np.loadtxt(framewise_displacement_jenkinson))}
# `relMeansRMSMotion` & `relMaxRMSMotion`
mot = np.genfromtxt(movement_parameters).T
# Relative RMS of translation
rms = np.sqrt(mot[3] ** 2 + mot[4] ** 2 + mot[5] ** 2)
rms_params = {"relMeansRMSMotion": [np.mean(rms)], "relMaxRMSMotion": [np.max(rms)]}
# `meanDVInit` & `meanDVFinal`
meanDV = {"meanDVInit": np.mean(np.loadtxt(dvars))}
try:
meanDV["motionDVCorrInit"] = dvcorr(dvars, framewise_displacement_jenkinson)
except ValueError as value_error:
meanDV["motionDVCorrInit"] = f"ValueError({value_error!s})"
meanDV["meanDVFinal"] = np.mean(np.loadtxt(dvars_after))
try:
meanDV["motionDVCorrFinal"] = dvcorr(
dvars_after, framewise_displacement_jenkinson
)
except ValueError as value_error:
meanDV["motionDVCorrFinal"] = f"ValueError({value_error!s})"
# Overlap
overlap_params = regisQ(
bold2t1w_mask=bold2t1w_mask,
t1w_mask=t1w_mask,
bold2template_mask=bold2template_mask,
template_mask=template_mask,
)
qc_dict = {
**from_bids,
**power_params,
**rms_params,
**shape_params,
**overlap_params,
**meanDV,
}
df = pd.DataFrame(qc_dict, columns=columns)
df.to_csv(qc_filepath, sep="\t", index=False)
return qc_filepath
[docs]
def get_bids_info(subject, scan, wf_name):
"""
Gather BIDS information from a strat_pool.
Parameters
----------
subject : str
subject ID
scan : str
scan ID
wf_name : str
workflow name
Returns
-------
subject : str
subject ID
session : str
session ID
task : str
task ID
run : str or int
run ID
Examples
--------
>>> subject, session, task, run = get_bids_info(
... subject='DavidBowman', scan='rest_acq-1_run-1',
... wf_name='cpac_DavidBowman_3')
>>> subject
'DavidBowman'
>>> session
'3'
>>> task
'rest'
>>> run
'1'
>>> get_bids_info(subject='sub-colornest035', scan='rest_run-01',
... wf_name='cpac_sub-colornest035_ses-1')
('colornest035', '1', 'rest', '01')
"""
returns = ("subject", "session", "task", "run")
ses = wf_name.split("_")[-1]
if not ses.startswith("ses-"):
ses = f"ses-{ses}"
if not subject.startswith("sub-"):
subject = f"sub-{subject}"
resource = "_".join([subject, ses, scan if "task" in scan else f"task-{scan}"])
entities = parse_file_entities(resource)
returns = {key: entities.get(key) for key in returns}
if any(value is None for value in returns.values()):
entity_parts = resource.split("_")
def get_entity_part(key):
key = f"{key}-"
matching_parts = [part for part in entity_parts if part.startswith(key)]
if matching_parts:
return matching_parts[0].replace(key, "")
return None
for key, value in returns.items():
if value is None:
if key == "task":
returns[key] = get_entity_part(key)
else:
returns[key] = get_entity_part(key[:3])
return tuple(str(returns.get(key)) for key in ["subject", "session", "task", "run"])
[docs]
@nodeblock(
name="qc_xcp",
config=["pipeline_setup", "output_directory", "quality_control"],
switch=["generate_xcpqc_files"],
inputs=[
(
"subject",
"scan",
"bold",
"desc-preproc_bold",
"space-T1w_sbref",
"space-T1w_desc-brain_mask",
"max-displacement",
"space-template_desc-preproc_bold",
"space-bold_desc-brain_mask",
["T1w-brain-template-mask", "EPI-template-mask"],
["space-template_desc-bold_mask", "space-EPItemplate_desc-bold_mask"],
"regressors",
["T1w-brain-template-funcreg", "EPI-brain-template-funcreg"],
[
"desc-movementParametersUnfiltered_motion",
"desc-movementParameters_motion",
],
"dvars",
"framewise-displacement-jenkinson",
)
],
outputs={
"space-template_desc-xcp_quality": {"Template": "T1w-brain-template-mask"}
},
)
def qc_xcp(wf, cfg, strat_pool, pipe_num, opt=None):
"""Generate an XCP-style quality-control file."""
# pylint: disable=invalid-name, unused-argument
if cfg[
"nuisance_corrections", "2-nuisance_regression", "run"
] and not strat_pool.check_rpool("regressors"):
return wf, {}
bids_info = pe.Node(
Function(
input_names=["subject", "scan", "wf_name"],
output_names=["subject", "session", "task", "run"],
imports=["from bids.layout import parse_file_entities"],
function=get_bids_info,
as_module=True,
),
name=f"bids_info_{pipe_num}",
)
bids_info.inputs.wf_name = wf.name
qc_file = pe.Node(
Function(
input_names=[
"sub",
"ses",
"task",
"run",
"desc",
"bold2t1w_mask",
"t1w_mask",
"bold2template_mask",
"template_mask",
"original_func",
"final_func",
"template",
"movement_parameters",
"dvars",
"censor_indices",
"regressors",
"framewise_displacement_jenkinson",
"dvars_after",
],
output_names=["qc_file"],
function=generate_xcp_qc,
as_module=True,
),
name=f"qcxcp_{pipe_num}",
)
qc_file.inputs.desc = "preproc"
qc_file.inputs.regressors = (
strat_pool.node_data("regressors")
.node.name.split("regressors_")[-1][::-1]
.split("_", 1)[-1][::-1]
)
bold_to_T1w_mask = pe.Node(
interface=fsl.ImageMaths(),
name=f"binarize_bold_to_T1w_mask_{pipe_num}",
op_string="-bin ",
)
nodes = {
key: strat_pool.node_data(key)
for key in [
"bold",
"desc-preproc_bold",
"max-displacement",
"scan",
"space-bold_desc-brain_mask",
"space-T1w_desc-brain_mask",
"space-T1w_sbref",
"space-template_desc-preproc_bold",
"subject",
*motion_params,
]
if strat_pool.check_rpool(key)
}
nodes["bold2template_mask"] = strat_pool.node_data(
["space-template_desc-bold_mask", "space-EPItemplate_desc-bold_mask"]
)
nodes["template_mask"] = strat_pool.node_data(
["T1w-brain-template-mask", "EPI-template-mask"]
)
nodes["template"] = strat_pool.node_data(
["T1w-brain-template-funcreg", "EPI-brain-template-funcreg"]
)
resample_bold_mask_to_template = pe.Node(
afni.Resample(),
name=f"resample_bold_mask_to_anat_res_{pipe_num}",
mem_gb=0,
mem_x=(0.0115, "in_file", "t"),
)
resample_bold_mask_to_template.inputs.outputtype = "NIFTI_GZ"
wf = _connect_motion(wf, nodes, strat_pool, qc_file, pipe_num=pipe_num)
wf.connect(
[
(nodes["subject"].node, bids_info, [(nodes["subject"].out, "subject")]),
(nodes["scan"].node, bids_info, [(nodes["scan"].out, "scan")]),
(
nodes["space-T1w_sbref"].node,
bold_to_T1w_mask,
[(nodes["space-T1w_sbref"].out, "in_file")],
),
(
nodes["space-T1w_desc-brain_mask"].node,
qc_file,
[(nodes["space-T1w_desc-brain_mask"].out, "t1w_mask")],
),
(bold_to_T1w_mask, qc_file, [("out_file", "bold2t1w_mask")]),
(
nodes["template_mask"].node,
qc_file,
[(nodes["template_mask"].out, "template_mask")],
),
(nodes["bold"].node, qc_file, [(nodes["bold"].out, "original_func")]),
(
nodes["space-template_desc-preproc_bold"].node,
qc_file,
[(nodes["space-template_desc-preproc_bold"].out, "final_func")],
),
(nodes["template"].node, qc_file, [(nodes["template"].out, "template")]),
(
nodes["template_mask"].node,
resample_bold_mask_to_template,
[(nodes["template_mask"].out, "master")],
),
(
nodes["bold2template_mask"].node,
resample_bold_mask_to_template,
[(nodes["bold2template_mask"].out, "in_file")],
),
(
resample_bold_mask_to_template,
qc_file,
[("out_file", "bold2template_mask")],
),
(
bids_info,
qc_file,
[
("subject", "sub"),
("session", "ses"),
("task", "task"),
("run", "run"),
],
),
]
)
return wf, {"space-template_desc-xcp_quality": (qc_file, "qc_file")}