diff --git a/crates/fgumi-consensus/src/caller.rs b/crates/fgumi-consensus/src/caller.rs index cb27ec635..16409b469 100644 --- a/crates/fgumi-consensus/src/caller.rs +++ b/crates/fgumi-consensus/src/caller.rs @@ -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 { diff --git a/crates/fgumi-consensus/src/codec_caller.rs b/crates/fgumi-consensus/src/codec_caller.rs index caeb6e15e..2e1bd2cd9 100644 --- a/crates/fgumi-consensus/src/codec_caller.rs +++ b/crates/fgumi-consensus/src/codec_caller.rs @@ -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; @@ -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 @@ -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) } }; @@ -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 diff --git a/crates/fgumi-consensus/src/duplex_caller.rs b/crates/fgumi-consensus/src/duplex_caller.rs index 34acf1323..44e8ccd78 100644 --- a/crates/fgumi-consensus/src/duplex_caller.rs +++ b/crates/fgumi-consensus/src/duplex_caller.rs @@ -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}; @@ -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]); @@ -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); @@ -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 @@ -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 = 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]; diff --git a/crates/fgumi-consensus/src/vanilla_caller.rs b/crates/fgumi-consensus/src/vanilla_caller.rs index 5da066e88..ea761a472 100644 --- a/crates/fgumi-consensus/src/vanilla_caller.rs +++ b/crates/fgumi-consensus/src/vanilla_caller.rs @@ -636,7 +636,8 @@ impl VanillaUmiConsensusCaller { // For EM-seq/TAPs: annotate methylation first (counts conversions), then normalize // source read bases before consensus scoring so that C↔T / G↔A conversion - // events at ref-C positions don't inflate error counts or depress quality. + // events at ref-C positions don't inflate error counts or depress quality. This runs + // over the full (uncapped) set, independent of the per-strand consensus cap below. let (methylation, source_reads) = if self.options.methylation_mode.is_enabled() { let (annot, normalized) = self.annotate_and_normalize(source_reads); (annot, normalized) @@ -644,9 +645,40 @@ impl VanillaUmiConsensusCaller { (None, source_reads) }; - // Build consensus from (possibly normalized) source reads + // Cap the reads contributing to the single-strand CONSENSUS, matching fgbio's + // `consensusCall` (`VanillaUmiConsensusCaller.scala:288`: `shuffle.take(maxReads)`). + // The duplex and codec callers reach the single-strand consensus through + // `consensus_call`, so without this cap `--max-reads-per-strand` is silently ignored + // on those paths (DUPLEX3-01). The vanilla path caps raw records earlier in + // `process_group` and does not route through here, so it is unaffected. + // + // The cap shapes ONLY the consensus bases/quals/depths: the FULL (uncapped) + // `source_reads` are retained on the output below. fgbio caps inside `consensusCall` + // but passes the *pre-cap* `filteredAbR1s ++ filteredBaR2s` to `duplexConsensus` + // (`DuplexConsensusCaller.scala:331,337`, verified at tag 4.1.0), so the duplex caller — + // which counts per-base errors and calls the consensus UMI over `output.source_reads` — + // must see the uncapped set. Capping it here would systematically undercount duplex + // errors under `--max-reads-per-strand`, biasing the `ce`/error-rate tags. + let capped_for_consensus = match self.options.max_reads { + Some(max_reads) if source_reads.len() > max_reads => { + Some(self.downsample_source_reads(&source_reads)) + } + _ => None, + }; + let consensus_reads: &[SourceRead] = + capped_for_consensus.as_deref().unwrap_or(source_reads.as_slice()); + + // A degenerate cap (`max_reads = 0`, or a cap below `min_reads`) can leave fewer reads + // than are required; return no consensus rather than reaching `lengths[min_reads - 1]` + // or the empty-source `bail!` in `create_consensus_from_source_reads`. This is a no-op + // for every sane configuration (`max_reads >= min_reads`, or `max_reads` unset). + if consensus_reads.len() < self.options.min_reads { + return Ok(None); + } + + // Build consensus from the (capped, possibly normalized) scoring reads. let (bases, quals, depths, errors) = - self.create_consensus_from_source_reads(&source_reads)?; + self.create_consensus_from_source_reads(consensus_reads)?; // Truncate methylation annotation to consensus length (the anchor read may be // longer than the consensus when min_reads > 1 trims to shorter coverage). @@ -767,6 +799,44 @@ impl VanillaUmiConsensusCaller { reads } + /// Downsamples already-filtered `SourceReads` if there are more than `max_reads`. + /// + /// This is the `SourceRead` analogue of [`Self::downsample_reads`], applied inside + /// [`Self::consensus_call`] so the per-strand cap reaches the duplex and codec callers + /// (which build consensus from pre-filtered `SourceReads` rather than raw records). It + /// mirrors fgbio's `consensusCall` (`shuffle.take(maxReads)`): the retained subset is + /// selected by shuffling with the caller's seeded RNG (default seed 42) so the selection is + /// deterministic and unbiased, then keeping `max_reads`. When `max_reads` is unset or the + /// strand is at or below the cap this simply clones the input, preserving byte-identical + /// output for the default configuration. + /// + /// To avoid deep-cloning the entire (potentially very large) family just to discard most of + /// it, this shuffles an *index* vector `[0, len)` and clones only the `max_reads` retained + /// reads. This is bit-identical to shuffling the reads themselves and truncating: + /// `SliceRandom::shuffle` is Fisher-Yates, and its RNG-call sequence depends only on the + /// slice LENGTH, not the element type. So shuffling `[0..len]` yields the SAME permutation as + /// shuffling the reads, hence `idxs[0..max_reads]` selects exactly the reads (in the same + /// order) that `reads.shuffle().truncate(max_reads)` would keep. Selection therefore stays + /// deterministic-per-seed and unchanged from the prior implementation. + /// + /// The *rule* matches fgbio, but the *selection* does not byte-match it when downsampling + /// actually triggers: fgumi shuffles with `StdRng` (`ChaCha`) while fgbio uses Scala's + /// `Random.shuffle` over `java.util.Random` (an LCG), so the same seed value yields a + /// different permutation and therefore a different surviving subset. Downsampled consensus + /// bases/quals/depths are thus deterministic within fgumi but not byte-identical to fgbio on + /// any strand that exceeds the cap. + fn downsample_source_reads(&mut self, source_reads: &[SourceRead]) -> Vec { + match self.options.max_reads { + Some(max_reads) if source_reads.len() > max_reads => { + let mut idxs: Vec = (0..source_reads.len()).collect(); + idxs.shuffle(&mut self.rng); + idxs.truncate(max_reads); + idxs.into_iter().map(|i| source_reads[i].clone()).collect() + } + _ => source_reads.to_vec(), + } + } + /// Implements phred-style quality trimming from the 3' end. /// /// This matches htsjdk's `TrimmingUtil.findQualityTrimPoint` algorithm: @@ -1556,6 +1626,8 @@ mod tests { use fgumi_raw_bam::{ ParsedBamRecord, SamBuilder, encode_op, num_bases_extending_past_mate_raw, }; + use proptest::prelude::*; + use rstest::rstest; /// Call `consensus_reads` with a batch of already-built [`RawRecord`]s. fn consensus_reads_from_raw( @@ -1863,6 +1935,219 @@ mod tests { assert_eq!(downsampled.len(), 3); } + #[test] + fn test_consensus_call_caps_reads_per_strand() { + // DUPLEX3-01: `consensus_call` is the single-strand consensus entry point used by the + // duplex and codec callers (the vanilla path downsamples raw records earlier and never + // reaches this method). fgbio caps the contributing reads inside `consensusCall` + // (`shuffle.take(maxReads)`), so `--max-reads-per-strand` must limit how many reads + // contribute here. Without the cap, every read contributes and the per-strand depth + // exceeds the requested maximum. + let options = VanillaUmiConsensusOptions { + min_reads: 1, + max_reads: Some(3), + seed: Some(42), + ..Default::default() + }; + let mut caller = + VanillaUmiConsensusCaller::new("consensus".to_string(), "A".to_string(), options); + + // Eight identical, full-length source reads: the deepest position sees all eight + // unless the per-strand cap is applied. + let source_reads: Vec = + (0..8).map(|_| create_source_read_with_cigar("4M")).collect(); + + let consensus = caller + .consensus_call("umi", source_reads) + .expect("consensus_call should succeed") + .expect("consensus should be produced"); + + let max_depth = consensus.depths.iter().copied().max().unwrap_or(0); + assert_eq!( + max_depth, 3, + "per-strand depth must be capped at max_reads (3), got {max_depth}" + ); + } + + /// The `--max-reads-per-strand` cap must shape ONLY the single-strand consensus + /// (bases/quals/depths), NOT the `source_reads` the output carries. The duplex caller + /// counts per-base errors and calls the consensus UMI over `output.source_reads`, and fgbio + /// uses the *pre-cap* `filteredAbR1s ++ filteredBaR2s` for `duplexConsensus` + /// (`DuplexConsensusCaller.scala:337`, verified at tag 4.1.0). Capping the retained source + /// reads too would shrink the duplex error denominator and systematically undercount duplex + /// errors, biasing the `ce`/error-rate tags a base/read is filtered on. So: the consensus + /// depth is capped, but the full uncapped read set survives in `source_reads`. + #[test] + fn consensus_call_caps_consensus_but_retains_full_source_reads() { + let options = VanillaUmiConsensusOptions { + min_reads: 1, + max_reads: Some(3), + seed: Some(42), + ..Default::default() + }; + let mut caller = + VanillaUmiConsensusCaller::new("consensus".to_string(), "A".to_string(), options); + + let source_reads: Vec = + (0..8).map(|_| create_source_read_with_cigar("4M")).collect(); + + let consensus = caller + .consensus_call("umi", source_reads) + .expect("consensus_call should succeed") + .expect("consensus should be produced"); + + // The consensus itself is capped at max_reads (DUPLEX3-01) ... + assert_eq!( + consensus.depths.iter().copied().max().unwrap_or(0), + 3, + "per-strand consensus depth must be capped at max_reads (3)" + ); + // ... but the retained source reads are the FULL uncapped set, so the duplex caller's + // downstream error counting / UMI calling sees every read, matching fgbio. + assert_eq!( + consensus.source_reads.as_ref().map(Vec::len), + Some(8), + "output.source_reads must retain all 8 uncapped reads for duplex error counting" + ); + } + + /// The per-strand downsampling is deterministic for a fixed seed: two `consensus_call`s over + /// the same reads from a caller seeded identically select the same subset and so produce a + /// byte-identical consensus. Reads carry a distinct base composition so *which* reads survive + /// the cap changes the consensus — the selection (not just the count) is what is pinned. + /// (Determinism only: the shuffle uses fgumi's `StdRng`, not fgbio's `java.util.Random`, so + /// the selection is not byte-parity with fgbio even at the same seed — see + /// `downsample_source_reads`.) + #[test] + fn consensus_call_downsampling_is_deterministic_for_a_fixed_seed() { + let make_reads = || -> Vec { + (0..8) + .map(|i| { + let mut sr = create_source_read_with_cigar("4M"); + // A third of the reads read 'C' instead of 'A', so the consensus base at each + // position depends on which reads the cap keeps. + if i % 3 == 0 { + sr.bases = vec![b'C'; sr.bases.len()]; + } + sr + }) + .collect() + }; + let call = || { + let options = VanillaUmiConsensusOptions { + min_reads: 1, + max_reads: Some(3), + seed: Some(42), + ..Default::default() + }; + let mut caller = + VanillaUmiConsensusCaller::new("consensus".to_string(), "A".to_string(), options); + caller + .consensus_call("umi", make_reads()) + .expect("consensus_call should succeed") + .expect("consensus should be produced") + }; + let first = call(); + let second = call(); + assert_eq!(first.bases, second.bases, "same seed must select the same reads → same bases"); + assert_eq!(first.quals, second.quals, "same seed → identical quals"); + assert_eq!(first.depths, second.depths, "same seed → identical depths"); + } + + #[rstest] + // Cap above the read count: no downsampling, every read contributes. + #[case::cap_above_count(Some(12), 8)] + // Cap equal to the read count: boundary, no downsampling, full depth. + #[case::cap_equals_count(Some(8), 8)] + // Cap below the read count: downsampled to the cap. + #[case::cap_below_count(Some(3), 3)] + // No cap: full depth (default configuration, byte-identical behavior). + #[case::no_cap(None, 8)] + // Degenerate zero cap: fewer than `min_reads` survive, so no consensus is produced + // (and, critically, no panic / empty-source bail leaks out). + #[case::zero_cap(Some(0), 0)] + fn test_consensus_call_max_reads_boundaries( + #[case] max_reads: Option, + #[case] expected_max_depth: u16, + ) { + let options = VanillaUmiConsensusOptions { + min_reads: 1, + max_reads, + seed: Some(42), + ..Default::default() + }; + let mut caller = + VanillaUmiConsensusCaller::new("consensus".to_string(), "A".to_string(), options); + + let source_reads: Vec = + (0..8).map(|_| create_source_read_with_cigar("4M")).collect(); + + let consensus = caller + .consensus_call("umi", source_reads) + .expect("consensus_call should never error on the cap boundaries"); + + let observed_max_depth = + consensus.map_or(0, |c| c.depths.iter().copied().max().unwrap_or(0)); + assert_eq!( + observed_max_depth, expected_max_depth, + "max_reads={max_reads:?} should yield max depth {expected_max_depth}" + ); + } + + proptest! { + #![proptest_config(ProptestConfig::with_cases(48))] + + /// DUPLEX3-01 cap invariants over arbitrary family sizes, caps, and `min_reads`. + /// Whenever a consensus is produced: (a) at most `max_reads` reads contribute to it — + /// observed as the maximum per-base depth, which equals the contributing-read count here + /// because every read is identical and full length — and (b) the FULL uncapped family is + /// always retained on `output.source_reads`, so the duplex caller's downstream per-base + /// error counting and UMI calling see every read (the invariant the PR establishes). + #[test] + fn consensus_call_cap_invariants( + n in 1usize..=12, + max_reads in 1usize..=12, + min_reads in 1usize..=4, + ) { + let options = VanillaUmiConsensusOptions { + min_reads, + max_reads: Some(max_reads), + seed: Some(42), + ..Default::default() + }; + let mut caller = VanillaUmiConsensusCaller::new( + "consensus".to_string(), + "A".to_string(), + options, + ); + + let source_reads: Vec = + (0..n).map(|_| create_source_read_with_cigar("4M")).collect(); + + let consensus = caller + .consensus_call("umi", source_reads) + .expect("consensus_call must not error on any cap/size/min_reads combination"); + + if let Some(consensus) = consensus { + // (a) The consensus is capped: the deepest position sees at most `max_reads` + // reads (the family exceeds the cap exactly when `n > max_reads`). + let observed_max_depth = consensus.depths.iter().copied().max().unwrap_or(0); + prop_assert!( + usize::from(observed_max_depth) <= max_reads, + "consensus depth {} exceeded cap {}", + observed_max_depth, + max_reads + ); + // (b) The full uncapped family survives on the output for duplex error counting. + prop_assert_eq!( + consensus.source_reads.as_ref().map(Vec::len), + Some(n), + "output.source_reads must retain all uncapped reads" + ); + } + } + } + /// Parse a CIGAR string into BAM-encoded `u32` ops for use with [`SamBuilder::cigar_ops`]. /// /// BAM op codes: 0=M, 1=I, 2=D, 3=N, 4=S, 5=H, 6=P, 7==, 8=X