Repository navigation
fix(consensus): clamp consensus depth at i16::MAX instead of wrapping u16 - #448
Conversation
|
Important Review skippedNo new commits to review since the last review. ⚙️ Run configurationConfiguration used: Path: .coderabbit.yaml Review profile: ASSERTIVE Plan: Pro Run ID: You can disable this status message by setting the Use the checkbox below for a quick retry:
WalkthroughPer-base observation counters in ChangesDeep UMI family overflow fix
Estimated code review effort🎯 2 (Simple) | ⏱️ ~10 minutes Possibly related PRs
Suggested labels
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✏️ Tip: You can configure your own custom pre-merge checks in the settings. ✨ Finishing Touches🧪 Generate unit tests (beta)
Comment |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #448 +/- ##
==========================================
- Coverage 91.12% 91.05% -0.08%
==========================================
Files 78 78
Lines 51414 51324 -90
==========================================
- Hits 46851 46732 -119
- Misses 4563 4592 +29 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
fe6f0c2 to
7f0d170
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
… u16 The per-base observation counter in `ConsensusBaseBuilder` was a `u16`, and the per-base consensus depth and error counts were stored as `u16`. On very deep UMI families — high-duplication cfDNA can exceed 65_535 reads at a single position — the observation counter wrapped modulo 65_536 (e.g. a 138_656-read family reported depth 7_584). The wrapped depth corrupts the cD/cM/cE consensus tags and, because cE = errors / depth, can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. fgbio's VanillaUmiConsensusCaller counts depth in `Int` and clamps the per-base depth/error to `Short.MaxValue` (32_767), deriving cD/cM/cE from the clamped arrays. Mirror that exactly: - Widen `ConsensusBaseBuilder::observations` (and `contributions()`, `observations_for_base()`, `all_observations()`) from `u16` to `u32` so the counter never wraps. - Clamp the per-base depth and error to `i16::MAX` when storing them, so cD/cM (max/min depth), cE (errors / depth) and the per-base cd/ce arrays all match fgbio. The unclamped depth is still used for the min-reads threshold, as fgbio does (VanillaUmiConsensusCaller.scala:329). Duplex needs no change: it builds per-strand consensus via the vanilla caller (now corrected) and sums the two clamped strand depths exactly as fgbio does (DuplexConsensusCaller.scala:480), without re-clamping the combined value. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux types. It previously handled only `c`/`s`/`i`, so a preceding tag written as `S` (which is how `cD`=32_767 now serializes) made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Tests: a deep family (70_000 observations) no longer wraps `contributions()`, and a 40_000-read family reports cD=cM=32_767 with cE computed from the clamped depth — matching fgbio.
7f0d170 to
4702609
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test.
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test.
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test.
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test.
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test.
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test.
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test.
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test. Cover the error half of the clamp. The depth test has one disagreeing observation, so errors.push(error_count.min(max_short)) was never exercised -- reverting just that .min() left it green. Saturating errors needs depth above 2 * i16::MAX, since the error count is depth minus the majority base's observations: 70,000 reads split 36,000 A / 34,000 C at position 0 puts 34,000 errors there. The oracle is fgbio's algorithm rather than fgumi's output -- summing the clamped arrays gives cE = 32,767 / (32,767 * 4) = 0.25 exactly, where an unclamped stored error gives 34,000 / 131,068 = 0.2594. The per-base ce array is clamped again when the tag is written, so the summary cE is the only tag that distinguishes the two; verified the test fails at 0.2594 with the error clamp reverted. Bounds-check the aux scan in testutil. find_i16_array_tag_in_aux advanced its cursor and indexed scalar, B-array, and Z/H payloads without checking them against the end of the block, so truncated or malformed aux bytes panicked instead of returning None. Every step now goes through a checked advance, an unknown B subtype fails closed rather than resyncing the cursor on a guessed element width, and an unterminated Z/H string is malformed rather than an implicit end-of-block. Adds a case table over the truncation shapes.
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test. Cover the error half of the clamp. The depth test has one disagreeing observation, so errors.push(error_count.min(max_short)) was never exercised -- reverting just that .min() left it green. Saturating errors needs depth above 2 * i16::MAX, since the error count is depth minus the majority base's observations: 70,000 reads split 36,000 A / 34,000 C at position 0 puts 34,000 errors there. The oracle is fgbio's algorithm rather than fgumi's output -- summing the clamped arrays gives cE = 32,767 / (32,767 * 4) = 0.25 exactly, where an unclamped stored error gives 34,000 / 131,068 = 0.2594. The per-base ce array is clamped again when the tag is written, so the summary cE is the only tag that distinguishes the two; verified the test fails at 0.2594 with the error clamp reverted. Bounds-check the aux scan in testutil. find_i16_array_tag_in_aux advanced its cursor and indexed scalar, B-array, and Z/H payloads without checking them against the end of the block, so truncated or malformed aux bytes panicked instead of returning None. Every step now goes through a checked advance, an unknown B subtype fails closed rather than resyncing the cursor on a guessed element width, and an unterminated Z/H string is malformed rather than an implicit end-of-block. Adds a case table over the truncation shapes. Assert the saturation tests against a programmatic fgbio oracle rather than hand-computed constants. fgbio_summary_tags_baseline reimplements fgbio's derivation -- store each per-position depth and error as a Short, then take cD as the max depth, cM as the min, and cE as sum(errors)/sum(depths) over those clamped arrays -- with no fgbio dependency, the same shape as fgbio_edit_baseline_partition in fgumi-umi. Both saturation tests now compare against the oracle applied to their fixture's known pileup. Re-verified non-vacuous: reverting the error clamp still fails at 0.2594 against the oracle's 0.25. Reject --min-reads above the depth ceiling. The threshold is compared against a u16 per-position depth that saturates at u16::MAX, so any larger value can never be met and would silently no-call every position of an arbitrarily deep family as NotEnoughReads. Simplex is the only reachable path -- duplex and codec pin their single-strand caller's min_reads to 1, so their own min_reads is a strand-count filter rather than this threshold. u16::MAX itself stays accepted; it is exactly reachable. Add proptest coverage for the aux scanner. The case table pins the truncation shapes worth naming; the property tests cover the ones nobody thought of -- arbitrary byte sequences must never panic and must not yield an array larger than the block it came from, and well-formed B:s arrays must still round-trip, so the bounds checks cannot drift into rejecting valid input.
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test. Cover the error half of the clamp. The depth test has one disagreeing observation, so errors.push(error_count.min(max_short)) was never exercised -- reverting just that .min() left it green. Saturating errors needs depth above 2 * i16::MAX, since the error count is depth minus the majority base's observations: 70,000 reads split 36,000 A / 34,000 C at position 0 puts 34,000 errors there. The oracle is fgbio's algorithm rather than fgumi's output -- summing the clamped arrays gives cE = 32,767 / (32,767 * 4) = 0.25 exactly, where an unclamped stored error gives 34,000 / 131,068 = 0.2594. The per-base ce array is clamped again when the tag is written, so the summary cE is the only tag that distinguishes the two; verified the test fails at 0.2594 with the error clamp reverted. Bounds-check the aux scan in testutil. find_i16_array_tag_in_aux advanced its cursor and indexed scalar, B-array, and Z/H payloads without checking them against the end of the block, so truncated or malformed aux bytes panicked instead of returning None. Every step now goes through a checked advance, an unknown B subtype fails closed rather than resyncing the cursor on a guessed element width, and an unterminated Z/H string is malformed rather than an implicit end-of-block. Adds a case table over the truncation shapes. Assert the saturation tests against a programmatic fgbio oracle rather than hand-computed constants. fgbio_summary_tags_baseline reimplements fgbio's derivation -- store each per-position depth and error as a Short, then take cD as the max depth, cM as the min, and cE as sum(errors)/sum(depths) over those clamped arrays -- with no fgbio dependency, the same shape as fgbio_edit_baseline_partition in fgumi-umi. Both saturation tests now compare against the oracle applied to their fixture's known pileup. Re-verified non-vacuous: reverting the error clamp still fails at 0.2594 against the oracle's 0.25. Reject --min-reads above the depth ceiling. The threshold is compared against a u16 per-position depth that saturates at u16::MAX, so any larger value can never be met and would silently no-call every position of an arbitrarily deep family as NotEnoughReads. Simplex is the only reachable path -- duplex and codec pin their single-strand caller's min_reads to 1, so their own min_reads is a strand-count filter rather than this threshold. u16::MAX itself stays accepted; it is exactly reachable. Add proptest coverage for the aux scanner. The case table pins the truncation shapes worth naming; the property tests cover the ones nobody thought of -- arbitrary byte sequences must never panic and must not yield an array larger than the block it came from, and well-formed B:s arrays must still round-trip, so the bounds checks cannot drift into rejecting valid input. # Conflicts: # src/lib/commands/simplex.rs
…den #448) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test. Cover the error half of the clamp. The depth test has one disagreeing observation, so errors.push(error_count.min(max_short)) was never exercised -- reverting just that .min() left it green. Saturating errors needs depth above 2 * i16::MAX, since the error count is depth minus the majority base's observations: 70,000 reads split 36,000 A / 34,000 C at position 0 puts 34,000 errors there. The oracle is fgbio's algorithm rather than fgumi's output -- summing the clamped arrays gives cE = 32,767 / (32,767 * 4) = 0.25 exactly, where an unclamped stored error gives 34,000 / 131,068 = 0.2594. The per-base ce array is clamped again when the tag is written, so the summary cE is the only tag that distinguishes the two; verified the test fails at 0.2594 with the error clamp reverted. Bounds-check the aux scan in testutil. find_i16_array_tag_in_aux advanced its cursor and indexed scalar, B-array, and Z/H payloads without checking them against the end of the block, so truncated or malformed aux bytes panicked instead of returning None. Every step now goes through a checked advance, an unknown B subtype fails closed rather than resyncing the cursor on a guessed element width, and an unterminated Z/H string is malformed rather than an implicit end-of-block. Adds a case table over the truncation shapes. Assert the saturation tests against a programmatic fgbio oracle rather than hand-computed constants. fgbio_summary_tags_baseline reimplements fgbio's derivation -- store each per-position depth and error as a Short, then take cD as the max depth, cM as the min, and cE as sum(errors)/sum(depths) over those clamped arrays -- with no fgbio dependency, the same shape as fgbio_edit_baseline_partition in fgumi-umi. Both saturation tests now compare against the oracle applied to their fixture's known pileup. Re-verified non-vacuous: reverting the error clamp still fails at 0.2594 against the oracle's 0.25. Reject --min-reads above the depth ceiling. The threshold is compared against a u16 per-position depth that saturates at u16::MAX, so any larger value can never be met and would silently no-call every position of an arbitrarily deep family as NotEnoughReads. Simplex is the only reachable path -- duplex and codec pin their single-strand caller's min_reads to 1, so their own min_reads is a strand-count filter rather than this threshold. u16::MAX itself stays accepted; it is exactly reachable. Add proptest coverage for the aux scanner. The case table pins the truncation shapes worth naming; the property tests cover the ones nobody thought of -- arbitrary byte sequences must never panic and must not yield an array larger than the block it came from, and well-formed B:s arrays must still round-trip, so the bounds checks cannot drift into rejecting valid input. # Conflicts: # src/lib/commands/simplex.rs
…den #448) (#555) feat-runall already avoids the u16 wrap by saturating the observation counter at u16::MAX (rather than #448's u32 widening) and clamps the per-base cd/ce arrays to i16::MAX. But it computes the SUMMARY tags cD/cM (max/min depth) and cE (errors/depth) from the UNCLAMPED depths/errors vectors, so on a UMI family deeper than i16::MAX (32_767) at a position the summary tags diverge from fgbio (which stores depths/errors as `Short` and derives cD/cM/cE from those clamped arrays). High-duplication cfDNA can exceed 32_767 reads at a single position, and because cE = errors / depth a divergent denominator can push a read's error rate past the filter threshold and silently change a keep/drop decision relative to fgbio. Clamp at push so the depths/errors vectors — and everything derived from them (the summary cD/cM/cE tags, the per-base cd/ce tags, and duplex's per-strand tags + combined strand-sum) — match fgbio. The unclamped `depth` local is still used for the min-reads threshold (as fgbio does, VanillaUmiConsensusCaller.scala:329); min-reads keys off source_reads.len(), so it is unaffected either way. Duplex needs no separate change: it derives its per-strand aD/aM/bD/bM and cd/ce arrays from the vanilla `ConsensusResult` depths (now clamped) and sums the two clamped strand depths into the combined cD/cM without re-clamping, exactly as fgbio's DuplexConsensusCaller.scala:480. Also completes the test-only aux-tag scanner in `testutil` (`find_i16_array_tag_in_aux`) to skip the unsigned `C`/`S`/`I` (and `A`/`H`) BAM aux scalar types — it previously handled only `c`/`s`/`i`/`f`, so a preceding tag written as an unsigned scalar made the scanner stop before reaching the per-base array tag. Surfaced by the new test. Adds a deep-family (40_000-read) regression test asserting cD=cM=32_767 with cE computed from the clamped depth — verified to fail before the fix ("cD must saturate at i16::MAX"). Ported from the merged #448 test. Cover the error half of the clamp. The depth test has one disagreeing observation, so errors.push(error_count.min(max_short)) was never exercised -- reverting just that .min() left it green. Saturating errors needs depth above 2 * i16::MAX, since the error count is depth minus the majority base's observations: 70,000 reads split 36,000 A / 34,000 C at position 0 puts 34,000 errors there. The oracle is fgbio's algorithm rather than fgumi's output -- summing the clamped arrays gives cE = 32,767 / (32,767 * 4) = 0.25 exactly, where an unclamped stored error gives 34,000 / 131,068 = 0.2594. The per-base ce array is clamped again when the tag is written, so the summary cE is the only tag that distinguishes the two; verified the test fails at 0.2594 with the error clamp reverted. Bounds-check the aux scan in testutil. find_i16_array_tag_in_aux advanced its cursor and indexed scalar, B-array, and Z/H payloads without checking them against the end of the block, so truncated or malformed aux bytes panicked instead of returning None. Every step now goes through a checked advance, an unknown B subtype fails closed rather than resyncing the cursor on a guessed element width, and an unterminated Z/H string is malformed rather than an implicit end-of-block. Adds a case table over the truncation shapes. Assert the saturation tests against a programmatic fgbio oracle rather than hand-computed constants. fgbio_summary_tags_baseline reimplements fgbio's derivation -- store each per-position depth and error as a Short, then take cD as the max depth, cM as the min, and cE as sum(errors)/sum(depths) over those clamped arrays -- with no fgbio dependency, the same shape as fgbio_edit_baseline_partition in fgumi-umi. Both saturation tests now compare against the oracle applied to their fixture's known pileup. Re-verified non-vacuous: reverting the error clamp still fails at 0.2594 against the oracle's 0.25. Reject --min-reads above the depth ceiling. The threshold is compared against a u16 per-position depth that saturates at u16::MAX, so any larger value can never be met and would silently no-call every position of an arbitrarily deep family as NotEnoughReads. Simplex is the only reachable path -- duplex and codec pin their single-strand caller's min_reads to 1, so their own min_reads is a strand-count filter rather than this threshold. u16::MAX itself stays accepted; it is exactly reachable. Add proptest coverage for the aux scanner. The case table pins the truncation shapes worth naming; the property tests cover the ones nobody thought of -- arbitrary byte sequences must never panic and must not yield an array larger than the block it came from, and well-formed B:s arrays must still round-trip, so the bounds checks cannot drift into rejecting valid input. # Conflicts: # src/lib/commands/simplex.rs
Summary
On very deep UMI families, fgumi silently overflowed the consensus depth counter, corrupting the
cD/cM/cEconsensus tags and — because the error rate iserrors / depth— changing a real filter keep/drop decision relative to fgbio. This makes fgumi clamp consensus depth ati16::MAX(32_767) exactly as fgbio does.How it was found
Benchmarking fgumi against fgbio 4.0.0 on a deep cfDNA panel (
idt-cfdna), the filter-stage BAM equivalency checks differed for every consensus strategy (simplex edit/identity/adjacency and duplex). The moleculeMI:40915has 277,312 raw reads (~138,656 per read-end). fgumi reportedcD = 7_584, which is exactly138_656 mod 65_536— au16wraparound. fgbio reportedcD = 32_767(Short.MaxValue).Root cause
ConsensusBaseBuilder::observationswas[u16; 4], incremented with+= 1, and the per-base depth/error vectors wereVec<u16>. A single position deeper than 65_535 wrapped the counter (a debug build panics with "attempt to add with overflow"; the release build that runs in production wraps silently). The wrapped depth then:cD/cM(max/min depth) and the per-basecdarray;cE(errors / depth), which in theMI:40915case inflated the error rate past--max-read-error-rate, so fgumi dropped a consensus read that fgbio kept (the 52_756 vs 52_758 record-count difference observed in the benchmark).This is not a v0.3.1 regression — the
u16types predate v0.3.0; the path was simply never exercised until a family exceeded 65_535 reads.Fix
fgbio's
VanillaUmiConsensusCallercounts depth inIntand clamps the per-base depth/error toShort.MaxValuewhen storing, derivingcD/cM/cEfrom the clamped arrays (VanillaUmiConsensusCaller.scala:324-340,380-382). This mirrors it:observations/contributions()/observations_for_base()/all_observations()fromu16tou32so the counter never wraps.i16::MAXwhen storing, socD/cM/cEand the per-basecd/cearrays all match fgbio. The unclamped depth is still used for the min-reads threshold, exactly as fgbio does (VanillaUmiConsensusCaller.scala:329).Duplex needs no change. It builds per-strand consensus through the (now-corrected) vanilla caller and sums the two clamped strand depths as
i32without re-clamping — which is exactly what fgbio does (DuplexConsensusCaller.scala:480-489). Its error accumulation already clamps toi16::MAX.Test-helper fix
The change also completes the test-only aux-tag scanner
find_i16_array_tag_in_aux(infgumi-raw-bam/testutil) to skip the unsignedC/S/I(andA/H) BAM aux types. It previously handled onlyc/s/i/f/Z, so whencDserializes asS(which a value of 32_767 now does) the scanner hit_ => breakand never reached the per-base array tag. The production BAM reader handles all aux types; this only affected the test helper, and the new test surfaced it.Tests
base_builder: 70_000 single-base observations no longer wrapcontributions()(was7_584viamod 65_536).vanilla_caller: a 40_000-read family (one disagreeing observation) reportscD = cM = 32_767andcE = 1 / (32_767 * read_len)— i.e. computed from the clamped depth, matching fgbio.cargo ci-test(2,214 tests),cargo ci-fmt, andcargo ci-lint(clippy pedantic,-D warnings) all pass.Summary by CodeRabbit
u16tou32.