Fix bugs

hcpdiffpy/interfaces/data.py

@@ -34,12 +34,13 @@ class InitData(SimpleInterface):
             for phase in self.inputs.config["phases"]:
                 for file_type in [".nii.gz", ".bval", ".bvec"]:
                     key = f"dir{ndir}_{phase}{file_type}"
-                    self._results["d_files"][key] = Path(d_dir, f"{self.inputs.subject}_dMRI_{key}")
+                    self._results["d_files"][key] = Path(
+                        d_dir, f"{self.inputs.config['subject']}_dMRI_{key}")
 
         anat_dir = Path(sub_dir, "T1w")
         self._results["t1_file"] = Path(anat_dir, 'T1w_acpc_dc.nii.gz')
         self._results["t1_restore_file"] = Path(anat_dir, 'T1w_acpc_dc_restore.nii.gz')
-        self._results["t1_restore_brain"] = Path(anat_dir, 'T1w_acpc_dc_restore_brain.nii.gz')
+        self._results["t1_brain_file"] = Path(anat_dir, 'T1w_acpc_dc_restore_brain.nii.gz')
         self._results["bias_file"] = Path(anat_dir, 'BiasField_acpc_dc.nii.gz')
         self._results["mask_file"] = Path(anat_dir, 'brainmask_fs.nii.gz')
 
@@ -63,7 +64,6 @@ class SaveData(SimpleInterface):
 
     def _run_interface(self, runtime):
         out_dir = self.inputs.config["output_dir"]
-        out_dir.mkdir(parents=True, exist_ok=True)
         copyfile(self.inputs.data_file, Path(out_dir, "data.nii.gz"))
         copyfile(self.inputs.bval_file, Path(out_dir, "bvals"))
         copyfile(self.inputs.bvec_file, Path(out_dir, "bvecs"))

hcpdiffpy/interfaces/preproc.py

