fix: recompute the BAM bin field after raw POS/CIGAR mutations (clip, zipper) - #591
Conversation
Add a bin module with the canonical SAM spec §5.3 reg2bin, the UNMAPPED_BIN constant, and bin_from_raw_bam/set_bin_from_raw_bam that compute a record's index bin from its POS and CIGAR reference span (unmapped records get 4680). Expose RawRecord::recompute_bin as an ergonomic wrapper. The raw BAM pipelines mutate encoded records in place and emit the bytes verbatim, so unlike the noodles/htsjdk encoders there is no serialization step to refresh the bin field after a POS or CIGAR change. These helpers are the raw-byte equivalent, to be called at each mutation site. Point the builder's private UNMAPPED_BIN at the shared constant.
…ps a read The raw clipper mutated POS (start-clip), the CIGAR (start/end-clip), and unmapped reads whose bases were fully clipped, but never refreshed the index bin (bytes 10-11) — and clip writes records verbatim, so the stale pre-clip bin shipped in the output. fgbio recomputes it (htsjdk refreshes the indexing bin on write), so fgumi diverged for any clip that moved a read into a different bin. Recompute the bin at the three raw mutation primitives (clip_start_of_alignment, clip_end_of_alignment, make_read_unmapped_raw) so a clipped record is always self-consistent, and in set_mate_info_raw when a clipped-unmapped read is relocated to its mapped mate's coordinate or a pair is cleared to unmapped. upgrade_clipping only converts clip type without moving POS or the aligned span, so it needs no recompute.
…unmaps a read Template::fix_mate_info (used by zipper) places an unmapped read at its mapped mate's coordinate and clears both reads to unmapped without touching the index bin, so a relocated read kept the stale unmapped bin (4680) instead of its placed position's bin. fgbio emits the placed bin. Recompute the bin after set_mate_info_one_unmapped relocates the read and after set_mate_info_both_unmapped clears the pair.
region_to_bin duplicated the SAM spec reg2bin reference code. Delegate to fgumi_raw_bam::reg2bin so the simulate encoder and the raw clip/zipper pipelines share one implementation of the binning scheme, keeping the same 1-based-inclusive input contract.
|
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 (8)
WalkthroughChangesBAM bin calculation is centralized in BAM bin recomputation
Estimated code review effort: 3 (Moderate) | ~25 minutes Sequence Diagram(s)sequenceDiagram
participant RawMutation
participant RawRecord
participant BinAPI
participant BAMBytes
RawMutation->>RawRecord: change POS, CIGAR, or mapping state
RawRecord->>BinAPI: recompute current raw record bin
BinAPI->>BAMBytes: write little-endian bin bytes
BAMBytes-->>RawRecord: updated BAM bin
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Comment |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #591 +/- ##
==========================================
- Coverage 93.03% 93.02% -0.01%
==========================================
Files 167 168 +1
Lines 103266 103429 +163
==========================================
+ Hits 96070 96217 +147
- Misses 7196 7212 +16 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
Summary
fgumi clip(and, in a narrower case,fgumi zipper) emitted a stale BAMbinfield. The raw pipelines mutate an already-encoded BAM record in place and write the bytes verbatim, so — unlike the noodles/htsjdk encoders, which recompute the indexing bin on write — there is no serialization step to refreshbin(bytes 10–11) afterPOSor the CIGAR changes. Any clip that moved a read into a different bin shipped the pre-clip value, diverging from fgbio.Reproduction
A read placed at 16300 with
200Mstraddles the 16 kb bin boundary (level-4 bin 585). A 116-base 5′ soft-clip moves it topos 16416,116S84M, whose correct bin is 4682:binpos 16417,116S84Mpos 16417,116S84MThe same class of bug affected
zipper's mate reconciliation: a read relocated to its mapped mate's coordinate kept the stale unmapped bin (4680) instead of its placed position's bin.Note: fgumi's own BAI/CSI writer recomputes start/end from
POS+CIGAR and ignores the record'sbin, so its own index was unaffected. The stale bin hurt external consumers (htsjdk/samtools-based tools reading fgumi output) and broke byte-parity with fgbio.What changed
fgumi-raw-bam::bin(new module): the canonical SAM spec §5.3reg2bin, theUNMAPPED_BINconstant, andbin_from_raw_bam/set_bin_from_raw_bam, which compute a record's bin from itsPOSand CIGAR reference span (unmappedPOS < 0→ 4680, keying off the alignment start like htsjdkcomputeIndexingBin, not the unmapped flag). Plus aRawRecord::recompute_bin()wrapper.clip: recompute the bin at the three raw mutation primitives (clip_start_of_alignment,clip_end_of_alignment,make_read_unmapped_raw) so a clipped record is always self-consistent, and inset_mate_info_rawwhen a clipped-unmapped read is relocated or a pair is cleared to unmapped.upgrade_clippingonly converts clip type (noPOS/span move), so it correctly needs no recompute.zipper(Template::fix_mate_info): recompute after relocating an unmapped mate to its mate's coordinate and after clearing a pair to unmapped.simulate:region_to_binnow delegates to the sharedreg2bininstead of duplicating the spec code; the privateUNMAPPED_BINinbuilder.rsalso shares the constant.Testing
set_mate_info_rawrelocation,fix_mate_inforelocation), each asserting the emitted bin equalsreg2binof the post-op coordinates — TDD, each written failing first.reg2binis covered by the full authoritative vector table ported from htsjdkGenomicIndexUtilTest.testRegionToBinDataProvider, exercising every binning level (16 kb through the whole-window bin 0).reg2binof the post-op coordinates, and directly matches fgbio ClipBam / ZipperBams output bins.cargo ci-fmt, andcargo ci-lintall pass.Reading order
crates/fgumi-raw-bam/src/bin.rs— the shared helper and its semantics.crates/fgumi-sam/src/clipper.rsandsrc/lib/commands/clip.rs— the clip fix.src/lib/template.rs— the zipper fix.src/lib/commands/simulate/mod.rs— the delegation refactor.Summary by CodeRabbit
New Features
Bug Fixes
Tests