Retain orthogonalized confounds (PennLINC#975)

xcp_d/interfaces/censoring.py

@@ -425,6 +425,7 @@ class _GenerateConfoundsOutputSpec(TraitedSpec):
             "It will also always have the linear trend and a constant column."
         ),
     )
+    confounds_metadata = traits.Dict(desc="Metadata associated with the confounds_file output.")
     motion_file = File(
         exists=True,
         desc="The filtered motion parameters.",
@@ -497,14 +498,65 @@ def _run_interface(self, runtime):
         )
 
         # Load nuisance regressors, but use filtered motion parameters.
-        confounds_df = load_confound_matrix(
+        confounds_df, confounds_metadata = load_confound_matrix(
             params=self.inputs.params,
             img_file=self.inputs.in_file,
             confounds_file=self._results["filtered_confounds_file"],
             confounds_json_file=self.inputs.fmriprep_confounds_json,
             custom_confounds=self.inputs.custom_confounds_file,
         )
 
+        if confounds_df is not None:
+            # Orthogonalize full nuisance regressors w.r.t. any signal regressors
+            signal_columns = [c for c in confounds_df.columns if c.startswith("signal__")]
+            if signal_columns:
+                LOGGER.warning(
+                    "Signal columns detected. "
+                    "Orthogonalizing nuisance columns w.r.t. the following signal columns: "
+                    f"{', '.join(signal_columns)}"
+                )
+                noise_columns = [c for c in confounds_df.columns if not c.startswith("signal__")]
+
+                # Don't orthogonalize the intercept or linear trend regressors
+                untouched_cols = ["linear_trend", "intercept"]
+                cols_to_orth = [c for c in noise_columns if c not in untouched_cols]
+                orth_confounds_df = confounds_df[noise_columns].copy()
+                orth_cols = [f"{c}_orth" for c in cols_to_orth]
+                orth_confounds_df = pd.DataFrame(
+                    index=confounds_df.index,
+                    columns=orth_cols + untouched_cols,
+                )
+                orth_confounds_df.loc[:, untouched_cols] = confounds_df[untouched_cols]
+
+                # Do the orthogonalization
+                signal_regressors = confounds_df[signal_columns].to_numpy()
+                noise_regressors = confounds_df[cols_to_orth].to_numpy()
+                signal_betas = np.linalg.lstsq(signal_regressors, noise_regressors, rcond=None)[0]
+                pred_noise_regressors = np.dot(signal_regressors, signal_betas)
+                orth_noise_regressors = noise_regressors - pred_noise_regressors
+
+                # Replace the old data
+                orth_confounds_df.loc[:, orth_cols] = orth_noise_regressors
+                confounds_df = orth_confounds_df
+
+                for col in cols_to_orth:
+                    desc_str = (
+                        "This regressor is orthogonalized with respect to the 'signal' regressors "
+                        f"({', '.join(signal_columns)}) after dummy scan removal, "
+                        "but prior to any censoring."
+                    )
+
+                    col_metadata = {}
+                    if col in confounds_metadata.keys():
+                        col_metadata = confounds_metadata.pop(col)
+                        if "Description" in col_metadata.keys():
+                            desc_str = f"{col_metadata['Description']} {desc_str}"
+
+                    col_metadata["Description"] = desc_str
+                    confounds_metadata[f"{col}_orth"] = col_metadata
+
+        self._results["confounds_metadata"] = confounds_metadata
+
         # get the output
         self._results["motion_file"] = fname_presuffix(
             self.inputs.fmriprep_confounds_file,

xcp_d/tests/data/test_ds001419_cifti_outputs.txt

@@ -57,6 +57,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-dcan_qc.hdf5
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_design.tsv
+xcp_d/sub-01/func/sub-01_task-imagery_run-01_design.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_outliers.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_space-fsLR_atlas-4S1052Parcels_coverage.tsv
@@ -169,6 +170,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-dcan_qc.hdf5
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_design.tsv
+xcp_d/sub-01/func/sub-01_task-imagery_run-02_design.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_outliers.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_space-fsLR_atlas-4S1052Parcels_coverage.tsv

xcp_d/tests/data/test_ds001419_nifti_outputs.txt

@@ -55,6 +55,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_outliers.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_space-MNI152NLin2009cAsym_atlas-4S1052Parcels_coverage.tsv
@@ -151,6 +152,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_run-01_space-MNI152NLin2009cAsym_res-2_reh
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_outliers.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_space-MNI152NLin2009cAsym_atlas-4S1052Parcels_coverage.tsv
@@ -265,6 +267,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_space-MNI152NLin2009cAsym_res-2_desc-denoi
 xcp_d/sub-01/func/sub-01_task-rest_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-rest_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-rest_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-rest_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-rest_outliers.json
 xcp_d/sub-01/func/sub-01_task-rest_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-rest_space-MNI152NLin2009cAsym_atlas-4S1052Parcels_coverage.tsv

xcp_d/tests/data/test_fmriprep_without_freesurfer_outputs.txt

@@ -55,6 +55,7 @@ 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-mixedgamblestask_run-1_desc-dcan_qc.hdf5
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_motion.json
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_motion.tsv
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_outliers.json
@@ -126,6 +127,7 @@ xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_space-MNI152NLin2009cAsym_r
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_space-MNI152NLin2009cAsym_reho.nii.gz
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_desc-dcan_qc.hdf5
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_motion.json
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_motion.tsv
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_outliers.json

xcp_d/tests/data/test_nibabies_outputs.txt

@@ -53,6 +53,7 @@ xcp_d/sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_space-MNIInfant_cohort-1_dseg.n
 xcp_d/sub-01/ses-1mo/func
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-dcan_qc.hdf5
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-preproc_design.tsv
+xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-preproc_design.json
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_motion.json
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_motion.tsv
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_outliers.json

xcp_d/tests/data/test_pnc_cifti_outputs.txt

@@ -104,6 +104,7 @@ xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleb
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-filtered_motion.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-filtered_motion.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_design.tsv
+xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_design.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_outliers.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_outliers.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_space-fsLR_atlas-4S1052Parcels_coverage.tsv

xcp_d/tests/data/test_pnc_nifti_outputs.txt

@@ -54,6 +54,7 @@ xcp_d/sub-1648798153/ses-PNC1/func
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_desc-filtered_motion.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_desc-filtered_motion.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_desc-preproc_design.tsv
+xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_desc-preproc_design.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_outliers.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_outliers.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_space-MNI152NLin6Asym_atlas-4S1052Parcels_coverage.tsv
@@ -150,6 +151,7 @@ xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_space-
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-filtered_motion.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-filtered_motion.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-preproc_design.tsv
+xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-preproc_design.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_outliers.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_outliers.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_space-MNI152NLin6Asym_atlas-4S1052Parcels_coverage.tsv

xcp_d/tests/test_confounds.py

@@ -44,7 +44,7 @@ def test_custom_confounds(ds001419_data, tmp_path_factory):
     )
     custom_confounds.to_csv(custom_confounds_file, sep="\t", index=False)
 
-    combined_confounds = load_confound_matrix(
+    combined_confounds, _ = load_confound_matrix(
         params="24P",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -56,7 +56,7 @@ def test_custom_confounds(ds001419_data, tmp_path_factory):
     assert "condition01" in combined_confounds.columns
     assert "condition02" in combined_confounds.columns
 
-    custom_confounds = load_confound_matrix(
+    custom_confounds, _ = load_confound_matrix(
         params="custom",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -126,63 +126,63 @@ def test_load_confounds(ds001419_data):
 
     N_VOLUMES = 60
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="24P",
         img_file=bold_file,
         confounds_file=confounds_file,
         confounds_json_file=confounds_json,
     )
     assert confounds_df.shape == (N_VOLUMES, 26)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="27P",
         img_file=bold_file,
         confounds_file=confounds_file,
         confounds_json_file=confounds_json,
     )
     assert confounds_df.shape == (N_VOLUMES, 29)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="36P",
         img_file=bold_file,
         confounds_file=confounds_file,
         confounds_json_file=confounds_json,
     )
     assert confounds_df.shape == (N_VOLUMES, 38)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="acompcor",
         img_file=bold_file,
         confounds_file=confounds_file,
         confounds_json_file=confounds_json,
     )
     assert confounds_df.shape == (N_VOLUMES, 30)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="acompcor_gsr",
         img_file=bold_file,
         confounds_file=confounds_file,
         confounds_json_file=confounds_json,
     )
     assert confounds_df.shape == (N_VOLUMES, 31)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="aroma",
         img_file=bold_file,
         confounds_file=confounds_file,
         confounds_json_file=confounds_json,
     )
     assert confounds_df.shape == (N_VOLUMES, 50)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="aroma_gsr",
         img_file=bold_file,
         confounds_file=confounds_file,
         confounds_json_file=confounds_json,
     )
     assert confounds_df.shape == (N_VOLUMES, 51)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="none",
         img_file=bold_file,
         confounds_file=confounds_file,

