hippb500 workflow, output_spaces

.gitignore

@@ -0,0 +1,3 @@
+hippunfold.egg-info
+hippunfold/.snakemake
+hippunfold/__pycache__

hippunfold/config/nnunet_model_urls.yml

@@ -1,5 +1,7 @@
+
 T1w: 'https://zenodo.org/record/4508747/files/trained_model.3d_fullres.Task101_hcp1200_T1w.nnUNetTrainerV2.model_best.tar'
 T2w: 'https://zenodo.org/record/4508747/files/trained_model.3d_fullres.Task102_hcp1200_T2w.nnUNetTrainerV2.model_best.tar'
 T1T2w: 'https://zenodo.org/record/4508747/files/trained_model.3d_fullres.Task103_hcp1200_T1T2w.nnUNetTrainerV2.model_best.tar'
 b1000: 'https://zenodo.org/record/4508747/files/trained_model.3d_fullres.Task104_hcp1200_b1000.nnUNetTrainerV2.model_best.tar'
 trimodal: 'https://zenodo.org/record/4508747/files/trained_model.3d_fullres.Task105_hcp1200_trimodal.nnUNetTrainerV2.model_best.tar'
+hippb500: 'https://www.dropbox.com/s/y827qqkuj1ey1jb/trained_model.3d_fullres.Task110_hcp1200_b1000crop.nnUNetTrainerV2.model_best.tar'

hippunfold/config/snakebids.yml

@@ -54,6 +54,16 @@ pybids_inputs:
       - run
       - dir
 
+  hippb500:
+    filters:
+      suffix: 'b500'
+      extension: '.nii.gz'
+      invalid_filters: 'allow'
+      datatype: 'dwi'
+    wildcards:
+      - subject
+      - session
+
   T1w:
     filters:
       suffix: 'T1w'
@@ -122,7 +132,8 @@ parse_args:
     choices:
       - T1w
       - T2w
-      - b500
+      - hippb500
+      - dwi
       - segT1w
       - segT2w
     default: 
@@ -159,7 +170,15 @@ parse_args:
     default: False
     action: 'store_true'
 
-
+  --output_spaces:
+    choices: 
+      - 'cropT1w'
+      - 'corobl'
+    default:
+      - 'cropT1w'
+    nargs: '+'
+    help: 'Sets the output space for results (default: %(default)s)'
+
 #--- workflow specific configuration -- 
 
 singularity:
@@ -189,11 +208,12 @@ modality:
 
 #these will be downloaded to ~/.cache/hippunfold
 nnunet_model:
-    T1w: trained_model.3d_fullres.Task101_hcp1200_T1w.nnUNetTrainerV2.model_best.tar
-    T2w: trained_model.3d_fullres.Task102_hcp1200_T2w.nnUNetTrainerV2.model_best.tar
-    T1T2w: trained_model.3d_fullres.Task103_hcp1200_T1T2w.nnUNetTrainerV2.model_best.tar
-    b1000: trained_model.3d_fullres.Task104_hcp1200_b1000.nnUNetTrainerV2.model_best.tar
-    trimodal: trained_model.3d_fullres.Task105_hcp1200_trimodal.nnUNetTrainerV2.model_best.tar
+  T1w: trained_model.3d_fullres.Task101_hcp1200_T1w.nnUNetTrainerV2.model_best.tar
+  T2w: trained_model.3d_fullres.Task102_hcp1200_T2w.nnUNetTrainerV2.model_best.tar
+  T1T2w: trained_model.3d_fullres.Task103_hcp1200_T1T2w.nnUNetTrainerV2.model_best.tar
+  b1000: trained_model.3d_fullres.Task104_hcp1200_b1000.nnUNetTrainerV2.model_best.tar
+  trimodal: trained_model.3d_fullres.Task105_hcp1200_trimodal.nnUNetTrainerV2.model_best.tar
+  hippb500: trained_model.3d_fullres.Task110_hcp1200_b1000crop.nnUNetTrainerV2.model_best.tar
 
 
 use_mcr: True
