From aff6d77ccda1c379364349ac0955f78ed0d26452 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Mon, 20 Jul 2026 13:22:19 -0700 Subject: [PATCH] fix(simulate): emit faithful duplex consensus tags MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two defects made `fgumi simulate consensus-reads --duplex` produce records that no real consensus caller would emit, so simulated duplex BAMs were unsound fixtures for validating `fgumi filter` and for fgbio-parity harnesses. Combined per-base arrays on duplex records. `CD_BASES`/`CE_BASES` were written unconditionally, but `duplex_caller` emits only the per-strand `AD_BASES`/ `AE_BASES`/`BD_BASES`/`BE_BASES` arrays alongside the scalar cD/cM/cE. Emit the combined arrays only in simplex mode. This is safe for the documented "suitable for input to `fgumi filter`" contract: filter branches on `is_duplex_consensus` and masks duplex reads through `mask_duplex_bases`, which reads the four per-strand arrays and never the combined ones, while simplex reads still get the `CD_BASES`/`CE_BASES` that `mask_bases` requires. Derive the duplex cD/cM/cE scalars from the per-strand sums rather than from the independently sampled simplex values, matching `duplex_caller`, which computes them over the combined per-base depth (ab_i + ba_i). Strand minimum depths. `aM` and `bM` were each computed as `strand_depth.min(cM)`, so both could equal `cM` and `aM + bM` routinely exceeded it. A duplex position's combined depth is the sum of its two strand depths, so the truth TSV and the BAM tags disagreed with each other and anyone validating a real run against `--truth` got impossible expectations. Split `cM` by the same strand fraction already used for `cD`, clamping the A share to `[cM - bD, min(aD, cM)]` — feasible because `cM <= cD == aD + bD` — which guarantees `aM + bM == cM`, `aM <= aD`, and `bM <= bD`. The combined scalars are summed from the per-strand scalars rather than reduced over the per-base arrays: the two agree only when the arrays carry the min anchor, which `per_base_arrays` adds solely for `read_len >= 2`, so a `--read-length 1` run would have reported `cM == cD` and drifted from the truth TSV again. A proptest over generated seeds (the strand fraction is sampled per read) asserts the summation and ordering invariants, an `rstest` case table pins the truth/tag agreement at read lengths 1, 2, and 50, and further tests assert duplex records carry no `CD_BASES`/`CE_BASES` while retaining all four per-strand arrays, and that simplex records still carry the combined arrays. --- src/lib/commands/simulate/consensus_reads.rs | 235 ++++++++++++++++++- 1 file changed, 223 insertions(+), 12 deletions(-) diff --git a/src/lib/commands/simulate/consensus_reads.rs b/src/lib/commands/simulate/consensus_reads.rs index 15a15b4ee..977df872e 100644 --- a/src/lib/commands/simulate/consensus_reads.rs +++ b/src/lib/commands/simulate/consensus_reads.rs @@ -378,9 +378,16 @@ fn generate_consensus_pair( let ad = ((cd as f64) * a_frac).round() as i32; let bd = cd - ad; - // Min depths for each strand - let am = (ad.min(cm)).max(0); - let bm = (bd.min(cm)).max(0); + // Min depths for each strand. These must satisfy aM + bM == cM, because a + // duplex consensus position's combined depth is the sum of its two strand + // depths — the real duplex caller derives cM as the per-base minimum of + // (ab_i + ba_i). Splitting cM by the same strand fraction as cD keeps the + // truth file and the emitted tags mutually consistent. + // + // The split is feasible because cM <= cD == aD + bD; clamping the A share to + // [cM - bD, min(aD, cM)] guarantees both aM <= aD and bM = cM - aM <= bD. + let am = (((cm as f64) * a_frac).round() as i32).clamp((cm - bd).max(0), ad.min(cm)); + let bm = cm - am; // Errors distributed proportionally let ae = ((ce as f64) * a_frac).round() as i32; @@ -729,18 +736,40 @@ fn build_consensus_record( if total_depth > 0 { (total_errors as f64 / total_depth as f64) as f32 } else { 0.0 } }; - // Single-strand consensus summary + per-base arrays. - let (cd_bases, ce_bases) = per_base_arrays(cd, cm, ce); - b.add_int_tag(SamTag::CD, cd) - .add_int_tag(SamTag::CM, cm) - .add_float_tag(SamTag::CE, error_rate(&ce_bases, &cd_bases)); - b.add_array_i16(SamTag::CD_BASES, &cd_bases).add_array_i16(SamTag::CE_BASES, &ce_bases); - - // Add duplex-specific tags if present: per-strand summary (aD/bD/aM/bM int, - // aE/bE float rate) plus the per-base strand depth/error arrays filter needs. if let Some((ad, bd, am, bm, ae, be)) = duplex_tags { + // Duplex: emit per-strand summary (aD/bD/aM/bM int, aE/bE float rate) plus + // the per-base strand arrays `fgumi filter` masks with (`mask_duplex_bases` + // reads AD/AE/BD/BE only). + // + // The combined cD/cM/cE scalars are derived from the strand sums rather than + // the independently sampled values, matching `duplex_caller`, which computes + // them over the per-base combined depth (ab_i + ba_i). + // + // No CD_BASES/CE_BASES here: the real duplex caller never emits the combined + // arrays, so writing them would make simulated duplex BAMs diverge from + // anything production produces. let (ad_bases, ae_bases) = per_base_arrays(ad, am, ae); let (bd_bases, be_bases) = per_base_arrays(bd, bm, be); + + let combined_depths: Vec = ad_bases + .iter() + .zip(&bd_bases) + .map(|(&a, &b)| i16::try_from(i32::from(a) + i32::from(b)).unwrap_or(i16::MAX)) + .collect(); + let combined_errors: Vec = ae_bases + .iter() + .zip(&be_bases) + .map(|(&a, &b)| i16::try_from(i32::from(a) + i32::from(b)).unwrap_or(i16::MAX)) + .collect(); + + // Sum the scalars rather than taking max/min over `combined_depths`: the two + // agree whenever the arrays carry the min anchor, but `per_base_arrays` only + // anchors it when `read_len >= 2`, so a 1-base read would report cM == cD (and + // a 0-base read cD == cM == 0) and drift from the truth TSV again. + b.add_int_tag(SamTag::CD, ad + bd) + .add_int_tag(SamTag::CM, am + bm) + .add_float_tag(SamTag::CE, error_rate(&combined_errors, &combined_depths)); + b.add_int_tag(SamTag::AD, ad) .add_int_tag(SamTag::BD, bd) .add_int_tag(SamTag::AM, am) @@ -751,6 +780,13 @@ fn build_consensus_record( .add_array_i16(SamTag::AE_BASES, &ae_bases) .add_array_i16(SamTag::BD_BASES, &bd_bases) .add_array_i16(SamTag::BE_BASES, &be_bases); + } else { + // Simplex: combined summary + the per-base arrays `mask_bases` reads. + let (cd_bases, ce_bases) = per_base_arrays(cd, cm, ce); + b.add_int_tag(SamTag::CD, cd) + .add_int_tag(SamTag::CM, cm) + .add_float_tag(SamTag::CE, error_rate(&ce_bases, &cd_bases)); + b.add_array_i16(SamTag::CD_BASES, &cd_bases).add_array_i16(SamTag::CE_BASES, &ce_bases); } // Add methylation tags if enabled. @@ -806,6 +842,7 @@ mod tests { use super::*; use crate::commands::simulate::common::generate_random_sequence; use crate::simulate::create_rng; + use rstest::rstest; /// Decode a raw BAM record into a noodles `RecordBuf` for higher-level test /// assertions. @@ -1929,6 +1966,180 @@ mod tests { ); } + /// Build `GenerationParams` over a small synthetic reference, duplex-configurable. + fn strand_test_params(fasta: &tempfile::NamedTempFile, duplex: bool) -> Arc { + strand_test_params_with_read_length(fasta, duplex, 50) + } + + /// As [`strand_test_params`], with an explicit read length so tests can reach the + /// short-read boundaries (`--read-length 1` skips the per-base min anchor). + fn strand_test_params_with_read_length( + fasta: &tempfile::NamedTempFile, + duplex: bool, + read_length: usize, + ) -> Arc { + let ref_genome = Arc::new(ReferenceGenome::load(fasta.path()).unwrap()); + Arc::new(GenerationParams { + read_length, + min_depth: 1, + max_depth: 30, + depth_mean: 12.0, + depth_stddev: 6.0, + error_rate_mean: 0.02, + error_rate_stddev: 0.01, + duplex, + consensus_quality: 40, + methylation_mode: MethylationMode::Disabled, + methylation_depth_dist: create_depth_distribution(5.0, 2.5), + cpg_methylation_rate: 0.75, + conversion_rate: 0.98, + ref_genome, + }) + } + + fn small_test_fasta() -> tempfile::NamedTempFile { + use std::io::Write as IoWrite; + let mut fasta = tempfile::NamedTempFile::new().unwrap(); + writeln!(fasta, ">chr1").unwrap(); + fasta.write_all(&b"ACGT".repeat(500)).unwrap(); + writeln!(fasta).unwrap(); + fasta.flush().unwrap(); + fasta + } + + /// Duplex strand minimum depths must sum to the combined minimum, and each + /// strand's minimum must not exceed its own maximum. + /// + /// Previously `aM` and `bM` were each computed as `strand_depth.min(cM)`, so both + /// could equal `cM` and `aM + bM` routinely exceeded it — making the truth file + /// and the emitted tags mutually inconsistent. The strand fraction is sampled per + /// read, so the invariants are checked as a property over generated seeds rather + /// than at a handful of fixed points; the reference and params are built once and + /// shared across cases because loading them is the expensive part. + #[test] + fn test_duplex_strand_minimums_sum_to_combined_minimum() { + use proptest::prelude::*; + + let fasta = small_test_fasta(); + let params = strand_test_params(&fasta, true); + let strand_bias = StrandBiasModel::new(5.0, 5.0); + + proptest!(|(seed in any::())| { + let pair = generate_consensus_pair(0, seed, ¶ms, &strand_bias); + let aux = fgumi_raw_bam::aux_data_slice(&pair.r1_record); + + let get = |tag: SamTag| -> i64 { + fgumi_raw_bam::find_int_tag(aux, tag) + .unwrap_or_else(|| panic!("tag {tag:?} missing for seed {seed}")) + }; + let (cd, cm, ad, bd, am, bm) = ( + get(SamTag::CD), + get(SamTag::CM), + get(SamTag::AD), + get(SamTag::BD), + get(SamTag::AM), + get(SamTag::BM), + ); + + // The truth TSV and the emitted tags must not drift apart: the truth + // tuple carries the sampled values while duplex CD/CM are derived from + // the strand sums, and only the aM+bM==cM invariant keeps them equal. + let (truth_cd, truth_cm) = (pair.truth.0, pair.truth.1); + prop_assert_eq!( + i64::from(truth_cd), + cd, + "truth cD ({}) disagrees with emitted CD ({}) for seed {}", + truth_cd, + cd, + seed + ); + prop_assert_eq!( + i64::from(truth_cm), + cm, + "truth cM ({}) disagrees with emitted CM ({}) for seed {}", + truth_cm, + cm, + seed + ); + + prop_assert_eq!(am + bm, cm, "aM + bM != cM (aM={} bM={} cM={})", am, bm, cm); + prop_assert_eq!(ad + bd, cd, "aD + bD != cD (aD={} bD={} cD={})", ad, bd, cd); + prop_assert!(am <= ad, "aM ({}) > aD ({})", am, ad); + prop_assert!(bm <= bd, "bM ({}) > bD ({})", bm, bd); + prop_assert!(am >= 0 && bm >= 0, "negative strand minimum (aM={} bM={})", am, bm); + }); + } + + /// The duplex cD/cM tags must match the truth TSV at short read lengths too. + /// + /// `per_base_arrays` only anchors the per-base minimum when `read_len >= 2`, so + /// deriving the scalars as max/min over the combined per-base depths would report + /// `cM == cD` for a 1-base read — the same truth-vs-tag drift this change fixes. + #[rstest] + #[case::single_base(1)] + #[case::two_bases(2)] + #[case::typical(50)] + fn test_duplex_combined_scalars_match_truth_at_short_read_lengths(#[case] read_length: usize) { + let fasta = small_test_fasta(); + let params = strand_test_params_with_read_length(&fasta, true, read_length); + let strand_bias = StrandBiasModel::new(5.0, 5.0); + + for seed in 0..8u64 { + let pair = generate_consensus_pair(0, seed, ¶ms, &strand_bias); + let aux = fgumi_raw_bam::aux_data_slice(&pair.r1_record); + let cd = fgumi_raw_bam::find_int_tag(aux, SamTag::CD).expect("CD missing"); + let cm = fgumi_raw_bam::find_int_tag(aux, SamTag::CM).expect("CM missing"); + let (truth_cd, truth_cm) = (pair.truth.0, pair.truth.1); + assert_eq!( + i64::from(truth_cd), + cd, + "truth cD ({truth_cd}) != emitted CD ({cd}) at read_length {read_length}, seed {seed}" + ); + assert_eq!( + i64::from(truth_cm), + cm, + "truth cM ({truth_cm}) != emitted CM ({cm}) at read_length {read_length}, seed {seed}" + ); + } + } + + /// Duplex records must not carry the combined per-base arrays: the real + /// `duplex_caller` emits only the per-strand AD/AE/BD/BE arrays, so emitting + /// `CD_BASES`/`CE_BASES` would make simulated duplex BAMs unfaithful fixtures. + /// Simplex records must still carry them (`fgumi filter`'s `mask_bases` reads them). + #[test] + fn test_duplex_omits_combined_per_base_arrays_simplex_keeps_them() { + let fasta = small_test_fasta(); + let strand_bias = StrandBiasModel::new(5.0, 5.0); + + let duplex_params = strand_test_params(&fasta, true); + let duplex_pair = generate_consensus_pair(0, 7, &duplex_params, &strand_bias); + let duplex_aux = fgumi_raw_bam::aux_data_slice(&duplex_pair.r1_record); + assert!( + fgumi_raw_bam::find_array_tag(duplex_aux, SamTag::CD_BASES).is_none(), + "duplex record must not carry CD_BASES (no real caller emits it)" + ); + assert!( + fgumi_raw_bam::find_array_tag(duplex_aux, SamTag::CE_BASES).is_none(), + "duplex record must not carry CE_BASES (no real caller emits it)" + ); + // The per-strand arrays filter masks with must still be present. + for tag in [SamTag::AD_BASES, SamTag::AE_BASES, SamTag::BD_BASES, SamTag::BE_BASES] { + assert!( + fgumi_raw_bam::find_array_tag(duplex_aux, tag).is_some(), + "duplex record missing {tag:?}" + ); + } + + let simplex_params = strand_test_params(&fasta, false); + let simplex_pair = generate_consensus_pair(0, 7, &simplex_params, &strand_bias); + let simplex_aux = fgumi_raw_bam::aux_data_slice(&simplex_pair.r1_record); + assert!( + fgumi_raw_bam::find_array_tag(simplex_aux, SamTag::CD_BASES).is_some(), + "simplex record must keep CD_BASES for filter's mask_bases" + ); + } + #[test] fn test_consensus_reads_produces_mapped_records() { use std::io::Write as IoWrite;