Freesurfer update

dbsc/workflow/rules/freesurfer.smk

@@ -1,69 +1,73 @@
-from os.path import join 
+from os.path import join
+
+# Directories
+freesurfer_dir = join(config["out_dir"], "freesurfer")
+
+# Freesurfer references (with additional in rules as necessary)
+# B. Fischl, A. van der Kouwe, C. Destrieux, E. Halgren, F. Ségonne, D.H. Salat, E. Busa, L.J. Seidman, J. Goldstein, D. Kennedy, V. Caviness, N. Makris, B. Rosen, A.M. Dale. Automatically parcellating the human cerebral cortex. Cereb. Cortex, 14 (2004), pp. 11-22, 10.1093/cercor/bhg087
+
 
 rule thalamic_segmentation:
+    """Perform thalamus segmentation
+
+    Reference: J.E. Iglesias, R. Insausti, G. Lerma-Usabiaga, M. Bocchetta, K. Van Leemput, D.N. Greve, A. van der Kouwe, B. Fischl, C. Caballero-Gaudes, P.M. Paz-Alonso. A probabilistic atlas of the human thalamic nuclei combining ex vivo MRI and histology. NeuroImage, 183 (2018), pp. 314-326, 10.1016/j.neuroimage.2018.08.012
+
+    NOTE: Output file is defined by Freesurfer script
+    """
     input:
-        freesurfer_dir=join(config['bids_dir'], 'derivatives/freesurfer'),
+        freesurfer_dir=freesurfer_dir,
     params:
         fs_license=config["fs_license"],
-        threads=workflow.cores,
     output:
-        thal_seg=bids(
-            root=join(config["output_dir"], "freesurfer"),
-            datatype="anat",
-            suffix="T1.mgz"
-            space="Freesurfer",
-            **config["subj_wildcards"],
-        ),
-    container: 
-        config['singularity']['freesurfer']
+        thal_seg=join(freesurfer_dir, "{subject}/mri/ThalamicNuclei.v12.T1.mgz"),
+    threads: workflow.cores
+    container:
+        config["singularity"]["freesurfer"]
     shell:
         "FS_LICENSE={params.fs_license} "
         "ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS={params.threads} "
         "SUBJECTS_DIR={input.freesurfer_dir} "
         "segmentThalamicNuclei.sh {{subject}}"
 
-rule mgz_to_nii:
-    input: 
+
+fs_out = bids(
+    root=freesurfer_dir,
+    datatype="anat",
+    suffix="{suffix}",
+    space="{space}",
+    **config["subj_wildcards"],
+)
+
+
+rule mgz2nii:
+    """
+    Convert from .mgz to .nii.gz
+
+    NOTE: During conversion, files are renamed to BIDS-esque formatting
+    """
+    input:
         thal=rules.thalamic_segmentation.output.thal_seg,
-        aparcaseg=join(config["output_dir"], "freesurfer/{subject}/mri/aparc+aseg.mgz")
-        lRibbon=join(config["output_dir"], "freesurfer/{subject}/mri/lh.ribbon.mgz")
-        rRibbon=join(config["output_dir"], "freesurfer/{subject}/mri/rh.ribbon.mgz")
+        aparcaseg=join(freesurfer_dir, "{subject}/mri/aparc+aseg.mgz"),
+        lRibbon=join(freesurfer_dir, "{subject}/mri/lh.ribbon.mgz"),
+        rRibbon=join(freesurfer_dir, "{subject}/mri/rh.ribbon.mgz"),
     params:
         fs_license=config["fs_license"],
-        threads=workflow.cores,
     output:
