带有 expand() 的 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

此方法是将文件重命名为不同的扩展名。如果这样的文件不存在，它会打印一个错误，但你可能认为这是一个特性...:[=13=]

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 2

似乎是 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"

Answer 3

也许你把事情弄得比必要的更复杂了。 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}
        """

带有 expand() 的 Snakemake 可选输出

Snakemake optional output with expand()

python

bioinformatics

snakemake