Temporal averaging¶

The next step is to average together, over time, all the motion-corrected functional images for a given participant.

Outputs¶

We want to produce a single NIFTI file per participant that has a single volume — the temporal average of all the motion-corrected data. We will give the desc entity the label tmean.

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

workflow/rules/tmean.smk¶

rule tmean:
    "Averages the motion-corrected images over time"
    output:
        img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_bold.nii.gz",

Inputs¶

Each output file requires two input files — one motion-corrected image for each of the two tasks. The mot_correct rule was applied separately to participants and tasks, so we cannot simply refer to the output of the rule. Instead, we can use the collect helper function (from Snakemake) to produce a sequence of inputs based on a rule’s output.

workflow/rules/tmean.smk¶

rule tmean:
    "Averages the motion-corrected images over time"
    input:
        imgs=lambda wildcards: collect(
            rules.mot_correct.output.img,
            sub_num=wildcards.sub_num,
            task=TASKS,
        ),
    output:
        img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_bold.nii.gz",

We need to use an input function here as we want to use the active subject number with the fixed set of tasks. Note that imgs here in the input is now a sequence rather than just a single string as it has been previously.

Parameters¶

We don’t need any particular parameters for this rule, so we can skip the params directive.

Mechanism¶

Container¶

We will use the same AFNI container that we used in the motion correction rule:

workflow/rules/tmean.smk¶

rule tmean:
    "Averages the motion-corrected images over time"
    input:
        imgs=lambda wildcards: collect(
            rules.mot_correct.output.img,
            sub_num=wildcards.sub_num,
            task=TASKS,
        ),
    output:
        img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_bold.nii.gz",
    container:
        CONTAINER_SOURCES["AFNI"]

Logging¶

As usual, we need to specify the file to store logging information:

workflow/rules/tmean.smk¶

rule tmean:
    "Averages the motion-corrected images over time"
    input:
        imgs=lambda wildcards: collect(
            rules.mot_correct.output.img,
            sub_num=wildcards.sub_num,
            task=TASKS,
        ),
    output:
        img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_bold.nii.gz",
    container:
        CONTAINER_SOURCES["AFNI"]
    log:
        "logs/tmean/tmean_{sub_num}.txt"

Shell command¶

Similarly to the first rule, we will directly specify shell commands to run rather than using a script. First, we use 3dTcat to combine the separate files into one file (tcat.nii.gz). Then, we use 3dTstat to average that file and produce our desired output.

workflow/rules/tmean.smk¶

rule tmean:
    "Averages the motion-corrected images over time"
    input:
        imgs=lambda wildcards: collect(
            rules.mot_correct.output.img,
            sub_num=wildcards.sub_num,
            task=TASKS,
        ),
    output:
        img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_bold.nii.gz",
    container:
        CONTAINER_SOURCES["AFNI"]
    log:
        "logs/tmean/tmean_{sub_num}.txt"
    shell:
        """
# concatenate the input images over time into one output image
3dTcat -prefix tcat.nii.gz {input.imgs} > {log} 2>&1
# average the concatenated image over time
3dTstat -prefix {output.img} -mean tcat.nii.gz >> {log} 2>&1
        """

Note

We use the >> operator, rather than >, to redirect to the log file in the second command. This is because the >> appends to an existing file, whereas > would overwrite the output from the first command.

Shadowing¶

There are two important to things to note about the above command that would cause problems if executed within the current state of the rule:

The intermediate file tcat.nii.gz would remain present after the job has completed.
Each job (i.e., across participants) would be writing to the same tcat.nii.gz file.

We can resolve both problems by using shadowing. When rules are shadowed, their jobs are executed in an isolated temporary directory and only the listed outputs are copied to their correct destination after the job completes. For our purposes here, that means that any intermediate files are discarded (resolving the first issue listed above) and that each job gets its own directory (resolving the second issue listed above). There are a few options for how Snakemake implements the shadowing; here we are using shallow, but see the Snakemake documentation for other options.

workflow/rules/tmean.smk¶

rule tmean:
    "Averages the motion-corrected images over time"
    input:
        imgs=lambda wildcards: collect(
            rules.mot_correct.output.img,
            sub_num=wildcards.sub_num,
            task=TASKS,
        ),
    output:
        img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_bold.nii.gz",
    container:
        CONTAINER_SOURCES["AFNI"]
    log:
        "logs/tmean/tmean_{sub_num}.txt"
    shadow:
        "shallow"
    shell:
        """
# concatenate the input images over time into one output image
3dTcat -prefix tcat.nii.gz {input.imgs} > {log} 2>&1
# average the concatenated image over time
3dTstat -prefix {output.img} -mean tcat.nii.gz >> {log} 2>&1
        """

Resources¶

We will again specify that Snakemake should budget for the job using up to 1 GB of RAM.

workflow/rules/tmean.smk¶

rule tmean:
    "Averages the motion-corrected images over time"
    input:
        imgs=lambda wildcards: collect(
            rules.mot_correct.output.img,
            sub_num=wildcards.sub_num,
            task=TASKS,
        ),
    output:
        img="results/sub-{sub_num}/func/sub-{sub_num}_desc-tmean_bold.nii.gz",
    container:
        CONTAINER_SOURCES["AFNI"]
    log:
        "logs/tmean/tmean_{sub_num}.txt"
    resources:
        mem="1GB",
    shadow:
        "shallow"
    shell:
        """
# concatenate the input images over time into one output image
3dTcat -prefix tcat.nii.gz {input.imgs} > {log} 2>&1
# average the concatenated image over time
3dTstat -prefix {output.img} -mean tcat.nii.gz >> {log} 2>&1
        """

Preparing for execution¶

As usual, our next step is to add the new rule file to the Snakefile:

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"

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

rule all:
    input:
        expand(rules.acquire_anat.output.img, sub_num=SUB_NUMS),
        expand(rules.mot_correct.output.img, sub_num=SUB_NUMS, task=TASKS),

The output from this rule is now one of our ‘highest’ level of output, in the sense of not being required for any other rules, so we need to update 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"

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

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

Executing the workflow¶

Finally, you can run Snakemake and execute the workflow:

$ uv run snakemake