diff --git a/CHANGES.md b/CHANGES.md index acfd58e..9852536 100644 --- a/CHANGES.md +++ b/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 ======= diff --git a/README.md b/README.md index 52e707f..54a1e2a 100644 --- a/README.md +++ b/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 diff --git a/src/somalierpkg/relate.nim b/src/somalierpkg/relate.nim index c286717..6e83257 100644 --- a/src/somalierpkg/relate.nim +++ b/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