Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nh13 merged 2 commits into from
Jul 17, 2026
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9f51b4d
fix(duplex): apply --max-reads-per-strand cap; document mi5252 C-vs-N…
nh13 Jul 10, 2026
40f8450
fix(consensus): share a Short-ceiling cap for combined duplex errors
nh13 Jul 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 22 additions & 0 deletions crates/fgumi-consensus/src/caller.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -359,6 +359,28 @@ pub fn clamp_per_base_to_fgbio_short(depth: u16) -> i32 {
i32::from(depth.min(FGBIO_SHORT_DEPTH_MAX))
}

/// Clamp a combined per-base duplex error count to fgbio's `Short` ceiling, matching
/// fgbio's per-base `errors` `Array[Short]` (capped at 32767 at storage).
///
/// The duplex and codec callers combine two single-strand consensuses per base, summing
/// each strand's per-base error (`ea + eb`, and the disagreement `ea + (db - eb)` /
/// `eb + (da - ea)` variants). Both callers build those strands through
/// `VanillaUmiConsensusCaller::consensus_call`, which already caps each per-base depth and
/// error at the `Short` ceiling at push, so in the real pipeline every term is `<= 32767`
/// and the combined sum is `<= 65534` (fits `u16`). This is therefore a **defensive**
/// bound: it is computed in a wider signed type and saturated so a hypothetical
/// above-ceiling strand value can never overflow the `u16` per-base error store or wrap
/// negative, exactly as fgbio's `Short` array would saturate. The downstream `cE`
/// numerator re-clamps to the same ceiling, so this changes only the stored intermediate,
/// never the emitted tag.
#[must_use]
pub fn clamp_combined_error_to_fgbio_short(error_sum: i64) -> u16 {
// The clamp bounds the value to `[0, 32767]`, which always fits `u16`; the `unwrap_or`
// fallback (the ceiling itself) is unreachable and keeps this panic-free.
u16::try_from(error_sum.clamp(0, i64::from(FGBIO_SHORT_DEPTH_MAX)))
.unwrap_or(FGBIO_SHORT_DEPTH_MAX)
}