@@ -206,10 +226,18 @@ java_home: '/cvmfs/soft.computecanada.ca/easybuild/software/2017/Core/java/1.8.0
 cnn_model:
   T2w: 'HCP1200-T2'
   T1w: 'HCP1200-T1'
-  b500: 'HCP1200-b1000'
+  hippb500: 'HCP1200-b1000'
   b1000: 'HCP1200-b1000'
 
 
+hippdwi_opts:
+  resample_dim: '734x720x67' # from 220x216x20 @ 1x1x1mm -> 0.3mm
+  bbox_x:
+    R: '383 510'
+    L: '224 351'
+  bbox_y: '198 453'
+
+
 #---- dwi options below, not used in the default workflow ----
 
 #eddy options only used for experimental dwi workflow

hippunfold/tests/hippunfold_test.py

@@ -15,17 +15,19 @@ def test_dry_runs(script_runner):
     with tempfile.TemporaryDirectory() as output_dir:
         ret = script_runner.run('hippunfold', 'test_data/bids_T1w',output_dir,'participant','-np','--modality','T1w')
         assert ret.success
-
+
+    with tempfile.TemporaryDirectory() as output_dir:
+        ret = script_runner.run('hippunfold', 'test_data/bids_hippb500',output_dir,'participant','-np','--modality','hippb500')
+        assert ret.success
+
     with tempfile.TemporaryDirectory() as output_dir:
         ret = script_runner.run('hippunfold', 'test_data/bids_T1w_longitudinal',output_dir,'participant','-np','--modality','T1w')
         assert ret.success
 
     with tempfile.TemporaryDirectory() as output_dir:
         ret = script_runner.run('hippunfold', 'test_data/bids_singleT2w_longitudinal',output_dir,'participant','-np')
-        assert ret.success
 
     with tempfile.TemporaryDirectory() as output_dir:
         ret = script_runner.run('hippunfold', 'test_data/bids_segT2w',output_dir,'participant','-np','--modality','segT2w')
         assert ret.success
 
-

hippunfold/workflow/Snakefile

@@ -5,9 +5,14 @@ from snakebids import bids
 
 configfile: 'config/snakebids.yml'
 
+#make adjustments to input params:
+if 'hippb500' in config['modality']: #right now hippb500 cannot use T1w output space
+    config['output_spaces'] = ['corobl']
+
 #get list of inputs to limit to (T1w + modalities)
 limit_to_list = set()
-limit_to_list.add('T1w') #always need T1w
+if 'cropT1w' in config['output_spaces']:
+    limit_to_list.add('T1w') 
 
 
 for modality in config['modality']:
@@ -48,16 +53,17 @@ def get_modality_key(modality):
 # case not all subjects possess the same modalities
 def get_final_specs():
     specs = []
-    for modality in config['modality']:
-        mod_index = get_modality_key(modality)
-        specs = specs + \
-            expand(bids(root='results',datatype='surf_{modality}',suffix='hippunfold.spec', **config['subj_wildcards']),
-                modality=modality,
-                subject=config['input_lists'][mod_index]['subject'],
-                session=config['sessions'])
+    if 'cropT1w' in  config['output_spaces']:
+        for modality in config['modality']:
+            mod_index = get_modality_key(modality)
+            specs = specs + \
+                expand(bids(root='results',datatype='surf_{modality}',suffix='hippunfold.spec', **config['subj_wildcards']),
+                    modality=modality,
+                    subject=config['input_lists'][mod_index]['subject'],
+                    session=config['sessions'])
 
     return specs
-          
+
 def get_final_subfields():
     subfields = []
     for modality in config['modality']:
@@ -66,7 +72,7 @@ def get_final_subfields():
             expand(bids(root='results',datatype='seg_{modality}',desc='subfields',suffix='dseg.nii.gz', space='{space}',hemi='{hemi}', **config['subj_wildcards']),
                 modality=modality,
                 hemi=['L','R'],
-                space=['cropT1w','unfold'],
+                space=config['output_spaces'],
                 subject=config['input_lists'][mod_index]['subject'],
                 session=config['sessions'])
 
@@ -77,23 +83,31 @@ def get_final_coords():
     for modality in config['modality']:
         mod_index = get_modality_key(modality)
         coords = coords + \
