Fix correctness bugs in pileup alignment metrics - #13
Draft
micah-olivas wants to merge 1 commit into
Draft
Conversation
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.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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:The cleaner-looking read was the more broken one. Now
PILEUP_DEL(*) andPILEUP_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: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 frameshiftor-3 bp in-frame indelso a shifted translation is legible rather than looking like garbage.3. Read counts contradicted their own percentage
n_readswas the post-filter aligned-row count whilefraccame from the pre-filter group, so a header could read0 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 carryn_reads(group size) andn_alignedseparately; headers read12 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
_rows_from_aligned_bam), so the symbol encoding cannot drift._build_pileup_from_bamwas dead code still carrying the old encoding — exactly how this bug would return._run_pileup_alignmentnow closes the parent's copy of minimap2's stdout and checks both exit codes. Previously, if samtools exited early minimap2 never received EPIPE andwait()could block forever; a failed toolchain also degraded silently to an empty pileup. Verified: the samtools-missing case now returns in 0.01s._make_sectionskipped the group.</so a payload cannot close the<script>block.[]is truthy, so it never ran).Verification
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.PILEUP_DELpositions and 96.67% identity (was 100%); perfect read 100%, 1-substitution read 99.89%.