xcp_d/tests/test_interfaces_censoring.py

@@ -1,4 +1,5 @@
 """Tests for framewise displacement calculation."""
+import json
 import os
 
 import nibabel as nb
@@ -16,16 +17,28 @@ def test_generate_confounds(ds001419_data, tmp_path_factory):
     confounds_json = ds001419_data["confounds_json"]
 
     df = pd.read_table(confounds_file)
+    with open(confounds_json, "r") as fo:
+        metadata = json.load(fo)
 
     # Replace confounds tsv values with values that should be omitted
     df.loc[1:3, "trans_x"] = [6, 8, 9]
     df.loc[4:6, "trans_y"] = [7, 8, 9]
     df.loc[7:9, "trans_z"] = [12, 8, 9]
 
+    # Modify JSON file
+    metadata["trans_x"] = {"test": "hello"}
+    confounds_json = os.path.join(tmpdir, "edited_confounds.json")
+    with open(confounds_json, "w") as fo:
+        json.dump(metadata, fo)
+
     # Rename with same convention as initial confounds tsv
     confounds_tsv = os.path.join(tmpdir, "edited_confounds.tsv")
     df.to_csv(confounds_tsv, sep="\t", index=False, header=True)
 
+    custom_confounds_file = os.path.join(tmpdir, "custom_confounds.tsv")
+    df2 = pd.DataFrame(columns=["signal__test"], data=np.random.random((df.shape[0], 1)))
+    df2.to_csv(custom_confounds_file, sep="\t", index=False, header=True)
+
     # Run workflow
     interface = censoring.GenerateConfounds(
         in_file=in_file,
@@ -35,7 +48,7 @@ def test_generate_confounds(ds001419_data, tmp_path_factory):
         head_radius=50,
         fmriprep_confounds_file=confounds_tsv,
         fmriprep_confounds_json=confounds_json,
-        custom_confounds_file=None,
+        custom_confounds_file=custom_confounds_file,
         motion_filter_type=None,
         motion_filter_order=4,
         band_stop_min=0,
@@ -47,6 +60,10 @@ def test_generate_confounds(ds001419_data, tmp_path_factory):
     assert os.path.isfile(results.outputs.confounds_file)
     assert os.path.isfile(results.outputs.motion_file)
     assert os.path.isfile(results.outputs.temporal_mask)
+    out_confounds_file = results.outputs.confounds_file
+    out_df = pd.read_table(out_confounds_file)
+    assert out_df.shape[1] == 26  # 24(P) + linear trend + constant
+    assert sum(out_df.columns.str.endswith("_orth")) == 24  # all 24(P), not linear trend/constant
 
 
 def test_random_censor(tmp_path_factory):

xcp_d/tests/test_utils_utils.py

@@ -112,30 +112,6 @@ def test_denoise_with_nilearn(ds001419_data, tmp_path_factory):
     assert uncensored_denoised_bold.shape == (n_volumes, n_voxels)
     assert interpolated_filtered_bold.shape == (n_volumes, n_voxels)
 
-    # Next, do the orthogonalization
-    reduced_confounds_df["signal__test"] = confounds_df["global_signal"]
-
-    # Move intercept to end of dataframe
-    reduced_confounds_df = reduced_confounds_df[
-        [c for c in reduced_confounds_df.columns if c not in ["intercept"]] + ["intercept"]
-    ]
-    orth_confounds_file = os.path.join(tmpdir, "orth_confounds.tsv")
-    reduced_confounds_df.to_csv(orth_confounds_file, sep="\t", index=False)
-    (
-        uncensored_denoised_bold,
-        interpolated_filtered_bold,
-    ) = utils.denoise_with_nilearn(
-        preprocessed_bold=preprocessed_bold_arr,
-        confounds_file=orth_confounds_file,
-        temporal_mask=temporal_mask,
-        low_pass=low_pass,
-        high_pass=high_pass,
-        filter_order=filter_order,
-        TR=TR,
-    )
-    assert uncensored_denoised_bold.shape == (n_volumes, n_voxels)
-    assert interpolated_filtered_bold.shape == (n_volumes, n_voxels)
-
     # Finally, run without denoising
     (
         uncensored_denoised_bold,

xcp_d/utils/confounds.py

@@ -1,6 +1,7 @@
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """Confound matrix selection based on Ciric et al. 2007."""
+import json
 import os
 import warnings
 
@@ -371,6 +372,8 @@ def load_confound_matrix(
         If "AROMA" is requested, then this DataFrame will include signal components as well.
         These will be named something like "signal_[XX]".
         If ``params`` is "none", ``confounds_df`` will be None.
+    confounds_metadata : :obj:`dict`
+        Metadata for the columns in the confounds file.
     """
     PARAM_KWARGS = {
         # Get rot and trans values, as well as derivatives and square
@@ -427,7 +430,7 @@ def load_confound_matrix(
     }
 
     if params == "none":
-        return None
+        return None, {}
 
     if params in PARAM_KWARGS:
         kwargs = PARAM_KWARGS[params]
@@ -462,10 +465,24 @@ def load_confound_matrix(
         custom_confounds_df = pd.read_table(custom_confounds, sep="\t")
         confounds_df = pd.concat([custom_confounds_df, confounds_df], axis=1)
 
+    with open(confounds_json_file, "r") as fo:
+        full_confounds_metadata = json.load(fo)
+
+    confounds_metadata = {
+        k: v for k, v in full_confounds_metadata.items() if k in confounds_df.columns
+    }
+    confounds_metadata["linear_trend"] = {
+        "Description": (
+            "Linear trend regressor, generated after dummy scan removal, "
+            "but before any censoring."
+        ),
+    }
+    confounds_metadata["intercept"] = {"Description": "Intercept regressor."}
+
     confounds_df["linear_trend"] = np.arange(confounds_df.shape[0])
     confounds_df["intercept"] = np.ones(confounds_df.shape[0])
 
-    return confounds_df
+    return confounds_df, confounds_metadata
 
 
 def _get_mixing_matrix(img_file):