synthSR dwi-t1 reg, motioncorr, new test datasets (#53)

- adds synthSR to images used for T1 to dwi registration
- ensures motion-correction (moco) is used before b0s are averaged at any point in the workflow 
   - moco rules moved to motioncorr.smk
- apply topup only on b0s (we were previously applying to whole dwi, then pulling out b0s) so this is just for efficiency and to ensure we use the moco'd b0s
- added new test datasets with only 14 dirs (and a single b0 in each scan) to test this edge case (they conveniently run faster than the previous test data too)
   - one for topup (LR+RL scan), and one for synthsrsdc workflows (single AP scan)
   - added DAGS for these

snakedwi/config/snakebids.yml

@@ -1,5 +1,7 @@
-bids_dir: 'test_data/bids_multishell_downsampled'
-output_dir: 'test_output/snakedwi_multishell_downsampled'
+bids_dir: 'test_data/bids_downsampled_1bzero' #for testing topup workflow
+#bids_dir: 'test_data/bids_downsampled_1bzero_1scan' #for testing synthsrsdc workflow
+
+output_dir: 'test_output/snakedwi_downsampled_1bzero'
 
 #enable printing debug statements during parsing -- disable if generating dag visualization
 debug: False

snakedwi/test_data/bids_downsampled_1bzero/dataset_description.json

@@ -0,0 +1,23 @@
+{
+    "Acknowledgements": "TODO: whom you want to acknowledge",
+    "Authors": [
+        "TODO:",
+        "First1 Last1",
+        "First2 Last2",
+        "..."
+    ],
+    "BIDSVersion": "1.0.1",
+    "DatasetDOI": "TODO: eventually a DOI for the dataset",
+    "Funding": [
+        "TODO",
+        "GRANT #1",
+        "GRANT #2"
+    ],
+    "HowToAcknowledge": "TODO: describe how to acknowledge -- either cite a corresponding paper, or just in acknowledgement section",
+    "License": "TODO: choose a license, e.g. PDDL (http://opendatacommons.org/licenses/pddl/)",
+    "Name": "TODO: name of the dataset",
+    "ReferencesAndLinks": [
+        "TODO",
+        "List of papers or websites"
+    ]
+}

snakedwi/test_data/bids_downsampled_1bzero/participants.tsv

@@ -0,0 +1,2 @@
+participant_id
+sub-001

snakedwi/test_data/bids_downsampled_1bzero/sub-001/dwi/sub-001_acq-LR_dwi.bval

@@ -0,0 +1 @@
+5 2590 2590 1305 2585 2600 2600 2605 1300 2600 2590 2595 2605 1295

snakedwi/test_data/bids_downsampled_1bzero/sub-001/dwi/sub-001_acq-LR_dwi.bvec

@@ -0,0 +1,3 @@
+0.560676 0.625298 0.737581 0.25836 -0.0507769 0.865369 0.162918 0.0447011 -0.959064 -0.756343 -0.563563 -0.622254 -0.474609 0.071442 
+-0.516028 0.780298 -0.596422 -0.189986 0.19533 -0.0487741 -0.938974 0.640698 -0.130358 -0.208598 0.574387 -0.54356 0.276886 -0.972928 
+-0.647579 -0.0116906 -0.316631 0.947183 -0.979422 0.498756 0.30296 0.766491 0.251401 0.620026 -0.593697 -0.563332 0.835512 -0.219789

snakedwi/test_data/bids_downsampled_1bzero/sub-001/dwi/sub-001_acq-LR_dwi.json

@@ -0,0 +1,4 @@
+{
+  "EffectiveEchoSpacing": 0.0014,
+  "PhaseEncodingDirection": "i"
+}

snakedwi/test_data/bids_downsampled_1bzero/sub-001/dwi/sub-001_acq-RL_dwi.bval

@@ -0,0 +1 @@
+5 2590 2590 1305 2585 2600 2600 2605 1300 2600 2590 2595 2605 1295 

snakedwi/test_data/bids_downsampled_1bzero/sub-001/dwi/sub-001_acq-RL_dwi.bvec

@@ -0,0 +1,3 @@
+-0.590593 0.625298 0.737581 0.25836 -0.0507769 0.865369 0.162918 0.0447011 -0.959064 -0.756343 -0.563563 -0.622254 -0.474609 0.071442 
+0.632757 0.780298 -0.596422 -0.189986 0.19533 -0.0487741 -0.938974 0.640698 -0.130358 -0.208598 0.574387 -0.54356 0.276886 -0.972928
+-0.500817 -0.0116906 -0.316631 0.947183 -0.979422 0.498756 0.30296 0.766491 0.251401 0.620026 -0.593697 -0.563332 0.835512 -0.219789 

snakedwi/test_data/bids_downsampled_1bzero/sub-001/dwi/sub-001_acq-RL_dwi.json

@@ -0,0 +1,4 @@
+{
+  "EffectiveEchoSpacing": 0.0014,
+  "PhaseEncodingDirection": "i-"
+}

snakedwi/test_data/bids_downsampled_1bzero_1scan/dataset_description.json

@@ -0,0 +1,23 @@
+{
+    "Acknowledgements": "TODO: whom you want to acknowledge",
+    "Authors": [
+        "TODO:",
+        "First1 Last1",
+        "First2 Last2",
+        "..."
+    ],
+    "BIDSVersion": "1.0.1",
+    "DatasetDOI": "TODO: eventually a DOI for the dataset",
+    "Funding": [
+        "TODO",
+        "GRANT #1",
+        "GRANT #2"
+    ],
+    "HowToAcknowledge": "TODO: describe how to acknowledge -- either cite a corresponding paper, or just in acknowledgement section",
+    "License": "TODO: choose a license, e.g. PDDL (http://opendatacommons.org/licenses/pddl/)",
+    "Name": "TODO: name of the dataset",
+    "ReferencesAndLinks": [
+        "TODO",
+        "List of papers or websites"
+    ]
+}

snakedwi/test_data/bids_downsampled_1bzero_1scan/participants.tsv

@@ -0,0 +1,2 @@
+participant_id
+sub-001

snakedwi/test_data/bids_downsampled_1bzero_1scan/sub-001/dwi/sub-001_acq-AP_dwi.bval

@@ -0,0 +1 @@
+5 2590 2590 1305 2585 2600 2600 2605 1300 2600 2590 2595 2605 1295

snakedwi/test_data/bids_downsampled_1bzero_1scan/sub-001/dwi/sub-001_acq-AP_dwi.bvec

@@ -0,0 +1,3 @@
+0.560676 0.625298 0.737581 0.25836 -0.0507769 0.865369 0.162918 0.0447011 -0.959064 -0.756343 -0.563563 -0.622254 -0.474609 0.071442 
+-0.516028 0.780298 -0.596422 -0.189986 0.19533 -0.0487741 -0.938974 0.640698 -0.130358 -0.208598 0.574387 -0.54356 0.276886 -0.972928 
+-0.647579 -0.0116906 -0.316631 0.947183 -0.979422 0.498756 0.30296 0.766491 0.251401 0.620026 -0.593697 -0.563332 0.835512 -0.219789

snakedwi/test_data/bids_downsampled_1bzero_1scan/sub-001/dwi/sub-001_acq-AP_dwi.json

@@ -0,0 +1,4 @@
+{
+  "EffectiveEchoSpacing": 0.0014,
+  "PhaseEncodingDirection": "j"
+}

snakedwi/workflow/Snakefile

@@ -5,6 +5,7 @@ include: "rules/setup.smk"
 include: "rules/common.smk"
 include: "rules/workflowopts.smk"
 include: "rules/prepdwi.smk"
+include: "rules/motioncorr.smk"
 include: "rules/sdc.smk"
 include: "rules/eddy.smk"
 include: "rules/bedpost.smk"

snakedwi/workflow/rules/motioncorr.smk

@@ -0,0 +1,189 @@
+
+rule moco_subj_bzeros_4d:
+    """ run motion-correction (rigid reg to init volume) on the bzeros 
+    from a concatenated b0 volume, per subj/ses"""
+    input:
+        nii_4d=bids(
+            root=work,
+            suffix="b0s.nii.gz",
+            datatype="dwi",
+            desc="degibbs",
+            **subj_wildcards
+        ),
+    output:
+        affine_dir=directory(
+            bids(
+                root=work,
+                suffix="transforms",
+                desc="moco",
+                datatype="dwi",
+                **subj_wildcards
+            )
+        ),
+        nii_4d=bids(
+            root=work,
+            suffix="b0s.nii.gz",
+            desc="moco",
+            datatype="dwi",
+            **subj_wildcards
+        ),
+        nii_avg3d=bids(
+            root=work,
+            suffix="b0.nii.gz",
+            desc="moco",
+            datatype="dwi",
+            **subj_wildcards
+        ),
+    threads: 8  #doesn't need to be more than the number of bzeros 
+    resources:
+        mem_mb=32000,
+    shadow:
+        "minimal"
+    container:
+        config["singularity"]["prepdwi"]  #-- this rule needs niftyreg, c3d and mrtrix
+    group:
+        "subj"
+    shell:
+        "c4d {input.nii_4d} -slice w 0:-1 -oo dwi_%03d.nii.gz && "
+        "if [ ! -e 'dwi_001.nii.gz' ]; then "
+        "    mkdir -p {output.affine_dir} &&  "
+        "    cp dwi_000.nii.gz {output.nii_4d} && "
+        "    mrmath {output.nii_4d} mean {output.nii_avg3d} -axis 3; "
+        "else "
+        "parallel --eta --jobs {threads} "
+        "reg_aladin -flo dwi_{{1}}.nii.gz  -ref dwi_000.nii.gz -res warped_{{1}}.nii.gz -aff affine_xfm_ras_{{1}}.txt "
+        " ::: `ls dwi_???.nii.gz | tail -n +2 | grep -Po '(?<=dwi_)[0-9]+'` && "
+        " mkdir -p {output.affine_dir} && cp affine_xfm_ras_*.txt {output.affine_dir} && "
+        " echo -e '1 0 0 0\n0 1 0 0\n0 0 1 0\n0 0 0 1' > {output.affine_dir}/affine_xfm_ras_000.txt && "
+        " mrcat dwi_000.nii.gz warped_*.nii.gz {output.nii_4d} && "
+        " mrmath {output.nii_4d} mean {output.nii_avg3d} -axis 3; "
+        "fi"
+
+
+rule moco_scan_bzeros_4d:
+    """ run motion-correction (rigid reg to init volume) on the bzeros from 
+    a single acquisition (ie before concatenating across acq/dir entities)"""
+    input:
+        nii_4d=bids(
+            root=work,
+            suffix="b0s.nii.gz",
+            datatype="dwi",
+            desc="degibbs",
+            **input_wildcards["dwi"]
+        ),
+    output:
+        affine_dir=directory(
+            bids(
+                root=work,
+                suffix="transforms",
+                desc="moco",
+                datatype="dwi",
+                **input_wildcards["dwi"]
+            )
+        ),
+        nii_avg3d=bids(
+            root=work,
+            suffix="b0.nii.gz",
+            desc="moco",
+            datatype="dwi",
+            **input_wildcards["dwi"]
+        ),
+    threads: 8  #doesn't need to be more than the number of bzeros 
+    resources:
+        mem_mb=32000,
+    shadow:
+        "minimal"
+    container:
+        config["singularity"]["prepdwi"]  #-- this rule needs niftyreg, c3d and mrtrix
+    group:
+        "subj"
+    shell:
+        "c4d {input.nii_4d} -slice w 0:-1 -oo dwi_%03d.nii.gz && "
+        "if [ ! -e 'dwi_001.nii.gz' ]; then "
+        "    mkdir -p {output.affine_dir} &&  "
+        "    cp dwi_000.nii.gz {output.nii_avg3d}; "
+        "else "
+        " parallel --eta --jobs {threads} "
+        " reg_aladin -flo dwi_{{1}}.nii.gz  -ref dwi_000.nii.gz -res warped_{{1}}.nii.gz -aff affine_xfm_ras_{{1}}.txt "
+        " ::: `ls dwi_???.nii.gz | tail -n +2 | grep -Po '(?<=dwi_)[0-9]+'` && "
+        " mkdir -p {output.affine_dir} && cp affine_xfm_ras_*.txt {output.affine_dir} && "
+        " echo -e '1 0 0 0\n0 1 0 0\n0 0 1 0\n0 0 0 1' > {output.affine_dir}/affine_xfm_ras_000.txt && "
+        " mrcat dwi_000.nii.gz warped_*.nii.gz merged_4d.nii.gz && "
+        " mrmath merged_4d.nii.gz mean {output.nii_avg3d} -axis 3;  "
+        "fi"
+
+
+def get_cmd_moco_bzeros_3d(wildcards, input, output, threads):
+
+    if len(input.b0s) == 1:
+        # skip registration
+        return f"cp {input.b0s} {output.nii_avg3d}"
+    else:
+        flo_indices = " ".join([f"{i}" for i in range(1, len(input.b0s))])
+        flo_imgs = " ".join(input.b0s[1:])
+
+        return (
+            f"parallel --eta --jobs {threads} --link "
+            f"reg_aladin -flo {{2}}  -ref {input.b0s[0]} -res warped_{{1}}.nii -aff affine_xfm_ras_{{1}}.txt "
+            f" ::: {flo_indices} ::: {flo_imgs}  && "
+            f" mkdir -p {output.affine_dir} && cp affine_xfm_ras_*.txt {output.affine_dir} && "
+            f" echo -e '1 0 0 0\n0 1 0 0\n0 0 1 0\n0 0 0 1' > {output.affine_dir}/affine_xfm_ras_000.txt && "
+            f" mrcat {input.b0s[0]} warped_*.nii {output.nii_4d} && "
+            f" mrmath {output.nii_4d} mean {output.nii_avg3d} -axis 3"
+        )
+
+
+rule moco_bzeros_3d:
+    """ motion-correct the avgb0 scans from each dwi acquisition"""
+    input:
+        b0s=lambda wildcards: get_dwi_indices(
+            expand(
+                bids(
+                    root=work,
+                    suffix="b0.nii.gz",
+                    datatype="dwi",
+                    desc="moco",
+                    **input_wildcards["dwi"]
+                ),
+                zip,
+                **filter_list(input_zip_lists["dwi"], wildcards)
+            ),
+            wildcards,
+        ),
+    params:
+        cmd=get_cmd_moco_bzeros_3d,
+    output:
+        affine_dir=directory(
+            bids(
+                root=work,
+                suffix="transforms",
+                desc="mocoavgb0",
+                datatype="dwi",
+                **subj_wildcards
+            )
+        ),
+        nii_4d=bids(
+            root=work,
+            suffix="b0s.nii.gz",
+            desc="mocoavgb0",
+            datatype="dwi",
+            **subj_wildcards
+        ),
+        nii_avg3d=bids(
+            root=work,
+            suffix="b0.nii.gz",
+            desc="mocoavgb0",
+            datatype="dwi",
+            **subj_wildcards
+        ),
+    threads: 8  #doesn't need to be more than the number of bzeros 
+    resources:
+        mem_mb=32000,
+    shadow:
+        "minimal"
+    container:
+        config["singularity"]["prepdwi"]  #-- this rule needs niftyreg, c3d and mrtrix
+    group:
+        "subj"
+    shell:
+        "{params.cmd}"

snakedwi/workflow/rules/prepdwi.smk

@@ -101,60 +101,6 @@ rule mrdegibbs:
         "cp {input[3]} {output[3]}"
 
 
-rule moco_bzeros:
-    """ run motion-correction (rigid reg to init volume) on the bzeros """
-    input:
-        nii_4d=bids(
-            root=work,
-            suffix="b0s.nii.gz",
-            datatype="dwi",
-            desc="degibbs",
-            **subj_wildcards
-        ),
-    output:
-        affine_dir=directory(
-            bids(
-                root=work,
-                suffix="transforms",
-                desc="moco",
-                datatype="dwi",
-                **subj_wildcards
-            )
-        ),
-        nii_4d=bids(
-            root=work,
-            suffix="b0s.nii.gz",
-            desc="moco",
-            datatype="dwi",
-            **subj_wildcards
-        ),
-        nii_avg3d=bids(
-            root=work,
-            suffix="b0.nii.gz",
-            desc="moco",
-            datatype="dwi",
-            **subj_wildcards
-        ),
-    threads: 8  #doesn't need to be more than the number of bzeros 
-    resources:
-        mem_mb=32000,
-    shadow:
-        "minimal"
-    container:
-        config["singularity"]["prepdwi"]  #-- this rule needs niftyreg, c3d and mrtrix
-    group:
-        "subj"
-    shell:
-        "c4d {input.nii_4d} -slice w 0:-1 -oo dwi_%03d.nii && "
-        "parallel --eta --jobs {threads} "
-        "reg_aladin -flo dwi_{{1}}.nii  -ref dwi_000.nii -res warped_{{1}}.nii -aff affine_xfm_ras_{{1}}.txt "
-        " ::: `ls dwi_???.nii | tail -n +2 | grep -Po '(?<=dwi_)[0-9]+'` && "
-        " mkdir -p {output.affine_dir} && cp affine_xfm_ras_*.txt {output.affine_dir} && "
-        " echo -e '1 0 0 0\n0 1 0 0\n0 0 1 0\n0 0 0 1' > {output.affine_dir}/affine_xfm_ras_000.txt && "
-        " mrcat dwi_000.nii warped_*.nii {output.nii_4d} && "
-        " mrmath {output.nii_4d} mean {output.nii_avg3d} -axis 3"
-
-
 # now have nii with just the b0's, want to create the topup phase-encoding text files for each one:
 rule get_phase_encode_txt:
     input:

snakedwi/workflow/rules/reg_dwi_to_t1.smk

@@ -105,13 +105,20 @@ rule n4_t1_withmask:
 
 rule reg_dwi_to_t1:
     input:
-        t1w=bids(
+        t1wsynth=bids(
             root=root,
-            suffix="T1w.nii.gz",
+            suffix="T1wSynthSR.nii.gz",
             desc="preproc",
             datatype="anat",
             **subj_wildcards
         ),
+        avgb0synth=bids(
+            root=work,
+            suffix="b0SynthSR.nii.gz",
+            desc="dwiref",
+            datatype="dwi",
+            **subj_wildcards
+        ),
         avgb0=bids(
             root=work,
             suffix="b0.nii.gz",
@@ -121,7 +128,7 @@ rule reg_dwi_to_t1:
         ),
     params:
         general_opts="-d 3",
-        rigid_opts="-m NMI -a -dof 6 -ia-{rigid_dwi_t1_init} -n {rigid_dwi_t1_iters}".format(
+        rigid_opts="-m SSD -a -dof 6 -ia-{rigid_dwi_t1_init} -n {rigid_dwi_t1_iters}".format(
             rigid_dwi_t1_init=config["rigid_dwi_t1_init"],
             rigid_dwi_t1_iters=config["rigid_dwi_t1_iters"],
         ),
@@ -151,8 +158,8 @@ rule reg_dwi_to_t1:
         bids(root="logs", suffix="reg_b0_to_t1.txt", datatype="dwi", **subj_wildcards),
     threads: 8
     shell:
-        "greedy -threads {threads} {params.general_opts} {params.rigid_opts}  -i {input.t1w} {input.avgb0} -o {output.xfm_ras}  &> {log}  && "
-        "greedy -threads {threads} {params.general_opts} -rf {input.t1w} -rm {input.avgb0} {output.warped_avgb0}  -r {output.xfm_ras} &>> {log}"
+        "greedy -threads {threads} {params.general_opts} {params.rigid_opts}  -i {input.t1wsynth} {input.avgb0synth} -o {output.xfm_ras}  &> {log}  && "
+        "greedy -threads {threads} {params.general_opts} -rf {input.t1wsynth} -rm {input.avgb0} {output.warped_avgb0}  -r {output.xfm_ras} &>> {log}"
 
 
 rule qc_reg_dwi_t1: