Merge pull request #8 from PennLINC/enh/ukb

Fix UKB ingression

ingress2qsirecon/cli/main.py

@@ -56,7 +56,9 @@ def _ingress2qsirecon(**kwargs):
     layouts = create_layout(input_dir, output_dir, input_pipeline, participant_label)
 
     # Create and run overall workflow, which will be broken down to single subject workflows
-    ingress2qsirecon_wf = create_ingress2qsirecon_wf(layouts, base_dir=work_dir, skip_mni2009c_norm=skip_mni2009c_norm)
+    ingress2qsirecon_wf = create_ingress2qsirecon_wf(
+        layouts, input_pipeline, base_dir=work_dir, skip_mni2009c_norm=skip_mni2009c_norm
+    )
     ingress2qsirecon_wf.run()
 
 

ingress2qsirecon/utils/functions.py

@@ -37,7 +37,7 @@
             "bvals": ["DTI", "dMRI", "dMRI", "bvals"],
             "bvecs": ["DTI", "dMRI", "dMRI", "bvecs"],
             "dwi": ["DTI", "dMRI", "dMRI", "data_ud.nii.gz"],
-            "dwiref": ["DTI", "dMRI", "dMRI", "dti_FA.nii.gz"],
+            # "dwiref": ["DTI", "dMRI", "dMRI", "dti_FA.nii.gz"],
             "t1w_brain": ["T1", "T1_unbiased_brain.nii.gz"],
             "brain_mask": ["T1", "T1_brain_mask.nii.gz"],
             # TODO: Add UKB XFM path
