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

Question

I'm using snakemake for the first time in order to build a basic pipeline using cutadapt, bwa and GATK (trimming ; mapping ; calling). I would like to run this pipeline on every fastq file contained in a directory, without having to specify their name or whatever in the snakefile or in the config file. I would like to succeed in doing this.

The first two steps (cutadapt and bwa / trimming and mapping) are running fine, but I'm encountering some problems with GATK.

First, I have to generate g.vcf files from bam files. I'm doing this using these rules:

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

These rules work fine. For example, if I have the files : Cla001d_S281_L001_R1_001.fastq.gz Cla001d_S281_L001_R2_001.fastq.gz

I can create one bam file (Cla001d_S281_L001_001.bam) and from that bam file create a GVCF file (Cla001d_S281_L001_001.g.vcf). I have a lot of sample like this one, and I need to create one GVCF file for each, and then merge these GVCF files into one file. The problem is that I'm unable to give the list of the file to merge to the following rule:

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

I tried several things in order to do that, but cannot succeed. The problem is the link between the two rules (GATK_raw_calling and GATK_merge that is supposed to merge the output of GATK_raw_calling). I can't output one single file if I'm specifying the output of GATK_raw_calling as the input of the following rule (Wildcards in input files cannot be determined from output files), and I'm unable to make a link between the two rules if I'm not specifying these files as an input...

Is there a way to succeed in doing that? The difficulty is that I'm not defining a list of names or whatever, I think.

Thanks you in advance for your help.

Answer 1

You can try to generate a list of sample IDs using glob_wildcards on the initial fastq.gz files:

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

Then, you can use this to expand the input of GATK_merge:

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

If the same run ID always come with the same sample ID, you will need to zip instead of expanding, in order to avoid non-existing combinations:

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

Answer 2

You can achieve this by using a python function as an input for your rule, as described in the snakemake documentation here.

Could look like this for example:

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

Hope this helps.