Merge pull request #365 from khanlab/create-template

New participant_create_template analysis_level

.dryrun_test_all.sh

@@ -9,8 +9,8 @@ hippunfold test_data/bids_T1w test_out participant -np --modality T1w
 hippunfold test_data/bids_hippb500 test_out participant -np --modality hippb500 
 hippunfold test_data/bids_T1w_longitudinal test_out participant -np --modality T1w 
 hippunfold test_data/bids_singleT2w_longitudinal test_out participant -np --modality T2w 
-hippunfold test_data/bids_manualseg test_out participant -np --modality manualseg --path_manualseg test_data/bids_manualseg/sub-{subject}_hemi-{hemi}_dseg.nii.gz
-hippunfold test_data/bids_manualseg_1hemi test_out participant -np --modality manualseg --hemi L
+hippunfold test_data/bids_dsegtissue test_out participant -np --modality dsegtissue --path_dsegtissue test_data/bids_dsegtissue/sub-{subject}_hemi-{hemi}_dseg.nii.gz
+hippunfold test_data/bids_dsegtissue_1hemi test_out participant -np --modality dsegtissue --hemi L
 hippunfold test_data/bids_singleT2w test_out participant -np --modality T2w --t1_reg_template 
 hippunfold test_data/bids_singleT2w test_out participant -np --modality T2w --output_space T1w 
 hippunfold test_data/bids_T1w test_out participant -np --modality T1w --use-template-seg 

.github/workflows/python-testing.yml

@@ -96,11 +96,11 @@ jobs:
 
       - name: Test manualseg bids, with path override
         run: |
-          poetry run hippunfold test_data/bids_manualseg test_out participant -np --modality dsegtissue --path_dsegtissue test_data/bids_manualseg/sub-{subject}_hemi-{hemi}_dseg.nii.gz
+          poetry run hippunfold test_data/bids_dsegtissue test_out participant -np --modality dsegtissue --path_dsegtissue test_data/bids_dsegtissue/sub-{subject}_hemi-{hemi}_dseg.nii.gz
 
       - name: Test manualseg bids, one hemisphere hemi
         run: |
-          poetry run hippunfold test_data/bids_manualseg_1hemi test_out participant -np --modality dsegtissue --hemi L
+          poetry run hippunfold test_data/bids_dsegtissue_1hemi test_out participant -np --modality dsegtissue --hemi L
 
       - name: Test T2w with T1w template registration
         run: |

hippunfold/config/snakebids.yml

@@ -14,13 +14,16 @@ derivatives: False #will search in bids/derivatives if True; can also be path(s)
 #list of analysis levels in the bids app 
 analysis_levels: &analysis_levels
  - participant
+ - participant_create_template
  - group
 
 
 #mapping from analysis_level to set of target rules or files
 targets_by_analysis_level:
   participant:
     - ''  # if '', then the first rule is run
+  participant_create_template:
+    - 'all_create_template'
   group:
     - 'all_group_tsv'
 
@@ -62,7 +65,7 @@ pybids_inputs:
       - acquisition
       - run
 
-  manualseg:
+  dsegtissue:
     filters:
       suffix: 'dseg'
       extension: '.nii.gz'
@@ -75,20 +78,31 @@ pybids_inputs:
       - run
 
 
+  dsegsubfields:
+    filters:
+      suffix: 'dseg'
+      extension: '.nii.gz'
+      datatype: 'anat'
+    wildcards:
+      - subject
+      - session
+      - hemi
+      - acquisition
+      - run
 
 
 #configuration for the command-line parameters to make available
 # passed on the argparse add_argument()
 parse_args:
 
   --modality:
-    help: 'Type of image to run hippunfold on. For manualseg and be sure to match the label scheme applied in the look up table (lut) online. For best performance, please also make sure the input is approximately aligned to the template space, with separate hemi-L|R in the name(s). (default: %(default)s)'
+    help: 'Type of image to run hippunfold on. For dsegtissue and be sure to match the label scheme applied in the look up table (lut) online. For best performance, please also make sure the input is approximately aligned to the template space, with separate hemi-L|R in the name(s). (default: %(default)s)'
     required: True
     choices:
       - T1w
       - T2w
       - hippb500
-      - manualseg
+      - dsegtissue
 
 
   --template:
@@ -171,8 +185,8 @@ parse_args:
     help: 'Sets the bounding box size for the crop native (e.g. cropT1w space). Make this larger if your hippocampi in crop{T1w,T2w} space are getting cut-off. This must be in voxels (vox) not millimeters (mm). (default: %(default)s)'
     default: '256x256x256vox'
 
-  --resample_manualseg:
-    help: 'Resample the manual segmentation (for --modality manualseg), keeping the bounding box the same, but changing the number of voxels in the image. This can be specified as voxels ("300x300x300"), or a percentage ("20%%") or percentage in each direction ("20x20x40%%"). Note: the segmentation will always be subsequently cropped around the segmentation, with a padding of 5 voxels.  (default: %(default)s)'
+  --resample_dsegtissue:
+    help: 'Resample the manual segmentation (for --modality dsegtissue), keeping the bounding box the same, but changing the number of voxels in the image. This can be specified as voxels ("300x300x300"), or a percentage ("20%%") or percentage in each direction ("20x20x40%%"). Note: the segmentation will always be subsequently cropped around the segmentation, with a padding of 5 voxels.  (default: %(default)s)'
     default: null
 
 