@@ -97,22 +97,47 @@ def make_bids_file_paths(subject_layout: dict) -> dict:
     bids_file_paths : :obj:`dict`
         A dictionary of BIDS-ified file paths.
     """
+    import nibabel as nb
+    import numpy as np
+
     bids_base = str(subject_layout["bids_base"])
     subject = str(subject_layout["subject"])
     session = subject_layout["session"]
     if session == None:
         sub_session_string = f"sub-{subject}"
+        bids_base_session = bids_base
     else:
         sub_session_string = f"sub-{subject}_ses-{session}"
+        bids_base_session = os.path.join(bids_base, f"ses-{session}")
     mni_template = str(subject_layout["MNI_template"])
 
+    # Check for DWI obliquity
+    dwi_img = nb.load(subject_layout["dwi"])
+    dwi_obliquity = np.any(nb.affines.obliquity(dwi_img.affine) > 1e-4)
+    if dwi_obliquity:
+        dwi_oblique_string = "_acq-VARIANTOblique"
+    else:
+        dwi_oblique_string = ""
+
     # BIDS-ify required files
-    bids_dwi_file = os.path.join(bids_base, "dwi", sub_session_string + "_space-T1w_desc-preproc_dwi.nii.gz")
-    bids_bval_file = os.path.join(bids_base, "dwi", sub_session_string + "_space-T1w_desc-preproc_dwi.bval")
-    bids_bvec_file = os.path.join(bids_base, "dwi", sub_session_string + "_space-T1w_desc-preproc_dwi.bvec")
-    bids_b_file = os.path.join(bids_base, "dwi", sub_session_string + "_space-T1w_desc-preproc_dwi.b")
-    bids_bmtxt_file = os.path.join(bids_base, "dwi", sub_session_string + "_space-T1w_desc-preproc_dwi.bmtxt")
-    bids_dwiref_file = os.path.join(bids_base, "dwi", sub_session_string + "_space-T1w_dwiref.nii.gz")
+    bids_dwi_file = os.path.join(
+        bids_base_session, "dwi", sub_session_string + dwi_oblique_string + "_space-T1w_desc-preproc_dwi.nii.gz"
+    )
+    bids_bval_file = os.path.join(
+        bids_base_session, "dwi", sub_session_string + dwi_oblique_string + "_space-T1w_desc-preproc_dwi.bval"
+    )
+    bids_bvec_file = os.path.join(
+        bids_base_session, "dwi", sub_session_string + dwi_oblique_string + "_space-T1w_desc-preproc_dwi.bvec"
+    )
+    bids_b_file = os.path.join(
+        bids_base_session, "dwi", sub_session_string + dwi_oblique_string + "_space-T1w_desc-preproc_dwi.b"
+    )
+    bids_bmtxt_file = os.path.join(
+        bids_base_session, "dwi", sub_session_string + dwi_oblique_string + "_space-T1w_desc-preproc_dwi.bmtxt"
+    )
+    bids_dwiref_file = os.path.join(
+        bids_base_session, "dwi", sub_session_string + dwi_oblique_string + "_space-T1w_dwiref.nii.gz"
+    )
 
     bids_file_paths = {
         "bids_dwi": Path(bids_dwi_file),
@@ -124,20 +149,42 @@ def make_bids_file_paths(subject_layout: dict) -> dict:
     }
 
     # Now for optional files
-    if subject_layout['t1w_brain']:
-        bids_t1w_brain = os.path.join(bids_base, "anat", sub_session_string + "_desc-preproc_T1w.nii.gz")
+    if 't1w_brain' in subject_layout:
+        # Check for T1w obliquity
+        t1_img = nb.load(subject_layout["t1w_brain"])
+        t1_obliquity = np.any(nb.affines.obliquity(t1_img.affine) > 1e-4)
+        if t1_obliquity:
+            t1_oblique_string = "_acq-VARIANTOblique"
+        else:
+            t1_oblique_string = ""
+
+        bids_t1w_brain = os.path.join(
+            bids_base, "anat", f"sub-{subject}" + t1_oblique_string + "_desc-preproc_T1w.nii.gz"
+        )
         bids_file_paths.update({"bids_t1w_brain": Path(bids_t1w_brain)})
-    if subject_layout['brain_mask']:
-        bids_brain_mask = os.path.join(bids_base, "anat", sub_session_string + "_desc-brain_mask.nii.gz")
+    if "brain_mask" in subject_layout:
+        bids_brain_mask = os.path.join(
+            bids_base, "anat", f"sub-{subject}" + t1_oblique_string + "_desc-brain_mask.nii.gz"
+        )
         bids_file_paths.update({"bids_brain_mask": Path(bids_brain_mask)})
-    if subject_layout['subject2MNI']:
+    if "subject2MNI" in subject_layout:
         bids_subject2MNI = os.path.join(
-            bids_base, "anat", sub_session_string + f"_from-T1w_to-{mni_template}_mode-image_xfm.h5"
+            bids_base, "anat", f"sub-{subject}" + f"_from-T1w_to-{mni_template}_mode-image_xfm.h5"
         )
         bids_file_paths.update({"bids_subject2MNI": Path(bids_subject2MNI)})
-    if subject_layout['MNI2subject']:
+    else:
+        bids_subject2MNI = os.path.join(
+            bids_base, "anat", f"sub-{subject}" + f"_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5"
+        )
+        bids_file_paths.update({"bids_subject2MNI": Path(bids_subject2MNI)})
+    if "MNI2subject" in subject_layout:
+        bids_MNI2subject = os.path.join(
+            bids_base, "anat", f"sub-{subject}" + f"_from-{mni_template}_to-T1w_mode-image_xfm.h5"
+        )
+        bids_file_paths.update({"bids_MNI2subject": Path(bids_MNI2subject)})
+    else:
         bids_MNI2subject = os.path.join(
-            bids_base, "anat", sub_session_string + f"_from-{mni_template}_to-T1w_mode-image_xfm.h5"
+            bids_base, "anat", f"sub-{subject}" + f"_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5"
         )
         bids_file_paths.update({"bids_MNI2subject": Path(bids_MNI2subject)})
 
@@ -186,8 +233,8 @@ def create_layout(input_dir: Path, output_dir: Path, input_pipeline: str, partic
 
         # Make BIDS base organization
         bids_base = output_dir / f"sub-{subject}"
-        if renamed_ses:
-            bids_base = bids_base / f"ses-{renamed_ses}"
+        # if renamed_ses:
+        #    bids_base = bids_base / f"ses-{renamed_ses}"
 
         file_paths = get_file_paths(potential_dir, input_pipeline)
         # check if any required files do not exist
@@ -204,6 +251,7 @@ def create_layout(input_dir: Path, output_dir: Path, input_pipeline: str, partic
             continue
 
         subject_layout = {
+            "original_name": potential_dir.name,
             "subject": subject,
             "session": renamed_ses,
             "path": Path(potential_dir),
@@ -221,7 +269,7 @@ def create_layout(input_dir: Path, output_dir: Path, input_pipeline: str, partic
     layout = sorted(layout, key=lambda x: x["subject"])
 
     # Raise warnings for requested subjects not in layout
-    missing_particpants = sorted(set(participant_label) - set([subject["subject"] for subject in layout]))
+    missing_particpants = sorted(set(participant_label) - set([subject["original_name"] for subject in layout]))
     if missing_particpants:
         warn(
             f"Requested participant(s) {missing_particpants} not found in layout, please confirm their data exists and are properly organized."

ingress2qsirecon/utils/interfaces.py

@@ -252,6 +252,7 @@ def _run_interface(self, runtime):
 
         else:
             LOGGER.info("Not applying reorientation to %s: already in %s", dwi_in_file, orientation)
+            input_img.to_filename(dwi_out_file)
             self._results["dwi_out_file"] = dwi_out_file
             # Copy and rename bvecs
             if not os.path.exists(bvec_out_file):

ingress2qsirecon/utils/workflows.py

@@ -3,6 +3,7 @@
 """
 
 import os
+import shutil
 from pathlib import Path
 
 from nipype.pipeline.engine import Workflow
@@ -26,7 +27,7 @@ def parse_layout(subject_layout):
     return tuple(subject_layout.values())
 
 
-def create_single_subject_wf(subject_layout, skip_mni2009c_norm=False):
+def create_single_subject_wf(subject_layout, input_pipeline, skip_mni2009c_norm=False):
     """
     Create a nipype workflow to ingest a single subject.
 
@@ -57,6 +58,9 @@ def create_single_subject_wf(subject_layout, skip_mni2009c_norm=False):
         A nipype workflow that operates on a single subject.
     """
     #### WHY DO I HAVE TO REIMPORT THIS STUFF??
+
+    import nibabel as nb
+    import numpy as np
     from nipype import (
         Node,
         Workflow,
@@ -75,8 +79,12 @@ def create_single_subject_wf(subject_layout, skip_mni2009c_norm=False):
     # Make BIDS subject output folder
     bids_base = subject_layout['bids_base']
     if not os.path.exists(bids_base):
-        os.makedirs(Path(bids_base / "anat").resolve())
-        os.makedirs(Path(bids_base / "dwi").resolve())
+        os.makedirs(Path(bids_base / "anat").resolve(), exist_ok=True)
+    session = subject_layout['session']
+    if session == None:
+        os.makedirs(Path(bids_base / "dwi").resolve(), exist_ok=True)
+    else:
+        os.makedirs(Path(bids_base / f"ses-{session}" / "dwi").resolve(), exist_ok=True)
 
     # Create single subject workflow
     wf_name = f"ingress2qsirecon_single_subject_{subject_name}_wf"
@@ -140,12 +148,14 @@ def create_single_subject_wf(subject_layout, skip_mni2009c_norm=False):
             (parse_layout_node, create_bfile_node, [("bids_b", "b_file_out")]),
         ]
     )
+    # if input_pipeline == "ukb":
+    #    conform_dwi_node.inputs.orientation = "LAS"
 
     # Create nodes to conform anatomicals and save to BIDS layout
-    # TMP If false because does not work yet
     if "t1w_brain" in subject_layout.keys():
         template_dimensions_node = Node(TemplateDimensions(), name="template_dimensions")
         conform_t1w_node = Node(Conform(), name="conform_t1w")
+
         wf.connect(
             [
                 (
@@ -194,61 +204,64 @@ def create_single_subject_wf(subject_layout, skip_mni2009c_norm=False):
                 )
             ]
         )
+    else:
+        shutil.copy(subject_layout["dwiref"], subject_layout["bids_dwiref"])
 
     # Convert FNIRT nii warps to ITK nii, then ITK nii to ITK H5
     # Start with subject2MNI
-    if "subject2MNI" in subject_layout.keys():
-        convert_warpfield_node_subject2MNI = Node(ConvertWarpfield(), name="convert_warpfield_subject2MNI")
-        convert_warpfield_node_subject2MNI.inputs.itk_out_xfm = str(subject_layout["bids_subject2MNI"]).replace(
-            ".h5", ".nii.gz"
-        )
-        nii_to_h5_node_subject2MNI = Node(NIFTItoH5(), name="nii_to_h5_subject2MNI")
-        wf.connect(
-            [
-                (
-                    parse_layout_node,
-                    convert_warpfield_node_subject2MNI,
-                    [("subject2MNI", "fnirt_in_xfm"), ("MNI_ref", "fnirt_ref_file")],
-                ),
-                (
-                    convert_warpfield_node_subject2MNI,
-                    nii_to_h5_node_subject2MNI,
-                    [("itk_out_xfm", "xfm_nifti_in")],
-                ),
-                (
-                    parse_layout_node,
-                    nii_to_h5_node_subject2MNI,
-                    [("bids_subject2MNI", "xfm_h5_out")],
-                ),
-            ]
-        )
+    if False: # We're going to skip this because it doesn't work great, remove workbench dependency
+        if "subject2MNI" in subject_layout.keys():
+            convert_warpfield_node_subject2MNI = Node(ConvertWarpfield(), name="convert_warpfield_subject2MNI")
+            convert_warpfield_node_subject2MNI.inputs.itk_out_xfm = str(subject_layout["bids_subject2MNI"]).replace(
+                ".h5", ".nii.gz"
+            )
+            nii_to_h5_node_subject2MNI = Node(NIFTItoH5(), name="nii_to_h5_subject2MNI")
+            wf.connect(
+                [
+                    (
+                        parse_layout_node,
+                        convert_warpfield_node_subject2MNI,
+                        [("subject2MNI", "fnirt_in_xfm"), ("MNI_ref", "fnirt_ref_file")],
+                    ),
+                    (
+                        convert_warpfield_node_subject2MNI,
+                        nii_to_h5_node_subject2MNI,
+                        [("itk_out_xfm", "xfm_nifti_in")],
+                    ),
+                    (
+                        parse_layout_node,
+                        nii_to_h5_node_subject2MNI,
+                        [("bids_subject2MNI", "xfm_h5_out")],
+                    ),
+                ]
+            )
 
