fix(consensus): apply --max-reads-per-strand cap, share combined-error Short cap, diagnose mi5252 (DUPLEX3-01/04)

[nh13]

Summary

Three consensus findings from the final-audit burn-down (W10d):

• DUPLEX3-01 — --max-reads-per-strand was a silent no-op on the duplex/codec single-strand path.
• DUPLEX3-04 — a disabled real-data test hid an mi5252 C-vs-N output divergence; investigated, root-caused as a near-tie floating-point ordering artifact (no consensus-logic change), and replaced with a precise writeup.
• Combined-error Short cap — the CODEC caller's per-base combined duplex error was summed in u16 (a defensive overflow concern); introduce one shared clamp_combined_error_to_fgbio_short helper and route every combine site (codec + duplex) through it.

Merge status: rebased onto main — the base consensus PR #552 (nh/fix-consensus-depth-clamp-parity) is now merged, so this branch targets main directly and inherits #552's fgbio_oracle_saturation_tests.

DUPLEX3-01 — --max-reads-per-strand was a silent no-op on the duplex path

fgbio caps the reads contributing to a single-strand consensus inside consensusCall (shuffle.take(maxReads), VanillaUmiConsensusCaller.scala:288). fgumi's duplex (and codec) callers reach the single-strand consensus through VanillaUmiConsensusCaller::consensus_call, which never downsampled — so with the flag set fgumi used every read and emitted higher per-strand depth/quality than fgbio while the cap was silently ignored. (The vanilla/simplex path caps raw records earlier in process_group and does not route through consensus_call, so it was never affected.)

The fix applies the cap inside consensus_call via a new downsample_source_reads helper (the SourceRead analogue of the existing downsample_reads), shuffling with the caller's seeded RNG and truncating to max_reads — exactly where fgbio caps. Two invariants preserved:

• Default output stays byte-identical — when the flag is unset the cap is a no-op.
• The full uncapped source_reads are retained on the output — the cap shapes only the single-strand consensus bases/quals/depths; fgbio passes the pre-cap filteredAbR1s ++ filteredBaR2s to duplexConsensus, so the duplex caller (which counts per-base errors and calls the consensus UMI over output.source_reads) must still see every read. Capping it would undercount duplex errors and bias the cE/error-rate tags.

A degenerate cap (max_reads = 0, or below min_reads) now returns no consensus instead of panicking / hitting the empty-source bail.

Reproduction (before → after)

On a synthetic deep duplex grouped BAM (per-strand read pairs up to 19), reading the per-strand depth tags aD/bD:

Run	max aD	max bD
fgumi, no flag	15	19
fgumi --max-reads-per-strand 3	3	3
fgbio --max-reads-per-strand 3	3	3

Before the fix the flag left aD/bD at full depth (the added unit test failed asserting depth 8 with max_reads=3); after the fix per-strand depth is capped, matching fgbio's contract. fgbio uses a random shuffle.take over java.util.Random (an LCG) while fgumi uses StdRng (ChaCha), so exact base/qual parity under downsampling is a documented, intentional divergence (see downsample_source_reads) — the contract is "≤ N reads per strand contribute", and determinism-per-seed is pinned in tests.

Combined-error Short cap — single-source the combine-step error clamp

The CODEC caller's per-base combined duplex error was ea + eb (and the disagreement-branch sums) on u16. This is a defensive overflow concern, not a reachable production bug: both the codec and duplex callers build their single-strand consensuses through VanillaUmiConsensusCaller::consensus_call, which already caps each per-base depth and error at fgbio's Short ceiling (32767) at push (vanilla_caller.rs:1413/1417). So every strand term feeding a combine step is ≤ 32767 and the combined sum is ≤ 65534 (fits u16) for real pipeline data. Still, a raw ea + eb on u16 would panic in debug (and wrap in release) on any above-ceiling strand value, so the combine step should saturate rather than silently assume the upstream invariant.

Introduce one shared clamp_combined_error_to_fgbio_short helper in caller.rs (beside clamp_per_base_to_fgbio_short) — sum the two strands' per-base errors in a wide signed type, then saturate to [0, 32767], matching fgbio's per-base errors Array[Short] — and route every combine site through it:

• CODEC build_duplex_consensus_from_padded — agreement, both disagreement branches, and the neither-has-data branch (previously an uncapped u16 sum).
• Duplex build_duplex_consensus (source-read + approximate methods) and the call_duplex_from_ss_pair test helper — previously an inline .clamp(0, i16::MAX), now the shared helper (behavior-identical, confirmed by the fgbio_oracle_saturation_tests).

The simplex/vanilla caller has no combine step; it caps each per-base error at push, which is the upstream invariant this helper defends. The downstream cE numerator re-clamps to the same ceiling, so this changes only the stored intermediate, never the emitted tag.

DUPLEX3-04 — disabled test hid an output divergence: investigated, root-caused, documented (no consensus-logic change)

test_mi5252_real_data was #[ignore]/commented ("bam::Reader issue") and recorded fgbio→C vs fgumi→N at R1 position 111 on real data. Investigation outcome: the divergence is a near-tie floating-point ordering artifact, not a consensus bug, so the consensus logic is deliberately left unchanged.

• fgumi's tie rule is faithful to fgbio's: a base is a no-call (N) when a competing likelihood is within one machine epsilon of the maximum. fgumi uses abs_diff_eq!(ll, max, epsilon = f64::EPSILON) (base_builder.rs::call); fgbio uses MathUtil.maxWithIndex(requireUniqueMaximum = true, epsilon = 1/2^52) (ConsensusCaller.call). Both epsilons are 2^-52 and both let a strictly-greater value clear the tie — same semantics, same order-sensitivity.
• At position 111 the two contending bases' likelihoods are equal to within ~epsilon. Whether one wins uniquely (C) or is flagged a within-epsilon tie (N) depends on the order in which per-base likelihoods are accumulated. fgumi accumulates with SIMD-vectorized Kahan summation over reads ordered by its port of fgbio's filterToMostCommonAlignment; fgbio accumulates with scalar Kahan summation over its own ordering. The two differ in the last ULP — enough to flip a unique max into a tie at this locus.
• This is not a correctness contract: fgbio PR #1120 (Kahan summation, shipped in fgbio ≥ 3.1.1 and therefore in the 4.1.0 used for the one-time manual comparison) exists precisely because these near-ties are numerically unstable. fgumi's N (no-call at a genuine tie) is the conservative, defensible outcome. Changing the SIMD summation order, the epsilon, or the read ordering to match fgbio at this one locus would be a speculative edit to shared consensus scoring that risks the unit tests pinning exact fgbio numbers elsewhere.

The dead test is removed (obsolete bam::Reader API, a committed-BAM dependency that no longer exists, and it violated the generate-test-data-programmatically convention) and replaced with a precise root-cause writeup. The live synthetic coverage of the same near-tie ordering mechanism is test_tie_breaking_for_simplex_consensus, whose stale, self-contradictory doc comment (claimed N, asserted A) is corrected to match reality.

Tests

• RED-first test_consensus_call_caps_reads_per_strand — fails at depth 8 before the DUPLEX3-01 fix, passes at depth 3 after.
• test_duplex_per_base_error_caps_before_summing — a defensive regression feeding an above-ceiling (errors = [40000, 40000]) single-strand pair (a contract-level input the real pipeline does not produce); panics on debug overflow with a raw u16 sum, asserts [32767, 32767] with the cap (verified by temporarily reverting it).
• test_consensus_call_max_reads_boundaries (rstest): cap-above / cap-equals / cap-below / no-cap / zero-cap boundaries, including the degenerate max_reads = 0 no-crash path.
• consensus_call_downsampling_is_deterministic_for_a_fixed_seed + consensus_call_cap_invariants (proptest): determinism-per-seed and the cap/full-source-retention invariants across arbitrary sizes/caps/min_reads.
• Inherited from fix(consensus)!: clamp codec/duplex scalar depth+error tags to fgbio's Short ceiling #552 via the rebase: fgbio_oracle_saturation_tests (codec + duplex) pin depth/error tags against values captured from a real fgbio 4.0.1 run at the 32767/65534 saturation boundary and an open-interval mixed-strand case.
• Full workspace suite green: 4609 tests pass, 26 skipped. cargo ci-fmt + cargo ci-lint clean.

Summary by CodeRabbit

• Bug Fixes
  • Improved handling of very high per-base coverage to prevent overflow and ensure depth and error-rate metrics remain accurate.
  • Aligned depth and error calculations with fgbio’s saturation limits for duplex and single-strand consensus results.
  • Applied read-count limits during consensus scoring while preserving complete source-read information for downstream analysis.
  • Improved deterministic downsampling behavior and handling of minimum and maximum read-count boundaries.
• Tests
  • Added regression coverage for high-depth saturation, overflow prevention, and deterministic consensus behavior.

[coderabbitai]

Walkthrough

The change aligns depth and error aggregation with fgbio Short saturation semantics across consensus callers, replaces overflowing per-base casts, and applies deterministic max_reads downsampling only to vanilla consensus scoring while preserving full source reads.

Changes

Consensus calling updates

Layer / File(s)	Summary
Short clamp contract crates/fgumi-consensus/src/caller.rs, crates/fgumi-consensus/Cargo.toml	Adds and tests the shared fgbio Short depth clamp, with development dependencies reordered.
Duplex metric saturation crates/fgumi-consensus/src/codec_caller.rs, crates/fgumi-consensus/src/duplex_caller.rs	Caps per-base depths and errors before duplex construction, scalar tag aggregation, and optional per-base SAM tag encoding; adds saturation and overflow regression tests.
Vanilla consensus read cap crates/fgumi-consensus/src/vanilla_caller.rs	Downsamples the consensus-scoring subset deterministically after annotation and normalization, retains uncapped source reads, and tests boundary, determinism, and property-based invariants.

Estimated code review effort: 4 (Complex) | ~45 minutes

Sequence Diagram(s)

sequenceDiagram
  participant consensus_call
  participant downsample_source_reads
  participant create_consensus_from_source_reads
  participant VanillaConsensusRead
  consensus_call->>downsample_source_reads: select max_reads scoring subset
  downsample_source_reads-->>consensus_call: deterministic retained reads
  consensus_call->>create_consensus_from_source_reads: compute consensus from subset
  create_consensus_from_source_reads-->>consensus_call: consensus metrics
  consensus_call->>VanillaConsensusRead: retain full source_reads with capped consensus

Loading

Possibly related PRs

• fulcrumgenomics/fgumi#552: Related fgbio Short ceiling handling for codec and duplex depth and error tags.

🚥 Pre-merge checks | ✅ 5

✅ Passed checks (5 passed)

Check name	Status	Explanation
Docstring Coverage	✅ Passed	No functions found in the changed files to evaluate docstring coverage. Skipping docstring coverage check.
Linked Issues check	✅ Passed	Check skipped because no linked issues were found for this pull request.
Out of Scope Changes check	✅ Passed	Check skipped because no linked issues were found for this pull request.
Description Check	✅ Passed	Check skipped - CodeRabbit’s high-level summary is enabled.
Title check	✅ Passed	The title matches the main changes: max-reads-per-strand downsampling, Short-cap saturation, and mi5252 investigation.

✨ Finishing Touches

📝 Generate docstrings

• Create stacked PR
• Commit on current branch

🧪 Generate unit tests (beta)

• Create PR with unit tests
• Commit unit tests in branch nh/fix-duplex-maxreads-and-mi5252

---

_{Comment @coderabbitai help to get the list of available commands.}

[codecov]

Codecov Report

❌ Patch coverage is 97.61905% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 92.87%. Comparing base (efdceef) to head (40f8450).
⚠️ Report is 2 commits behind head on main.

Files with missing lines	Patch %	Lines
crates/fgumi-consensus/src/codec_caller.rs	92.59%	2 Missing ⚠️
crates/fgumi-consensus/src/vanilla_caller.rs	98.91%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #563      +/-   ##
==========================================
- Coverage   92.89%   92.87%   -0.03%     
==========================================
  Files         167      167              
  Lines      102777   102895     +118     
==========================================
+ Hits        95479    95563      +84     
- Misses       7298     7332      +34     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[nh13]

@coderabbitai review

[coderabbitai]

✅ Action performed

Review finished.

Note: CodeRabbit is an incremental review system and does not re-review already reviewed commits. This command is applicable only when automatic reviews are paused.

[coderabbitai]

Actionable comments posted: 3

🤖 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 `@crates/fgumi-consensus/src/codec_caller.rs`:
- Around line 3224-3330: Add programmatically generated fgbio-baseline
regression coverage at crates/fgumi-consensus/src/codec_caller.rs lines
3224-3330 and crates/fgumi-consensus/src/duplex_caller.rs lines 4991-5139.
Extend the existing codec and duplex tests to compare scalar depth/error tags
and capped per-base arrays against generated fgbio output for values 32767,
32768, and 65535, asserting identity or documenting any intentional divergence;
the requested changes apply at both sites.
- Around line 1179-1188: Update the duplex_error construction in the caller
logic so ea + eb and disagreement-branch additions occur in a wider integer
type, then saturate the intermediate result to fgbio’s Short ceiling before
converting to the final type. Preserve the existing cE/error-cap behavior while
preventing debug overflow and release wrapping for large u16 inputs; leave the
duplex_depth calculation unchanged.

In `@crates/fgumi-consensus/src/vanilla_caller.rs`:
- Around line 822-838: Add coverage for downsampling in downsample_source_reads
by either implementing an fgbio-compatible shuffle or creating a programmatic
differential test that compares capped-family consensus bases, qualities, and
depths against generated fgbio output. Ensure the test explicitly documents and
asserts the intended divergence if the existing StdRng behavior remains,
including cases where max_reads truncates the source reads.

🪄 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: 71635712-9e61-42a5-b049-2f43fa3f5bb6

📥 Commits

Reviewing files that changed from the base of the PR and between efdceef and 9630ad9.

📒 Files selected for processing (5)

• crates/fgumi-consensus/Cargo.toml
• crates/fgumi-consensus/src/caller.rs
• crates/fgumi-consensus/src/codec_caller.rs
• crates/fgumi-consensus/src/duplex_caller.rs
• crates/fgumi-consensus/src/vanilla_caller.rs