-        thal=bids(
-            root=join(config["output_dir"], "freesurfer"),
-            datatype="anat",
-            suffix="T1.nii.gz"
-            space="Freesurfer",
-            **config["subj_wildcards"],
+        thal=expand(
+            fs_out, suffix="thalamus.nii.gz", space="Freesurfer", allow_missing=True
         ),
-        aparcaseg=bids(
-            root=join(config["output_dir"], "freesurfer"),
-            datatype="anat",
-            suffix="aparcaseg.nii.gz"
-            space="Freesurfer",
-            **config["subj_wildcards"],
+        aparcaseg=expand(
+            fs_out, suffix="aparcaseg.nii.gz", space="Freesurfer", allow_missing=True
         ),
-        ribbon_mgz=temp(
-            bids(
-                root=join(config['output_dir'], 'freesurfer'),
-                datatype='anat',
-                suffix='ribbon.mgz'
-                space='Freesurfer',
-                **config["subj_wildcards"],
-            )
+        ribbon_mgz=expand(
+            fs_out, suffix="ribbon.mgz", space="Freesurfer", allow_missing=True
         ),
-        ribbon=bids(
-            root=join(config['output_dir'], 'freesurfer'),
-            datatype='anat',
-            suffix='ribbon.nii.gz'
-            space='Freesurfer',
-            **config["subj_wildcards"],
-        )
+        ribbon=expand(
+            fs_out, suffix="ribbon.nii.gz", space="Freesurfer", allow_missing=True
+        ),
+    threads: workflow.cores
     container:
-        config['singularity']['freesurfer']
+        config["singularity"]["freesurfer"]
     shell:
         "FS_LICENSE={params.fs_license} "
         "ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS={params.threads} "
@@ -72,42 +76,28 @@ rule mgz_to_nii:
         "mergeseg --src {input.lRibbon} --merge {input.rRibbon} --o {output.ribbon_mgz}"
         "mri_convert {output.ribbon_mgz} {output.ribbon}"
 
+
 rule fs_xfm_to_native:
+    """Transform from Freesurfer space to T1w space"""
     input:
-        thal=rules.mgz_to_nii.output.thal,
-        aparcaseg=rules.mgz_to_nii.output.aparcaseg,
-        ribbon=rules.mgz_to_nii.output.ribbon,
+        thal=rules.mgz2nii.output.thal,
+        aparcaseg=rules.mgz2nii.output.aparcaseg,
+        ribbon=rules.mgz2nii.output.ribbon,
         ref=bids(
             root=config["bids_dir"],
             datatype="anat",
             suffix="T1w.nii.gz",
             **config["subj_wildcards"],
         ),
     output:
-        thal=bids(
-            root=join(config["output_dir"], "freesurfer"),
-            datatype="anat",
-            suffix="T1.nii.gz"
-            space="T1w",
-            **config["subj_wildcards"],
+        thal=expand(fs_out, suffix="thalamus.nii.gz", space="T1w", allow_missing=True),
+        aparcaseg=expand(
+            fs_out, suffix="aparcaseg.nii.gz", space="T1w", allow_missing=True
         ),
-        aparcaseg=bids(
-            root=join(config["output_dir"], 'freesurfer'),
-            datatype='anat',
-            suffix="aparcaseg.nii.gz",
-            space='T1w',
-            **config["subj_wildcards"],
-        )
-        ribbon=bids(
-            root=join(config["output_dir"], 'freesurfer'),
-            datatype='anat',
-            suffix='ribbon.nii.gz',
-            space='T1w',
-            **config["subj_wildcards"],
-        )
+        ribbon=expand(fs_out, suffix="ribbon.nii.gz", space="T1w", allow_missing=True),
     container:
-        config['singularity']['neuroglia-core'],
+        config["singularity"]["neuroglia-core"]
     shell:
         "antsApplyTransforms -d 3 -n MultiLabel -i {input.thal} -r {input.ref} -o {output.thal} "
         "antsApplyTransforms -d 3 -n MultiLabel -i {input.aparcaseg} -r {input.ref} -o {output.aparcaseg} "
-        "antsApplyTransforms -d 3 -n MultiLabel -i {input.ribbon} -r {input.ref} -o {output.ribbon} "
+        "antsApplyTransforms -d 3 -n MultiLabel -i {input.ribbon} -r {input.ref} -o {output.ribbon} "