@@ -44,7 +44,7 @@ class ExtractB0(SimpleInterface):
 
         for b in bvals:
             roi_file = Path(
-                self.inputs.config["work_dir"],
+                self.inputs.config["tmp_dir"],
                 f"roi{vol_count}_{Path(self.inputs.data_file).name}")
             if b < b0maxbval and b0dist is None:
                 roi = fsl.ExtractROI(
@@ -85,7 +85,7 @@ class ExtractB0(SimpleInterface):
 
 class _RescaleInputSpec(BaseInterfaceInputSpec):
     config = traits.Dict(mandatory=True, desc="Workflow configurations")
-    scale_files = traits.List(mandatory=True, dtype=str, desc="filenames of scale files")
+    scale_files = traits.List(mandatory=True, desc="filenames of scale files")
     d_files = traits.Dict(mandatory=True, dtype=Path, desc="filenames of diffusion data")
     fsl_cmd = traits.Any(mandatory=True, desc="FSL command for using singularity image (or not)")
 
@@ -161,6 +161,32 @@ class PrepareTopup(SimpleInterface):
         return runtime
 
 
+class _BETMaskInputSpec(BaseInterfaceInputSpec):
+    in_file = traits.File(mandatory=True, exists=True, desc="Input file to skull strip")
+    fsl_cmd = traits.Any(mandatory=True, desc="FSL command for using singularity image (or not)")
+    config = traits.Dict(mandatory=True, desc="Workflow configurations")
+
+
+class _BETMaskOutputSpec(TraitedSpec):
+    mask_file = traits.File(exists=True, desc="Binary brain mask file")
+
+
+class BETMask(SimpleInterface):
+    """Run FSL BET to get a mask with subprocess"""
+    input_spec = _BETMaskInputSpec
+    output_spec = _BETMaskOutputSpec
+
+    def _run_interface(self, runtime):
+        out_file = Path(self.inputs.config["tmp_dir"], f"{Path(self.inputs.in_file).stem}.nii.gz")
+        self._results["mask_file"] = Path(
+            self.inputs.config["tmp_dir"], f"{Path(self.inputs.in_file).stem}_mask.nii.gz")
+        subprocess.run(
+            self.inputs.fsl_cmd.cmd("bet").split() + [
+                self.inputs.in_file, out_file, "-f", "0.20", "-m"], check=True)
+
+        return runtime
+
+
 class _MergeBFilesInputSpec(BaseInterfaceInputSpec):
     config = traits.Dict(mandatory=True, desc="Workflow configurations")
     bval_files = traits.List(mandatory=True, dtype=Path, desc="list of bval files to merge")
@@ -189,8 +215,8 @@ class MergeBFiles(SimpleInterface):
             bvecs = pd.concat(
                 [bvecs, pd.read_csv(bvec_file[0], delim_whitespace=True, header=None)], axis=1)
 
-        self._results["bval_merged"] = Path(self.inputs.config["work_dir"], "merged.bval")
-        self._results["bvec_merged"] = Path(self.inputs.config["work_dir"], "merged.bvec")
+        self._results["bval_merged"] = Path(self.inputs.config["tmp_dir"], "merged.bval")
+        self._results["bvec_merged"] = Path(self.inputs.config["tmp_dir"], "merged.bvec")
         bvals.to_csv(self._results["bval_merged"], sep='\t', header=False, index=False)
         bvecs.to_csv(self._results["bvec_merged"], sep='\t', header=False, index=False)
 
@@ -248,7 +274,7 @@ class EddyIndex(SimpleInterface):
                         indices.append(neg_count)
             vol_prev = vol_curr
 
-        self._results["index_file"] = Path(self.inputs.config["work_dir"], "index.txt")
+        self._results["index_file"] = Path(self.inputs.config["tmp_dir"], "index.txt")
         pd.DataFrame(indices).to_csv(
             self._results["index_file"], sep="\t", header=False, index=False)
 
@@ -298,11 +324,11 @@ class EddyPostProc(SimpleInterface):
             corrvols.append([dim4, dim4])
             tsizes.append(bval.shape[1])
 
-        bval_merged = Path(self.inputs.config["work_dir"], f"{dirs}.bval")
+        bval_merged = Path(self.inputs.config["tmp_dir"], f"{dirs}.bval")
         bvals.to_csv(bval_merged, sep='\t', header=False, index=False)
-        bvec_merged = Path(self.inputs.config["work_dir"], f"{dirs}.bvec")
+        bvec_merged = Path(self.inputs.config["tmp_dir"], f"{dirs}.bvec")
         bvecs.to_csv(bvec_merged, sep='\t', header=False, index=False)
-        corrvols_file = Path(self.inputs.config["work_dir"], f"{dirs}_volnum.txt")
+        corrvols_file = Path(self.inputs.config["tmp_dir"], f"{dirs}_volnum.txt")
         pd.DataFrame(corrvols).to_csv(corrvols_file, sep='\t', header=False, index=False)
 
         bval_tsize = bvals.shape[1]
@@ -344,8 +370,8 @@ class EddyPostProc(SimpleInterface):
             avg_bvals[i] = np.rint(eigvals[eigvalmax] ** 0.5)
             avg_bvecs[:, i] = eigvecs[:, eigvalmax]
 
-        self._results["rot_bvals"] = Path(self.inputs.config["work_dir"], "rotated.bval")
-        self._results["rot_bvecs"] = Path(self.inputs.config["work_dir"], "rotated.bvec")
+        self._results["rot_bvals"] = Path(self.inputs.config["tmp_dir"], "rotated.bval")
+        self._results["rot_bvecs"] = Path(self.inputs.config["tmp_dir"], "rotated.bvec")
         pd.DataFrame(avg_bvals).T.to_csv(
             self._results["rot_bvals"], sep=' ', header=False, index=False)
         pd.DataFrame(avg_bvecs).to_csv(
@@ -363,9 +389,9 @@ class EddyPostProc(SimpleInterface):
         subprocess.run(
             self.inputs.fsl_cmd.cmd("eddy_combine").split() + [
                 pos_dwi, pos_bval, pos_bvec, pos_corrvols, neg_dwi, neg_bval, neg_bvec,
-                neg_corrvols, self.inputs.config["work_dir"], "1"],
+                neg_corrvols, self.inputs.config["tmp_dir"], "1"],
             check=True)
-        self._results["combined_dwi_file"] = Path(self.inputs.config["work_dir"], "data.nii.gz")
+        self._results["combined_dwi_file"] = Path(self.inputs.config["tmp_dir"], "data.nii.gz")
         self._rotate_b(pos_tsize, neg_tsize, pos_bvals, neg_bvals)
 
         return runtime
@@ -388,7 +414,7 @@ class WBDilate(SimpleInterface):
     output_spec = _WBDilateOutputSpec
 
     def _run_interface(self, runtime):
-        self._results["out_file"] = Path(self.inputs.config["work_dir"], "data_dilated.nii.gz")
+        self._results["out_file"] = Path(self.inputs.config["tmp_dir"], "data_dilated.nii.gz")
         args = (
             f"-volume-dilate {self.inputs.data_file} {self.inputs.dilate} NEAREST "
             f"{self._results['out_file']}")
@@ -423,7 +449,7 @@ class DilateMask(SimpleInterface):
 
         for _ in range(6):
             resamp_curr = fsl.ImageMaths(
-                command=self.inputs.fsl_cmd.run_cmd("fslmaths"), in_file=mask_prev, args=args)
+                command=self.inputs.fsl_cmd.cmd("fslmaths"), in_file=mask_prev, args=args)
             resamp_curr.run()
             mask_prev = resamp_curr.aggregate_outputs().out_file
         self._results["out_file"] = mask_prev
