1

我正在编写一个管道，其中参考基因组上的读取的第一个映射导致一个 sam 文件（规则map_on_genome），然后是一个排序和索引的 bam 文件（规则sam2indexedbam）。处理完这个 bam 文件（规则small_RNA_seq_annotate）后，我提取了一些我想用相同参数重新映射到同一个基因组上的读数。我还从第一次映射（规则）期间未映射的读数中修剪额外的核苷酸extract_nomap_siRNAs，并且我想尝试以类似方式再次映射生成的修剪读数。

理想情况下，rules map_on_genomeand因此将在这些规则的输出sam2indexedbam之前small_RNA_seq_annotate和之后再次使用。extract_nomap_siRNAs

在我尝试这样做时，我命名的通配符read_type应该可以帮助 snakemake 区分第一个映射和另一个映射，但显然我失败了，因为我最终得到以下错误：

Cyclic dependency on rule sam2indexedbam.

实现我想要的目标的推荐方法是什么（即避免代码重复）？

例子

我试图在我的真实管道的简化版本中重现该错误。这是一个导致CyclicGraphException：

rule all:
    input:
        expand("results/{read_type}_on_genome_sorted.bam", read_type=["reads", "nomap_siRNA", "piRNA", "siRNA", "miRNA"]),


def source_fastq(wildcards):
    if wildcards.read_type == "nomap_siRNA":
        return rules.extract_nomap_siRNA.output.fastq
    elif wildcards.read_type == "piRNA":
        return rules.small_RNA_seq_annotate.output.pi_fastq
    elif wildcards.read_type == "siRNA":
        return rules.small_RNA_seq_annotate.output.si_fastq
    elif wildcards.read_type == "miRNA":
        return rules.small_RNA_seq_annotate.output.mi_fastq
    elif wildcards.read_type == "reads":
        return "data/reads.fastq"
    else:
        raise NotImplementedError("Unknown read type: %s" % wildcards.read_type)


rule map_on_genome:
    input:
        source_fastq
    output:
        sam = "results/{read_type}_on_genome.sam",
        nomap = "results/{read_type}_nomap.fastq"
    shell:
        """
        touch {output}
        """

rule extract_nomap_siRNA:
    input:
        rules.map_on_genome.output.sam
        #"results/reads_nomap.fastq"
    output:
        fastq = "results/{read_type}_unmapped_siRNA.fastq"
    shell:
        """
        touch {output}
        """

rule sam2indexedbam:
    input:
        rules.map_on_genome.output.sam
    output:
        "results/{read_type}_on_genome_sorted.bam"
    shell:
        """
        touch {output}
        """

rule small_RNA_seq_annotate:
    input:
        "results/reads_on_genome_sorted.bam"
    output:
        pi_fastq = "results/piRNA.fastq",
        si_fastq = "results/siRNA.fastq",
        mi_fastq = "results/miRNA.fastq"
    shell:
        """
        touch {output.pi_fastq}
        touch {output.si_fastq}
        touch {output.mi_fastq}
        """

解决方案

如果我将规则的输入更改为extract_nomap_siRNA硬编码文件路径而不是map_on_genome规则的输出，错误就会消失。

我想这种硬编码打破了循环依赖。

4

0 回答 0