CODEWITHSUNDEEP
python-3.xsnakemake
njBernstein
njBernstein

Reputation: 115

how to propapagate snakemake wildcards through a lambda function

I'm trying to merge vcf files together using snakemake, but get the error:

Building DAG of jobs...
MissingInputException in line 21 of 

Missing input files for rule all:
outputs/****.g.vcf.gz

The goal is to call variants on separate chromosomes and then merge them back together. It seems like the sample wildcard isn't getting propagated through my lambda function. I've tried a few different iterations and haven't been able to crack it. I'm sure that the rest of the code is okay because when I remove the merge function and just call variants across all chromosomes the file works fine.

Any help would be much appreciated.

import glob

configfile: "config.json"

chroms = [1, 2, 3, 4, 5]
str_chroms = ["chr{}".format(chr) for chr in chroms]


def get_fq1(wildcards):
    # code that returns a list of fastq files for read 1 based on
    # *wildcards.sample* e.g.
    return sorted(glob.glob(wildcards.sample + '*_R1_001.fastq.gz'))


def get_fq2(wildcards):
    # code that returns a list of fastq files for read 2 based
    # on *wildcards.sample*, e.g.
    return sorted(glob.glob(wildcards.sample + '*_R2_001.fastq.gz'))


rule all:
    input:
        "outputs/" + config['sample'] + "_picard_alignment_metrics_output.txt",
        "outputs/" + config['sample'] + "_fastqc",
        "outputs/" + config['sample'] + "_analyze_covariates.pdf",
        "outputs/" + config['sample'] + ".g.vcf.gz",
        "outputs/" + config['sample'] + ".coverage"

rule bwa_map:
    input:
        config['reference_file'],
        get_fq1,
        get_fq2
    output:
        "outputs/{sample}_sorted.bam"
    shell:
        "bwa mem -t 16 {input} |  samtools view -bS - | \
        samtools sort -@ 16 -m 7G - -o {output}"

#a bunch of intermediate steps that are not the issue

rule variant_calling:
    input:
        bam = "outputs/{sample}_recal_reads.bam",
        bai = "outputs/{sample}_recal_reads.bam.bai",
        reference_file = config['reference_file']
    output:
       "outputs/{sample}_{chr}.g.vcf.gz"
    shell:
        """gatk --java-options "-Xmx128g" HaplotypeCaller \
       -R {reference_file} -I {input.bam} -L {wildcards.chr}\
       -O {output} -ERC GVCF"""

rule merge_vcfs:
    input:
        lambda wildcards: expand("outputs/{sample}_{chr}.g.vcf.gz",
                             chr=str_chroms, 
                      sample=wildcards.sample)
     output:
         "output/{sample}.g.vcf.gz"
     shell:
        "vcf-merge {input} | bgzip -c > {output}"

Upvotes: 1

Views: 2421

Answers (2)

JohnnyBD
JohnnyBD

Reputation: 161

Yes, as @JeeYem mentiones you have a typo in output file for rule merge.

Also I do not see the need for the lambda in rule merge? You are passing the same set of chromosome irrespective of sample? The str_chroms is independent of sample in your setup so you can rewrite it as:

rule merge_vcfs:
    input: expand("outputs/{{sample}}_{chr}.g.vcf.gz",chr=str_chroms)
    output: "output/{sample}.g.vcf.gz"
    shell: "vcf-merge {input} | bgzip -c > {output}"

Upvotes: 2

Manavalan Gajapathy
Manavalan Gajapathy

Reputation: 4089

Typo in output of rule merge_vcfs. It should be outputs/{sample}.g.vcf.gz (ie. outputs instead of output)

Upvotes: 3

Related Questions