/// Reasons why reads might be rejected and not used in consensus calling
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum RejectionReason {
Expand Down
63 changes: 57 additions & 6 deletions crates/fgumi-consensus/src/codec_caller.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,8 @@

use crate::caller::{
ConsensusCaller, ConsensusCallingStats, ConsensusOutput,
RejectionReason as CallerRejectionReason, clamp_per_base_to_fgbio_short,
RejectionReason as CallerRejectionReason, clamp_combined_error_to_fgbio_short,
clamp_per_base_to_fgbio_short,
};
use crate::phred::{MIN_PHRED, NO_CALL_BASE, NO_CALL_BASE_LOWER, PhredScore};
use crate::simple_umi::consensus_umis;
Expand Down Expand Up @@ -1164,17 +1165,22 @@ impl CodecConsensusCaller {
// - When bases agree: sum of both errors
// - When bases disagree: count errors from the chosen strand + non-errors from other strand
// (because the non-errors from the other strand now disagree with the consensus)
// Sum in `i64` (each term is a `u16` per-base count) and saturate at the
// `Short` ceiling via the shared combine-step cap: single-strand errors are
// already Short-capped at push, so this is a defensive bound (a raw `ea + eb`
// on `u16` would overflow if a strand ever exceeded the ceiling).
let duplex_error = if ba == bb {
// Agreement: sum errors from both strands
ea + eb
i64::from(ea) + i64::from(eb)
} else if ba == raw_base {
// We chose A's base (including equal quality disagreements where fgbio uses aBase)
// Count A's errors + B's non-errors (which now disagree with consensus)
ea + db.saturating_sub(eb)
i64::from(ea) + i64::from(db.saturating_sub(eb))
} else {
// We chose B's base: count B's errors + A's non-errors
eb + da.saturating_sub(ea)
i64::from(eb) + i64::from(da.saturating_sub(ea))
};
let duplex_error = clamp_combined_error_to_fgbio_short(duplex_error);

// Per-base duplex depth = sum of each strand's per-base depth,
// each CAPPED at fgbio's `Short` ceiling first (fgbio's
Expand Down Expand Up @@ -1208,8 +1214,10 @@ impl CodecConsensusCaller {
(false, false) => {
// Neither has data - N
// fgbio still calculates errors for these positions:
// since both bases are the same (N == N), errors = a.errors + b.errors
let duplex_error = ea + eb;
// since both bases are the same (N == N), errors = a.errors + b.errors.
// Same shared combine-step cap as the duplex branch above.
let duplex_error =
clamp_combined_error_to_fgbio_short(i64::from(ea) + i64::from(eb));
(NO_CALL_BASE, MIN_PHRED, 0, duplex_error)
}
};
Expand Down Expand Up @@ -3263,6 +3271,49 @@ mod tests {
);
}

/// Defensive regression: the per-base combined duplex ERROR count must sum the two
/// strands' per-base errors in a type wider than `u16` and saturate at fgbio's `Short`
/// ceiling, so a strand value above the ceiling can never overflow the `u16` error store.
/// The single-strand caller already Short-caps per-base errors at push, so this
/// above-ceiling input (`errors = [40000, …]`) is a contract-level state the real
/// pipeline does not produce — but a raw `ea + eb` on `u16` would still panic in debug
/// (and wrap in release) on it. Companion of
/// `test_duplex_per_base_depth_caps_each_strand_before_summing` for the error array.
#[test]
fn test_duplex_per_base_error_caps_before_summing() {
let options = CodecConsensusOptions::default();
let mut caller = CodecConsensusCaller::new("codec".to_string(), "RG1".to_string(), options);

// Both strands agree on both bases (drives the `ea + eb` agreement branch), each
// carrying a per-base error above the `Short` ceiling (a contract-level input; the
// single-strand caller would have capped these to 32767 at push).
let deep = SingleStrandConsensus {
bases: b"AC".to_vec(),
quals: vec![40, 40],
depths: vec![u16::MAX, u16::MAX],
errors: vec![40_000, 40_000], // 40000 + 40000 = 80000 > u16::MAX
raw_read_count: u16::MAX as usize,
ref_start: 0,
ref_end: 1,
is_negative_strand: false,
};
let ss_a = deep.clone();
let ss_b = SingleStrandConsensus { is_negative_strand: true, ..deep };

let duplex = caller
.build_duplex_consensus_from_padded(&ss_a, &ss_b)
.expect("Should build duplex consensus without overflowing the error sum");

// min(40000+40000, 32767) = 32767 per base — capped at the Short ceiling, not the
// (overflowing) raw u16 sum. The downstream cE numerator re-clamps to the same value.
assert_eq!(
duplex.errors,
vec![32_767, 32_767],
"per-base duplex error must be the Short-capped wide-int sum of the strand errors, \
not the raw (overflowing) ea + eb"
);
}

/// The emitted error-rate NUMERATORS (aE/bE/cE) must sum per-base errors capped at fgbio's
/// `Short` ceiling, matching fgbio's capped `errors` `Array[Short]` — not just the depth
/// denominators that #552 already capped. A per-base error above 32767 left uncapped inflates
Expand Down
135 changes: 51 additions & 84 deletions crates/fgumi-consensus/src/duplex_caller.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,8 @@ use noodles::sam::alignment::record_buf::RecordBuf;

use crate::caller::ConsensusOutput;
use crate::caller::{
ConsensusCaller, ConsensusCallingStats, RejectionReason, clamp_per_base_to_fgbio_short,
ConsensusCaller, ConsensusCallingStats, RejectionReason, clamp_combined_error_to_fgbio_short,
clamp_per_base_to_fgbio_short,
};
use crate::phred::MAX_PHRED;
use crate::phred::{MIN_PHRED, PhredScore};
Expand Down Expand Up @@ -955,7 +956,7 @@ impl DuplexConsensusCaller {
num_errors += 1;
}
}
num_errors.clamp(0, i32::from(i16::MAX)) as u16
clamp_combined_error_to_fgbio_short(i64::from(num_errors))
} else {
// Approximate method when source reads unavailable
let a_err = i32::from(a.errors[i]);
Expand All @@ -970,7 +971,7 @@ impl DuplexConsensusCaller {
} else {
b_err + (a_dep - a_err)
};
err.clamp(0, i32::from(i16::MAX)) as u16
clamp_combined_error_to_fgbio_short(i64::from(err))
};

errors.push(error_count);
Expand Down Expand Up @@ -1484,7 +1485,10 @@ impl DuplexConsensusCaller {
num_errors
};

duplex_errors.push(error_at_i.clamp(0, i32::from(i16::MAX)) as i16);
// Cap at fgbio's `Short` ceiling via the shared combine-step helper (the same one
// the codec caller uses). It returns `u16` in `[0, 32767]`, which casts losslessly
// to the `i16` per-base error store fgbio's `Array[Short]` uses.
duplex_errors.push(clamp_combined_error_to_fgbio_short(i64::from(error_at_i)) as i16);
}

// Create duplex consensus record using ss_a as template
Expand Down Expand Up @@ -4040,90 +4044,53 @@ mod tests {
Ok(())
}

