Coregistration

The next step is to coregister the anatomical and functional images.

Outputs

We want to produce a single NIFTI file per participant that is their anatomical image in alignment with the mean functional image. We will give the desc entity the label coreg.

We start by creating a new rule (workflow/rules/coreg.smk) that has this output information:

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",

For subsequent visualisation, it is also useful to have a copy of the mean functional image after it has been resampled into the grid of the coregistered anatomical image. We add that file as an additional output of the rule:

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",
        func_img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_grid-anat_bold.nii.gz",

Inputs

The coregistration process requires two inputs per participant: their anatomical image and their mean functional image. We enter these as separate items, with their paths given by the output from the relevant rule.

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    input:
        anat=rules.acquire_anat.output.img,
        func=rules.tmean.output.img,
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",
        func_img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_grid-anat_bold.nii.gz",

Parameters

Because we will want to rename an intermediate output during execution of this rule, it is helpful for us to be able to access a path’s ‘stem’ (i.e., the filename without its suffix). We can use the Snakemake helper function subpath to do that; the strip_suffix option allows us to remove the .nii.gz from the path and the basename=True removes the directory information.

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    input:
        anat=rules.acquire_anat.output.img,
        func=rules.tmean.output.img,
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",
        func_img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_grid-anat_bold.nii.gz",
    params:
        anat_stem=subpath(input.anat, strip_suffix=".nii.gz", basename=True),

Mechanism

Container

We will use the same AFNI container that we have used in previous rules:

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    input:
        anat=rules.acquire_anat.output.img,
        func=rules.tmean.output.img,
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",
        func_img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_grid-anat_bold.nii.gz",
    params:
        anat_stem=subpath(input.anat, strip_suffix=".nii.gz", basename=True),
    container:
        CONTAINER_SOURCES["AFNI"]

Logging

As usual, we need to specify the file to store logging information (particularly useful for this rule, given the importance of the coregistration process summary that is printed by AFNI):

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    input:
        anat=rules.acquire_anat.output.img,
        func=rules.tmean.output.img,
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",
        func_img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_grid-anat_bold.nii.gz",
    params:
        anat_stem=subpath(input.anat, strip_suffix=".nii.gz", basename=True),
    container:
        CONTAINER_SOURCES["AFNI"]
    log:
        "logs/coreg/coreg_{sub_num}.txt"

Shell command

We will use the align_epi_anat.py AFNI command to run the coregistration and resample the mean functional image:

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    input:
        anat=rules.acquire_anat.output.img,
        func=rules.tmean.output.img,
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",
        func_img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_grid-anat_bold.nii.gz",
    params:
        anat_stem=subpath(input.anat, strip_suffix=".nii.gz", basename=True),
    container:
        CONTAINER_SOURCES["AFNI"]
    log:
        "logs/coreg/coreg_{sub_num}.txt"
    shell:
        """
# run the coregistration
align_epi_anat.py \
-anat {input.anat} -epi {input.func} \
-epi_base 0 -giant_move \
> {log} 2>&1
# convert the resulting files into NIFTI format
3dcopy \
{params.anat_stem}_al+orig \
{output.anat_img} \
>> {log} 2>&1
# save a copy of the functional in the anatomical grid
3dresample \
-inset {input.func} \
-master {output.anat_img} \
-prefix {output.func_img} \
>> {log} 2>&1
        """

Shadowing

Because this command produces quite a few intermediate outputs, we will again use rule shadowing:

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    input:
        anat=rules.acquire_anat.output.img,
        func=rules.tmean.output.img,
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",
        func_img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_grid-anat_bold.nii.gz",
    params:
        anat_stem=subpath(input.anat, strip_suffix=".nii.gz", basename=True),
    container:
        CONTAINER_SOURCES["AFNI"]
    log:
        "logs/coreg/coreg_{sub_num}.txt"
    shadow:
        "shallow"
    shell:
        """
# run the coregistration
align_epi_anat.py \
-anat {input.anat} -epi {input.func} \
-epi_base 0 -giant_move \
> {log} 2>&1
# convert the resulting files into NIFTI format
3dcopy \
{params.anat_stem}_al+orig \
{output.anat_img} \
>> {log} 2>&1
# save a copy of the functional in the anatomical grid
3dresample \
-inset {input.func} \
-master {output.anat_img} \
-prefix {output.func_img} \
>> {log} 2>&1
        """

Resources

We will give Snakemake a bit more RAM allowance for this job:

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    input:
        anat=rules.acquire_anat.output.img,
        func=rules.tmean.output.img,
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",
        func_img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_grid-anat_bold.nii.gz",
    params:
        anat_stem=subpath(input.anat, strip_suffix=".nii.gz", basename=True),
    container:
        CONTAINER_SOURCES["AFNI"]
    log:
        "logs/coreg/coreg_{sub_num}.txt"
    shadow:
        "shallow"
    resources:
        mem="4GB",
    shell:
        """
# run the coregistration
align_epi_anat.py \
-anat {input.anat} -epi {input.func} \
-epi_base 0 -giant_move \
> {log} 2>&1
# convert the resulting files into NIFTI format
3dcopy \
{params.anat_stem}_al+orig \
{output.anat_img} \
>> {log} 2>&1
# save a copy of the functional in the anatomical grid
3dresample \
-inset {input.func} \
-master {output.anat_img} \
-prefix {output.func_img} \
>> {log} 2>&1
        """

We will also allow it to use two cores when executing each job using the threads directive:

workflow/rules/coreg.smk

rule coreg:
    "Coregisters the anatomical to the mean functional"
    input:
        anat=rules.acquire_anat.output.img,
        func=rules.tmean.output.img,
    output:
        anat_img="results/sub-{sub_num}/anat/sub-{sub_num}_desc-coreg_T1w.nii.gz",
        func_img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_grid-anat_bold.nii.gz",
    params:
        anat_stem=subpath(input.anat, strip_suffix=".nii.gz", basename=True),
    container:
        CONTAINER_SOURCES["AFNI"]
    log:
        "logs/coreg/coreg_{sub_num}.txt"
    shadow:
        "shallow"
    resources:
        mem="4GB",
    threads: 2
    shell:
        """
# run the coregistration
align_epi_anat.py \
-anat {input.anat} -epi {input.func} \
-epi_base 0 -giant_move \
> {log} 2>&1
# convert the resulting files into NIFTI format
3dcopy \
{params.anat_stem}_al+orig \
{output.anat_img} \
>> {log} 2>&1
# save a copy of the functional in the anatomical grid
3dresample \
-inset {input.func} \
-master {output.anat_img} \
-prefix {output.func_img} \
>> {log} 2>&1
        """

Preparing for execution

As usual, our next step is to add the new rule file to the Snakefile and adjust the output of the all rule:

workflow/Snakefile

include: "rules/common.smk"
include: "rules/acquire_anat.smk"
include: "rules/acquire_func.smk"
include: "rules/mot_correct.smk"
include: "rules/tmean.smk"
include: "rules/coreg.smk"

wildcard_constraints:
    task="|".join(TASKS)

rule all:
    input:
        expand(rules.coreg.output.img, sub_num=SUB_NUMS),

Note that we only need to expand this single rule now, given it depends on output from all of the other rules.

Executing the workflow

Finally, you can run Snakemake and execute the workflow:

$ uv run snakemake