Skip to content

Commit

Permalink
fix bug in --inference
Browse files Browse the repository at this point in the history
  • Loading branch information
brentp committed Apr 7, 2021
1 parent d11e061 commit df0e63c
Show file tree
Hide file tree
Showing 2 changed files with 5 additions and 2 deletions.
1 change: 1 addition & 0 deletions CHANGES.md
@@ -1,6 +1,7 @@
v0.2.13
=======
+ add "Heterozygosity rate" as a per-sample metric to the html output. (Thanks Irenaeus and Kelly for the suggestion)
+ fix inference for some cases. obvious parent-child pairs were sometimes missed.


v0.2.12
Expand Down
6 changes: 4 additions & 2 deletions src/somalierpkg/relate.nim
Expand Up @@ -235,6 +235,7 @@ proc relatedness(r: var relation_matrices, j: int, k:int): relation {.inline.} =
ibs0: r.ibs[j * r.n_samples + k],
shared_hets: r.ibs[k * r.n_samples + j],
shared_hom_alts: r.shared_hom_alts[j * r.n_samples + k],
het_ab: r.shared_hom_alts[k * r.n_samples + j],
ibs2: r.n[k * r.n_samples + j],
n: r.n[j * r.n_samples + k],
x_ibs0: r.x[j * r.n_samples + k],
Expand All @@ -254,7 +255,7 @@ proc relatedness(r: var relation_matrices, j: int, k:int): relation {.inline.} =
r.x[j * r.n_samples + k] = xir.IBS0.uint16
r.x[k * r.n_samples + j] = xir.IBS2.uint16

return relation(#sample_a: sample_names[j],
result = relation(#sample_a: sample_names[j],
#sample_b: sample_names[k],
hets_a: hets_j, hets_b: hets_k,
hom_alts_a: r.gt_counts[2][j], hom_alts_b: r.gt_counts[2][k],
Expand Down Expand Up @@ -453,7 +454,7 @@ const missing = [".", "0", "-9", ""]
proc high_quality(gt_counts: array[5, seq[uint16]], i:int): bool {.inline.} =
# less than 3% of sites with allele balance outside of 0.1 .. 0.9
result = gt_counts[4][i].float / (gt_counts[0][i] + gt_counts[1][i] + gt_counts[2][i] + gt_counts[3][i] + gt_counts[4][i]).float < 0.06
if not result:
if not result:
return false
result = gt_counts[2][i].float / gt_counts[1][i].float > 0.7

Expand Down Expand Up @@ -995,6 +996,7 @@ specified as comma-separated groups per line e.g.:
let expected_relatedness = if idx == -1: -1'f else: groups[idx].rel
let rr = rel.rel

# PARENT_CHILD
if rr > 0.4 and rr < 0.6 and rel.ibs0.float / rel.ibs2.float < 0.005:
parent_child_pair.mgetOrPut(rel.sample_a, @[]).add(rel.sample_b)
parent_child_pair.mgetOrPut(rel.sample_b, @[]).add(rel.sample_a)
Expand Down

0 comments on commit df0e63c

Please sign in to comment.