-    # Then MNI2Subject
-    if "MNI2subject" in subject_layout.keys():
-        convert_warpfield_node_MNI2subject = Node(ConvertWarpfield(), name="convert_warpfield_MNI2subject")
-        convert_warpfield_node_MNI2subject.inputs.itk_out_xfm = str(subject_layout["bids_MNI2subject"]).replace(
-            ".h5", ".nii.gz"
-        )
-        nii_to_h5_node_MNI2subject = Node(NIFTItoH5(), name="nii_to_h5_MNI2subject")
-        wf.connect(
-            [
-                (
-                    parse_layout_node,
-                    convert_warpfield_node_MNI2subject,
-                    [("MNI2subject", "fnirt_in_xfm"), ("MNI_ref", "fnirt_ref_file")],
-                ),
-                (
-                    convert_warpfield_node_MNI2subject,
-                    nii_to_h5_node_MNI2subject,
-                    [("itk_out_xfm", "xfm_nifti_in")],
-                ),
-                (
-                    parse_layout_node,
-                    nii_to_h5_node_MNI2subject,
-                    [("bids_MNI2subject", "xfm_h5_out")],
-                ),
-            ]
-        )
+        # Then MNI2Subject
+        if "MNI2subject" in subject_layout.keys():
+            convert_warpfield_node_MNI2subject = Node(ConvertWarpfield(), name="convert_warpfield_MNI2subject")
+            convert_warpfield_node_MNI2subject.inputs.itk_out_xfm = str(subject_layout["bids_MNI2subject"]).replace(
+                ".h5", ".nii.gz"
+            )
+            nii_to_h5_node_MNI2subject = Node(NIFTItoH5(), name="nii_to_h5_MNI2subject")
+            wf.connect(
+                [
+                    (
+                        parse_layout_node,
+                        convert_warpfield_node_MNI2subject,
+                        [("MNI2subject", "fnirt_in_xfm"), ("MNI_ref", "fnirt_ref_file")],
+                    ),
+                    (
+                        convert_warpfield_node_MNI2subject,
+                        nii_to_h5_node_MNI2subject,
+                        [("itk_out_xfm", "xfm_nifti_in")],
+                    ),
+                    (
+                        parse_layout_node,
+                        nii_to_h5_node_MNI2subject,
+                        [("bids_MNI2subject", "xfm_h5_out")],
+                    ),
+                ]
+            )
 
     # Now get transform to MNI2009cAsym
     MNI_template = subject_layout["MNI_template"]
@@ -333,7 +346,9 @@ def save_xfm_outputs(
     return wf
 
 
-def create_ingress2qsirecon_wf(layouts, name="ingress2qsirecon_wf", base_dir=os.getcwd(), skip_mni2009c_norm=False):
+def create_ingress2qsirecon_wf(
+    layouts, input_pipeline, name="ingress2qsirecon_wf", base_dir=os.getcwd(), skip_mni2009c_norm=False
+):
     """
     Creates the overall ingress2qsirecon workflow.
 
@@ -342,6 +357,9 @@ def create_ingress2qsirecon_wf(layouts, name="ingress2qsirecon_wf", base_dir=os.
     layouts : list of dict
         A list of dictionaries, one per subject, from the create_layout function.
 
+    input_pipeline : str
+        The name of the input pipeline (e.g. 'hcpya', 'ukb')
+
     name : str, optional
         The name of the workflow. Default is "ingress2qsirecon_wf".
 
@@ -361,7 +379,9 @@ def create_ingress2qsirecon_wf(layouts, name="ingress2qsirecon_wf", base_dir=os.
     print(f"Subject(s) to run: {subjects_to_run}")
 
     for subject_layout in layouts:
-        single_subject_wf = create_single_subject_wf(subject_layout, skip_mni2009c_norm=skip_mni2009c_norm)
+        single_subject_wf = create_single_subject_wf(
+            subject_layout, input_pipeline, skip_mni2009c_norm=skip_mni2009c_norm
+        )
         wf.add_nodes([single_subject_wf])
 
     return wf

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "Ingress2QSIRecon"
-version = "0.2.1"
+version = "0.2.2"
 description = "Tool to ingress data from other pipelines for use in QSIRecon"
 authors = ["Steven Meisler <[email protected]>"]
 readme = "README.md"