6

我第一次使用snakemake，以便使用cutadapt、bwa 和GATK（修剪；映射；调用）构建基本管道。我想在目录中包含的每个 fastq 文件上运行这个管道，而不必在蛇文件或配置文件中指定它们的名称或任何内容。我想成功地做到这一点。

前两个步骤（cutadapt 和 bwa / 修剪和映射）运行良好，但我遇到了 GATK 的一些问题。

首先，我必须从 bam 文件生成 g.vcf 文件。我正在使用以下规则执行此操作：

configfile: "config.yaml"

import os
import glob

rule all:
    input:
        "merge_calling.g.vcf"

rule cutadapt:
    input:
        read="data/Raw_reads/{sample}_R1_{run}.fastq.gz", 
        read2="data/Raw_reads/{sample}_R2_{run}.fastq.gz" 
    output:
        R1=temp("trimmed_reads/{sample}_R1_{run}.fastq.gz"),
        R2=temp("trimmed_reads/{sample}_R2_{run}.fastq.gz") 
    threads:
        10
    shell:
        "cutadapt -q {config[Cutadapt][Quality_value]} -m {config[Cutadapt][min_length]} -a {config[Cutadapt][forward_adapter]} -A  {config[Cutadapt][reverse_adapter]} -o {output.R1} -p '{output.R2}' {input.read} {input.read2}"

rule bwa_map:
    input:
        genome="data/genome.fasta",
        read=expand("trimmed_reads/{{sample}}_{pair}_{{run}}.fastq.gz", pair=["R1", "R2"]) 
    output:
        temp("mapped_bam/{sample}_{run}.bam")
    threads:
        10
    params:
        rg="@RG\\tID:{sample}\\tPL:ILLUMINA\\tSM:{sample}"
    shell:
        "bwa mem -t 2 -R '{params.rg}' {input.genome} {input.read} | samtools view -Sb - > {output}"

rule picard_sort:
    input:
        "mapped_bam/{sample}.bam"
    output:
        "sorted_reads/{sample}.bam"
    shell:
        "java -Xmx4g -jar /home/alexandre/picard-tools/picard.jar SortSam I={input} O={output} SO=coordinate VALIDATION_STRINGENCY=SILENT"

rule picard_rmdup:
    input:
        bam="sorted_reads/{sample}.bam"
    output:
        "rmduped_reads/{sample}.bam",
        "picard_stats/{sample}.bam"
    params:
        reads="rmduped_reads/{sample}.bam",
        stats="picard_stats/{sample}.bam",
    shell:
        "java -jar -Xmx2g /home/alexandre/picard-tools/picard.jar MarkDuplicates "
        "I={input.bam} "
        "O='{params.reads}' "
        "VALIDATION_STRINGENCY=SILENT "
        "MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000 "
        "REMOVE_DUPLICATES=TRUE "
        "M='{params.stats}'"

rule samtools_index:
    input:
        "rmduped_reads/{sample}.bam"
    output:
        "rmduped_reads/{sample}.bam.bai"
    shell:
        "samtools index {input}"

rule GATK_raw_calling:
    input:
        bam="rmduped_reads/{sample}.bam",
        bai="rmduped_reads/{sample}.bam.bai",
        genome="data/genome.fasta"
    output:
        "Raw_calling/{sample}.g.vcf",
    shell:
        "java -Xmx4g -jar /home/alexandre/GenomeAnalysisTK-3.7/GenomeAnalysisTK.jar -ploidy 2 --emitRefConfidence GVCF -T HaplotypeCaller -R {input.genome} -I {input.bam} --genotyping_mode DISCOVERY -o {output}"

这些规则运作良好。例如，如果我有文件： Cla001d_S281_L001_R1_001.fastq.gz Cla001d_S281_L001_R2_001.fastq.gz

我可以创建一个 bam 文件 ( Cla001d_S281_L001_001.bam)，然后从该 bam 文件创建一个 GVCF 文件 ( Cla001d_S281_L001_001.g.vcf)。我有很多这样的样本，我需要为每个样本创建一个 GVCF 文件，然后将这些 GVCF 文件合并到一个文件中。问题是我无法提供要合并到以下规则的文件列表：

rule GATK_merge:
    input:
        ???
    output:
        "merge_calling.g.vcf"
    shell:
        "java -Xmx4g -jar /home/alexandre/GenomeAnalysisTK-3.7/GenomeAnalysisTK.jar "
        "-T CombineGVCFs "
        "-R data/genome.fasta "
        "--variant {input} "
        "-o {output}"

为了做到这一点，我尝试了几件事，但无法成功。问题是两个规则之间的链接（GATK_raw_calling并且GATK_merge应该合并 的输出GATK_raw_calling）。如果我将输出指定GATK_raw_calling为以下规则的输入，我无法输出一个文件（输入文件中的通配符无法从输出文件中确定），并且如果我无法在这两个规则之间建立链接我没有将这些文件指定为输入...

有没有办法成功地做到这一点？我认为，困难在于我没有定义名称列表或其他任何东西。

提前感谢您的帮助。

4

2 回答 2

4

您可以尝试在初始 fastq.gz 文件中使用glob_wildcards生成示例 ID 列表：

sample_ids, run_ids = glob_wildcards("data/Raw_reads/{sample}_R1_{run}.fastq.gz")

然后，您可以将其用于expand以下输入GATK_merge：

rule GATK_merge:
    input:
        expand("Raw_calling/{sample}_{run}.g.vcf",
               sample=sample_ids, run=run_ids)

如果相同的运行 ID 总是带有相同的样本 ID，您将需要压缩而不是展开，以避免不存在的组合：

rule GATK_merge:
    input:
        ["Raw_calling/{sample}_{run}.g.vcf".format(
            sample=sample_id,
            run=run_id) for sample_id, run_id in zip(sample_ids, run_ids)]
于 2017-07-06T09:54:13.440 回答
2

您可以通过使用 python 函数作为规则的输入来实现此目的，如此处的蛇形文档中所述。

例如，可能看起来像这样：

# Define input files
def gatk_inputs(wildcards):
    files = expand("Raw_calling/{sample}.g.vcf", sample=<samples list>)
    return files

# Rule
rule gatk:
    input: gatk_inputs
    output: <output file name>
    run: ...

希望这可以帮助。

于 2017-07-04T15:48:21.403 回答