From 778c96f6eea253952bb2f6420aac833914613a25 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Mon, 25 Sep 2023 15:10:59 -0400 Subject: [PATCH] Load T1w-to-standard transform to same space as volumetric BOLD scan (#926) --- .../data/test_ds001419_cifti_outputs.txt | 4 +- xcp_d/tests/test_utils_bids.py | 11 +- xcp_d/tests/test_utils_utils.py | 18 +- xcp_d/utils/bids.py | 56 +++- xcp_d/workflows/anatomical.py | 8 + xcp_d/workflows/base.py | 6 +- xcp_d/workflows/execsummary.py | 291 +++++------------- 7 files changed, 154 insertions(+), 240 deletions(-) diff --git a/xcp_d/tests/data/test_ds001419_cifti_outputs.txt b/xcp_d/tests/data/test_ds001419_cifti_outputs.txt index a9dd8b71d..e4bd18a37 100644 --- a/xcp_d/tests/data/test_ds001419_cifti_outputs.txt +++ b/xcp_d/tests/data/test_ds001419_cifti_outputs.txt @@ -47,8 +47,8 @@ xcp_d/space-fsLR_atlas-HCP_den-91k_dseg.dlabel.nii xcp_d/sub-01 xcp_d/sub-01.html xcp_d/sub-01/anat -xcp_d/sub-01/anat/sub-01_space-MNI152NLin6Asym_desc-preproc_T1w.nii.gz -xcp_d/sub-01/anat/sub-01_space-MNI152NLin6Asym_dseg.nii.gz +xcp_d/sub-01/anat/sub-01_space-MNI152NLin2009cAsym_desc-preproc_T1w.nii.gz +xcp_d/sub-01/anat/sub-01_space-MNI152NLin2009cAsym_dseg.nii.gz xcp_d/sub-01/func xcp_d/sub-01/func/sub-01_task-imagery_desc-dcan_qc.hdf5 xcp_d/sub-01/func/sub-01_task-imagery_desc-filtered_motion.tsv diff --git a/xcp_d/tests/test_utils_bids.py b/xcp_d/tests/test_utils_bids.py index 1c31851e1..ff99b55ca 100644 --- a/xcp_d/tests/test_utils_bids.py +++ b/xcp_d/tests/test_utils_bids.py @@ -14,6 +14,11 @@ def test_collect_participants(datasets): This also covers BIDSError and BIDSWarning. """ bids_dir = datasets["ds001419"] + + # Pass in non-BIDS folder to get BIDSError. + with pytest.raises(xbids.BIDSError, match="Could not find participants"): + xbids.collect_participants(os.path.dirname(bids_dir), participant_label="fail") + with pytest.raises(xbids.BIDSError, match="Could not find participants"): xbids.collect_participants(bids_dir, participant_label="fail") @@ -30,7 +35,7 @@ def test_collect_participants(datasets): assert found_labels == ["01"] -def test_collect_data_pnc(datasets): +def test_collect_data_ds001419(datasets): """Test the collect_data function.""" bids_dir = datasets["ds001419"] @@ -69,8 +74,8 @@ def test_collect_data_pnc(datasets): assert "space-fsLR" in subj_data["bold"][0] assert "space-" not in subj_data["t1w"] assert os.path.basename(subj_data["t1w"]) == "sub-01_desc-preproc_T1w.nii.gz" - assert "to-MNI152NLin6Asym" in subj_data["anat_to_template_xfm"] - assert "from-MNI152NLin6Asym" in subj_data["template_to_anat_xfm"] + assert "to-MNI152NLin2009cAsym" in subj_data["anat_to_template_xfm"] + assert "from-MNI152NLin2009cAsym" in subj_data["template_to_anat_xfm"] def test_collect_data_nibabies(datasets): diff --git a/xcp_d/tests/test_utils_utils.py b/xcp_d/tests/test_utils_utils.py index ddac086d2..c2064a8ef 100644 --- a/xcp_d/tests/test_utils_utils.py +++ b/xcp_d/tests/test_utils_utils.py @@ -155,18 +155,18 @@ def test_denoise_with_nilearn(ds001419_data, tmp_path_factory): def test_list_to_str(): """Test the list_to_str function.""" - lst = ["a"] - string = utils.list_to_str(lst) + string = utils.list_to_str(["a"]) assert string == "a" - lst = ["a", "b"] - string = utils.list_to_str(lst) + string = utils.list_to_str(["a", "b"]) assert string == "a and b" - lst = ["a", "b", "c"] - string = utils.list_to_str(lst) + string = utils.list_to_str(["a", "b", "c"]) assert string == "a, b, and c" + with pytest.raises(ValueError, match="Zero-length list provided."): + utils.list_to_str([]) + def test_get_bold2std_and_t1w_xfms(ds001419_data): """Test get_bold2std_and_t1w_xfms.""" @@ -331,7 +331,7 @@ def test_get_std2bold_xfms(ds001419_data): """ bold_file_nlin2009c = ds001419_data["nifti_file"] - # MNI152NLin2009cAsym --> MNI152NLin6Asym + # MNI152NLin6Asym --> MNI152NLin2009cAsym xforms_to_mni = utils.get_std2bold_xfms(bold_file_nlin2009c) assert len(xforms_to_mni) == 1 @@ -343,7 +343,7 @@ def test_get_std2bold_xfms(ds001419_data): xforms_to_mni = utils.get_std2bold_xfms(bold_file_nlin6asym) assert len(xforms_to_mni) == 1 - # MNIInfant --> MNI152NLin6Asym + # MNI152NLin6Asym --> MNIInfant bold_file_infant = bold_file_nlin2009c.replace( "space-MNI152NLin2009cAsym_", "space-MNIInfant_cohort-1_", @@ -351,7 +351,7 @@ def test_get_std2bold_xfms(ds001419_data): xforms_to_mni = utils.get_std2bold_xfms(bold_file_infant) assert len(xforms_to_mni) == 2 - # tofail --> MNI152NLin6Asym + # MNI152NLin6Asym --> tofail bold_file_tofail = bold_file_nlin2009c.replace("space-MNI152NLin2009cAsym_", "space-tofail_") with pytest.raises(ValueError, match="Space 'tofail'"): utils.get_std2bold_xfms(bold_file_tofail) diff --git a/xcp_d/utils/bids.py b/xcp_d/utils/bids.py index 5e38fe999..1321ab536 100644 --- a/xcp_d/utils/bids.py +++ b/xcp_d/utils/bids.py @@ -32,8 +32,8 @@ "nibabies": { "cifti": ["fsLR"], "nifti": [ - "MNIInfant", "MNI152NLin6Asym", + "MNIInfant", "MNI152NLin2009cAsym", ], }, @@ -279,18 +279,34 @@ def collect_data( if cifti: # Select the appropriate volumetric space for the CIFTI template. # This space will be used in the executive summary and T1w/T2w workflows. - temp_query = queries["anat_to_template_xfm"].copy() - volumetric_space = ASSOCIATED_TEMPLATES[space] + allowed_spaces = INPUT_TYPE_ALLOWED_SPACES.get( + input_type, + DEFAULT_ALLOWED_SPACES, + )["nifti"] + + temp_bold_query = queries["bold"].copy() + temp_bold_query.pop("den", None) + temp_bold_query["extension"] = ".nii.gz" + + temp_xfm_query = queries["anat_to_template_xfm"].copy() - temp_query["to"] = volumetric_space - transform_files = layout.get(**temp_query) - if not transform_files: + for volspace in allowed_spaces: + temp_bold_query["space"] = volspace + bold_data = layout.get(**temp_bold_query) + temp_xfm_query["to"] = volspace + transform_files = layout.get(**temp_xfm_query) + + if bold_data and transform_files: + # will leave the best available space in the query + break + + if not bold_data or not transform_files: raise FileNotFoundError( - f"No nifti transforms found to allowed space ({volumetric_space})" + f"No BOLD NIfTI or transforms found to allowed space ({volspace})" ) - queries["anat_to_template_xfm"]["to"] = volumetric_space - queries["template_to_anat_xfm"]["from"] = volumetric_space + queries["anat_to_template_xfm"]["to"] = volspace + queries["template_to_anat_xfm"]["from"] = volspace else: # use the BOLD file's space if the BOLD file is a nifti. queries["anat_to_template_xfm"]["to"] = queries["bold"]["space"] @@ -542,19 +558,20 @@ def collect_morphometry_data(layout, participant_label): @fill_doc -def collect_run_data(layout, input_type, bold_file, cifti, primary_anat): +def collect_run_data(layout, bold_file, cifti, primary_anat, target_space): """Collect data associated with a given BOLD file. Parameters ---------- %(layout)s - %(input_type)s bold_file : :obj:`str` Path to the BOLD file. %(cifti)s Whether to collect files associated with a CIFTI image (True) or a NIFTI (False). primary_anat : {"T1w", "T2w"} The anatomical modality to use for the anat-to-native transform. + target_space + Used to find NIfTIs in the appropriate space if ``cifti`` is ``True``. Returns ------- @@ -598,24 +615,29 @@ def collect_run_data(layout, input_type, bold_file, cifti, primary_anat): suffix="xfm", ) else: - allowed_nifti_spaces = INPUT_TYPE_ALLOWED_SPACES.get( - input_type, - DEFAULT_ALLOWED_SPACES, - )["nifti"] + # Split cohort out of the space for MNIInfant templates. + cohort = None + if "+" in target_space: + target_space, cohort = target_space.split("+") + run_data["boldref"] = layout.get_nearest( bids_file.path, strict=False, - space=allowed_nifti_spaces, + space=target_space, + cohort=cohort, suffix="boldref", extension=[".nii", ".nii.gz"], + invalid_filters="allow", ) run_data["nifti_file"] = layout.get_nearest( bids_file.path, strict=False, - space=allowed_nifti_spaces, + space=target_space, + cohort=cohort, desc="preproc", suffix="bold", extension=[".nii", ".nii.gz"], + invalid_filters="allow", ) LOGGER.log( diff --git a/xcp_d/workflows/anatomical.py b/xcp_d/workflows/anatomical.py index e91feac5e..874b6e54b 100644 --- a/xcp_d/workflows/anatomical.py +++ b/xcp_d/workflows/anatomical.py @@ -30,6 +30,7 @@ ) from xcp_d.utils.bids import get_freesurfer_dir, get_freesurfer_sphere from xcp_d.utils.doc import fill_doc +from xcp_d.utils.utils import list_to_str from xcp_d.workflows.execsummary import ( init_brainsprite_figures_wf, init_execsummary_anatomical_plots_wf, @@ -204,6 +205,13 @@ def init_postprocess_anat_wf( workflow.connect([(inputnode, ds_t2w_std, [("t2w", "in_file")])]) else: + out = ( + ["T1w"] if t1w_available else [] + ["T2w"] if t2w_available else [] + ["segmentation"] + ) + workflow.__desc__ = f"""\ +Native-space {list_to_str(out)} images were transformed to {target_space} space at 1 mm3 +resolution. +""" warp_anat_dseg_to_template = pe.Node( ApplyTransforms( num_threads=2, diff --git a/xcp_d/workflows/base.py b/xcp_d/workflows/base.py index 7573d57ca..7620c3c94 100644 --- a/xcp_d/workflows/base.py +++ b/xcp_d/workflows/base.py @@ -669,11 +669,11 @@ def init_subject_wf( for j_run, bold_file in enumerate(task_files): run_data = collect_run_data( - layout, - input_type, - bold_file, + layout=layout, + bold_file=bold_file, cifti=cifti, primary_anat=primary_anat, + target_space=target_space, ) post_scrubbing_duration = flag_bad_run( diff --git a/xcp_d/workflows/execsummary.py b/xcp_d/workflows/execsummary.py index f98b66aab..ed053db67 100644 --- a/xcp_d/workflows/execsummary.py +++ b/xcp_d/workflows/execsummary.py @@ -268,9 +268,7 @@ def init_brainsprite_figures_wf( # fmt:off workflow.connect([ - (modify_pngs_template_scene, create_scenewise_pngs, [ - ("out_file", "scene_file"), - ]), + (modify_pngs_template_scene, create_scenewise_pngs, [("out_file", "scene_file")]), (get_png_scene_names, create_scenewise_pngs, [ ("scene_index", "scene_name_or_number"), ]), @@ -348,7 +346,9 @@ def init_execsummary_functional_plots_wf( Set from the parameter. %(boldref)s t1w + T1w image in a standard space, taken from the output of init_postprocess_anat_wf. t2w + T2w image in a standard space, taken from the output of init_postprocess_anat_wf. """ workflow = Workflow(name=name) @@ -356,68 +356,63 @@ def init_execsummary_functional_plots_wf( niu.IdentityInterface( fields=[ "preproc_nifti", - "boldref", + "boldref", # a nifti boldref "t1w", "t2w", # optional - ] - ), # a nifti boldref + ], + ), name="inputnode", ) - if preproc_nifti: - inputnode.inputs.preproc_nifti = preproc_nifti - - # Only grab the bb_registration_file if the preprocessed BOLD file is a parameter. - # Get bb_registration_file prefix from fmriprep - # TODO: Replace with interfaces. - current_bold_file = os.path.basename(preproc_nifti) - if "_space" in current_bold_file: - bb_register_prefix = current_bold_file.split("_space")[0] - else: - bb_register_prefix = current_bold_file.split("_desc")[0] - - bold_t1w_registration_files = layout.get( - desc=["bbregister", "coreg", "bbr"], - extension=".svg", - suffix="bold", - return_type="file", + if not preproc_nifti: + raise ValueError( + "No preprocessed NIfTI found. Executive summary figures cannot be generated." ) - bold_t1w_registration_file = fnmatch.filter( - bold_t1w_registration_files, - f"*/{bb_register_prefix}*", - )[0] - ds_registration_figure = pe.Node( - DerivativesDataSink( - base_directory=output_dir, - in_file=bold_t1w_registration_file, - dismiss_entities=["den"], - datatype="figures", - desc="bbregister", - ), - name="ds_registration_figure", - run_without_submitting=True, - ) + inputnode.inputs.preproc_nifti = preproc_nifti - workflow.connect([(inputnode, ds_registration_figure, [("preproc_nifti", "source_file")])]) + # Get bb_registration_file prefix from fmriprep + # TODO: Replace with interfaces. + current_bold_file = os.path.basename(preproc_nifti) + if "_space" in current_bold_file: + bb_register_prefix = current_bold_file.split("_space")[0] else: - LOGGER.warning( - "Preprocessed NIFTI file not provided as a parameter, " - "so the BBReg figure will not be extracted." - ) + bb_register_prefix = current_bold_file.split("_desc")[0] - # Plot the mean bold image + bold_t1w_registration_files = layout.get( + desc=["bbregister", "coreg", "bbr"], + extension=".svg", + suffix="bold", + return_type="file", + ) + bold_t1w_registration_file = fnmatch.filter( + bold_t1w_registration_files, + f"*/{bb_register_prefix}*", + )[0] + + ds_registration_figure = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + in_file=bold_t1w_registration_file, + dismiss_entities=["den"], + datatype="figures", + desc="bbregister", + ), + name="ds_registration_figure", + run_without_submitting=True, + ) + + workflow.connect([(inputnode, ds_registration_figure, [("preproc_nifti", "source_file")])]) + + # Calculate the mean bold image calculate_mean_bold = pe.Node( BinaryMath(expression="np.mean(img, axis=3)"), name="calculate_mean_bold", ) - plot_meanbold = pe.Node(AnatomicalPlot(), name="plot_meanbold") + workflow.connect([(inputnode, calculate_mean_bold, [("preproc_nifti", "in_file")])]) - # fmt:off - workflow.connect([ - (inputnode, calculate_mean_bold, [("preproc_nifti", "in_file")]), - (calculate_mean_bold, plot_meanbold, [("out_file", "in_file")]), - ]) - # fmt:on + # Plot the mean bold image + plot_meanbold = pe.Node(AnatomicalPlot(), name="plot_meanbold") + workflow.connect([(calculate_mean_bold, plot_meanbold, [("out_file", "in_file")])]) # Write out the figures. ds_meanbold_figure = pe.Node( @@ -440,7 +435,6 @@ def init_execsummary_functional_plots_wf( # Plot the reference bold image plot_boldref = pe.Node(AnatomicalPlot(), name="plot_boldref") - workflow.connect([(inputnode, plot_boldref, [("boldref", "in_file")])]) # Write out the figures. @@ -464,115 +458,54 @@ def init_execsummary_functional_plots_wf( # Start plotting the overlay figures # T1 in Task, Task in T1, Task in T2, T2 in Task - if t1w_available: - # Resample T1w to match resolution of task data - resample_t1w = pe.Node( + anatomicals = ["t1w"] if t1w_available else [] + ["t2w"] if t2w_available else [] + for anat in anatomicals: + # Resample BOLD to match resolution of T1w/T2w data + resample_bold_to_anat = pe.Node( ResampleToImage(), - name="resample_t1w", + name=f"resample_bold_to_{anat}", ) # fmt:off workflow.connect([ - (inputnode, resample_t1w, [ - ("t1w", "in_file"), - ]), - (calculate_mean_bold, resample_t1w, [ - ("out_file", "target_file"), - ]), + (inputnode, resample_bold_to_anat, [(anat, "target_file")]), + (calculate_mean_bold, resample_bold_to_anat, [("out_file", "in_file")]), ]) # fmt:on - plot_t1w_on_task_wf = init_plot_overlay_wf( + plot_anat_on_task_wf = init_plot_overlay_wf( output_dir=output_dir, - desc="T1wOnTask", - name="plot_t1w_on_task_wf", + desc=f"{anat[0].upper()}{anat[1:]}OnTask", + name=f"plot_{anat}_on_task_wf", ) # fmt:off workflow.connect([ - (inputnode, plot_t1w_on_task_wf, [ + (inputnode, plot_anat_on_task_wf, [ ("preproc_nifti", "inputnode.name_source"), + (anat, "inputnode.overlay_file"), ]), - (calculate_mean_bold, plot_t1w_on_task_wf, [ + (resample_bold_to_anat, plot_anat_on_task_wf, [ ("out_file", "inputnode.underlay_file"), ]), - (resample_t1w, plot_t1w_on_task_wf, [ - ("out_file", "inputnode.overlay_file"), - ]), ]) # fmt:on - plot_task_on_t1w_wf = init_plot_overlay_wf( + plot_task_on_anat_wf = init_plot_overlay_wf( output_dir=output_dir, - desc="TaskOnT1w", - name="plot_task_on_t1w_wf", + desc=f"TaskOn{anat[0].upper()}{anat[1:]}", + name=f"plot_task_on_{anat}_wf", ) # fmt:off workflow.connect([ - (inputnode, plot_task_on_t1w_wf, [ + (inputnode, plot_task_on_anat_wf, [ ("preproc_nifti", "inputnode.name_source"), + (anat, "inputnode.underlay_file"), ]), - (calculate_mean_bold, plot_task_on_t1w_wf, [ + (resample_bold_to_anat, plot_task_on_anat_wf, [ ("out_file", "inputnode.overlay_file"), ]), - (resample_t1w, plot_task_on_t1w_wf, [ - ("out_file", "inputnode.underlay_file"), - ]), - ]) - # fmt:on - - if t2w_available: - # Resample T2w to match resolution of task data - resample_t2w = pe.Node( - ResampleToImage(), - name="resample_t2w", - ) - - # fmt:off - workflow.connect([ - (inputnode, resample_t2w, [ - ("t2w", "in_file"), - ]), - (calculate_mean_bold, resample_t2w, [ - ("out_file", "target_file"), - ]), - ]) - # fmt:on - - plot_t2w_on_task_wf = init_plot_overlay_wf( - output_dir=output_dir, - desc="T2wOnTask", - name="plot_t2w_on_task_wf", - ) - - # fmt:off - workflow.connect([ - (inputnode, plot_t2w_on_task_wf, [ - ("preproc_nifti", "inputnode.underlay_file"), - ("preproc_nifti", "inputnode.name_source"), - ]), - (resample_t2w, plot_t2w_on_task_wf, [ - ("out_file", "inputnode.overlay_file"), - ]), - ]) - # fmt:on - - plot_task_on_t2w_wf = init_plot_overlay_wf( - output_dir=output_dir, - desc="TaskOnT2w", - name="plot_task_on_t2w_wf", - ) - - # fmt:off - workflow.connect([ - (inputnode, plot_task_on_t2w_wf, [ - ("preproc_nifti", "inputnode.overlay_file"), - ("preproc_nifti", "inputnode.name_source"), - ]), - (resample_t2w, plot_task_on_t2w_wf, [ - ("out_file", "inputnode.underlay_file"), - ]), ]) # fmt:on @@ -614,7 +547,9 @@ def init_execsummary_anatomical_plots_wf( Inputs ------ t1w + T1w image, after warping to standard space. t2w + T2w image, after warping to standard space. template """ workflow = Workflow(name=name) @@ -631,109 +566,53 @@ def init_execsummary_anatomical_plots_wf( ) # Start plotting the overlay figures - # Atlas in T1w, T1w in Atlas - if t1w_available: - # Resample T1w to match resolution of template data - resample_t1w = pe.Node( + # Atlas in T1w/T2w, T1w/T2w in Atlas + anatomicals = ["t1w"] if t1w_available else [] + ["t2w"] if t2w_available else [] + for anat in anatomicals: + # Resample anatomical to match resolution of template data + resample_anat = pe.Node( ResampleToImage(), - name="resample_t1w", + name=f"resample_{anat}", ) # fmt:off workflow.connect([ - (inputnode, resample_t1w, [ - ("t1w", "in_file"), + (inputnode, resample_anat, [ + (anat, "in_file"), ("template", "target_file"), - ]) - ]) - # fmt:on - - plot_t1w_on_atlas_wf = init_plot_overlay_wf( - output_dir=output_dir, - desc="AnatOnAtlas", - name="plot_t1w_on_atlas_wf", - ) - - # fmt:off - workflow.connect([ - (inputnode, plot_t1w_on_atlas_wf, [ - ("template", "inputnode.underlay_file"), - ("t1w", "inputnode.name_source"), - ]), - (resample_t1w, plot_t1w_on_atlas_wf, [ - ("out_file", "inputnode.overlay_file"), ]), ]) # fmt:on - plot_atlas_on_t1w_wf = init_plot_overlay_wf( - output_dir=output_dir, - desc="AtlasOnAnat", - name="plot_atlas_on_t1w_wf", - ) - - # fmt:off - workflow.connect([ - (inputnode, plot_atlas_on_t1w_wf, [ - ("template", "inputnode.overlay_file"), - ("t1w", "inputnode.name_source"), - ]), - (resample_t1w, plot_atlas_on_t1w_wf, [ - ("out_file", "inputnode.underlay_file"), - ]), - ]) - # fmt:on - - # Atlas in T2w, T2w in Atlas - if t2w_available: - # Resample T2w to match resolution of template data - resample_t2w = pe.Node( - ResampleToImage(), - name="resample_t2w", - ) - - # fmt:off - workflow.connect([ - (inputnode, resample_t2w, [ - ("t2w", "in_file"), - ("template", "target_file"), - ]) - ]) - # fmt:on - - plot_t2w_on_atlas_wf = init_plot_overlay_wf( + plot_anat_on_atlas_wf = init_plot_overlay_wf( output_dir=output_dir, desc="AnatOnAtlas", - name="plot_t2w_on_atlas_wf", + name=f"plot_{anat}_on_atlas_wf", ) # fmt:off workflow.connect([ - (inputnode, plot_t2w_on_atlas_wf, [ + (inputnode, plot_anat_on_atlas_wf, [ ("template", "inputnode.underlay_file"), - ("t2w", "inputnode.name_source"), - ]), - (resample_t2w, plot_t2w_on_atlas_wf, [ - ("out_file", "inputnode.overlay_file"), + (anat, "inputnode.name_source"), ]), + (resample_anat, plot_anat_on_atlas_wf, [("out_file", "inputnode.overlay_file")]), ]) # fmt:on - plot_atlas_on_t2w_wf = init_plot_overlay_wf( + plot_atlas_on_anat_wf = init_plot_overlay_wf( output_dir=output_dir, desc="AtlasOnAnat", - name="plot_atlas_on_t2w_wf", + name=f"plot_atlas_on_{anat}_wf", ) # fmt:off workflow.connect([ - (inputnode, plot_atlas_on_t2w_wf, [ + (inputnode, plot_atlas_on_anat_wf, [ ("template", "inputnode.overlay_file"), - ("t2w", "inputnode.name_source"), - ]), - (resample_t2w, plot_atlas_on_t2w_wf, [ - ("out_file", "inputnode.underlay_file"), + (anat, "inputnode.name_source"), ]), + (resample_anat, plot_atlas_on_anat_wf, [("out_file", "inputnode.underlay_file")]), ]) # fmt:on