@@ -449,7 +475,7 @@ class RotateBVec2Str(SimpleInterface):
     def _run_interface(self, runtime):
         bvecs = pd.read_csv(self.inputs.bvecs_file, delim_whitespace=True, header=None)
         rotated_bvecs = np.matmul(np.array(self.inputs.rot)[:3, :3], bvecs)
-        self._results["rotated_file"] = Path(self.inputs.config["work_dir"], "rotated2str.bvec")
+        self._results["rotated_file"] = Path(self.inputs.config["tmp_dir"], "rotated2str.bvec")
         pd.DataFrame(rotated_bvecs).to_csv(
             self._results["rotated_file"], sep=' ', header=False, index=False,
             float_format="%10.6f", quoting=3, escapechar=' ')

hcpdiffpy/interfaces/utilities.py

@@ -50,7 +50,7 @@ class UpdateDiffFiles(SimpleInterface):
         for key in self.inputs.config["keys"]:
             dwi_key = [
                 d_key for d_key in self.inputs.d_files if key in d_key and '.nii.gz' in d_key]
-            dwi_replace = [d_file for d_file in self.inputs.dwi_replacements if key in str(d_file)]
+            dwi_replace = [d_file for d_file in self.inputs.data_files if key in str(d_file)]
             self._results["d_files"][dwi_key[0]] = dwi_replace[0]
         return runtime
 
@@ -85,7 +85,7 @@ class _DiffResInputSpec(BaseInterfaceInputSpec):
 
 
 class _DiffResOutputSpec(TraitedSpec):
-    res = traits.Int(desc="Diffusion data resolution (assuming isotropic)")
+    res = traits.Float(desc="Diffusion data resolution (assuming isotropic)")
     dilate = traits.Int(desc="Dilation range")
 
 
@@ -136,8 +136,8 @@ class CreateList(SimpleInterface):
 class _CombineStringsInputSpec(BaseInterfaceInputSpec):
     input1 = traits.Str(mandatory=True, desc="Input string 1")
     input2 = traits.Str(mandatory=True, desc="Input string 2")
