使用snakemake touch使用bwa和gatk创建参考基因组的索引

如何解决使用snakemake touch使用bwa和gatk创建参考基因组的索引

我将读取与bwa对齐，并将变体与gatk调用。 gatk需要为参考基因组创建dict，而bwa需要创建索引。当我对他们两个都使用触摸时，会出现此错误：

AmbiguousRuleException:
Rules bwa_index and gatk_refdict are ambiguous for the file ref.
Expected input files:
        bwa_index: ref.fasta
        gatk_refdict: ref.fasta

这是代码：

rule bwa_index:
    input:
        database="ref.fasta"
    output:
        done =touch("ref")
    shell:
        """
        bwa index -p ref {input.database}
        """
rule bwa_mem:
    input:
        bwa_index_done = "ref",fastq1="{sample}_R1.trimmed.fastq.gz",fastq2="{sample}_R2.trimmed.fastq.gz"
    output:
        bam = temp("{sample}.bam")
    shell:
        """
        bwa mem ref {input.fastq1} {input.fastq2} -o {output.bam}
        """
rule_gatk_refdict:
    input:
        ref="ref.fasta"
    output:
        done =touch("ref")
    shell:
        """
        java -jar gatk-package-4.1.9.0-local.jar CreateSequenceDictionary -R {input.ref} -O {output.done}
        """
rule gatk:
    input:
        gatk_refdict_done = "ref",bam="bam_list"
    output:
        outf ="{chr}.vcf"
    shell:
        """
        java -jar gatk-package-4.1.9.0-local.jar HaplotypeCaller -L {wildcards.chr} -R ref -I {input.bam} --min-base-quality-score 20 -O {output.outf}
        """

或者，我指定了索引.dict，但是它也不起作用，因为gatk在创建dict之前就调用了变体，因此我得到一个错误，即没有dict文件：

rule_gatk_refdict:
    input:
        ref="ref.fasta"
    output:
        outf ="ref.dict"
    shell:
        """
        java -jar gatk-package-4.1.9.0-local.jar CreateSequenceDictionary -R {input.ref} -O {output.outf}
        """
rule gatk:
    input:
        ref = "ref.fasta",bam="bam_list"
    output:
        outf ="{chr}.vcf"
    shell:
        """
        java -jar gatk-package-4.1.9.0-local.jar HaplotypeCaller -L {wildcards.chr} -R {input.ref} -I {input.bam} --min-base-quality-score 20 -O {output.outf}
        """

如何解决这个问题？

解决方法

为什么不简单地将dict文件定义为gatk规则的输入，而将“ ndex”定义为bwa规则的输入呢？

rule bwa_index:
    input:
        database="ref.fasta"
    output:
        done =touch("ref")
    shell:
        """
        bwa index -p ref {input.database}
        """
rule bwa_mem:
    input:
        bwa_index_done = "ref",fastq1="{sample}_R1.trimmed.fastq.gz",fastq2="{sample}_R2.trimmed.fastq.gz"
    output:
        bam = temp("{sample}.bam")
    shell:
        """
        bwa mem ref {input.fastq1} {input.fastq2} -o {output.bam}
        """
rule gatk_refdict:
    input:
        ref="ref.fasta"
    output:
        done = "ref.dict"
    shell:
        """
        java -jar gatk-package-4.1.9.0-local.jar CreateSequenceDictionary -R {input.ref} -O {output.done}
        """
rule gatk:
    input:
        ref = "ref.fasta",dict = "ref.dict",bam="bam_list"
    output:
        outf ="{chr}.vcf"
    shell:
        """
        java -jar gatk-package-4.1.9.0-local.jar HaplotypeCaller -L {wildcards.chr} -R {input.ref} -I {input.bam} --min-base-quality-score 20 -O {output.outf}
        """

您得到的AmbiguousRuleException是因为snakemake不知道要运行哪个规则，因为两个规则具有相同的输出。不要忘了snakemake会尝试从规则开始构建DAG。在运行rule gatk时，您将"ref"定义为输入。由于可以使用两个规则来生成此文件，因此snakemake不知道必须使用rule gatk_refdict还是rule bwa_index。

您在这里有错字（“ _”不应在其中）：

----v
rule_gatk_refdict:
    input:
        ref="ref.fasta"
    ...