Skip to content

Fix correctness bugs in pileup alignment metrics - #13

Draft
micah-olivas wants to merge 1 commit into
mainfrom
worktree-fix-pileup-alignment-metrics
Draft

Fix correctness bugs in pileup alignment metrics#13
micah-olivas wants to merge 1 commit into
mainfrom
worktree-fix-pileup-alignment-metrics

Conversation

@micah-olivas

Copy link
Copy Markdown
Contributor

Fixes four correctness bugs in the per-well pileup viewer (src/usortm/demux/streakout.py). The canvas rendering itself was fine — the problems were all in the metrics layered on top of the alignment, and because this view is what decides whether a well looks recoverable, each could push that call the wrong way.

Each was confirmed against synthetic reads before and after, and the deletion encoding was additionally verified end-to-end through a real minimap2/samtools alignment.

1. Deletions were invisible to "Read identity"

Deletions and uncovered positions were both encoded ("-", True), and identity excluded every - from the denominator:

read before after
30 bp deletion, zero substitutions 100.0% 96.7%
5 substitutions, full coverage 92.4% 92.4%

The cleaner-looking read was the more broken one. Now PILEUP_DEL (*) and PILEUP_NOCOV (-) are distinct, and only uncovered positions are excluded — a deletion is evidence of disagreement, not missing data.

2. The protein track silently preserved frame across a frameshift

Deleted positions were padded with N, one per reference position, so the consensus stayed reference-length and never shifted. A 1 bp deletion rendered as one changed residue:

before   refAA : KKKKKKKKKKFFFFFFFFFF
         consAA: KKXKKKKKKKFFFFFFFFFF    ← frame preserved, reads as a point substitution
after    consAA is 19 aa vs 20 aa, and residues downstream diverge

The insert is now sliced in reference coordinates first, then deleted positions are dropped, so the frame genuinely shifts. The net indel is labelled -1 bp frameshift or -3 bp in-frame indel so a shifted translation is legible rather than looking like garbage.

3. Read counts contradicted their own percentage

n_reads was the post-filter aligned-row count while frac came from the pre-filter group, so a header could read 0 reads (75%). The empty-state message interpolated that same zeroed count and therefore always printed (0 reads unaligned) — the one number it existed to convey. Sections now carry n_reads (group size) and n_aligned separately; headers read 12 of 30 reads aligned (45%) when they differ. This matters most where the midpoint filter drops concatemers, which is common and deliberate.

4. No depth floor on flagged columns

The >10% disagreement threshold was a fraction of covered reads, so one read covering one position flagged it at 100% — spurious triangles at every read end. Flagging now requires PILEUP_MIN_FLAG_DEPTH (3) covering reads. Deletions count toward disagreement.

Also fixed

  • All three row builders now share one parser (_rows_from_aligned_bam), so the symbol encoding cannot drift. _build_pileup_from_bam was dead code still carrying the old encoding — exactly how this bug would return.
  • Pipeline could hang. _run_pileup_alignment now closes the parent's copy of minimap2's stdout and checks both exit codes. Previously, if samtools exited early minimap2 never received EPIPE and wait() could block forever; a failed toolchain also degraded silently to an empty pileup. Verified: the samtools-missing case now returns in 0.01s.
  • Consistent empty-group handling — the orient-ref branch emitted a section with no reference where _make_section skipped the group.
  • Deletions render in their own colour with legend and tooltip entries; the legend now states that insertions are not shown (inherent to one column per reference position, previously undocumented).
  • Canvas backing store clamped to a max area — large pileups exceeded iOS Safari's ~16.8M device-px cap and rendered blank.
  • Inline JSON escapes </ so a payload cannot close the <script> block.
  • Removed an unreachable JS fallback ([] is truthy, so it never ran).

Verification

  • New tests/test_pileup_html.py, 30 tests — the file had zero coverage before. The render-layer tests need no aligner; the alignment-failure tests skip when minimap2/samtools are absent.
  • Full suite: 202 passed.
  • Real-alignment check: a 30 bp deletion through actual minimap2/samtools yields exactly 30 PILEUP_DEL positions and 96.67% identity (was 100%); perfect read 100%, 1-substitution read 99.89%.

The pileup viewer's rendering was sound, but the metrics layered on top of the
alignment could misreport what the reads actually showed. Since this view is
used to judge whether a well is recoverable, each of these could push the call
the wrong way.

Deletions were encoded identically to uncovered positions, as ("-", True), so
they were excluded from the identity denominator: a read with a 30 bp deletion
and no substitutions reported 100% identity, scoring better than a read with
five substitutions. Deletions and uncovered positions are now distinct symbols
(PILEUP_DEL vs PILEUP_NOCOV) and only uncovered positions are excluded.

The consensus padded deleted positions with N, keeping the sequence at
reference length, so a frameshift translated in frame and looked like a single
point substitution. The insert is now sliced in reference coordinates and
deleted positions are then dropped, so the frame actually shifts, and the net
indel is labelled as a frameshift or an in-frame indel.

n_reads held the post-filter aligned-row count while frac came from the
pre-filter group, so a header could read "0 reads (75%)". The empty-state
message interpolated the same zeroed count and so always said "(0 reads
unaligned)" - the one number it existed to report. Sections now carry n_reads
(group size) and n_aligned separately.

Problem-position flagging applied a >10% threshold with no depth floor, so a
single read covering a single position flagged it at 100%. Flagging now
requires PILEUP_MIN_FLAG_DEPTH (3) covering reads.

Also:
- All three row builders now route through one _rows_from_aligned_bam, so the
  symbol encoding cannot drift between them. _build_pileup_from_bam was dead
  code still carrying the old encoding.
- _run_pileup_alignment closes the parent's copy of minimap2's stdout and
  checks both exit codes; previously, if samtools died early minimap2 never
  got EPIPE and wait() could block forever, and a failed pipeline silently
  looked like an empty pileup.
- Deletions render in their own colour, with legend and tooltip entries, and
  the legend states that insertions are not shown.
- Canvas backing store is clamped to a maximum area; large pileups exceeded
  iOS Safari's cap and rendered blank.
- Inline JSON escapes "</" so a payload cannot close the script block.
- Removed an unreachable JS fallback branch ([] is truthy).

Adds tests/test_pileup_html.py (30 tests) covering all of the above; the
render layer needs no aligner, and the alignment-failure tests skip when
minimap2/samtools are absent. Full suite: 202 passed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant