Repository navigation
fix(raw-bam): classify dovetail FR pairs with coincident 5' ends as FR - #1022
Conversation
htsjdk 5.0.0 SamPairUtil.getPairOrientation (htsjdk#1771), which fgbio pins, classifies a pair as FR when the positive-strand 5' position is <= the negative-strand 5' position. fgumi used a strict <, so a dovetail FR pair whose reverse read's aligned end equals the forward read's aligned start (the HEK293T geometry from #505) was classified RF from the reverse record. Because is_primary_fr_pair_raw evaluates the reverse record's arm, such pairs were rejected: - fgumi codec dropped them as NotPrimaryFrPair; - fgumi clip skipped overlap and past-mate clipping for them; - the simplex and duplex callers (num_bases_extending_past_mate_raw, via the MC-tag forward arm) did not clip read-through bases past the mate. Use the inclusive comparison in is_fr_pair_raw and the MC-based forward arm (is_fr_pair_with_mate_cigar_raw), and in the TLEN-based per-record orientation helpers (record_utils::get_pair_orientation and template::get_pair_orientation_raw). The forward TLEN arms now compute the mate 5' as htsjdk's CoordMath.getEnd(start, tlen) = start + tlen - 1, which with <= is identical to the previous start + tlen with <, so only the 5' tie changes. That arithmetic is now done in i64 in every copy. Compute alignment ends exactly as htsjdk does (start + refLength - 1, unclamped) in every orientation copy. The raw and typed copies clamped a zero reference span to end = start; with the inclusive comparison that would have made a mapped reverse read with no reference-consuming ops (e.g. 100S) at its forward mate's start FR and clipped the forward read down to one base. htsjdk ends such a read at start - 1, so the pair is RF and left untouched. Admitting tie pairs makes a degenerate overlap-clip case reachable: the midpoint equals the reverse read's alignment end, so there is no reference position after it to keep. fgbio's readPosAtRefPos returns 0 there and leaves the reverse read unclipped (no upgrade or masking in any mode); clip_overlapping_reads now does the same instead of clipping away the reverse read's whole alignment. This changes output for these pairs, matching fgbio: codec now emits consensus reads for them, clip now clips overlap and past-mate bases, and simplex/duplex now trim read-through bases (e.g. 7 bases per read on the HEK293T pair) before consensus.
|
Navigate logical layers of code changes, visualize relationships, and explore their blast radius. 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 configuration
📒 Files selected for processing (5)
Included review availability: This review used your included allowance. 2 included reviews remain after this review. Your included PR review attempts over the past 7 days set your current allowance at 3 reviews per hour.
WalkthroughFR orientation now treats coincident 5′ ends as FR across raw BAM, SAM, and template checks. Overlap clipping handles dovetail and zero-reference-span cases. Consensus tests cover how dovetail pairs are handled under two overlap-window modes. ChangesFR orientation and dovetail handling
Priority: ➖ Normal Estimated code review effort: 3 (Moderate) | ~25 minutes Change: Bug fix Suggested labels: Merge Risk: ⚪ Minimal · up to No actionable issue is established for this change; it is mergeable after normal checks. 🚥 Pre-merge checks | ✅ 3✅ Passed checks (3 passed)
Comment |
|
@coderabbitai pause |
|
@coderabbitai review |
✅ Action performedReviews paused. |
✅ Action performedReview finished.
|
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1022 +/- ##
==========================================
- Coverage 96.47% 96.46% -0.01%
==========================================
Files 299 299
Lines 152124 152297 +173
==========================================
+ Hits 146756 146920 +164
- Misses 5368 5377 +9 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
This branch was successfully deployed
1 participant
Changes output for
codec,clip(overlap and past-mate clipping),simplexandduplexon dovetail FR pairs whose 5' ends coincide. The new output matches fgbio.htsjdk 5.0.0
SamPairUtil.getPairOrientation(samtools/htsjdk#1771, pinned by fgbio) classifies a pair as FR when the positive-strand 5' position is<=the negative-strand 5' position. fgumi used a strict<. So a dovetail FR pair whose reverse read's aligned end equals the forward read's aligned start was classified RF from the reverse record. This is the HEK293T CODEC geometry #505 targeted; see fgbioCodecConsensusCallerTest.scala:210and:359. As a result:codecdropped these pairs asNotPrimaryFrPair.clipskipped overlap and past-mate clipping for them.simplex/duplex(num_bases_extending_past_mate_raw, MC arm) did not trim read-through bases past the mate.Changes
<=:is_fr_pair_raw,is_fr_pair_with_mate_cigar_raw,record_utils::get_pair_orientationandtemplate::get_pair_orientation_raw.CoordMath.getEnd(start, tlen) = start + tlen - 1. Combined with<=, this is identical to the oldstart + tlenwith<, so only the tie changes.start + refLength - 1, unclamped. The raw and typed copies clamped a zero reference span toend = start. With<=, that would have made a mapped100Sreverse read at its mate's start FR, and its forward mate would have been clipped to one base. htsjdk givesstart - 1, so the pair is RF and left untouched.readPosAtRefPosthen returns 0 and leaves the reverse read unclipped, with no upgrade or masking.clip_overlapping_readsnow does the same instead of clipping away the reverse read's whole alignment.Tests
Each test cites its fgbio or htsjdk source. A mutation check of each change made at least one test fail.
is_primary_fr_pair_rawon the fgbio pair, both argument orders.codec: the pair is no longer rejected as non-FR. With the default window it gives 1 consensus read. Under--legacy-overlap-windowit is rejected asDovetail.clip:(52, 0), and the reverse read is untouched.(99, 99).simplex/duplexread-through:num_bases_extending_past_mate_rawon tie pairs gives 7/7 (HEK293T pair) and 99/99 (htsjdk pair).Out of scope
The typed and template per-record helpers still use TLEN on the forward arm, where htsjdk 5.0.0 prefers MC when present. Their docs now say so. Consolidating the four orientation copies into one helper is left for a follow-up.