snakemake

Question

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

前两个步骤（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}"

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

有没有办法成功做到这一点？困难在于我没有定义名称列表或其他任何内容，我想。

预先感谢您的帮助。

Answer 1

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

例如：

# 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: ...

希望这对您有所帮助。

Answer 2

您可以尝试在初始 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)]

snakemake - 在一个规则中从多个输入文件输出一个文件

snakemake - output one only file from multiple input files in one rule