/* DISABLED DUE TO bam::Reader ISSUE
#[test]
#[ignore] // Run with: cargo test --release -- --ignored test_mi5252_real_data
fn test_mi5252_real_data() -> Result<()> {
use noodles::sam::alignment::RecordBuf;
use noodles::bam;
use std::fs::File;

// Read the test BAM file containing only MI 5252 reads
let test_bam_path = "/Users/nhomer/work/git/fgumi/debug/test_mi5252_only.bam";
let mut reader = File::open(test_bam_path)
.map(bam::Reader::new)
.expect("Failed to open test BAM file");

let _header = reader.read_header()?;

// Read all records from the BAM file
let mut records = Vec::new();
for result in reader.records() {
let record = result?;
let record_buf: RecordBuf = record.try_into()?;
records.push(record_buf);
}


println!("Read {} records from test BAM", records.len());
assert_eq!(records.len(), 18, "Should have exactly 18 reads for MI 5252");

// Create a duplex consensus caller with minimum thresholds
let mut caller = DuplexConsensusCaller::new(
"test".to_string(),
"RG1".to_string(),
vec![1], // min_reads = 1 (M=1 in fgbio)
20, // min_base_quality
false, // trim
false, // sort_order check
None, // no mask
None, // no error_rate_pre_umi
false, // not producing per-base tags yet
45, // error_rate_post_umi
40, // min_consensus_base_quality
)?;

// Call consensus on these reads
let consensus_records = caller.consensus_reads_from_sam_records(records)?;

println!("Generated {} consensus records", consensus_records.len());
assert_eq!(consensus_records.len(), 2, "Should generate 2 consensus records (R1 and R2)");

// Find the R1 record (flag should indicate first of pair)
let r1 = consensus_records
.iter()
.find(|r| r.flags().is_first_segment())
.expect("Should have R1 consensus");

// Extract the sequence at position 111 (0-based)
let seq: Vec<u8> = r1.sequence().as_ref().to_vec();
assert!(seq.len() > 111, "Sequence should be longer than 111 bases");

println!("R1 consensus base at position 111: {}", seq[111] as char);

// EXPECTED: fgbio outputs 'C' at position 111
// ACTUAL: fgumi currently outputs 'N' at position 111
// This test documents the discrepancy
assert_eq!(
seq[111] as char,
'C',
"Position 111 should be 'C' to match fgbio output (currently fgumi outputs 'N')"
);

Ok(())
}

*/
// DUPLEX3-04 root cause (removed real-data test `test_mi5252_real_data`).
//
// The removed test read a committed real-data BAM
// (`.../debug/test_mi5252_only.bam`, no longer present) through the stale
// `noodles::bam::Reader` API (the "bam::Reader issue" it was commented out for) and asserted
// that fgumi's R1 consensus base at position 111 equalled fgbio's `C`, documenting that fgumi
// instead emitted `N`. It was removed rather than re-enabled because it (a) violates the
// "generate test data programmatically, do not commit BAM fixtures" convention, (b) depended
// on an external absolute path that no longer exists, and (c) used an obsolete reader API.
//
// Investigation outcome: the `C`-vs-`N` 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)` from
// `ConsensusCaller.call`. Both epsilons are 2^-52 and both treat a strictly-greater value
// as clearing the tie, so the semantics (including their order-sensitivity) match.
// * 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 the per-base likelihoods are accumulated. fgumi accumulates with
// SIMD-vectorized Kahan summation (`base_builder.rs`) over reads ordered by its port of
// fgbio's `filterToMostCommonAlignment` (`filter_by_alignment` ->
// `select_most_common_alignment_group`); fgbio accumulates with scalar Kahan summation
// over its own ordering. The two orderings differ in the last ULP, which is 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 pinned fgbio 4.1.0 this divergence was compared against in
// a one-time manual investigation — NOT an automated/CI-enforced baseline) exists
// precisely because these near-ties are numerically unstable. fgumi's `N` (no-call at a genuine tie) is the conservative,
// defensible outcome. Changing 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 ~84 unit tests pinning exact fgbio numbers elsewhere.
//
// Live synthetic coverage of the identical near-tie ordering mechanism is
// `test_tie_breaking_for_simplex_consensus` below.
#[test]
fn test_tie_breaking_for_simplex_consensus() -> Result<()> {
use fgumi_raw_bam::SamBuilder;

// This test reproduces the tie-breaking scenario:
// 8 BA read pairs with 4 having 'A' at position 5 and 4 having 'C' at position 5
// With equal evidence (same quality, same count), this is a true tie.
// With Kahan summation for numeric stability (matching fgbio PR #1120),
// the likelihoods are exactly equal, so we correctly return 'N' (no-call).
// Previously, numerical instability would incorrectly break the tie.
// Near-tie base-calling scenario: 8 BA read pairs, 4 with 'A' and 4 with 'C' at
// position 5, all at equal quality. The two bases' likelihoods are equal to within
// machine epsilon, so the outcome is decided by floating-point accumulation ORDER
// (see the DUPLEX3-04 note above). With reads ordered by descending length (matching
// fgbio's `filterToMostCommonAlignment` behavior) fgumi's Kahan-summed likelihood for
// 'A' ends up strictly (by > epsilon) above 'C', so 'A' wins uniquely rather than
// no-calling. This pins fgumi's deterministic resolution of the near-tie.

// Quality array for 10 bases, all quality 38
let quals = vec![38u8; 10];
Expand Down
Loading
Loading