Snakemake optional output with expand()

user_17078546 0
Viewed 315

I am new to Snakemake, and I'm wondering if I'm able to put optional output files in a snakemake rule while using expand().

I'm using bowtie2-build to create an indexing of my reference genome, but depending on the genome size, bowtie2 creates indexing files with different extensions: .bt2 for small genomes, and .bt21 for big genomes.

I have the following rule:

rule bowtie2_build:
    """
    """
    input:
        "reference/"+config["reference_genome"]+".fa"
    output:
                                                                   # every possible extension in expand
        expand("reference/"+config["reference_genome"]+"{suffix}", suffix=[".1.bt2", ".2.bt2", ".3.bt2", ".4.bt2", ".rev.1.bt2", ".rev.2.bt2", ".1.bt21", ".2.bt21", ".3.bt21", ".4.bt21", ".rev.1.bt21", ".rev.2.bt21"])
    params:
        output_prefix=config["reference_genome"]
    shell:
        "bowtie2-build {input} reference/{params.output_prefix}"

But now, snakemake will always look for output with all extensions, which will give an error while running, because only files with either one of the two extensions .bt2 or .bt21 are actually created depending on genome size.

I have tried to use regex like so:

output:
        "reference/"+config["reference_genome"]+"{suffix, \.(1|2|3|4|rev)\.(bt2|bt21|1|2)\.?(bt2|bt21)?}"

And this works, but I feel like there should be an easier way for it.

snakemake bioinformatics python
3 Answers

Perhaps you are making things more complicated than necessary. Bowtie2 (the aligner) takes in input the prefix of the index files and it will find the actual index files by itself. So I wouldn't list the index files as output of any rule. Just use a flag file to indicate that indexing has been completed. For example:

rule all:
    input:
        expand('{sample}.sam', sample= ['A', 'B', 'C']),

rule bowtie_build:
    input:
        fa= 'genome.fa',
    output:
        touch('index.done'), # Flag file
    params:
        index_prefix= 'genome',
    shell:
        r"""
        bowtie2-build {input.fa} {params.index_prefix}
        """

rule bowtie_align:
    input:
        idx= 'index.done',
        fq1= '{sample}.R1.fastq.gz',
    output:
        sam= '{sample}.sam',
    params:
        index_prefix= 'genome'
    shell:
        r"""
        bowtie2 -x {params.index_prefix} -U {input.fq1} -S {output.sam}
        """
user_11144530

This approach is to rename the files with a different extension. It will print an error if such file does not exist, but you might consider this a feature...:

rule bowtie2_build:
    input:
        "reference/"+config["reference_genome"]+".fa"
    output:
        expand("reference/"+config["reference_genome"]+"{suffix}", suffix=[".1.bt2", ".2.bt2", ".3.bt2", ".4.bt2", ".rev.1.bt2", ".rev.2.bt2"])
    params:
        output_prefix=config["reference_genome"],
        file_name = "reference/"+config["reference_genome"],
    shell:
        """
        bowtie2-build {input} reference/{params.output_prefix}
        mv -v {params.file_name}.1.bt21 {params.file_name}.1.bt2
        mv -v {params.file_name}.2.bt21 {params.file_name}.2.bt2
        mv -v {params.file_name}.3.bt21 {params.file_name}.3.bt2
        # etc... it's also possible to put this in a bash loop
        """

In other examples, where the input file has information about the output, e.g. if the genome sequence size is known from the outset, then one option is to generate two rules, bowtie2_build and bowtie2_build_big, that will specify different extensions.

user_106935960

Seems like a use for checkpoints. With checkpoints, the DAG will be reevaluated after the checkpoint's execution. You could do something like this:

from glob import glob

def get_ref_index(wildcards):
    """Returns the reference index for wildcards"""
    checkpoints.bowtie2_build.get()  # Checks to see if the checkpoint rule has been run, might be an issue without wildcards
    suffix=[".1.bt2", ".2.bt2", ".3.bt2", ".4.bt2", ".rev.1.bt2", ".rev.2.bt2", ".1.bt21", ".2.bt21", ".3.bt21", ".4.bt21", ".rev.1.bt21", ".rev.2.bt21"]
    ref_files = glob.glob("reference/" + config["reference_genome"] + "*") # List files with reference genome pattern
    
    for ref in ref_files:
        suffix = ref.split(".", 1)[1]  # Split on first period to grab all suffixes
        if suffix in suffixes:
            return ref
    
checkpoint bowtie2_build:
    input:
        "reference/"+config["reference_genome"]+".fa"
    output:
        touch("reference/"+config["reference_genome"]+".done")  # Breadcrumb file to link to downstream rules, since we don't know what indexes we'll get
    params:
        output_prefix=config["reference_genome"]
    shell:
        "bowtie2-build {input} reference/{params.output_prefix}"
        
rule donwstream_rule:
    input:
        ref = "reference/"+config["reference_genome"]+".fa"
        ref_index = get_ref_index
        ref_index_done = "reference/"+config["reference_genome"]+".done"
    output:
        foo = "bar"
    shell:
        "touch bar"
user_142123400
Related