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

Question

我第一次使用snakemake来构建一个基本的pipeline，使用cutadapt、bwa和GATK（修剪；比对；调用）。我想在目录中运行此pipeline以处理每个fastq文件，而不必在snakefile或config文件中指定它们的名称或其他信息。我希望能够成功地做到这一点。

前两个步骤（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的输出指定为下一个规则的输入，则无法输出单个文件（无法从输出文件确定输入文件中的通配符），如果我不将这些文件指定为输入，则无法在两个规则之间建立链接...

有没有办法成功做到这一点？困难在于我没有定义名称列表或其他什么东西。

非常感谢您的帮助。

- Alexandre.S

---

如果你的工作流程不太长，或许你可以完整地发布它。另外，您可能会对此FAQ条目中提到的glob_wildcards函数感兴趣：https://snakemake.readthedocs.io/en/stable/project_info/faq.html#glob-wildcards（“如何在特定目录中运行我的规则以处理所有文件？”） - bli

这两个规则是可行的，但确实这个 {run} 语句有点麻烦。我尝试通过一个正则表达式来改进它，该正则表达式应该在开头匹配我的 fastq 文件。我在教程中看到了这个： SAMPLES, = glob_wildcards(join(FASTQ_DIR, '{sample,Samp[^/]+}.R1.fastq.gz')) 现在我试图为我的工作流程调整这种命令。想法是最终只用样本名称命名文件（例如：Cla001d.bam）。 - Alexandre.S

在哪个步骤中，您会“合并”同一样本的运行？此外，您展示的all的输入与GATK_merge的输出不同。 - bli

我的“全部”规则不正确，我会更改它。 在映射的第一步中，我合并R1和R2文件：规则bwa_map。 - Alexandre.S

是的，但是R1和R2并不像data/Raw_reads/{sample}_R1_{run}.fastq.gz中的run一样。如果这个run部分预期始终为“001”，那么在此处不使用通配符会更简单。我不知道这是否有保证。 - bli

Answer 1

您可以尝试使用 glob_wildcards 在初始的fastq.gz文件上生成样本ID列表：

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

然后，您可以使用此功能 扩展 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)]

- bli

---

对我来说不起作用。我做了完全相同的事情：barcodes, = glob_wildcards('/path/to/fastq_pass/{barcode}')，然后在我的聚合规则中，我有input: expand("nanopolish/{barcode}/{barcode}.consensus.fasta", barcode=barcodes)，但是我得到了AttributeError: 'Wildcards' object has no attribute 'barcode'。 - BCArg

@BCArg 通配符的存在取决于输出文件模式。也许您可以提出一个更详细描述问题的单独问题。 - bli

Answer 2

您可以通过在规则中使用Python函数作为输入来实现此目的，如 Snakemake 文档中所述 （这里）。

例如，可以像这样编写：

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

希望这有所帮助。

- rioualen

---

谢谢您的回答。问题是我没有文件列表，所以"samples =<样本列表>"不起作用。这个工作流程的唯一输入是包含fastq文件的目录。我还尝试在函数中使用"os.listdir(Raw_calling/)"或将其作为规则的输入。但是snakemake告诉我"Raw_calling"目录不存在。我认为这是因为Python代码在工作流运行之前执行，因此确实还没有创建该目录。或者可能无法将合并规则与上一个规则连接？ - Alexandre.S

{sample} 应该与您的规则 GATK_raw_calling 中的相同，是一个文件列表、标识符或其他，在某个地方定义。如果规则 GATK_raw_calling 正常工作，则这不应该成为问题... - rioualen

@Alexandre.S 你说你没有文件列表，但我猜你以某种方式拥有样本ID列表，这些ID用于输入顶级规则以驱动“GATK_raw_calling”规则的执行。如果您能发布更多工作流程的内容将会很有帮助。最好提供一个包含到“GATK_raw_calling”规则的工作示例。 - bli

你好。 难点在于我正在尝试创建一个只使用一个输入（即fastq文件目录）的流水线。 这意味着我没有ID列表，只有一个目录。我的想法是使这个流水线可用于大多数需要分析的数据，并在不到1分钟的时间内完成 -> 只需将fastq放入正确的目录并运行即可。 我会发布更详细的内容，感谢您的回答！ - Alexandre.S

有更新吗？我仍在尝试构建一个具有相同目标的流水线：将输入文件放入特定文件夹并运行流水线。 - fakechek