@@ -368,7 +382,7 @@ template_files:
     crop_ref: tpl-dHCP_cohort-1_res-1_space-corobl_hemi-{hemi}_T2w.nii.gz
     crop_refT1w: tpl-dHCP_cohort-1_res-1_space-corobl_hemi-{hemi}_T1w.nii.gz
     Mask_crop: tpl-dHCP_cohort-1_res-1_space-corobl_hemi-{hemi}_desc-hipp_mask.nii.gz
-    dseg: tpl-dHCP_space-corobl_hemi-{hemi}_desc-tissuemanualseg_dseg.nii.gz
+    dseg: tpl-dHCP_space-corobl_hemi-{hemi}_desc-tissuedsegtissue_dseg.nii.gz
     coords: tpl-dHCP_dir-{dir}_hemi-{hemi}_space-corobl_label-{label}_desc-laplace_coords.nii.gz
     hemi_wildcards:
       - L
@@ -722,4 +736,4 @@ root: results
 workdir: null
 use_template_seg: False
 template_seg_smoothing_factor: 2
-resample_manualseg: null
+resample_dsegtissue: null

hippunfold/workflow/Snakefile

@@ -31,7 +31,7 @@ if (
     limit_to_list.add("T2w")
 if (
     config["modality"] == "hippb500"
-    or config["modality"] == "manualseg"
+    or config["modality"] == "dsegtissue"
     or "corobl" in config["output_spaces"]
 ):
     ref_spaces.append("corobl")
@@ -54,6 +54,8 @@ if "T2w" in config["modality"]:
 elif "T1w" in config["modality"]:
     template_modality = "T1w"
 
+if config["analysis_level"] == "participant_create_template":
+    limit_to_list.add("dsegsubfields")
 
 limit_to_list.add(config["modality"])
 
@@ -100,9 +102,9 @@ if "T1w" in limit_to_list:
 
 
 # include rules only as they are needed..
-if config["modality"] == "manualseg":
+if config["modality"] == "dsegtissue":
 
-    include: "rules/preproc_manualseg.smk"
+    include: "rules/preproc_dsegtissue.smk"
 
 else:
 
@@ -154,6 +156,11 @@ rule all:
         get_final_output(),
 
 
+rule all_create_template:
+    input:
+        get_create_template_output(),
+
+
 rule all_group_tsv:
     input:
         tsv=expand(

hippunfold/workflow/rules/common.smk

@@ -262,3 +262,72 @@ def get_download_dir():
         dirs = AppDirs("hippunfold", "khanlab")
         download_dir = dirs.user_cache_dir
     return download_dir
+
+
+def get_cifti_metric_types(label):
+    types_list = config["cifti_metric_types"][label]
+    if config["generate_myelin_map"]:
+        types_list.append("myelin.dscalar")
+    return types_list
+
+
+def get_gifti_metric_types(label):
+    types_list = config["gifti_metric_types"][label]
+    if config["generate_myelin_map"]:
+        types_list.append("myelin.shape")
+    return types_list
+
+
+def get_create_template_output():
+
+    files = []
+    for label in config["autotop_labels"]:
+        files.extend(
+            inputs[config["modality"]].expand(
+                bids(
+                    root=root,
+                    datatype="surf",
+                    suffix="{metric}.gii",
+                    space="{space}",
+                    hemi="{hemi_}",
+                    label=label,
+                    **inputs.subj_wildcards,
+                ),
+                metric=get_gifti_metric_types(label),
+                space="corobl",
+                hemi_=config["hemi"],
+            )
+        )
+        files.extend(
+            inputs[config["modality"]].expand(
+                bids(
+                    root=root,
+                    datatype="surf",
+                    suffix="{surftype}.surf.gii",
+                    space="{space}",
+                    hemi="{hemi_}",
+                    label=label,
+                    **inputs.subj_wildcards,
+                ),
+                metric=get_gifti_metric_types(label),
+                space=["corobl", "unfold"],
+                surftype=["inner", "outer", "midthickness"],
+                hemi_=config["hemi"],
+            )
+        )
+        files.extend(
+            inputs[config["modality"]].expand(
+                bids(
+                    root=root,
+                    datatype="surf",
+                    suffix="subfields.label.gii",
+                    space="corobl",
+                    hemi="{hemi_}",
+                    label="hipp",
+                    **inputs.subj_wildcards,
+                ),
+                space="corobl",
+                hemi_=config["hemi"],
+            )
+        )
+    return files

hippunfold/workflow/rules/coords.smk

@@ -1,5 +1,8 @@
 def get_labels_for_laplace(wildcards):
-    if config["skip_inject_template_labels"]:
+    if (
+        config["skip_inject_template_labels"]
+        or config["analysis_level"] == "participant_create_template"
+    ):
         seg = get_input_for_shape_inject(wildcards)
     else:
         seg = bids(

hippunfold/workflow/rules/gifti.smk

@@ -297,20 +297,6 @@ def get_cmd_spec_file(wildcards, input, output):
     return " && ".join(cmds)
 
 
-def get_cifti_metric_types(label):
-    types_list = config["cifti_metric_types"][label]
-    if config["generate_myelin_map"]:
-        types_list.append("myelin.dscalar")
-    return types_list
-
-
-def get_gifti_metric_types(label):
-    types_list = config["gifti_metric_types"][label]
-    if config["generate_myelin_map"]:
-        types_list.append("myelin.shape")
-    return types_list
-
-
 rule create_spec_file_hipp:
     input:
         metrics=lambda wildcards: expand(

hippunfold/workflow/rules/native_surf.smk

@@ -549,10 +549,11 @@ rule register_midthickness:
         "subj"
     container:
         config["singularity"]["autotop"]
+    threads: 16
     conda:
         "../envs/greedy.yaml"
     shell:
-        "greedy -d 3 -i {input.fixed} {input.moving} -n 30 -o {output.warp}"
+        "greedy -threads {threads} -d 3 -i {input.fixed} {input.moving} -n 30x0 -o {output.warp}"
 
 
 rule apply_halfsurf_warp_to_img:

hippunfold/workflow/rules/preproc_manualseg.smk → hippunfold/workflow/rules/preproc_dsegtissue.smk

@@ -1,12 +1,12 @@
 
-rule import_manualseg_with_resample:
+rule import_dseg_tissue:
     input:
-        in_img=partial(get_single_bids_input, component="manualseg"),
+        in_img=partial(get_single_bids_input, component="dsegtissue"),
     params:
         resample_cmd=(
             ""
-            if config["resample_manualseg"] == None
-            else "-resample {res}".format(res=config["resample_manualseg"])
+            if config["resample_dsegtissue"] == None
+            else "-resample {res}".format(res=config["resample_dsegtissue"])
         ),
         crop_cmd="-trim 5vox",  #leave 5 voxel padding
     output:

hippunfold/workflow/rules/qc.smk

@@ -126,7 +126,7 @@ def get_bg_img_for_subfield_qc(wildcards):
             hemi="{hemi}",
             **inputs.subj_wildcards,
         )
-    elif config["modality"] == "manualseg":
+    elif config["modality"] == "dsegtissue":
         # blank image as bg
         return bids(
             root=work,

hippunfold/workflow/rules/shape_inject.smk

@@ -3,7 +3,7 @@
 
 
 def get_input_for_shape_inject(wildcards):
-    if config["modality"] == "manualseg":
+    if config["modality"] == "dsegtissue":
         seg = bids(
             root=work,
             datatype="anat",
@@ -26,7 +26,7 @@ def get_input_for_shape_inject(wildcards):
 
 
 def get_input_splitseg_for_shape_inject(wildcards):
-    if config["modality"] == "manualseg":
+    if config["modality"] == "dsegtissue":
         seg = bids(
             root=work,
             datatype="anat",

hippunfold/workflow/rules/subfields.smk

@@ -1,3 +1,64 @@
+rule import_dseg_subfields:
+    """read in volumetric subfield dseg, used for participant_create_template only"""
+    input:
+        vol_dseg=partial(get_single_bids_input, component="dsegsubfields"),
+        label_list=Path(workflow.basedir) / "../resources/atlas-v2/labellist_withdg.txt",
+    output:
+        label_dseg=bids(
+            root=work,
+            datatype="anat",
+            **inputs.subj_wildcards,
+            desc="subfields",
+            suffix="dseg.nii.gz",
+            space="corobl",
+            hemi="{hemi}",
+        ),
+    container:
+        config["singularity"]["autotop"]
+    conda:
+        "../envs/workbench.yaml"
+    group:
+        "subj"
+    shell:
+        "wb_command -volume-label-import {input.vol_dseg} {input.label_list} {output.label_dseg}"
+
+
+rule subfields_to_label_gifti:
+    input:
+        vol=bids(
+            root=work,
+            datatype="anat",
+            **inputs.subj_wildcards,
+            desc="subfields",
+            suffix="dseg.nii.gz",
+            space="corobl",
+            hemi="{hemi}",
+        ),
+        surf_gii=bids(
+            root=work,
+            datatype="surf",
+            suffix="midthickness.surf.gii",
+            space="corobl",
+            hemi="{hemi}",
+            label="hipp",
+            **inputs.subj_wildcards,
+        ),
+    output:
+        label_gii=bids(
+            root=root,
+            datatype="surf",
+            suffix="subfields.label.gii",
+            space="corobl",
+            hemi="{hemi}",
+            label="hipp",
+            **inputs.subj_wildcards,
+        ),
+    conda:
+        "../envs/workbench.yaml"
+    container:
+        config["singularity"]["autotop"]
+    shell:
+        "wb_command -volume-label-to-surface-mapping {input.vol} {input.surf_gii} {output.label_gii}"
 
 
 rule label_subfields_from_vol_coords_corobl: