closes #56

CHANGES.md

@@ -1,6 +1,7 @@
 v0.2.11 (dev)
 =============
 + more informative error message on bad sample name (#53)
++ allow setting SOMALIER_AB_HOM_CUTOFF to change which calls are considered hom-ref (#56)
 
 v0.2.10
 =======

README.md

@@ -166,6 +166,9 @@ all cases.
 
 By default `somalier` will only consider variants that have a "PASS" or "RefCall" FILTER. To extend this list, set
 the environment variable `SOMALIER_ALLOWED_FILTERS` to a comma-delimited list of additional filters to allow.
+
+by default sites with an allele balance < 0.01 will be considered homozygous reference. To adjust this, use e.g. : 
+`SOMALIER_AB_HOM_CUTOFF=0.04 somalier relate ...`
 
 ## Other Work
 

src/somalierpkg/relate.nim

@@ -293,7 +293,16 @@ proc ab*(c:allele_count, min_depth:int): float {.inline.} =
     return 0
   result = c.nalt.float / (c.nalt + c.nref).float
 
-proc alts*(ab:float, min_ab:float, ab_cutoff:float=0.01): int8 {.inline.} =
+var ab_cutoff: float = 0.01
+try:
+  ab_cutoff = parseFloat(getEnv("SOMALIER_AB_HOM_CUTOFF"))
+  if ab_cutoff > 0.5:
+    stderr.writeline("[somalier] error setting SOMALIER_AB_HOM_CUTOFF to:" & getEnv("SOMALIER_AB_HOM_CUTOFF"))
+    ab_cutoff = 0.01
+except:
+  discard
+
+proc alts*(ab:float, min_ab:float, ab_cutoff:float=ab_cutoff): int8 {.inline.} =
   if ab < 0: return -1
   if ab < ab_cutoff: return 0
   if ab > (1 - ab_cutoff): return 2