2
$\begingroup$

I have written a rule for CombineGVCFs in gatk4. The rule is as follow

all_gvcf = get_all_gvcf_list()

rule cohort:
  input:
    all_gvcf_list = all_gvcf,
    ref="/data/refgenome/hg38.fa",
    interval_list = prefix+"/bedfiles/hg38.interval_list",
  params:
    extra = "--variant",
  output:
    prefix+"/vcf/cohort.g.vcf",
  shell:
    "gatk CombineGVCFs -R {input.ref} {params.extra} {input.all_gvcf_list} -O {output} --tmp-dir=/data/tmp -L {input.interval_list}"

all_gvcf is the dataset for all gvcf files which will be combined. But the problem is I need to add --variant parameter before every input. Command I'm getting right now is as follows

gatk CombineGVCFs -R /data/refgenome/hg38.fa --variant /data/prjna644607/vcf/SRR12165216_HC.g.vcf /data/prjna644607/vcf/SRR12165217_HC.g.vcf /data/prjna644607/vcf/SRR12165218_HC.g.vcf /data/prjna644607/vcf/SRR12165219_HC.g.vcf -O /data/prjna644607/vcf/cohort.g.vcf --tmp-dir=/data/tmp -L /data/prjna644607/bedfiles/hg38.interval_list

The command I want to achieve is as follows

gatk CombineGVCFs -R /data/refgenome/hg38.fa --variant /data/prjna644607/vcf/SRR12165216_HC.g.vcf --variant /data/prjna644607/vcf/SRR12165217_HC.g.vcf --variant /data/prjna644607/vcf/SRR12165218_HC.g.vcf --variant /data/prjna644607/vcf/SRR12165219_HC.g.vcf -O /data/prjna644607/vcf/cohort.g.vcf --tmp-dir=/data/tmp -L /data/prjna644607/bedfiles/hg38.interval_list

How can I add this extra "--variant" tag before every input? I have added it in the get_all_gvcf_list() function. But then snakmake gives me inputfiles not found problem.

$\endgroup$

2 Answers 2

5
$\begingroup$

Found out the problem. Turns out I can write a lambda function as follows

params:
        extra=lambda wildcards, input: ' -V '.join(input.all_gvcf_list)

and add '-V' before {params.extra} in the "shell" sction. That solves the problem

$\endgroup$
2
  • $\begingroup$ i am facing a similar issue can you please paste the whole code for your solution? $\endgroup$
    – Ayam Gupta
    Commented Jul 15, 2022 at 8:54
  • $\begingroup$ That is the whole rule/code actually. get_all_gvcf_list() will give the absolute path of the files in a dataset. You just need to see if the problem really matches with you. If u still need any more help then maybe you can write the entire problem and tag the question with me. $\endgroup$ Commented Jul 16, 2022 at 10:51
-1
$\begingroup$

well I came up with the solution. I think it is a little unconventional.

    rule gatk_CombineGVCFs:                                                 
    
        input:
            vcf_dummy = expand("gvcf/{sample}.g.vcf.gz", sample = SAMPLES), # a dummy vcf to connect this rule to gatk_HaplotypeCaller
            fa = "data/genome.fa",
    
        output:
            gvcf = "gvcf_combined/all.g.vcf.gz",
            #index = directory("gvcf_index")                                 #index generated without seperate output
        params:
            list = "gvcfs.list"
            log: "logs/CombineGVCFs.log"
        shell:
            "ls gvcf/*.gz > gvcfs.list && gatk --java-options -Xmx7g CombineGVCFs -R {input.fa} --variant {params.list} -O {output.gvcf} "   
             #Two bash comands running sequentially to make a gvcf.list as during my testing the combineGVCFs tool takes lists as input. 
$\endgroup$
3
  • 1
    $\begingroup$ You came up with a solution, not the solution (by your own admission). $\endgroup$
    – Ram RS
    Commented Jul 26, 2022 at 14:30
  • $\begingroup$ The solution to your problem should be described elaborately. Either describe your answer elaborately else add your question separately and then answer $\endgroup$ Commented Jul 26, 2022 at 20:21
  • $\begingroup$ The comments in the code describe the solution pretty elaborately. The problem was in executing the the combineGVCF to which my solution was to add --variant in the shell and first make a list of vcfs with ls and pass it as a list to combineGVCF. $\endgroup$
    – Ayam Gupta
    Commented Jul 27, 2022 at 4:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.