Snakemake可选输出，带有扩展()

Question

我是第一次接触Snakemake，我想知道在使用expand()时是否可以将可选的输出文件放入snakemake规则中。

我正在使用bowtie2-build为我的参考基因组创建索引，但根据基因组大小的不同，bowtie2会创建具有不同扩展名的索引文件：.bt2用于小基因组，.bt21用于大基因组。

我有以下规则：

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}"

但现在，snakemake将始终查找所有扩展名的输出，这将在运行时给出错误，因为只有具有两个扩展名.bt2或.bt21之一的文件才会根据基因组大小实际创建。

我试着像这样使用正则表达式：

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

这是可行的，但我觉得应该有一种更简单的方法。

Answer 1

也许你把事情搞得太复杂了。Bowtie2 (校准器)接受输入的索引文件的前缀，它将自己找到实际的索引文件。所以我不会将索引文件作为任何规则的输出列出。只需使用一个标志文件来指示索引已完成。例如：

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}
        """

Answer 2

这种方法是用不同的扩展名来重命名文件。如果这样的文件不存在，它将打印一个错误，但您可能认为这是一个功能...：

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
        """

在其他示例中，在输入文件具有关于输出的信息的情况下，例如，如果基因组序列大小从一开始就知道，则一种选择是生成将指定不同扩展的两个规则bowtie2_build和bowtie2_build_big。

Answer 3

看起来像是checkpoints的用途。使用检查点时，将在执行检查点后重新评估DAG。你可以这样做：

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"