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}" 

我試過幾件事情，爲了做到這一點，但不能得逞。問題在於兩條規則之間的聯繫（GATK_raw_calling和GATK_merge，應該合併GATK_raw_calling的輸出）。如果我將GATK_raw_calling的輸出指定爲以下規則的輸入（輸入文件中的通配符無法從輸出文件中確定），則無法輸出單個文件，而且我無法在兩個規則之間建立鏈接如果我不指定這些文件作爲輸入...

有沒有辦法成功地做到這一點？我想，難點在於我沒有定義名單或其他名稱。

非常感謝您的幫助。

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，則需要使用zip而不是展開，以避免不存在的組合：

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)]