Repository navigation
fix(zipper,clip): compute supplementary TLEN instead of copying the mate primary's - #684
Conversation
|
Note Reviews pausedUse the following commands to manage reviews:
Use the checkboxes below for quick actions:
No actionable comments were generated in the recent review. 🎉 ℹ️ Recent review info⚙️ Run configurationConfiguration used: Path: .coderabbit.yaml Review profile: ASSERTIVE Plan: Pro Run ID: 📒 Files selected for processing (4)
WalkthroughSupplementary TLEN is now computed from each supplementary alignment and its mate primary. Insert-size calculation handles strand orientation, unmapped mates, cross-reference pairs, overflow, and stale TLEN values. ChangesSupplementary TLEN recalculation
Estimated code review effort: 4 (Complex) | ~45 minutes Possibly related issues
Possibly related PRs
Sequence Diagram(s)sequenceDiagram
participant Alignment as Alignment record
participant MateRepair as Template or clip mate repair
participant InsertSize as compute_insert_size_from_ends
participant MatePrimary as Mate primary
Alignment->>MateRepair: Provide supplementary coordinates
MateRepair->>MatePrimary: Read mate-primary coordinates
MateRepair->>InsertSize: Compute strand-aware TLEN
InsertSize-->>MateRepair: Return TLEN
MateRepair-->>Alignment: Store independently computed TLEN
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Comment |
|
@coderabbitai pause |
✅ Action performedReviews paused. |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #684 +/- ##
========================================
Coverage 93.94% 93.95%
========================================
Files 178 178
Lines 108058 108166 +108
========================================
+ Hits 101518 101630 +112
+ Misses 6540 6536 -4 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
There was a problem hiding this comment.
Actionable comments posted: 1
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@src/lib/template.rs`:
- Around line 823-880: Extract the shared htsjdk-style insert-size calculation
over ref_id, five_prime, and is_unmapped into a common function reusable by
template.rs and clip.rs. In src/lib/template.rs:823-880, have
compute_insert_size_from_ends delegate to it while preserving snapshot
conversion; in src/lib/commands/clip.rs:856-933, replace
compute_insert_size_raw’s duplicated formula with the same helper and retain
MateSnap only for mapq and cigar data.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: 18b2bceb-a307-4882-bbe5-aadd48246ca5
📒 Files selected for processing (4)
src/lib/commands/clip.rssrc/lib/commands/zipper.rssrc/lib/template.rstests/integration/test_clip_command.rs
…ate primary's `Template::fix_mate_info` and `clip`'s `set_supplemental_mate_info_raw` set a supplementary alignment's TLEN by negating its mate primary's TLEN, porting htsjdk's `SamPairUtil.setMateInformationOnSupplementalAlignment`. That assumes the supplementary sits where its own primary sits — same reference, same side of the mate — which is false by construction for the split alignments this path touches. Two consequences: a supplementary on a different reference than its mate got a non-zero TLEN, and one lying beyond its mate got both the wrong sign and the wrong magnitude (the primary pair's insert size, describing coordinates the supplementary does not occupy). Compute TLEN from the supplementary's own alignment against the mate primary instead. `compute_insert_size_raw` already returns 0 for unmapped and cross-reference pairs, so both cases fall out of routing through it; it is split into an `InsertSizeEnd` snapshot plus `compute_insert_size_from_ends` so the supplementary loops can use it without a second borrow of `self.records`. This matches bwa-mem and minibwa, which compute per record from its own 5' position and emit 0 across references, including for supplementary records. The existing supplementary-TLEN tests asserted the copied value and were written for fgbio parity; they now assert the computed value, so fgumi intentionally diverges from fgbio and Picard MergeBamAlignment on these records until the upstream fix lands. Adds case tables covering cross-reference, beyond-mate, before-mate and coincident-5'-end geometries for both zipper and clip. Closes #673
f0669a1 to
563875b
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
Closes #673.
Template::fix_mate_info(zipper) andset_supplemental_mate_info_raw(clip) set a supplementary alignment's TLEN by negating its mate primary's TLEN, a faithful port of htsjdk'sSamPairUtil.setMateInformationOnSupplementalAlignment. That encodes an unstated assumption — that the supplementary occupies the same place in the template geometry as its own primary, same reference and same side of the mate — which is false by construction for the split alignments this path exists to handle.Two failure modes
Cross-reference. A supplementary mapped to a different reference than its mate gets a non-zero TLEN, so the record carries
RNAME != RNEXTwithTLEN != 0.Same reference, supplementary beyond its mate. Both sign and magnitude are wrong. With R1 at chr1:1,000 forward, R2 at chr1:1,350 reverse and an R1 supplementary at chr1:5,000, the supplementary received
+450— the primary pair's insert size — when it is the rightmost segment of the template and so must be negative.Neither case had test coverage: every existing supplementary-TLEN test placed the supplementary on the same reference as its mate and asserted the copied value.
The fix
Compute TLEN from the supplementary's own alignment against the mate primary.
compute_insert_size_rawalready returns0for unmapped and cross-reference pairs, so both failure modes fall out of routing the supplementary path through it — the correct logic was already sitting in the same file, just not wired to this path.In
template.rsit is split into anInsertSizeEndsnapshot pluscompute_insert_size_from_ends, so the supplementary loops can compute against the mate primary without holding a second borrow ofself.records. No per-record allocation is added on either path.This also removes a latent ordering dependency: the old code required the primary pair to be fixed first so the supplementary would read an updated TLEN (fgbio carries the same constraint, noted at
Bams.scala:111). Computing from coordinates makes the two steps independent.Why this basis
The SAM specification is silent on supplementary TLEN rather than violated by the old behaviour — hts-specs #522 scoped the definition to primary reads and left non-primary records undefined, and #842 leaves the computation aligner-defined. The case rests on two other grounds:
compute_insert_size_rawguards cross-reference; the supplementary path did not.bwa-mem(bwamem.c:887-892) andminibwa(format.c:243-262) compute per emitted record from that record's own 5′ position and emit0across references — including for supplementary records.samtools fixmateleaves supplementaries untouched. Only the htsjdk lineage copies the value from a different record.The existing 5′-based pairwise basis is deliberately retained rather than htslib's chain-wide leftmost-to-rightmost basis: chain-wide extents would also change the primaries' TLEN whenever a supplementary falls outside the pair's span, diverging from bwa on records that all implementations currently agree on.
Scope and divergence
Three sites across two commands:
template.rs:533and:578(zipper),clip.rs:1004(clip). The issue as originally filed covered onlyzipper.The existing supplementary-TLEN tests were written for fgbio parity and asserted the copied value; they now assert the computed value. This is an intentional divergence from fgbio and Picard
MergeBamAlignmenton these records until the upstream fix lands — thecompareharness will report it.Root cause is
htsjdk SamPairUtil.java:354, introduced in 9e03608a (2014) with behaviour unchanged since. Filed upstream as samtools/htsjdk#1795, tracked for fgbio as fulcrumgenomics/fgbio#1165.Tests
Eight existing assertions updated (seven unit, one integration), each previously asserting the copied value. Eight new
rstestcases added acrosszipperandclipcovering cross-reference, beyond-mate, before-mate, and coincident-5′-end geometries; the case tables seed a sentinel TLEN so a regression cannot pass by coincidence.Worth noting that
test_clip_command_threads_mode_supplementary_mate_repairalready had exactly the geometry at issue — a supplementary 5 kb from its mate — and was asserting+298where the correct value is-4604. Nothing was checking supplementary TLEN correctly.cargo ci-fmt,cargo ci-lint, andcargo ci-testall pass (6898 tests).Investigation and reproduction assisted by Claude Code (Anthropic). All code references and the htsjdk 5.0.0 reproduction cited in #673 were verified by hand.
Summary by CodeRabbit