-            expand(bids(root='results',datatype='seg_{modality}',dir='{dir}',suffix='coords.nii.gz', space='cropT1w',hemi='{hemi}', **config['subj_wildcards']),
+            expand(bids(root='results',datatype='seg_{modality}',dir='{dir}',suffix='coords.nii.gz', space='{space}',hemi='{hemi}', **config['subj_wildcards']),
                 modality=modality,
                 dir=['AP','PD','IO'],
                 hemi=['L','R'],
+                space=config['output_spaces'],
                 subject=config['input_lists'][mod_index]['subject'],
                 session=config['sessions'])
 
     return coords
 
 def get_final_transforms():
     transforms = []
+
+    if 'cropT1w' in  config['output_spaces']: 
+        output_ref = 'T1w'
+    else:
+        output_ref = 'corobl'
+
     for modality in config['modality']:
         mod_index = get_modality_key(modality)
         transforms = transforms + \
-            expand(bids(root='results',datatype='seg_{modality}',**config['subj_wildcards'],suffix='xfm.nii.gz',hemi='{hemi}',from_='T1w',to='unfold',mode='image'),
+            expand(bids(root='results',datatype='seg_{modality}',**config['subj_wildcards'],suffix='xfm.nii.gz',hemi='{hemi}',from_='{space}',to='unfold',mode='image'),
                 modality=modality,
                 hemi=['L','R'],
+                space=output_ref,
                 subject=config['input_lists'][mod_index]['subject'],
                 session=config['sessions'])
 
@@ -102,27 +116,34 @@ def get_final_transforms():
 def get_final_anat():
     anat = []
     for modality in config['modality']:
+        if modality == 'hippb500':
+            continue
         mod_index = get_modality_key(modality)
         if mod_index == 'seg':
             modality_suffix = modality[3:]
         else:
             modality_suffix = modality
         anat = anat + \
-            expand(bids(root='results',datatype='seg_{modality}',desc='preproc',suffix='{modality_suffix}.nii.gz', space='cropT1w',hemi='{hemi}', **config['subj_wildcards']),
+            expand(bids(root='results',datatype='seg_{modality}',desc='preproc',suffix='{modality_suffix}.nii.gz', space='{space}',hemi='{hemi}', **config['subj_wildcards']),
                 modality=modality,
                 modality_suffix=modality_suffix,
                 hemi=['L','R'],
+                space=config['output_spaces'],
                 subject=config['input_lists'][mod_index]['subject'],
                 session=config['sessions'])
 
     return anat
 
+
 def get_final_qc():
     qc = []
-    qc = qc + expand(bids(root='work',datatype='qc',**config['subj_wildcards'],suffix='regqc.png',
+    if 'T1w' in config['modality']:
+        qc = qc + expand(bids(root='work',datatype='qc',**config['subj_wildcards'],suffix='regqc.png',
                 from_='subject', to=config['template']),**config['input_lists']['T1w'])
 
     for modality in config['modality']:
+        if modality == 'hippb500':
+            continue
         mod_index = get_modality_key(modality)
         qc = qc + \
             expand(bids(root='work',datatype='qc',suffix='dseg.png', desc='subfields',from_='{modality}',space='cropT1w',hemi='{hemi}', **config['subj_wildcards']),
@@ -138,6 +159,7 @@ def get_final_qc():
 
     return qc
 
+
 rule all:
     input: 
         get_final_specs(),
@@ -152,6 +174,8 @@ rule all_group_tsv:
         tsv = expand(bids(root='results',prefix='group',from_='{modality}',desc='subfields',suffix='volumes.tsv'),
                     modality=config['modality'])
 
+
+
 include: 'rules/common.smk'
 include: 'rules/preproc_t1.smk'
 
@@ -165,9 +189,12 @@ else:
 if 'segT2w' in config['modality'] or 'T2w' in config['modality']:
     include: 'rules/preproc_t2.smk'
 
-if 'b500' in config['modality'] or 'b1000' in config['modality']:
+if 'dwi' in config['modality']:
     include: 'rules/preproc_dwi.smk'
 
+if 'hippb500' in config['modality']:
+    include: 'rules/preproc_hippb500.smk'
+
 include: 'rules/autotop.smk'
 include: 'rules/gifti.smk'
 include: 'rules/subfields.smk'

hippunfold/workflow/rules/autotop.smk

@@ -26,8 +26,8 @@ def get_autotop_input (wildcards):
         nii = bids(root='work',datatype='anat',**config['subj_wildcards'],suffix='T2w.nii.gz',desc='cropped',space='corobl',hemi='{hemi}'),
     elif wildcards.modality == 'T1w':
         nii = bids(root='work',datatype='anat',**config['subj_wildcards'],suffix='T1w.nii.gz',desc='cropped',space='corobl',hemi='{hemi}'),
-    elif wildcards.modality == 'b500':
-        nii = bids(root='work/preproc_dwi',suffix='b500.nii.gz',desc='cropped',datatype='dwi',**config['subj_wildcards'],space='corobl',hemi='{hemi}'),
+    elif wildcards.modality == 'hippb500':
+        nii = bids(root='work',suffix='b500.nii.gz',desc='cropped',datatype='dwi',**config['subj_wildcards'],space='corobl',hemi='{hemi}'),
     elif wildcards.modality == 'b1000':
         nii = bids(root='work/preproc_dwi',suffix='b1000.nii.gz',desc='cropped',datatype='dwi',**config['subj_wildcards'],space='corobl',hemi='{hemi}'),
     else:
@@ -150,6 +150,17 @@ rule unflip_autotop_nii:
     shell: 'c3d {input} -flip x {output}'
 
 
+rule cp_coords_to_bids:
+    input:
+        nii = bids(root='work',**config['subj_wildcards'],suffix='autotop/coords-{dir}.nii.gz',desc='cropped',space='corobl',hemi='{hemi}',modality='{modality}'),
+    output:
+        nii = bids(root='work',datatype='seg_{modality}',dir='{dir}',suffix='coords.nii.gz', space='corobl',hemi='{hemi}', **config['subj_wildcards'])
+    container: config['singularity']['autotop']
+    group: 'subj'
+    shell: 'cp {input} {output}'
+
+
+
 rule compose_warps_corobl2unfold_rhemi:
     """ Compose corobl to unfold (unfold-template), for right hemi (ie no flip)"""
     input:
@@ -190,5 +201,4 @@ rule compose_warps_t1_to_unfold:
     group: 'subj'
     shell: 'ComposeMultiTransform 3 {output} -R {input.ref} {input.corobl2unfold} {input.t1w2corobl}'
 
-
 

hippunfold/workflow/rules/common.smk

@@ -7,3 +7,9 @@ def get_avg_or_cp_scans_cmd (wildcards, input, output):
     return cmd
 
 
+rule copy_to_results:
+    """ Generic rule for copying data from work to results"""
+    input: 'work/{file}'
+    output: 'results/{file}'
+    group: 'subj'
+    shell: 'cp {input} {output}'

hippunfold/workflow/rules/gifti.smk

@@ -72,12 +72,12 @@ def get_unfolded_surf_R_Lflip (wildcards):
         return  bids(root='work',datatype='surf_{modality}',suffix='{surfname}.surf.gii', space='unfolded',hemi='{hemi}flip', **config['subj_wildcards']).format(**wildcards)
 
 
-rule cp_unfolded_to_results:
-    """copy unfolded surf to results folder, from Lflip/R to L/R"""
+rule unflip_gii_unfolded:
+    """copy unfolded from Lflip to L"""
     input:
-        gii = get_unfolded_surf_R_Lflip
+        gii = bids(root='work',datatype='surf_{modality}',suffix='{surfname}.surf.gii', space='unfolded',hemi='{hemi}flip', **config['subj_wildcards'])
     output:
-        gii = bids(root='results',datatype='surf_{modality}',suffix='{surfname}.surf.gii', space='unfolded',hemi='{hemi,L|R}', **config['subj_wildcards'])
+        gii = bids(root='work',datatype='surf_{modality}',suffix='{surfname}.surf.gii', space='unfolded',hemi='{hemi,L}', **config['subj_wildcards'])
     container: config['singularity']['autotop']
     group: 'subj'
     shell:

hippunfold/workflow/rules/nnunet.smk

@@ -6,6 +6,8 @@ def get_nnunet_input (wildcards):
         nii = bids(root='work',datatype='anat',**config['subj_wildcards'],suffix='T2w.nii.gz',desc='cropped',space='corobl',hemi='{hemi}'),
     elif wildcards.modality == 'T1w':
         nii = bids(root='work',datatype='anat',**config['subj_wildcards'],suffix='T1w.nii.gz',desc='cropped',space='corobl',hemi='{hemi}'),
+    elif wildcards.modality == 'hippb500':
+        nii = bids(root='work',datatype='dwi',hemi='{hemi}',desc='cropped',space='corobl',suffix='b500.nii.gz',**config['subj_wildcards'] )
     else:
         raise ValueError('modality not supported for nnunet!')
     return nii

hippunfold/workflow/rules/preproc_hippb500.smk

@@ -0,0 +1,36 @@
+#get _b500.nii.gz from the bids dataset - this is a preprocessed avg b500 image
+
+rule resample_hippdwi_to_template:
+    """ Hipp DWI is already corobl, but just needs to be cropped
+        to get L and R subvolumes. We use predefined X and Y 
+        bounding boxes, and keep all Z slices (since is Z 
+        is smaller than in template). """
+    input:
+        b500 = config['input_path']['hippb500']
+    params:
+        resample_dim = config['hippdwi_opts']['resample_dim'],
+        bbox_x = lambda wildcards: config['hippdwi_opts']['bbox_x'][wildcards.hemi],
+        bbox_y = config['hippdwi_opts']['bbox_y']
+    output:
+        crop_b500 = bids(root='work',datatype='dwi',hemi='{hemi,L|R}',desc='cropped',space='corobl',suffix='b500.nii.gz',**config['subj_wildcards'] )
+    group: 'subj'
+    shell:
+        'c3d {input} -resample {params.resample_dim} -as UPSAMPLED '
+        ' -push UPSAMPLED -cmv -pop -popas COORDY -popas COORDX '
+        ' -push COORDX -thresh {params.bbox_x} 1 0 -as MASKX '
+        ' -push COORDY -thresh {params.bbox_y} 1 0 -as MASKY '
+        ' -push MASKX -push MASKY -multiply '
+        ' -push UPSAMPLED -multiply '
+        ' -trim 0vox -o {output}'
+
+
+rule lr_flip_b500:
+    input:
+        nii = bids(root='work',datatype='dwi',**config['subj_wildcards'],suffix='b500.nii.gz',desc='cropped',space='corobl',hemi='{hemi}'),
+    output:
+        nii = bids(root='work',datatype='dwi',**config['subj_wildcards'],suffix='b500.nii.gz',desc='cropped',space='corobl',hemi='{hemi,L}flip'),
+    container: config['singularity']['autotop']
+    group: 'subj'
+    shell:
+        'c3d {input} -flip x -o  {output}'
+

hippunfold/workflow/rules/resample_final_to_crop_t1.smk

@@ -59,7 +59,7 @@ rule resample_subfields_native_crop:
 
 rule resample_coords_native_crop:
     input:
-        nii = bids(root='work',**config['subj_wildcards'],suffix='autotop/coords-{dir}.nii.gz',desc='cropped',space='corobl',hemi='{hemi}',modality='{modality}'),
+        nii = bids(root='work',datatype='seg_{modality}',dir='{dir}',suffix='coords.nii.gz', space='corobl',hemi='{hemi}', **config['subj_wildcards']),
         xfm = bids(root='work',datatype='anat',**config['subj_wildcards'],suffix='xfm.txt',from_='T1w',to='corobl',desc='affine',type_='itk'),
         ref = bids(root='work',datatype='seg_{modality}',suffix='cropref.nii.gz', space='T1w',hemi='{hemi}', **config['subj_wildcards'])
     output:
@@ -94,6 +94,4 @@ rule resample_t2_to_crop:
     shell:
         'ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS={threads} '
         'antsApplyTransforms -d 3 --interpolation Linear -i {input.nii} -o {output.nii} -r {input.ref} -t {input.xfm}' 
-
-
-
+