-    input3 = traits.Str("", desc="Input string 3")
-    input4 = traits.Str("", desc="Input string 4")
+    input3 = traits.Str("", usedefault=True, desc="Input string 3")
+    input4 = traits.Str("", usedefault=True, desc="Input string 4")
 
 
 class _StringOutputSpec(TraitedSpec):

hcpdiffpy/main.py

@@ -7,7 +7,7 @@ import nipype.pipeline as pe
 from hcpdiffpy.container import SimgCmd
 from hcpdiffpy.interfaces.data import InitData, SaveData
 from hcpdiffpy.interfaces.preproc import (
-    EddyIndex, EddyPostProc, ExtractB0, DilateMask, MergeBFiles, Rescale, RotateBVec2Str,
+    BETMask, EddyIndex, EddyPostProc, ExtractB0, DilateMask, MergeBFiles, Rescale, RotateBVec2Str,
     PrepareTopup, WBDilate)
 from hcpdiffpy.interfaces.utilities import (
     CreateList, CombineStrings, DiffRes, FlattenList, ListItem, PickDiffFiles, SplitDiffFiles,
@@ -20,56 +20,53 @@ def main() -> None:
     parser = argparse.ArgumentParser(
         description="HCP Pipeline for diffusion preprocessing",
         formatter_class=lambda prog: argparse.ArgumentDefaultsHelpFormatter(prog, width=100))
-    parser.add_argument(
+    required = parser.add_argument_group("required arguments")
+    required.add_argument(
         "subject_dir", type=Path,
         help="Absolute path to the subject's data folder (organised in HCP-like structure)")
-    parser.add_argument("subject", type=str, help="Subject ID")
-    parser.add_argument("ndirs", nargs="+", help="List of numbers of directions")
-    parser.add_argument("phases", nargs="+", help="List of 2 phase encoding directions")
-    parser.add_argument("echo_spacing", type=float, help="Echo spacing used for acquisition in ms")
-    parser.add_argument(
-        "--workdir", type=Path, dest="work_dir", default=Path.cwd(),
-        help="Absolute path to work directory")
-    parser.add_argument(
-        "--output_dir", type=Path, dest="output_dir", default=Path.cwd(),
-        help="Absolute path to output directory")
-    parser.add_argument(
-        "--fsl_simg", type=Path, dest="fsl_simg", default=None,
-        help="singularity image to use for command line functions from FSL")
-    parser.add_argument(
-        "--fs_simg", type=Path, dest="fs_simg", default=None,
-        help="singularity image to use for command line functions from FreeSurfer")
-    parser.add_argument(
-        "--wb_simg", type=Path, dest="wb_simg", default=None,
-        help="singularity image to use for command line functions from Connectome Workbench")
-    parser.add_argument(
-        "--condordag", dest="condordag", action="store_true",
-        help="Submit workflow as DAG to HTCondor")
+    required.add_argument("subject", type=str, help="Subject ID")
+    required.add_argument(
+        "echo_spacing", type=float, help="Echo spacing used for acquisition in ms")
+    required.add_argument("--ndirs", required=True, nargs="+", help="List of numbers of directions")
+    required.add_argument(
+        "--phases", required=True, nargs="+", help="List of phase encoding directions")
+    optional = parser.add_argument_group("optional arguments")
+    optional.add_argument(
+        "--work_dir", type=Path, default=Path.cwd(), help="Absolute path to work directory")
+    optional.add_argument(
+        "--output_dir", type=Path, default=Path.cwd(), help="Absolute path to output directory")
+    optional.add_argument("--fsl_simg", type=Path, default=None, help="singularity image for FSL")
+    optional.add_argument(
+        "--fs_simg", type=Path, default=None, help="singularity image for FreeSurfer")
+    optional.add_argument(
+        "--wb_simg", type=Path, default=None, help="singularity image for Connectome Workbench")
+    optional.add_argument("--condordag", action="store_true", help="Submit as DAG to HTCondor")
     config = vars(parser.parse_args())
 
     # Set-up
-    config["tmp_dir"] = Path(config["work_dir"], "hcp_proc_tmp")
+    config["output_dir"].mkdir(parents=True, exist_ok=True)
+    config["tmp_dir"] = Path(config["work_dir"], f"hcpdiff_{config['subject']}_tmp")
     config["tmp_dir"].mkdir(parents=True, exist_ok=True)
     config["keys"] = [
         f"dir{ndir}_{phase}" for ndir in sorted(config["ndirs"])
         for phase in sorted(config["phases"])]
-    d_iterables = [("ndir", config["ndirs"]), ("phases", config["phases"])]
+    d_iterables = [("ndir", config["ndirs"]), ("phase", config["phases"])]
     fsl_cmd = SimgCmd(config, config["fsl_simg"])
     fs_cmd = SimgCmd(config, config["fs_simg"])
     wb_cmd = SimgCmd(config, config["wb_simg"])
     topup_config = Path(base_dir, "utilities", "b02b0.cnf")
     sch_file = Path(base_dir, "utilities", "bbr.sch")
-    fs_dir = Path(config["subject_dir"], "T1w", config["subject"])
+    fs_dir = Path(config["subject_dir"], "T1w")
 
     # Workflow set-up
     hcpdiff_wf = pe.Workflow(f"hcpdiff_{config['subject']}_wf", base_dir=config["work_dir"])
-    hcpdiff_wf.config["execution"]["try_hard_link_datasink"] = "false"
+    hcpdiff_wf.config["execution"]["remove_node_directories"] = "true"
     hcpdiff_wf.config["execution"]["crashfile_format"] = "txt"
     hcpdiff_wf.config["execution"]["stop_on_first_crash"] = "true"
 
     # Get data
     init_data = pe.Node(InitData(config=config), "init_data")
-    d_files = pe.Node(PickDiffFiles(), "split_diff_files", ietrables=d_iterables)
+    d_files = pe.Node(PickDiffFiles(), "d_files", iterables=d_iterables)
 
     hcpdiff_wf.connect([(init_data, d_files, [("d_files", "d_files")])])
 
@@ -83,7 +80,7 @@ def main() -> None:
     scale = pe.Node(fsl.ImageMeants(command=fsl_cmd.cmd("fslmeants")), "scale")
     rescale = pe.JoinNode(
         Rescale(fsl_cmd=fsl_cmd, config=config), "rescale",
-        joinfield=["scale_files"], joinsource="split_diff_files")
+        joinfield=["scale_files"], joinsource="d_files")
 
     hcpdiff_wf.connect([
         (d_files, mean_dwi, [("data_file", "in_file")]),
@@ -93,7 +90,7 @@ def main() -> None:
         (merge_b0s, mean_b0, [("merged_file", "in_file")]),
         (mean_b0, scale, [("out_file", "in_file")]),
         (init_data, rescale, [("d_files", "d_files")]),
-        (scale, rescale, [("out_file", "scale_file")])])
+        (scale, rescale, [("out_file", "scale_files")])])
 
     # 1.2. Prepare b0s and index files for topup
     update_d_files = pe.Node(UpdateDiffFiles(config=config), "update_d_files")
