From 89fcb12fdb18017fd73c0b36e0e9d9e1fead3bab Mon Sep 17 00:00:00 2001 From: Brent Pedersen Date: Tue, 30 Jun 2020 10:56:42 -0600 Subject: [PATCH] closes #56 --- CHANGES.md | 1 + README.md | 3 +++ src/somalierpkg/relate.nim | 11 ++++++++++- 3 files changed, 14 insertions(+), 1 deletion(-) 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