fix (regression) for gvcf with emtpy alts

closes #46

CHANGES.md

@@ -2,6 +2,7 @@ v0.2.10
 =======
 + fix extra output column in pairs.tsv (#47)
 + change output file name of ancestry to include
++ fix for gvcf with empty alts (#46)
 
 v0.2.9
 ======

src/somalier.nim

@@ -36,7 +36,7 @@ proc looks_like_gvcf_variant(v:Variant): bool {.inline.} =
 
 proc get_variant(ivcf:VCF, site:Site): Variant =
   for v in ivcf.query(&"{site.chrom}:{site.position+1}-{site.position+2}"):
-    if v.start == site.position and (
+    if v.start == site.position and v.ALT.len > 0 and (
       (v.REF == site.A_allele and v.ALT[0] == site.B_allele) or
       (v.REF == site.B_allele and v.ALT[0] == site.A_allele)):
       return(v.copy())
@@ -126,7 +126,12 @@ proc get_ref_alt_counts(ivcf:VCF, sites:seq[Site], fai:Fai=nil): seq[counts] =
 
     var mult = int(AD.len / vcf_samples.len)
     for j, s in vcf_samples:
-      var ac = allele_count(nref: max(0, AD[mult * j]).uint32, nalt: max(0, AD[mult * j + 1]).uint32, nother: 0)
+      var ac:allele_count
+      if mult >= 2:
+        ac = allele_count(nref: max(0, AD[mult * j]).uint32, nalt: max(0, AD[mult * j + 1]).uint32, nother: 0)
+      else:
+        ac = allele_count(nref: 0'u32, nalt: 0'u32, nother: 0)
+
       doAssert site.A_allele < site.B_allele
       # we always store a site by alphabetical order, so sometimes we have to
       # flip the alts.

tests/functional-tests.sh

@@ -46,3 +46,8 @@ assert_exit_code 0
 # checks that final column, expected relatedness is 1.
 assert_equal "1" $(awk 'NR == 1 && NR == 1.0' test_prefix_A/out.pairs.tsv | wc -l) 
 rm -rf test_prefix_A test_prefix_B
+
+
+run check_gvcf_no_alt $exe extract -s tests/x.gvcf.gz -f tests/test.fa tests/x.gvcf.gz
+assert_exit_code 0
+