@@ -131,7 +128,7 @@ def main() -> None:
     b01_files = pe.Node(CreateList(), "b01_files")
     apply_topup = pe.Node(
         fsl.ApplyTOPUP(command=fsl_cmd.cmd("applytopup"), method="jac"), "apply_topup")
-    nodiff_mask = pe.Node(fsl.BET(command=fsl_cmd.cmd("bet"), frac=0.2, mask=True), "nodiff_mask")
+    nodiff_mask = pe.Node(BETMask(config=config, fsl_cmd=fsl_cmd), "nodiff_mask")
 
     hcpdiff_wf.connect([
         (update_d_files, prepare_topup, [("d_files", "d_files")]),
@@ -178,18 +175,20 @@ def main() -> None:
     fov_mask = pe.Node(
         fsl.ImageMaths(command=fsl_cmd.cmd("fslmaths"), args="-abs -Tmin -bin -fillh"), "fov_mask")
     mask_args = pe.Node(CombineStrings(input1="-mas "), "mask_args")
-    mask_data = pe.Node(fsl.ImageMaths(command=fsl_cmd.cmd("fslmaths")), "mask_data")
+    fmask_data = pe.Node(fsl.ImageMaths(command=fsl_cmd.cmd("fslmaths")), "fmask_data")
     thr_data = pe.Node(fsl.ImageMaths(command=fsl_cmd.cmd("fslmaths"), args="-thr 0"), "thr_data")
 
     hcpdiff_wf.connect([
         (d_filetype, postproc, [
             ("bval_files", "bval_files"), ("bvec_files", "bvec_files"),
             ("data_files", "rescaled_files")]),
+        (eddy, postproc, [
+            ("out_corrected", "eddy_corrected_file"), ("out_rotated_bvecs", "eddy_bvecs_file")]),
         (postproc, fov_mask, [("combined_dwi_file", "in_file")]),
         (fov_mask, mask_args, [("out_file", "input2")]),
-        (postproc, mask_data, [("combined_dwi_file", "in_file")]),
-        (mask_args, mask_data, [("output", "args")]),
-        (mask_data, thr_data, [("out_file", "in_file")])])
+        (postproc, fmask_data, [("combined_dwi_file", "in_file")]),
+        (mask_args, fmask_data, [("output", "args")]),
+        (fmask_data, thr_data, [("out_file", "in_file")])])
 
     # 3.2. DiffusionToStructural
     # 3.2.1. nodiff-to-T1
