Interval list as a wildcard in snakemake

Question

[User000]

Dear all,

I'm doing variant calling with GATK. My reference genome is a fasta file with 200K scaffolds and I prepared a list called interval.list:

scaffold1
scaffold18
scaffold23

I want to use this interval list as (-L) wildcards, to accelerate/parallelise the process. How to do this?

 INTERVAL= "interval.list"
     rule all:
         input:
             expand("{int}.vcf", int=INTERVAL)
     rule gatk:
         input:
             ref_done="ref.fasta",
            bam="bam.list"
        output:
            outf ="{int}.vcf"
        shell:
             """
             /Tools/gatk/gatk --java-options "-Xmx16g -XX:ParallelGCThreads=10" HaplotypeCaller -L {wildcards.int} -R {input.ref_done} -I {input.bam} -O {output.outf}
             """
     rule merge:
         input:
             int="{int}.vcf"
         output:
             outf = "results/res.vcf"
         shell:
             "java -jar /Tools/picard/build/libs/picard.jar GatherVcfs {wildcards.int} O={output.outf}"

[Jeremy Leipzig]

you are almost there. actually read in the intervals file into a python list (what you call INTERVAL), set a target rule that uses that list as input (you did this already). Your wildcard rule is set up properly. The merge rule will take some modification because it can't be a wildcard loop but the input is the whole list of outputs.

[User000]

So I do not know how to read the interval.list file..... as a normal python list?

[Jeremy Leipzig]

yes open the file and append each line to the list in normal python

[User000]

I added what I did so far, does it look good? seems to be working

[User000]

no, this merges only scaffold1.... :/

Answer 1

[User000]

Something like this? Please if you have suggestions how to improve the code let me know since I am learning snakemake it could be very important for me

f = open('Interval2.list')
my_interval = f.read().split("\n")
f.close()

rule gatk:
    input:
        gatk_ref_done="ref.fasta",
        bam="bam.list",
    output:
        outf ="scaffolds/{int}.vcf"
    shell:
        """
        /Tools/gatk/gatk --java-options "-Xmx16g -XX:ParallelGCThreads=10" HaplotypeCaller -L {wildcards.int} -R {input.gatk_ref_done} -I {input.bam} -O {output.outf}
        """

rule merge:
    input:
        lambda w: expand("scaffolds/{int}.vcf", int=my_interval)
    output:
        outf = "rh.vcf"
    params:
       lambda w: "I=" + " I=".join(expand("scaffolds/{int}.vcf", int=my_interval))
    shell:
        "java -jar /Tools/picard/build/libs/picard.jar GatherVcfs {params} O={output.outf}"

[Jeremy Leipzig]

this should work without your lambda w:. Just get rid of those and type snakemake merge