Fix variant calling with deepvariant that outputs non-allowed variant…

…s for hapcut2

src/blr/Snakefile

@@ -935,6 +935,7 @@ rule call_variants:
                 " --output_vcf={output.vcf}"
                 " --reads={input.bam}"
                 " --ref={params.reference}"
+                " 2> {log}" 
                 }
         shell(commands[config["variant_caller"]])
 
@@ -995,7 +996,14 @@ rule filter_vcfs:
     input:
         vcf = "chunks/{base}.variants.called.{type}.filtered.vcf"
     shell:
-        "bcftools view -f 'PASS' {input.vcf} > {output.vcf}"
+        "bcftools view"
+        # Only heterozygous
+        " -g het"
+        # Only keep genotypes that hapcut2 accepts. I.e. in {0, 1, 2}
+        # Mainly filters variants from DeepVariant
+        " -i 'GT == \"0/1\" || GT == \"0/2\" || GT == \"1/2\"'"
+        " -f 'PASS'"
+        " {input.vcf} > {output.vcf}"
 
 
 rule merge_filtered_snps_indels:

src/blr/rules/phasing.smk

@@ -31,6 +31,7 @@ rule hapcut2_extracthairs:
         "extractHAIRS"
         " --10X 1"
         " --indels {params.indels}"
+        # TODO - test to include option `--triallelic` to include GT 1/2 
         " --realign_variants 1"  # Improves overall error-rate
         " --ref {params.reference}"
         " --bam {input.bam}"