@@ -234,7 +233,8 @@ def main() -> None:
             subject_id=config["subject"]),
         "bbr_epi2t1")
     tkr_diff2str = pe.Node(
-        freesurfer.Tkregister2(command=fs_cmd.cmd("tkregister2"), noedit=True), "tkr_diff2str")
+        freesurfer.Tkregister2(command=fs_cmd.cmd("tkregister2"), noedit=True, fsl_out=True),
+        "tkr_diff2str")
     diff2str = pe.Node(
         fsl.ConvertXFM(command=fsl_cmd.cmd("convert_xfm"), concat_xfm=True), "diff2str")
 
@@ -282,8 +282,8 @@ def main() -> None:
     dilate_mask = pe.Node(DilateMask(fsl_cmd=fsl_cmd), "dilate_mask")
     thr_fmask = pe.Node(
         fsl.ImageMaths(command=fsl_cmd.cmd("fslmaths"), args="-thr 0.999 -bin"), "thr_fmask")
-    masks_args = pe.Node(CombineStrings(input1="-mas ", input3="-mas "), "masks_args")
-    mask_data = pe.Node(fsl.ImageMaths(command=fsl_cmd.cmd("fslmaths")), "fmask_data")
+    masks_args = pe.Node(CombineStrings(input1="-mas ", input3=" -mas "), "masks_args")
+    mask_data = pe.Node(fsl.ImageMaths(command=fsl_cmd.cmd("fslmaths")), "mask_data")
     nonneg_data = pe.Node(
         fsl.ImageMaths(command=fsl_cmd.cmd("fslmaths"), args="-thr 0"), "nonneg_data")
     mean_mask = pe.Node(fsl.ImageMaths(command=fsl_cmd.cmd("fslmaths"), args="-Tmean"), "mean_mask")
@@ -323,7 +323,7 @@ def main() -> None:
         hcpdiff_wf.run(
             plugin="CondorDAGMan",
             plugin_args={
-                "dagman_args": f"-outfile_dir {config['work_dir']} -import_env",
+                "dagman_args": f"-outfile_dir {config['tmp_dir']} -import_env",
                 "wrapper_cmd": Path(base_dir, "utilities", "venv_wrapper.sh"),
                 "override_specs": "request_memory = 5 GB\nrequest_cpus = 1"})
     else: