1
$\begingroup$

I am trying to generate several snakemake rules. However, I would like to include some restrictions due to storage limitations

I first would like to copy the the bamfiles, 1 for each sample, to a staging area. However, while I have over 50 samples, I can only process three samples at a time.

Then I want to compare the bam file to chunked candidate files, 1 for each chromosome. This however, I would like to run on 23 cores per sample so that once each bam file is copied 69 cores are working.

finally, if all chunks for a given sample have been processed I want to clean up the staging area, copy the bamfile for the next sample and process the chunks for that sample.

My current implementation, see below, just continues copying bamfiles and crashes when I run out of strorage;

rule copy_bam:
    input:
        bamfile=os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam"),
        bamindex=os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam.bai")
    output:
        staging_bamfile="{staging_area}/{sample}.bam",
        staging_bamindex="{staging_area}/{sample}.bam.bai"
    wildcard_constraints:
        sample="[A-Za-z0-9_]+"
    resources:
        bam_copy=1
    shell:
        """
        cp {input.bamfile} {output.staging_bamfile}
        cp {input.bamindex} {output.staging_bamindex}
        """

rule process_chunks:
    input:
        chunkfile="{staging_area}/chunks/{chunk}_with_header.bcf.gz",
        bamfile=rules.copy_bam.output.staging_bamfile
    output:
        dynamic("{staging_area}/calls/{sample}/{chunk}.bcf")

    threads: 23
    shell:
        """
        process chunks
        """

rule cleanup:
    input:
        bamfile="{staging_area}/{sample}.bam",
        bamindex="{staging_area}/{sample}.bam.bai",
        processed=lambda wildcards: expand("{staging_area}/calls/{sample}/{chunk}.bcf", chunk=get_chunks()[wildcards.sample], staging_area=config['staging_area'], sample=wildcards.sample)
    output:
        touch("{staging_area}/{sample}/cleanup_done.txt")
    shell:
        """
        rm {input.bamfile}
        rm {input.bamindex}
        rm {input.processed}
        """

run with: snakemake --resources bam_copy=3 --cores 70

$\endgroup$
1
  • $\begingroup$ One thing to note about any approach here that uses Snakemake's resources is that they're defined per running job, versus with the "in use until something happens in a later job" idea here. Jame's Hawley's suggestion of temp() should help, just note that I don't think the resource tracking even matters in that case (the resource is only "in use" while copy_bam is actively running.) $\endgroup$
    – Jesse
    May 25, 2023 at 16:16

1 Answer 1

2
$\begingroup$

One Snakemake feature you can use is temporary outputs.

If you make your outputs from the copy_bam rule temporary, they'll be cleared up after no other rule requires that file as an input. Your Snakefile would look something like this:

rule copy_bam:
    input:
        bam = os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam"),
        idx = os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam.bai"),
    output:
        bam = temp("{staging_area}/{sample}.bam"),
        idx = temp("{staging_area}/{sample}.bam.bai"),
    wildcard_constraints:
        sample="[A-Za-z0-9_]+"
    resources:
        bam_copy=1
    shell:
        """
        cp {input.bam} {output.bam}
        cp {input.idx} {output.idx}
        """
# omitted for brevity

You can then specify how many jobs run concurrently with the --cores/--jobs options to keep the space usage small.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.