Add g.vcf file to default output of the pipeline

Als include some tests for the content of the g.vcf file

Snakefile

@@ -174,6 +174,7 @@ rule all:
         report="multiqc_report/multiqc_report.html",
         bais=expand("{sample}/bams/{sample}.markdup.bam.bai",sample=SAMPLES),
         vcfs=expand("{sample}/vcf/{sample}.vcf.gz",sample=SAMPLES),
+        gvcfs=expand("{sample}/vcf/{sample}.g.vcf.gz",sample=SAMPLES),
         stats=metrics()
 
 
@@ -352,6 +353,26 @@ rule gvcf_scatter:
            "-variant_index_type LINEAR -variant_index_parameter 128000 "
            "-BQSR {input.bqsr}"
 
+def aggregate_gvcf(wildcards):
+    checkpoint_output = checkpoints.scatterregions.get(**wildcards).output[0]
+    return expand("{{sample}}/vcf/{{sample}}.{i}.g.vcf.gz",
+           i=glob_wildcards(os.path.join(checkpoint_output, 'scatter-{i}.bed')).i)
+
+def aggregate_gvcf_tbi(wildcards):
+    checkpoint_output = checkpoints.scatterregions.get(**wildcards).output[0]
+    return expand("{{sample}}/vcf/{{sample}}.{i}.g.vcf.gz.tbi",
+           i=glob_wildcards(os.path.join(checkpoint_output, 'scatter-{i}.bed')).i)
+
+rule gvcf_gather:
+    """ Join the gvcf files together """
+    input:
+        gvcfs = aggregate_gvcf,
+        tbis = aggregate_gvcf_tbi,
+    output:
+        gvcf = "{sample}/vcf/{sample}.g.vcf.gz"
+    singularity: containers["bcftools"]
+    shell: "bcftools concat {input.gvcfs} --allow-overlaps --output {output.gvcf} "
+           "--output-type z"
 
 rule genotype_scatter:
     """Run GATK's GenotypeGVCFs by chunk"""
@@ -370,7 +391,7 @@ rule genotype_scatter:
            "-V {input.gvcf} -o '{output.vcf}'"
 
 
-def aggregate_input(wildcards):
+def aggregate_vcf(wildcards):
     checkpoint_output = checkpoints.scatterregions.get(**wildcards).output[0]
     return expand("{{sample}}/vcf/{{sample}}.{i}.vcf.gz",
            i=glob_wildcards(os.path.join(checkpoint_output, 'scatter-{i}.bed')).i)
@@ -379,7 +400,7 @@ def aggregate_input(wildcards):
 rule genotype_gather:
     """Gather all genotyping VCFs"""
     input:
-        vcfs = aggregate_input
+        vcfs = aggregate_vcf
     output:
         vcf = "{sample}/vcf/{sample}.vcf.gz"
     singularity: containers["bcftools"]

tests/test_dry_run.yml

@@ -13,6 +13,8 @@
     contains:
       - "Job counts"
       - "localrule all:"
+      - "output: micro/vcf/micro.vcf.gz"
+      - "output: micro/vcf/micro.g.vcf.gz"
     must_not_contain:
       - "rror"
   stderr:

tests/test_integration_run.yml

@@ -31,8 +31,21 @@
 
              - "chrM\t4745\t.\tA\tG\t5655.77\t."
              - "GT:AD:DP:GQ:PGT:PID:PL\t1/1:1,133:134:99:1|1:4745_A_G:5684,404,0"
+
              - "chrM\t4769\t.\tA\tG\t5182.77\t."
              - "GT:AD:DP:GQ:PGT:PID:PL\t1/1:1,120:121:99:1|1:4745_A_G:5211,363,0"
 
              - "chrM\t16023\t.\tG\tA\t1880.77\t."
              - "GT:AD:DP:GQ:PL\t0/1:72,73:145:99:1909,0,1872"
+
+   - path: "micro/vcf/micro.g.vcf.gz"
+     contains:
+             - "chrM\t1\t.\tG\t<NON_REF>\t.\t.\tEND=151\tGT:DP:GQ:MIN_DP:PL\t0/0:165:99:137:0,120,1800"
+
+             - "chrM\t152\t.\tT\tC,<NON_REF>\t3960.77"
+             - "GT:AD:DP:GQ:PL:SB\t1/1:0,130,0:130:99:3989,388,0,3989,388,3989:0,0,47,83"
+
+             - "chrM\t16023\t.\tG\tA,<NON_REF>\t1880.77\t."
+             - "GT:AD:DP:GQ:PL:SB\t0/1:72,73,0:145:99:1909,0,1872,2125,2089,4214:35,37,36,37"
+
+             - "chrM\t16560\t.\tC\t<NON_REF>\t.\t.\tEND=16569\tGT:DP:GQ:MIN_DP:PL\t0/0:188:0:180:0,0,0"