Repository navigation
fix(group,dedup): pair duplex strands sharing an unclipped 5' coordinate - #485
Conversation
|
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)
WalkthroughThis PR makes raw-record R1 ordering strand-aware across dedup and group, and adds duplex molecule-base mismatch counting to grouping comparisons in ChangesStrand-aware tie-break and duplex validation
Estimated code review effort: 4 (Complex) | ~50 minutes Possibly related issues
Possibly related PRs
Suggested labels: 🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Comment |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #485 +/- ##
==========================================
+ Coverage 91.06% 91.12% +0.06%
==========================================
Files 78 78
Lines 51376 51540 +164
==========================================
+ Hits 46783 46964 +181
+ Misses 4593 4576 -17 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
Caution Failed to replace (edit) comment. This is likely due to insufficient permissions or the comment being deleted. Error details |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
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/commands/dedup.rs`:
- Around line 469-489: The R1 tie-break ordering logic is duplicated in both
is_r1_genomically_earlier_raw and the matching helper in group, so it should be
centralized to avoid behavior drift. Extract the shared strand-aware raw-BAM
comparison into a common helper that both dedup and group call, preserving the
existing ref-id, unclipped 5' position, and reverse-strand tie-break behavior.
Update the callers to use the shared helper and keep the symbols
is_r1_genomically_earlier_raw and the corresponding group-side function aligned
through that shared implementation.
🪄 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: f3fbd800-bd12-4dae-b47e-383b50e1a1b8
📒 Files selected for processing (3)
src/lib/commands/compare/bams.rssrc/lib/commands/dedup.rssrc/lib/commands/group.rs
fgumi's `paired` grouping strategy could split a single duplex molecule into two when a template's two mates share an unclipped 5' coordinate (fully overlapping / short-insert fragments). `is_r1_genomically_earlier_raw` broke the position tie with `r1_pos <= r2_pos`, which returns true for both the top (R1 forward) and bottom (R1 reverse) strand. That assigns the lower/higher paired-UMI prefix inconsistently between the two strands, so their prefixed UMIs are no longer reverses of each other and fail to pair. Mirror fgbio `GroupReadsByUmi.umiForRead`'s tie-break (`pos1 == pos2 && r1.positiveStrand`): on a tie, R1 is earlier iff it is on the forward strand. Applied to both `group` and `dedup`, which shared the buggy helper. Also fix `compare bams` grouping mode: it keyed only on the full MI (base + `/A`|`/B`) and so reported such a split as EQUIVALENT (the per-MiKey mapping stays a consistent bijection). Add a base-level, strand-suffix- stripped molecule-membership check so a strand-pairing difference makes the groupings DIFFER; a molecule relabel or an `/A`<->`/B` swap stays EQUIVALENT. Tests: add a regression test for the overlapping-fragment tie case; repair `test_paired_assigner_explicit_ab_ba_symmetry` (the port used an FR-only builder and asserted two groups, inverting fgbio's contract) to build a real reverse-strand bottom read and assert one molecule with `/A` + `/B`; add unit tests for the compare base-level check.
5db2e1a to
d9e623f
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
This branch was previously deployed
1 participant
Summary
fgumi's
pairedgrouping strategy can split a single duplex molecule into two when a template's two mates share an unclipped 5' coordinate — i.e. fully-overlapping / short-insert fragments. This diverges from fgbioGroupReadsByUmi, which pairs the two strands into one molecule.Root cause
Duplex strand pairing prefixes each half of the paired UMI with a lower- vs higher-coordinate-read marker (
lowerReadUmiPrefix/higherReadUmiPrefix) so the top strand's prefixed UMI is the exact reverse of the bottom strand's. That prefix is chosen byis_r1_genomically_earlier_raw, which broke the position tie withr1_pos <= r2_pos. On a tie (r1_pos == r2_pos) that returnstruefor both strands (top = R1 forward, bottom = R1 reverse), so both get the lower prefix on R1's half, their prefixed UMIs are no longer reverses of each other, they fail to reverse-match in the assigner, and the molecule is split in two.fgbio breaks the same tie on strand (
GroupReadsByUmi.umiForRead):pos1 < pos2 || (pos1 == pos2 && r1.positiveStrand). This PR mirrors that — on a tie, R1 is "earlier" iff it is on the forward strand. The buggy helper was shared bygroupanddedup; both are fixed.compare bamswas blind to itThe grouping-mode comparison keyed only on the full MI (
base+/A|/B), treating the two strands as independent groups. A strand-pairing split keeps the per-MiKey mapping a consistent bijection (X/A↔X/A,X/B↔Y/A), so it reported the split asEQUIVALENTwith zero mismatches — which is why the benchmark's owngroup.pairedequivalence check never flagged this. This PR adds a base-level (strand-suffix-stripped) molecule-membership check: reads that share a molecule in one BAM must share a molecule in the other. A genuine strand-pairing split now makes the groupingsDIFFER; a molecule relabel or an/A↔/Bswap staysEQUIVALENT.Tests
test_paired_assigner_pairs_overlapping_duplex_strands: builds the two strands of one duplex molecule at equal unclipped 5' positions and asserts they group into one molecule. Verified RED on the pre-fix code (["0/A","0/A","1/A","1/A"], two bases), GREEN after.test_paired_assigner_explicit_ab_ba_symmetry: the port of fgbio's "correctly group reads with the paired assigner when the two UMIs are the same" used an FR-only builder (so it never built the reverse-strand bottom read) and asserted two separate groups, inverting fgbio's contract. It now builds a real reverse-strand bottom read and asserts one molecule with/A+/B, matching fgbio./A↔/Bswap, and single-strand relabel not flagged).Full suite green (2215 tests); fmt and clippy (pedantic) clean.
Validation on real data
Reproduced end-to-end on the agilent-hs2 vendor sample (SRR30485713): fgumi's own
group → duplex → filterproduced 1,025,610 filtered records vs fgbio's 1,025,612 — a single degenerate duplex molecule (chr11:119116712; 8 top-strand + 3 bottom-strand templates of a fully-overlapping fragment) that fgumi split and fgbio kept. After the fix fgumi produces 1,025,612, byte-identical to fgbio. The behavior is present onmain(0.4.0), so it is not a recent regression. The fixedcompare bamsflags exactly that one molecule out of 10,220,662 (Grouping mismatches: 1,DIFFER) where the old tool reportedEQUIVALENT.Summary by CodeRabbit