From b0fd52c7c5d7bc2158ce22c48f943b8b5d98c204 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Wed, 8 Jul 2026 17:53:51 -0700 Subject: [PATCH 1/4] fix(metrics): align TSV metric output with fgbio Metric format (#498) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Make fgumi's metrics `.txt` outputs readable by fgbio's `Metric.read` and match fgbio's rows, column names/order, and value tokens on a shared input. Verified end-to-end against fgbio 4.1.0's real `Metric.read` (see fixture below). - Non-finite floats (MET-02): add a shared serde helper (`float` module) that writes f64 fields as fgbio's `Infinity`/`-Infinity`/`NaN` tokens (integral values without a trailing `.0`), and apply it to every serialized f64 field across correct/shared/consensus/simplex/duplex/group. Previously only `correct.rs` was inf-safe; the rest emitted ryu `inf`/`-inf`, which fgbio's `Double.parseDouble` cannot read. - Empty-metrics header (MET-03): `write_metrics` now always emits the column header, even for a zero-row slice (the csv writer's lazy headers otherwise produce a 0-byte file that fgbio rejects with "No header found"). Route the duplex-metrics and simplex-metrics per-family/UMI/yield writes through it. - UmiGroupingMetrics (MET-04): serialize exactly fgbio's 5-column `UmiGroupingMetric` (`accepted_sam_records`, `discarded_non_pf`, `discarded_poor_alignment`, `discarded_ns_in_umi`, `discarded_umis_to_short`). The fgumi-internal summary fields (`total_records`, family aggregates) are retained in memory for logging but `#[serde(skip)]`-ed; they're derivable and the family distribution already lives in `.family_sizes.txt`. - Duplex family sizes sparse: emit one row per observed family size (matching fgbio `familySizeMetrics`) instead of a dense `1..=max` range. - Duplex UMI counts determinism: add `umi` as a secondary sort key so ties no longer fall back to nondeterministic HashMap order. - Consensus KV rejection text (MET-05, partial): match fgbio's exact rejection descriptions ("Reads has a different, and minority, set of indels"; "Only Generating One Strand of Duplex Consensus"). The caller-type-aware row-set seeding (non_paired_reads / potential_umi_collision, etc.) is deferred to the consensus PR (fgbio's usedByVanilla/Duplex/Codec seeding + R2-UCC-01), as it requires new RejectionReason variants wired into the duplex/codec callers. - format_float doc (MET-01): correct the doc-comment that falsely claimed it matches fgbio; numeric equality (not byte format) is the metric contract. CorrectUmis zero-count seeding (#498 item 5) already reproduces as fixed on main: both correct execution paths pre-seed a row per expected UMI plus the all-N sentinel. No change needed; kept as-is. fgbio-oracle evidence (reports/fgbio-parity-fixtures/met-02-03-04-fgbio-read.sh): fgbio 4.1.0 `Metric.read[UmiGroupingMetric]` parses fgumi's 5-column file and its header-only empty file (0 rows, no error); `Metric.read[UmiCorrectionMetrics]` parses `representation=Infinity` — while the old ryu `inf` token is rejected. CI: cargo ci-fmt / ci-lint clean; full suite 2215 passed. --- crates/fgumi-metrics/src/consensus.rs | 2 + crates/fgumi-metrics/src/correct.rs | 58 +----------- crates/fgumi-metrics/src/duplex.rs | 94 +++++++++++++++---- crates/fgumi-metrics/src/float.rs | 120 ++++++++++++++++++++++++ crates/fgumi-metrics/src/group.rs | 93 +++++++++++++++--- crates/fgumi-metrics/src/lib.rs | 16 +++- crates/fgumi-metrics/src/rejection.rs | 13 ++- crates/fgumi-metrics/src/shared.rs | 2 + crates/fgumi-metrics/src/simplex.rs | 7 ++ crates/fgumi-metrics/src/writer.rs | 49 ++++++++-- src/lib/commands/clip.rs | 2 +- src/lib/commands/duplex_metrics.rs | 31 +++--- src/lib/commands/group.rs | 32 +++++-- src/lib/commands/simplex_metrics.rs | 18 ++-- src/lib/metrics/mod.rs | 2 +- tests/integration/test_group_command.rs | 7 +- 16 files changed, 400 insertions(+), 146 deletions(-) create mode 100644 crates/fgumi-metrics/src/float.rs diff --git a/crates/fgumi-metrics/src/consensus.rs b/crates/fgumi-metrics/src/consensus.rs index eef802606..bf787bbd0 100644 --- a/crates/fgumi-metrics/src/consensus.rs +++ b/crates/fgumi-metrics/src/consensus.rs @@ -57,9 +57,11 @@ pub struct ConsensusMetrics { pub umi_groups_failed: u64, /// Average input reads per consensus read + #[serde(with = "crate::float")] pub avg_input_reads_per_consensus: f64, /// Average raw read depth per consensus read + #[serde(with = "crate::float")] pub avg_raw_read_depth: f64, /// Minimum raw read depth diff --git a/crates/fgumi-metrics/src/correct.rs b/crates/fgumi-metrics/src/correct.rs index 5e1353658..ba6a43cc9 100644 --- a/crates/fgumi-metrics/src/correct.rs +++ b/crates/fgumi-metrics/src/correct.rs @@ -3,57 +3,11 @@ //! This module provides metrics for tracking UMI correction operations. use anyhow::Result; -use serde::{Deserialize, Deserializer, Serialize, Serializer}; +use serde::{Deserialize, Serialize}; use std::path::Path; use crate::Metric; -/// Serializes an f64 as an integer string if it's a whole number. -/// -/// This matches fgbio's output format where 0.0 is output as "0". -#[expect(clippy::trivially_copy_pass_by_ref, reason = "serde requires &T signature")] -fn serialize_f64_as_int(value: &f64, serializer: S) -> Result -where - S: Serializer, -{ - let value = *value; - // Check for NaN and Infinity first - if value.is_nan() { - return serializer.serialize_str("NaN"); - } - if value.is_infinite() { - return serializer.serialize_str(if value > 0.0 { "Infinity" } else { "-Infinity" }); - } - - // If the value is a whole number, serialize without decimal point - #[expect( - clippy::cast_precision_loss, - reason = "i64::MAX boundary check, exact precision not needed" - )] - let max_safe = i64::MAX as f64; - if value.fract() == 0.0 && value.abs() < max_safe { - #[expect(clippy::cast_possible_truncation, reason = "guarded by abs check above")] - let int_value = value as i64; - serializer.serialize_str(&format!("{int_value}")) - } else { - serializer.serialize_str(&format!("{value}")) - } -} - -/// Deserializes an f64 from a string, handling special values like NaN and Infinity. -fn deserialize_f64_from_str<'de, D>(deserializer: D) -> Result -where - D: Deserializer<'de>, -{ - let s: String = String::deserialize(deserializer)?; - match s.as_str() { - "NaN" => Ok(f64::NAN), - "Infinity" => Ok(f64::INFINITY), - "-Infinity" => Ok(f64::NEG_INFINITY), - _ => s.parse().map_err(serde::de::Error::custom), - } -} - /// Metrics tracking how well observed UMIs match expected UMI sequences. /// /// These metrics are generated per-UMI and track the distribution of match types @@ -90,17 +44,11 @@ pub struct UmiCorrectionMetrics { pub other_matches: u64, /// The fraction of all UMIs that matched or were corrected to this UMI. - #[serde( - serialize_with = "serialize_f64_as_int", - deserialize_with = "deserialize_f64_from_str" - )] + #[serde(with = "crate::float")] pub fraction_of_matches: f64, /// The `total_matches` for this UMI divided by the mean `total_matches` for all UMIs. - #[serde( - serialize_with = "serialize_f64_as_int", - deserialize_with = "deserialize_f64_from_str" - )] + #[serde(with = "crate::float")] pub representation: f64, } diff --git a/crates/fgumi-metrics/src/duplex.rs b/crates/fgumi-metrics/src/duplex.rs index e32025a00..619e52757 100644 --- a/crates/fgumi-metrics/src/duplex.rs +++ b/crates/fgumi-metrics/src/duplex.rs @@ -23,20 +23,26 @@ pub struct FamilySizeMetric { /// Count of CS families with this size pub cs_count: usize, /// Fraction of all CS families with this size + #[serde(with = "crate::float")] pub cs_fraction: f64, /// Fraction of CS families with size >= `family_size` + #[serde(with = "crate::float")] pub cs_fraction_gt_or_eq_size: f64, /// Count of SS families with this size pub ss_count: usize, /// Fraction of all SS families with this size + #[serde(with = "crate::float")] pub ss_fraction: f64, /// Fraction of SS families with size >= `family_size` + #[serde(with = "crate::float")] pub ss_fraction_gt_or_eq_size: f64, /// Count of DS families with this size pub ds_count: usize, /// Fraction of all DS families with this size + #[serde(with = "crate::float")] pub ds_fraction: f64, /// Fraction of DS families with size >= `family_size` + #[serde(with = "crate::float")] pub ds_fraction_gt_or_eq_size: f64, } @@ -83,8 +89,10 @@ pub struct DuplexFamilySizeMetric { /// Count of families with these AB/BA sizes pub count: usize, /// Fraction of all duplex families with these sizes + #[serde(with = "crate::float")] pub fraction: f64, /// Fraction of duplex families with AB >= `ab_size` and BA >= `ba_size` + #[serde(with = "crate::float")] pub fraction_gt_or_eq_size: f64, } @@ -132,6 +140,7 @@ impl PartialEq for DuplexFamilySizeMetric { #[derive(Debug, Clone, Serialize, Deserialize)] pub struct DuplexYieldMetric { /// Approximate fraction of full dataset used + #[serde(with = "crate::float")] pub fraction: f64, /// Number of read pairs upon which metrics are based pub read_pairs: usize, @@ -144,8 +153,10 @@ pub struct DuplexYieldMetric { /// Number of DS families that are duplexes (min reads on both strands) pub ds_duplexes: usize, /// Fraction of DS families that are duplexes + #[serde(with = "crate::float")] pub ds_fraction_duplexes: f64, /// Expected fraction of DS families that should be duplexes under ideal model + #[serde(with = "crate::float")] pub ds_fraction_duplexes_ideal: f64, } @@ -192,10 +203,13 @@ pub struct DuplexUmiMetric { /// Number of double-stranded tag families observing this duplex UMI pub unique_observations: usize, /// Fraction of all raw observations + #[serde(with = "crate::float")] pub fraction_raw_observations: f64, /// Fraction of all unique observations + #[serde(with = "crate::float")] pub fraction_unique_observations: f64, /// Expected fraction based on individual UMI frequencies + #[serde(with = "crate::float")] pub fraction_unique_observations_expected: f64, } @@ -313,30 +327,30 @@ impl DuplexMetricsCollector { self.duplex_umi_counts.record(umi, raw_count, error_count, is_unique); } - /// Generates family size metrics - /// - /// # Panics - /// - /// Panics if the internal array-of-three maximum computation fails, which cannot - /// happen since the array is always non-empty. + /// Generates family size metrics, one row per observed family size (sparse), sorted + /// ascending by `family_size`. #[must_use] pub fn family_size_metrics(&self) -> Vec { - // Find max family size across all types - let max_size = *[ - self.cs_family_sizes.keys().max().unwrap_or(&0), - self.ss_family_sizes.keys().max().unwrap_or(&0), - self.ds_family_sizes.keys().max().unwrap_or(&0), - ] - .iter() - .max() - .expect("array of three elements always has a maximum"); - let coord_strand_total: usize = self.cs_family_sizes.values().sum(); let single_strand_total: usize = self.ss_family_sizes.values().sum(); let double_strand_total: usize = self.ds_family_sizes.values().sum(); - let mut metrics = Vec::new(); - for size in 1..=*max_size { + // fgbio emits one row per *observed* family size (sparse), not a dense + // `1..=max` range: `CollectDuplexSeqMetrics.familySizeMetrics` builds a map + // keyed only by sizes that occur in the CS/SS/DS counters, then sorts + // ascending by `family_size`. Take the sorted union of observed sizes. + let mut sizes: Vec = self + .cs_family_sizes + .keys() + .chain(self.ss_family_sizes.keys()) + .chain(self.ds_family_sizes.keys()) + .copied() + .collect(); + sizes.sort_unstable(); + sizes.dedup(); + + let mut metrics = Vec::with_capacity(sizes.len()); + for size in sizes { let mut metric = FamilySizeMetric::new(size); metric.cs_count = *self.cs_family_sizes.get(&size).unwrap_or(&0); @@ -469,8 +483,13 @@ impl DuplexMetricsCollector { }) .collect(); - // Sort by unique observations descending (matching Scala's sort order) - metrics.sort_by(|a, b| b.unique_observations.cmp(&a.unique_observations)); + // Sort by unique observations descending (matching fgbio's `sortBy(-unique_observations)`), + // with the UMI sequence as a deterministic secondary key. fgbio's own tie order comes from + // nondeterministic map iteration; keying ties on `umi` makes fgumi's output stable across + // runs and thread counts (an acceptance-criterion requirement) without changing the row set. + metrics.sort_by(|a, b| { + b.unique_observations.cmp(&a.unique_observations).then_with(|| a.umi.cmp(&b.umi)) + }); metrics } } @@ -775,6 +794,41 @@ mod tests { assert!((duplex_metrics[0].fraction_unique_observations_expected - 0.25).abs() < 0.001); } + #[test] + fn test_family_size_metrics_are_sparse() { + let mut collector = DuplexMetricsCollector::new(false); + // Observe only sizes 1 and 5 (a gap at 2, 3, 4). fgbio's `familySizeMetrics` + // emits a row per *observed* size; the old dense `1..=max` range emitted + // phantom all-zero rows for the unobserved sizes. + collector.record_cs_family(1); + collector.record_cs_family(5); + + let metrics = collector.family_size_metrics(); + let sizes: Vec = metrics.iter().map(|m| m.family_size).collect(); + assert_eq!(sizes, vec![1, 5], "only observed family sizes should be emitted"); + + // The reverse-cumulative `>= size` fraction stays correct across the gap: + // each size holds half the CS families, so size 1 → 1.0, size 5 → 0.5. + assert!((metrics[0].cs_fraction_gt_or_eq_size - 1.0).abs() < 1e-9); + assert!((metrics[1].cs_fraction_gt_or_eq_size - 0.5).abs() < 1e-9); + } + + #[test] + fn test_duplex_umi_metrics_deterministic_tie_order() { + // Three duplex UMIs each observed exactly once → equal `unique_observations` + // (a tie). Insertion order is reverse-sorted; the output must be ordered by + // `umi` ascending (the deterministic secondary sort key) rather than by the + // nondeterministic HashMap iteration order the primary key alone would leave. + let mut collector = DuplexMetricsCollector::new(true); + for umi in ["GGGG-CCCC", "AAAA-TTTT", "CCCC-GGGG"] { + collector.record_duplex_umi(umi, 1, 0, true); + } + + let metrics = collector.duplex_umi_metrics(&[]); + let order: Vec<&str> = metrics.iter().map(|m| m.umi.as_str()).collect(); + assert_eq!(order, vec!["AAAA-TTTT", "CCCC-GGGG", "GGGG-CCCC"]); + } + #[test] fn test_empty_collector() { let collector = DuplexMetricsCollector::new(false); diff --git a/crates/fgumi-metrics/src/float.rs b/crates/fgumi-metrics/src/float.rs new file mode 100644 index 000000000..96f3eca94 --- /dev/null +++ b/crates/fgumi-metrics/src/float.rs @@ -0,0 +1,120 @@ +//! Shared serde helpers for writing `f64` metric fields in fgbio's `Metric` TSV format. +//! +//! fgbio serializes doubles with `Double.toString`, which renders non-finite values +//! as `Infinity`, `-Infinity`, and `NaN` — the exact tokens `Double.parseDouble` +//! accepts on read. Rust's default `f64` serialization (via ryu) instead emits +//! `inf`/`-inf`/`NaN`, and `inf`/`-inf` are **unparseable** by fgbio's `Metric.read`. +//! +//! Every `f64` field in an fgumi metric that is meant to be read back by fgbio must +//! carry `#[serde(with = "crate::float")]` so non-finite values round-trip through +//! both tools. Finite whole numbers are written without a trailing `.0` (e.g. `0` +//! rather than `0.0`), matching fgbio's `DecimalFormat` output for integral values; +//! other finite values keep full precision (numeric equality — not byte-for-byte +//! formatting — is the metric contract, see the crate `format_float` docs). + +use serde::{Deserialize, Deserializer, Serializer}; + +/// Serializes an `f64` in fgbio's `Metric` TSV format. +/// +/// Non-finite values become `Infinity`/`-Infinity`/`NaN` (readable by fgbio's +/// `Double.parseDouble`); finite whole numbers drop the fractional part; all other +/// finite values keep full precision. +#[expect(clippy::trivially_copy_pass_by_ref, reason = "serde requires &T signature")] +pub(crate) fn serialize(value: &f64, serializer: S) -> Result +where + S: Serializer, +{ + let value = *value; + if value.is_nan() { + return serializer.serialize_str("NaN"); + } + if value.is_infinite() { + return serializer.serialize_str(if value > 0.0 { "Infinity" } else { "-Infinity" }); + } + + // Render an integral value without a decimal point (fgbio's DecimalFormat does the same). + #[expect( + clippy::cast_precision_loss, + reason = "i64::MAX boundary check, exact precision not needed" + )] + let max_safe = i64::MAX as f64; + if value.fract() == 0.0 && value.abs() < max_safe { + #[expect(clippy::cast_possible_truncation, reason = "guarded by abs check above")] + let int_value = value as i64; + serializer.serialize_str(&format!("{int_value}")) + } else { + serializer.serialize_str(&format!("{value}")) + } +} + +/// Deserializes an `f64` written by [`serialize`], accepting the non-finite +/// tokens `NaN`, `Infinity`, and `-Infinity` in addition to ordinary numbers. +pub(crate) fn deserialize<'de, D>(deserializer: D) -> Result +where + D: Deserializer<'de>, +{ + let s: String = String::deserialize(deserializer)?; + match s.as_str() { + "NaN" => Ok(f64::NAN), + "Infinity" => Ok(f64::INFINITY), + "-Infinity" => Ok(f64::NEG_INFINITY), + _ => s.parse().map_err(serde::de::Error::custom), + } +} + +#[cfg(test)] +mod tests { + use fgoxide::io::DelimFile; + use serde::{Deserialize, Serialize}; + use tempfile::NamedTempFile; + + #[derive(Debug, Serialize, Deserialize)] + struct Row { + #[serde(with = "super")] + value: f64, + } + + /// Writes a one-column TSV via the real `DelimFile` path and returns the value cell + /// (the second physical line, after the `value` header). + fn to_cell(value: f64) -> String { + let file = NamedTempFile::new().expect("temp file"); + DelimFile::default().write_tsv(file.path(), [Row { value }]).expect("write"); + let content = std::fs::read_to_string(file.path()).expect("read"); + content.lines().nth(1).expect("value row").to_string() + } + + #[test] + fn writes_fgbio_readable_non_finite_tokens() { + assert_eq!(to_cell(f64::INFINITY), "Infinity"); + assert_eq!(to_cell(f64::NEG_INFINITY), "-Infinity"); + assert_eq!(to_cell(f64::NAN), "NaN"); + } + + #[test] + fn writes_integral_values_without_decimal_point() { + assert_eq!(to_cell(0.0), "0"); + assert_eq!(to_cell(42.0), "42"); + assert_eq!(to_cell(-7.0), "-7"); + } + + #[test] + fn keeps_full_precision_for_fractional_values() { + assert_eq!(to_cell(0.5), "0.5"); + assert_eq!(to_cell(1.25), "1.25"); + } + + #[test] + fn round_trips_non_finite_and_finite() { + for &v in &[0.0_f64, 1.5, -3.25, 1000.0, f64::INFINITY, f64::NEG_INFINITY] { + let file = NamedTempFile::new().expect("temp file"); + DelimFile::default().write_tsv(file.path(), [Row { value: v }]).expect("write"); + let rows: Vec = DelimFile::default().read_tsv(file.path()).expect("read"); + assert_eq!(rows.len(), 1); + if v.is_infinite() { + assert!(rows[0].value.is_infinite() && rows[0].value.signum() == v.signum()); + } else { + assert!((rows[0].value - v).abs() < f64::EPSILON, "round-trip {v}"); + } + } + } +} diff --git a/crates/fgumi-metrics/src/group.rs b/crates/fgumi-metrics/src/group.rs index 0c4e6b9fb..8179977d9 100644 --- a/crates/fgumi-metrics/src/group.rs +++ b/crates/fgumi-metrics/src/group.rs @@ -39,14 +39,27 @@ fn build_size_distribution( /// Metrics for UMI grouping operations. /// -/// These metrics track how reads are grouped by UMI and provide insight into -/// data quality and molecule representation. +/// The **serialized** form matches fgbio's `UmiGroupingMetric` exactly — the same +/// five columns, in the same order, with the same (fgbio-spelled) names: +/// `accepted_sam_records`, `discarded_non_pf`, `discarded_poor_alignment`, +/// `discarded_ns_in_umi`, `discarded_umis_to_short`. This keeps the +/// `.grouping_metrics.txt` output readable by fgbio's `Metric.read`, which hard-fails +/// on any unexpected column. +/// +/// The remaining fields are `#[serde(skip)]` — they are computed and retained for +/// fgumi's own summary logging ([`crate::ProcessingMetrics`] and +/// `log_umi_grouping_summary`) but are **not** written to the metrics file. They are +/// all derivable from the accepted count and the family-size distribution, which is +/// emitted separately in `.family_sizes.txt`. #[derive(Debug, Clone, Default, Serialize, Deserialize)] pub struct UmiGroupingMetrics { - /// Total SAM records processed + /// Total SAM records processed (fgumi-internal; not part of fgbio's schema — + /// derivable as accepted + discarded). + #[serde(skip)] pub total_records: u64, - /// Records accepted for grouping + /// Records accepted for grouping (fgbio column `accepted_sam_records`). + #[serde(rename = "accepted_sam_records")] pub accepted_records: u64, /// Records discarded (not passing filter) @@ -58,25 +71,33 @@ pub struct UmiGroupingMetrics { /// Records discarded (Ns in UMI) pub discarded_ns_in_umi: u64, - /// Records discarded (UMI too short) + /// Records discarded (UMI too short). fgbio spells this column + /// `discarded_umis_to_short` (sic); the serialized name matches it exactly. + #[serde(rename = "discarded_umis_to_short")] pub discarded_umi_too_short: u64, - /// Number of unique molecule IDs assigned + /// Number of unique molecule IDs assigned (fgumi-internal summary only). + #[serde(skip)] pub unique_molecule_ids: u64, - /// Total number of UMI families/groups + /// Total number of UMI families/groups (fgumi-internal summary only). + #[serde(skip)] pub total_families: u64, - /// Average reads per molecule + /// Average reads per molecule (fgumi-internal summary only). + #[serde(skip)] pub avg_reads_per_molecule: f64, - /// Median reads per molecule + /// Median reads per molecule (fgumi-internal summary only). + #[serde(skip)] pub median_reads_per_molecule: u64, - /// Minimum reads per molecule + /// Minimum reads per molecule (fgumi-internal summary only). + #[serde(skip)] pub min_reads_per_molecule: u64, - /// Maximum reads per molecule + /// Maximum reads per molecule (fgumi-internal summary only). + #[serde(skip)] pub max_reads_per_molecule: u64, } @@ -123,9 +144,11 @@ pub struct FamilySizeMetrics { pub count: u64, /// Fraction of all families with this size + #[serde(with = "crate::float")] pub fraction: f64, /// Cumulative fraction (families with size >= this value) + #[serde(with = "crate::float")] pub fraction_gt_or_eq_family_size: f64, } @@ -176,9 +199,11 @@ pub struct PositionGroupSizeMetrics { pub count: u64, /// Fraction of all position groups with this size + #[serde(with = "crate::float")] pub fraction: f64, /// Cumulative fraction (position groups with size >= this value) + #[serde(with = "crate::float")] pub fraction_gt_or_eq_position_group_size: f64, } @@ -233,6 +258,52 @@ mod tests { assert_eq!(metrics.unique_molecule_ids, 0); } + #[test] + fn test_umi_grouping_metrics_serializes_fgbio_five_columns() { + use fgoxide::io::DelimFile; + use tempfile::NamedTempFile; + + let metrics = vec![UmiGroupingMetrics { + total_records: 100, + accepted_records: 90, + discarded_non_pf: 4, + discarded_poor_alignment: 3, + discarded_ns_in_umi: 2, + discarded_umi_too_short: 1, + // fgumi-internal fields that must NOT be serialized: + unique_molecule_ids: 42, + total_families: 42, + avg_reads_per_molecule: 2.5, + median_reads_per_molecule: 2, + min_reads_per_molecule: 1, + max_reads_per_molecule: 9, + }]; + + let file = NamedTempFile::new().expect("temp file"); + DelimFile::default().write_tsv(file.path(), &metrics).expect("write"); + let content = std::fs::read_to_string(file.path()).expect("read"); + let mut lines = content.lines(); + + // Header and row must be exactly fgbio's 5-column `UmiGroupingMetric`, in order. + assert_eq!( + lines.next(), + Some( + "accepted_sam_records\tdiscarded_non_pf\tdiscarded_poor_alignment\t\ + discarded_ns_in_umi\tdiscarded_umis_to_short" + ) + ); + assert_eq!(lines.next(), Some("90\t4\t3\t2\t1")); + assert_eq!(lines.next(), None); + + // And it reads back (fgumi-internal fields default to zero, fgbio columns intact). + let read_back: Vec = + DelimFile::default().read_tsv(file.path()).expect("read back"); + assert_eq!(read_back.len(), 1); + assert_eq!(read_back[0].accepted_records, 90); + assert_eq!(read_back[0].discarded_umi_too_short, 1); + assert_eq!(read_back[0].total_records, 0); + } + #[test] fn test_umi_grouping_metrics_default() { let metrics = UmiGroupingMetrics::default(); diff --git a/crates/fgumi-metrics/src/lib.rs b/crates/fgumi-metrics/src/lib.rs index 43d80f508..742623489 100644 --- a/crates/fgumi-metrics/src/lib.rs +++ b/crates/fgumi-metrics/src/lib.rs @@ -13,6 +13,7 @@ pub mod clip; pub mod consensus; pub mod correct; pub mod duplex; +pub(crate) mod float; pub mod group; pub mod rejection; pub mod shared; @@ -21,13 +22,20 @@ pub mod writer; use serde::{Deserialize, Serialize}; -/// Number of decimal places used for float metrics (matches fgbio). +/// Number of decimal places used by [`format_float`] for fixed-precision float metrics. pub const FLOAT_PRECISION: usize = 6; -/// Formats a float value with the standard precision for metrics. +/// Formats a float with a fixed 6 decimal places (e.g. `0.9` → `"0.900000"`). /// -/// This ensures consistent float formatting across all metrics output, -/// matching fgbio's 6 decimal place precision. +/// This does **not** reproduce fgbio's textual format. fgbio uses +/// `DecimalFormat("0.######")`, which trims trailing zeros (`0.9` → `"0.9"`, +/// `0.0` → `"0"`) and can switch to scientific notation for very small values. The +/// metric contract between fgumi and fgbio is *numeric* equality — `fgumi compare +/// metrics` parses both sides and compares with a tolerance — so this +/// representation difference is accepted (see the issue's "Accepted (not a bug)" +/// note). f64 fields that must be *read back* by fgbio use +/// [`float::serialize_f64`] instead, which additionally renders non-finite values +/// as fgbio-parseable `Infinity`/`-Infinity`/`NaN` tokens. /// /// # Example /// ``` diff --git a/crates/fgumi-metrics/src/rejection.rs b/crates/fgumi-metrics/src/rejection.rs index 97e59907b..673cdb889 100644 --- a/crates/fgumi-metrics/src/rejection.rs +++ b/crates/fgumi-metrics/src/rejection.rs @@ -14,7 +14,7 @@ use std::fmt; pub enum RejectionReason { /// Insufficient reads to generate a consensus InsufficientSupport, - /// Read has a different, and minority, set of indels + /// Reads has a different, and minority, set of indels MinorityAlignment, /// Too few reads agreed on the strand orientation InsufficientStrandSupport, @@ -56,7 +56,8 @@ impl RejectionReason { pub fn description(&self) -> &'static str { match self { Self::InsufficientSupport => "Insufficient reads to generate a consensus", - Self::MinorityAlignment => "Read has a different, and minority, set of indels", + // fgbio spells this "Reads has …" (sic) — matched verbatim for metric parity. + Self::MinorityAlignment => "Reads has a different, and minority, set of indels", Self::InsufficientStrandSupport => "Too few reads agreed on the strand orientation", Self::LowBaseQuality => "Base quality scores were below threshold", Self::ExcessiveNBases => "Read group had too many N bases", @@ -106,7 +107,7 @@ impl RejectionReason { pub fn kv_description(&self) -> &'static str { match self { Self::InsufficientSupport => "Insufficient reads to generate a consensus", - Self::MinorityAlignment => "Read has a different, and minority, set of indels", + Self::MinorityAlignment => "Reads has a different, and minority, set of indels", Self::InsufficientStrandSupport => "Insufficient strand support for consensus", Self::LowBaseQuality => "Low base quality", Self::ExcessiveNBases => "Excessive N bases in read", @@ -119,7 +120,8 @@ impl RejectionReason { Self::InsufficientMinDepth => "Insufficient minimum read depth", Self::ExcessiveErrorRate => "Excessive error rate", Self::UmiTooShort => "UMI sequence too short", - Self::SameStrandOnly => "Only generating one strand of duplex consensus", + // fgbio's exact title-cased string for `single_strand_only`, matched for parity. + Self::SameStrandOnly => "Only Generating One Strand of Duplex Consensus", Self::DuplicateUmi => "Duplicate UMI detected", Self::OrphanConsensus => "Only one of R1 or R2 consensus generated", Self::ZeroBasesPostTrimming => "Read or mate had zero bases post trimming", @@ -167,9 +169,10 @@ mod tests { fn test_rejection_reason_description() { assert!(RejectionReason::LowBaseQuality.description().contains("quality")); assert!(RejectionReason::InsufficientSupport.description().contains("Insufficient")); + // Matches fgbio's exact (grammatically-odd) wording for metric parity. assert_eq!( RejectionReason::MinorityAlignment.to_string(), - "Read has a different, and minority, set of indels" + "Reads has a different, and minority, set of indels" ); } diff --git a/crates/fgumi-metrics/src/shared.rs b/crates/fgumi-metrics/src/shared.rs index ec3a8bdd7..e1be58ed5 100644 --- a/crates/fgumi-metrics/src/shared.rs +++ b/crates/fgumi-metrics/src/shared.rs @@ -24,8 +24,10 @@ pub struct UmiMetric { /// Number of double-stranded tag families observing this UMI pub unique_observations: usize, /// Fraction of all raw observations + #[serde(with = "crate::float")] pub fraction_raw_observations: f64, /// Fraction of all unique observations + #[serde(with = "crate::float")] pub fraction_unique_observations: f64, } diff --git a/crates/fgumi-metrics/src/simplex.rs b/crates/fgumi-metrics/src/simplex.rs index e7ae1ffc9..726f382a7 100644 --- a/crates/fgumi-metrics/src/simplex.rs +++ b/crates/fgumi-metrics/src/simplex.rs @@ -22,14 +22,18 @@ pub struct SimplexFamilySizeMetric { /// Count of CS families with this size pub cs_count: usize, /// Fraction of all CS families with this size + #[serde(with = "crate::float")] pub cs_fraction: f64, /// Fraction of CS families with size >= `family_size` + #[serde(with = "crate::float")] pub cs_fraction_gt_or_eq_size: f64, /// Count of SS families with this size pub ss_count: usize, /// Fraction of all SS families with this size + #[serde(with = "crate::float")] pub ss_fraction: f64, /// Fraction of SS families with size >= `family_size` + #[serde(with = "crate::float")] pub ss_fraction_gt_or_eq_size: f64, } @@ -65,6 +69,7 @@ impl Metric for SimplexFamilySizeMetric { #[derive(Debug, Clone, Serialize, Deserialize)] pub struct SimplexYieldMetric { /// Approximate fraction of full dataset used + #[serde(with = "crate::float")] pub fraction: f64, /// Number of read pairs upon which metrics are based pub read_pairs: usize, @@ -73,10 +78,12 @@ pub struct SimplexYieldMetric { /// Number of SS (Single-Strand by UMI) families pub ss_families: usize, /// Mean SS family size + #[serde(with = "crate::float")] pub mean_ss_family_size: f64, /// Number of SS singleton families (size 1) pub ss_singletons: usize, /// Fraction of SS families that are singletons + #[serde(with = "crate::float")] pub ss_singleton_fraction: f64, /// Number of SS families with size >= consensus minimum pub ss_consensus_families: usize, diff --git a/crates/fgumi-metrics/src/writer.rs b/crates/fgumi-metrics/src/writer.rs index 852e567f4..af7afc62e 100644 --- a/crates/fgumi-metrics/src/writer.rs +++ b/crates/fgumi-metrics/src/writer.rs @@ -13,11 +13,19 @@ use crate::Metric; /// Write metrics to a TSV file with consistent error handling. /// /// This is a convenience wrapper around `DelimFile::write_tsv` that provides -/// consistent error messages across all commands. +/// consistent error messages across all commands, and — unlike a bare +/// `DelimFile::write_tsv` — **always writes the column header**, even for an empty +/// slice. +/// +/// The underlying csv writer emits headers lazily (on the first record), so writing +/// a zero-row slice through it produces a 0-byte file. fgbio writes the header +/// eagerly in its `Metric` writer, so fgbio's `Metric.read` rejects such a file with +/// "No header found". Emitting a header-only file keeps an empty metrics output +/// re-readable by both fgbio and fgumi (it parses back to zero rows). /// /// # Arguments /// * `path` - Path to the output TSV file -/// * `metrics` - The metrics to write (must implement Serialize) +/// * `metrics` - The metrics to write (must implement `Serialize` and `Default`) /// * `description` - Human-readable description of the metrics for error messages /// /// # Errors @@ -29,7 +37,7 @@ use crate::Metric; /// use serde::Serialize; /// use std::path::Path; /// -/// #[derive(Serialize)] +/// #[derive(Serialize, Default)] /// struct MyMetrics { /// count: usize, /// value: f64, @@ -42,17 +50,36 @@ use crate::Metric; /// /// write_metrics(Path::new("metrics.txt"), &metrics, "processing").unwrap(); /// ``` -pub fn write_metrics, T: Serialize>( +pub fn write_metrics, T: Serialize + Default>( path: P, metrics: &[T], description: &str, ) -> Result<()> { let path_ref = path.as_ref(); + if metrics.is_empty() { + return write_header_only::<_, T>(path_ref).with_context(|| { + format!("Failed to write {} metrics header: {}", description, path_ref.display()) + }); + } DelimFile::default() .write_tsv(path_ref, metrics) .with_context(|| format!("Failed to write {} metrics: {}", description, path_ref.display())) } +/// Writes a header-only metrics file for a metric type with no rows. +/// +/// Serializes a single `T::default()` row through the normal `DelimFile` path so the +/// header is formatted identically to a populated file, then truncates the file to +/// just that header line. The result parses back to zero metrics. +fn write_header_only, T: Serialize + Default>(path: P) -> Result<()> { + let path_ref = path.as_ref(); + DelimFile::default().write_tsv(path_ref, [T::default()])?; + let content = std::fs::read_to_string(path_ref)?; + let header = content.lines().next().unwrap_or_default(); + std::fs::write(path_ref, format!("{header}\n"))?; + Ok(()) +} + /// Write metrics implementing the Metric trait to a TSV file. /// /// This version uses the metric's own name for error messages, providing @@ -163,14 +190,22 @@ mod tests { } #[test] - fn test_write_metrics_empty() -> Result<()> { + fn test_write_metrics_empty_is_header_only_and_round_trips() -> Result<()> { let temp_file = NamedTempFile::new()?; let metrics: Vec = vec![]; write_metrics(temp_file.path(), &metrics, "empty")?; - // Verify the file was created (even if empty/header-only) - assert!(temp_file.path().exists()); + // An empty metrics slice must still produce the column header (fgbio's + // Metric.read rejects a 0-byte file with "No header found"), and no data rows. + let content = fs::read_to_string(temp_file.path())?; + let lines: Vec<&str> = content.lines().collect(); + assert_eq!(lines.len(), 1, "expected header-only file, got: {content:?}"); + assert_eq!(lines[0], "name\tcount\tvalue"); + + // And it round-trips back to zero rows. + let read_back: Vec = read_metrics(temp_file.path(), "empty")?; + assert!(read_back.is_empty()); Ok(()) } diff --git a/src/lib/commands/clip.rs b/src/lib/commands/clip.rs index 674c2a53d..da2bc4e9f 100644 --- a/src/lib/commands/clip.rs +++ b/src/lib/commands/clip.rs @@ -367,7 +367,7 @@ impl Clip { if let Some(mut metrics) = metrics_collection { info!("Writing metrics to {}", metrics_path.display()); metrics.finalize(); - let all_metrics = metrics.all_metrics(); + let all_metrics = metrics.all_metrics().map(Clone::clone); write_metrics_tsv(metrics_path, &all_metrics, "clipping")?; info!("Metrics written successfully"); } diff --git a/src/lib/commands/duplex_metrics.rs b/src/lib/commands/duplex_metrics.rs index 21062b174..a5f6d0f4d 100644 --- a/src/lib/commands/duplex_metrics.rs +++ b/src/lib/commands/duplex_metrics.rs @@ -12,9 +12,8 @@ use crate::metrics::duplex::{DuplexMetricsCollector, DuplexYieldMetric, FamilySi use crate::simple_umi_consensus::SimpleUmiConsensusCaller; use crate::umi::extract_mi_base; use crate::validation::validate_file_exists; -use anyhow::{Context, Result}; +use anyhow::Result; use clap::Parser; -use fgoxide::io::DelimFile; use log::info; use statrs::distribution::{Binomial, DiscreteCDF}; use std::path::PathBuf; @@ -200,30 +199,26 @@ impl Command for DuplexMetrics { // Generate and write metrics info!("Writing metrics..."); - // Family size metrics + // Family size metrics. Use `write_metrics_auto` (not a bare `DelimFile::write_tsv`) + // so an empty distribution still produces a header-only, fgbio-readable file. let family_size_metrics = main_collector.family_size_metrics(); let family_size_path = format!("{}.family_sizes.txt", self.output.display()); - DelimFile::default() - .write_tsv(&family_size_path, family_size_metrics) - .with_context(|| format!("Failed to write family size metrics: {family_size_path}"))?; + crate::metrics::writer::write_metrics_auto(&family_size_path, &family_size_metrics)?; info!("Wrote family size metrics to {family_size_path}"); // Duplex family size metrics let duplex_family_size_metrics = main_collector.duplex_family_size_metrics(); let duplex_family_size_path = format!("{}.duplex_family_sizes.txt", self.output.display()); - DelimFile::default() - .write_tsv(&duplex_family_size_path, duplex_family_size_metrics) - .with_context(|| { - format!("Failed to write duplex family size metrics: {duplex_family_size_path}") - })?; + crate::metrics::writer::write_metrics_auto( + &duplex_family_size_path, + &duplex_family_size_metrics, + )?; info!("Wrote duplex family size metrics to {duplex_family_size_path}"); // UMI metrics let umi_metrics = main_collector.umi_metrics(); let umi_path = format!("{}.umi_counts.txt", self.output.display()); - DelimFile::default() - .write_tsv(&umi_path, umi_metrics) - .with_context(|| format!("Failed to write UMI metrics: {umi_path}"))?; + crate::metrics::writer::write_metrics_auto(&umi_path, &umi_metrics)?; info!("Wrote UMI metrics to {umi_path}"); // Duplex UMI metrics (if enabled) @@ -231,17 +226,13 @@ impl Command for DuplexMetrics { let duplex_umi_metrics = main_collector.duplex_umi_metrics(&main_collector.umi_metrics()); let duplex_umi_path = format!("{}.duplex_umi_counts.txt", self.output.display()); - DelimFile::default().write_tsv(&duplex_umi_path, duplex_umi_metrics).with_context( - || format!("Failed to write duplex UMI metrics: {duplex_umi_path}"), - )?; + crate::metrics::writer::write_metrics_auto(&duplex_umi_path, &duplex_umi_metrics)?; info!("Wrote duplex UMI metrics to {duplex_umi_path}"); } // Yield metrics let yield_path = format!("{}.duplex_yield_metrics.txt", self.output.display()); - DelimFile::default() - .write_tsv(&yield_path, yield_metrics) - .with_context(|| format!("Failed to write yield metrics: {yield_path}"))?; + crate::metrics::writer::write_metrics_auto(&yield_path, &yield_metrics)?; info!("Wrote yield metrics to {yield_path}"); // Generate PDF plots using R script (optional) diff --git a/src/lib/commands/group.rs b/src/lib/commands/group.rs index 941e69c94..32b59adab 100644 --- a/src/lib/commands/group.rs +++ b/src/lib/commands/group.rs @@ -2771,12 +2771,17 @@ mod tests { ); // ---- Grouping metrics ---- - // 9 read pairs = 18 records accepted, 6 unique families + // 9 read pairs = 18 records accepted, 6 unique families. + // The `.grouping_metrics.txt` matches fgbio's 5-column `UmiGroupingMetric`; + // `total_records`/`total_families` are fgumi-internal and no longer serialized + // (family count is covered by the `.family_sizes.txt` assertions above). let grouping: Vec = DelimFile::default().read_tsv(&grouping_path)?; assert_eq!(grouping.len(), 1); - assert_eq!(grouping[0].total_records, 18); assert_eq!(grouping[0].accepted_records, 18); - assert_eq!(grouping[0].total_families, 6); + assert_eq!(grouping[0].discarded_non_pf, 0); + assert_eq!(grouping[0].discarded_poor_alignment, 0); + assert_eq!(grouping[0].discarded_ns_in_umi, 0); + assert_eq!(grouping[0].discarded_umi_too_short, 0); // ---- Position group size histogram ---- // 3 position groups: one with 1 family, one with 2, one with 3 @@ -2909,11 +2914,11 @@ mod tests { // Grouping counters: UmiGroupingMetrics has no PartialEq, so compare // the integer fields that the per-thread accumulators feed. + // The serialized metrics match fgbio's 5-column `UmiGroupingMetric`; the input + // total (20 = 18 accepted + 2 discarded) is fgumi-internal and no longer a column. let grouping: Vec = DelimFile::default().read_tsv(&grouping_path)?; assert_eq!(grouping.len(), 1); - assert_eq!(grouping[0].total_records, 20); assert_eq!(grouping[0].accepted_records, 18); - assert_eq!(grouping[0].total_families, 6); assert_eq!(grouping[0].discarded_non_pf, 0); assert_eq!(grouping[0].discarded_poor_alignment, 0); assert_eq!(grouping[0].discarded_ns_in_umi, 2); @@ -2999,9 +3004,22 @@ mod tests { DelimFile::default().read_tsv(&prefix_grouping)?; assert_eq!(individual_grouping.len(), 1); assert_eq!(prefix_grouping.len(), 1); - assert_eq!(individual_grouping[0].total_records, prefix_grouping[0].total_records); + // Compare the serialized fgbio columns (total_records/total_families are + // fgumi-internal and no longer written). assert_eq!(individual_grouping[0].accepted_records, prefix_grouping[0].accepted_records); - assert_eq!(individual_grouping[0].total_families, prefix_grouping[0].total_families); + assert_eq!(individual_grouping[0].discarded_non_pf, prefix_grouping[0].discarded_non_pf); + assert_eq!( + individual_grouping[0].discarded_poor_alignment, + prefix_grouping[0].discarded_poor_alignment + ); + assert_eq!( + individual_grouping[0].discarded_ns_in_umi, + prefix_grouping[0].discarded_ns_in_umi + ); + assert_eq!( + individual_grouping[0].discarded_umi_too_short, + prefix_grouping[0].discarded_umi_too_short + ); // Position group sizes: only written via --metrics prefix let position: Vec = diff --git a/src/lib/commands/simplex_metrics.rs b/src/lib/commands/simplex_metrics.rs index 5e21cc6e2..71bd2fa0c 100644 --- a/src/lib/commands/simplex_metrics.rs +++ b/src/lib/commands/simplex_metrics.rs @@ -11,9 +11,8 @@ use crate::logging::OperationTimer; use crate::metrics::simplex::{SimplexMetricsCollector, SimplexYieldMetric}; use crate::simple_umi_consensus::SimpleUmiConsensusCaller; use crate::validation::validate_file_exists; -use anyhow::{Context, Result}; +use anyhow::Result; use clap::Parser; -use fgoxide::io::DelimFile; use log::info; use std::path::PathBuf; @@ -157,27 +156,22 @@ impl Command for SimplexMetrics { // Generate and write metrics info!("Writing metrics..."); - // Family size metrics + // Family size metrics. Use `write_metrics_auto` (not a bare `DelimFile::write_tsv`) + // so an empty distribution still produces a header-only, fgbio-readable file. let family_size_metrics = main_collector.family_size_metrics(); let family_size_path = format!("{}.family_sizes.txt", self.output.display()); - DelimFile::default() - .write_tsv(&family_size_path, family_size_metrics) - .with_context(|| format!("Failed to write family size metrics: {family_size_path}"))?; + crate::metrics::writer::write_metrics_auto(&family_size_path, &family_size_metrics)?; info!("Wrote family size metrics to {family_size_path}"); // UMI metrics let umi_metrics = main_collector.umi_metrics(); let umi_path = format!("{}.umi_counts.txt", self.output.display()); - DelimFile::default() - .write_tsv(&umi_path, umi_metrics) - .with_context(|| format!("Failed to write UMI metrics: {umi_path}"))?; + crate::metrics::writer::write_metrics_auto(&umi_path, &umi_metrics)?; info!("Wrote UMI metrics to {umi_path}"); // Yield metrics let yield_path = format!("{}.simplex_yield_metrics.txt", self.output.display()); - DelimFile::default() - .write_tsv(&yield_path, yield_metrics) - .with_context(|| format!("Failed to write yield metrics: {yield_path}"))?; + crate::metrics::writer::write_metrics_auto(&yield_path, &yield_metrics)?; info!("Wrote yield metrics to {yield_path}"); // Generate PDF plots using R script (optional) diff --git a/src/lib/metrics/mod.rs b/src/lib/metrics/mod.rs index 09bbbb91f..633462104 100644 --- a/src/lib/metrics/mod.rs +++ b/src/lib/metrics/mod.rs @@ -36,4 +36,4 @@ pub use duplex::{ }; pub use group::{FamilySizeMetrics, PositionGroupSizeMetrics, UmiGroupingMetrics}; pub use shared::UmiMetric; -pub use writer::{read_metrics, write_metrics}; +pub use writer::{read_metrics, write_metrics, write_metrics_auto}; diff --git a/tests/integration/test_group_command.rs b/tests/integration/test_group_command.rs index 3870f01cd..d0d642f6b 100644 --- a/tests/integration/test_group_command.rs +++ b/tests/integration/test_group_command.rs @@ -55,10 +55,11 @@ fn test_group_command_writes_new_metrics() { let metric = &metrics[0]; - // Verify basic fields are populated - assert!(metric.total_records > 0, "total_records should be positive"); + // Verify the serialized fgbio columns are populated. The `.grouping_metrics.txt` + // matches fgbio's 5-column `UmiGroupingMetric`; `total_records`/`unique_molecule_ids` + // are fgumi-internal (`#[serde(skip)]`) and read back as zero, so they are not asserted + // here. assert!(metric.accepted_records > 0, "accepted_records should be positive"); - assert!(metric.unique_molecule_ids > 0, "unique_molecule_ids should be positive"); // Rejection fields are validated by their presence in the struct (no runtime check needed // since they are unsigned integers that are always >= 0) From 732b6c3b757d2b00ac04fa8e75c3ab662ddecf93 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Wed, 8 Jul 2026 22:13:17 -0700 Subject: [PATCH 2/4] fix(consensus): map fragment/collision rejections to fgbio non_paired_reads/potential_umi_collision fgbio's UmiConsensusCaller rejects unpaired/fragment reads as `non_paired_reads` and potential MI collisions between independent duplex molecules as `potential_umi_collision` (UmiConsensusCaller.scala:81,83). fgumi's caller-side `RejectionReason::FragmentRead` and `PotentialCollision` were mapped by `to_centralized()` to the wrong centralized reasons (`single_strand_only` and `duplicate_umi`), so codec fragment rejections and duplex collision rejections were emitted under the wrong `raw_reads_rejected_for_*` metric keys. Add `NonPairedReads` and `PotentialUmiCollision` to the metrics `RejectionReason` (fgbio-exact code + description), wire them through `ConsensusMetrics` (count/add/emit), and correct the `to_centralized()` mapping so codec now emits `non_paired_reads` and duplex now emits `potential_umi_collision`. Also adds `RejectionReason::used_by_duplex()` (mirroring fgbio's `usedByDuplex` flags) for the follow-on per-caller rejection-row seeding (R2-MET-05). Part of R2-UCC-01. Branch stacked on #504 (nh/fix-metrics-fgbio-parity). --- crates/fgumi-consensus/src/caller.rs | 14 ++++++----- crates/fgumi-metrics/src/consensus.rs | 20 +++++++++++++++- crates/fgumi-metrics/src/rejection.rs | 34 +++++++++++++++++++++++++++ 3 files changed, 61 insertions(+), 7 deletions(-) diff --git a/crates/fgumi-consensus/src/caller.rs b/crates/fgumi-consensus/src/caller.rs index ac3ecb62e..af5ac4fea 100644 --- a/crates/fgumi-consensus/src/caller.rs +++ b/crates/fgumi-consensus/src/caller.rs @@ -448,16 +448,18 @@ impl RejectionReason { Self::QualityTooLow | Self::QualityTrimmed => CentralReason::LowBaseQuality, Self::FailedQC => CentralReason::NotPassingFilter, Self::MissingUmi => CentralReason::MissingUmi, - // For reasons that don't have a direct centralized equivalent, - // we use the closest semantic match - Self::FragmentRead => CentralReason::SameStrandOnly, + // fgbio maps unpaired/fragment reads to `non_paired_reads` and MI + // collisions to `potential_umi_collision` + // (UmiConsensusCaller.scala:81,83), not to `single_strand_only` / + // `duplicate_umi` — match fgbio's rejection accounting. + Self::FragmentRead => CentralReason::NonPairedReads, Self::Unmapped | Self::Mapped | Self::SecondaryOrSupplementary | Self::IndelErrorBetweenStrands => CentralReason::NoValidAlignment, Self::ZeroLengthAfterTrimming => CentralReason::ZeroBasesPostTrimming, Self::OrphanConsensus => CentralReason::OrphanConsensus, - Self::PotentialCollision => CentralReason::DuplicateUmi, + Self::PotentialCollision => CentralReason::PotentialUmiCollision, } } } @@ -1067,7 +1069,7 @@ mod tests { // Test every variant maps to a centralized reason let mappings = [ - (RejectionReason::FragmentRead, CentralReason::SameStrandOnly), + (RejectionReason::FragmentRead, CentralReason::NonPairedReads), (RejectionReason::InsufficientReads, CentralReason::InsufficientSupport), (RejectionReason::QualityTooLow, CentralReason::LowBaseQuality), (RejectionReason::Unmapped, CentralReason::NoValidAlignment), @@ -1082,7 +1084,7 @@ mod tests { (RejectionReason::InsufficientOverlap, CentralReason::InsufficientSupport), (RejectionReason::OrphanConsensus, CentralReason::OrphanConsensus), (RejectionReason::IndelErrorBetweenStrands, CentralReason::NoValidAlignment), - (RejectionReason::PotentialCollision, CentralReason::DuplicateUmi), + (RejectionReason::PotentialCollision, CentralReason::PotentialUmiCollision), (RejectionReason::Other, CentralReason::InsufficientSupport), ]; diff --git a/crates/fgumi-metrics/src/consensus.rs b/crates/fgumi-metrics/src/consensus.rs index bf787bbd0..ae28e52f2 100644 --- a/crates/fgumi-metrics/src/consensus.rs +++ b/crates/fgumi-metrics/src/consensus.rs @@ -123,6 +123,14 @@ pub struct ConsensusMetrics { /// Reads rejected due to zero bases after trimming pub rejected_zero_bases_post_trimming: u64, + + /// Reads rejected due to being unpaired/fragment reads (fgbio + /// `non_paired_reads`; duplex/codec) + pub rejected_non_paired_reads: u64, + + /// Reads rejected due to a potential collision between independent duplex + /// molecules (fgbio `potential_umi_collision`; duplex/codec) + pub rejected_potential_umi_collision: u64, } /// Core rejection reasons always included in KV metrics output. @@ -134,7 +142,7 @@ const CORE_REJECTIONS: [RejectionReason; 4] = [ ]; /// Optional rejection reasons included in KV metrics output only when non-zero. -const OPTIONAL_REJECTIONS: [RejectionReason; 14] = [ +const OPTIONAL_REJECTIONS: [RejectionReason; 16] = [ RejectionReason::InsufficientStrandSupport, RejectionReason::LowBaseQuality, RejectionReason::ExcessiveNBases, @@ -149,6 +157,8 @@ const OPTIONAL_REJECTIONS: [RejectionReason; 14] = [ RejectionReason::UmiTooShort, RejectionReason::SameStrandOnly, RejectionReason::DuplicateUmi, + RejectionReason::NonPairedReads, + RejectionReason::PotentialUmiCollision, ]; /// Returns an iterator over all rejection reasons (core + optional). @@ -189,6 +199,8 @@ impl ConsensusMetrics { rejected_duplicate_umi: 0, rejected_orphan_consensus: 0, rejected_zero_bases_post_trimming: 0, + rejected_non_paired_reads: 0, + rejected_potential_umi_collision: 0, } } @@ -214,6 +226,8 @@ impl ConsensusMetrics { RejectionReason::DuplicateUmi => self.rejected_duplicate_umi, RejectionReason::OrphanConsensus => self.rejected_orphan_consensus, RejectionReason::ZeroBasesPostTrimming => self.rejected_zero_bases_post_trimming, + RejectionReason::NonPairedReads => self.rejected_non_paired_reads, + RejectionReason::PotentialUmiCollision => self.rejected_potential_umi_collision, } } @@ -244,6 +258,10 @@ impl ConsensusMetrics { RejectionReason::ZeroBasesPostTrimming => { self.rejected_zero_bases_post_trimming += count; } + RejectionReason::NonPairedReads => self.rejected_non_paired_reads += count, + RejectionReason::PotentialUmiCollision => { + self.rejected_potential_umi_collision += count; + } } } diff --git a/crates/fgumi-metrics/src/rejection.rs b/crates/fgumi-metrics/src/rejection.rs index 673cdb889..e4752ed8c 100644 --- a/crates/fgumi-metrics/src/rejection.rs +++ b/crates/fgumi-metrics/src/rejection.rs @@ -48,6 +48,12 @@ pub enum RejectionReason { OrphanConsensus, /// Read had zero bases after quality trimming ZeroBasesPostTrimming, + /// Unpaired/fragment reads given to the duplex/codec caller (fgbio + /// `non_paired_reads`; used by duplex and codec) + NonPairedReads, + /// Potential collision between independent duplex molecules (fgbio + /// `potential_umi_collision`; used by duplex and codec) + PotentialUmiCollision, } impl RejectionReason { @@ -74,6 +80,10 @@ impl RejectionReason { Self::DuplicateUmi => "Duplicate UMI at same genomic position", Self::OrphanConsensus => "Only one of R1 or R2 consensus generated", Self::ZeroBasesPostTrimming => "Read or mate had zero bases post trimming", + Self::NonPairedReads => "Unpaired/fragment reads not supported by Duplex caller", + Self::PotentialUmiCollision => { + "Potential collision between independent duplex molecules" + } } } @@ -99,6 +109,8 @@ impl RejectionReason { Self::DuplicateUmi => "raw_reads_rejected_for_duplicate_umi", Self::OrphanConsensus => "raw_reads_rejected_for_orphan_consensus", Self::ZeroBasesPostTrimming => "raw_reads_rejected_for_zero_bases_post_trimming", + Self::NonPairedReads => "raw_reads_rejected_for_non_paired_reads", + Self::PotentialUmiCollision => "raw_reads_rejected_for_potential_umi_collision", } } @@ -125,8 +137,30 @@ impl RejectionReason { Self::DuplicateUmi => "Duplicate UMI detected", Self::OrphanConsensus => "Only one of R1 or R2 consensus generated", Self::ZeroBasesPostTrimming => "Read or mate had zero bases post trimming", + Self::NonPairedReads => "Unpaired/fragment reads not supported by Duplex caller", + Self::PotentialUmiCollision => { + "Potential collision between independent duplex molecules" + } } } + + /// Whether the fgbio duplex caller emits (seeds) this rejection reason. + /// Mirrors fgbio's `usedByDuplex` flags on `RejectionReason` + /// (`UmiConsensusCaller.scala:60-88`) for the reasons fgumi shares with + /// fgbio; used to always-emit (seed to zero) the corresponding KV rows. + #[must_use] + pub fn used_by_duplex(&self) -> bool { + matches!( + self, + Self::InsufficientSupport + | Self::MinorityAlignment + | Self::OrphanConsensus + | Self::ZeroBasesPostTrimming + | Self::NonPairedReads + | Self::SameStrandOnly + | Self::PotentialUmiCollision + ) + } } impl fmt::Display for RejectionReason { From d04e3d9e0ecf6c3aaf2caa641a17bb99b5ff4a96 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Wed, 8 Jul 2026 22:37:03 -0700 Subject: [PATCH 3/4] fix(consensus): count duplex fragments as non_paired_reads, seed duplex rejection rows MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two fgbio parity fixes in the duplex caller and its rejection metrics. R2-UCC-01 — duplex fragment accounting. fgbio partitions unpaired/fragment reads out and rejects them as `non_paired_reads` (DuplexConsensusCaller.scala:204-205). fgumi dropped them from the R1/R2 strand groups but folded them into the input total; they are now counted as `FragmentRead` (→ `non_paired_reads`). R2-MET-05 — always-emit duplex rejection rows. fgbio seeds each caller's `usedBy*` rejection reasons to zero so their rows always appear (initializeRejectCounts). Add `ConsensusMetrics::to_kv_metrics_seeded` + `DUPLEX_SEEDED_REJECTIONS`; the duplex command seeds `non_paired_reads`, `single_strand_only`, and `potential_umi_collision` so those KV rows are always present. R2-UCC-02 (fgbio's distinct cell-barcode `require`) is deliberately NOT ported: fgumi's grouper keys on `(MI, cell)`, so a molecule spanning multiple cell barcodes is split into per-cell groups and consensus-called per cell (lossless), whereas fgbio's `require(barcodes.length <= 1)` fails fast. The caller therefore never sees more than one cell barcode via the CLI (verified: a same-MI/two-CB molecule yields two per-cell consensus reads), so a caller-level guard would be unreachable dead code. This is a deliberate, lossless architectural divergence. Tests: `to_kv_metrics_seeded` always-emits-at-zero and is not duplicated when non-zero. Full suite passes; ci-fmt/ci-lint clean. Branch stacked on #504. --- crates/fgumi-consensus/src/duplex_caller.rs | 20 ++++- crates/fgumi-metrics/src/consensus.rs | 85 ++++++++++++++++++++- src/lib/commands/duplex.rs | 6 +- 3 files changed, 106 insertions(+), 5 deletions(-) diff --git a/crates/fgumi-consensus/src/duplex_caller.rs b/crates/fgumi-consensus/src/duplex_caller.rs index 9a4a26b85..fa7af35f8 100644 --- a/crates/fgumi-consensus/src/duplex_caller.rs +++ b/crates/fgumi-consensus/src/duplex_caller.rs @@ -1771,7 +1771,11 @@ impl DuplexConsensusCaller { return Ok((ConsensusOutput::default(), stats, rejected_raw)); } - // Extract cell barcode from source reads (matching fgbio: recs.head.get[String](cellTag)) + // Extract the cell barcode from the source reads. fgumi's grouper keys + // on (MI, cell), so every record in this molecule already shares one + // cell barcode — taking the first is correct. (fgbio instead computes + // `distinct` + `require(<=1)`; fgumi's upstream per-cell partition is a + // deliberate lossless divergence, so the caller never sees >1 cell.) let cell_barcode: Option = cell_tag.and_then(|tag| { let tag_bytes: [u8; 2] = <[u8; 2]>::from(tag); a_records.first().or_else(|| b_records.first()).and_then(|r| { @@ -1782,6 +1786,20 @@ impl DuplexConsensusCaller { }) }); + // fgbio partitions out unpaired/fragment reads and rejects them as + // `non_paired_reads` (DuplexConsensusCaller.scala:204-205). fgumi drops + // them from the R1/R2 strand groups below; count them here as + // `FragmentRead` so they are accounted for in the rejection metrics + // rather than silently folded into the input total. + let fragment_count = a_records + .iter() + .chain(b_records.iter()) + .filter(|r| RawRecordView::new(r.as_ref()).flags() & flags::PAIRED == 0) + .count(); + if fragment_count > 0 { + stats.record_rejection(RejectionReason::FragmentRead, fragment_count); + } + // Split reads into R1/R2 groups for AB and BA strands (done once, reused below). // Require PAIRED for both: otherwise unpaired/fragment reads (which have both // segment bits clear) would be misclassified as R2. diff --git a/crates/fgumi-metrics/src/consensus.rs b/crates/fgumi-metrics/src/consensus.rs index ae28e52f2..8f5a61321 100644 --- a/crates/fgumi-metrics/src/consensus.rs +++ b/crates/fgumi-metrics/src/consensus.rs @@ -166,6 +166,15 @@ fn all_rejections() -> impl Iterator { CORE_REJECTIONS.into_iter().chain(OPTIONAL_REJECTIONS) } +/// Rejection reasons the duplex caller always emits (seeds to zero even when the +/// count is zero), mirroring fgbio's `usedByDuplex` reasons that are not already +/// in [`CORE_REJECTIONS`] (fgbio `UmiConsensusCaller.initializeRejectCounts`). +pub const DUPLEX_SEEDED_REJECTIONS: [RejectionReason; 3] = [ + RejectionReason::NonPairedReads, + RejectionReason::SameStrandOnly, + RejectionReason::PotentialUmiCollision, +]; + impl ConsensusMetrics { /// Creates a new consensus metrics struct with all counts initialized to zero. #[must_use] @@ -290,6 +299,17 @@ impl ConsensusMetrics { /// metrics output format used by `CallMolecularConsensusReads`. #[must_use] pub fn to_kv_metrics(&self) -> Vec { + self.to_kv_metrics_seeded(&[]) + } + + /// Like [`Self::to_kv_metrics`], but additionally always-emits (seeds to + /// zero) the given rejection reasons, matching fgbio's per-caller + /// `initializeRejectCounts` seeding. The duplex caller passes + /// [`DUPLEX_SEEDED_REJECTIONS`] so that `non_paired_reads`, + /// `single_strand_only`, and `potential_umi_collision` rows are always + /// present even when their counts are zero. + #[must_use] + pub fn to_kv_metrics_seeded(&self, seeded: &[RejectionReason]) -> Vec { let mut metrics = Vec::new(); let raw_reads_used = self.total_input_reads.saturating_sub(self.filtered_reads); @@ -326,8 +346,22 @@ impl ConsensusMetrics { )); } - // Optional rejection reasons (only included when non-zero) - self.push_optional_rejections(&mut metrics); + // Caller-seeded rejection reasons (always included, even when zero), + // skipping any already covered by CORE_REJECTIONS. + for reason in seeded { + if CORE_REJECTIONS.contains(reason) { + continue; + } + metrics.push(ConsensusKvMetric::new( + reason.tsv_key(), + self.rejection_count(*reason).to_string(), + reason.kv_description(), + )); + } + + // Optional rejection reasons (only included when non-zero, and not + // already emitted as a seeded reason above). + self.push_optional_rejections(&mut metrics, seeded); // Final output metric metrics.push(ConsensusKvMetric::new( @@ -340,8 +374,15 @@ impl ConsensusMetrics { } /// Appends fgumi-specific rejection metrics with non-zero counts. - fn push_optional_rejections(&self, metrics: &mut Vec) { + fn push_optional_rejections( + &self, + metrics: &mut Vec, + seeded: &[RejectionReason], + ) { for reason in &OPTIONAL_REJECTIONS { + if seeded.contains(reason) { + continue; + } let count = self.rejection_count(*reason); if count > 0 { metrics.push(ConsensusKvMetric::new( @@ -583,6 +624,44 @@ mod tests { assert!(en.is_none()); } + #[test] + fn test_to_kv_metrics_seeded_always_emits_duplex_reasons() { + // With all counts zero, the plain KV output omits the duplex-seeded + // reasons (they are optional/non-zero), but the seeded output always + // emits exactly one row per seeded reason. + let metrics = ConsensusMetrics::new(); + + let plain = metrics.to_kv_metrics(); + for reason in DUPLEX_SEEDED_REJECTIONS { + assert!( + plain.iter().all(|m| m.key != reason.tsv_key()), + "{} should be omitted from the un-seeded output at zero", + reason.tsv_key() + ); + } + + let seeded = metrics.to_kv_metrics_seeded(&DUPLEX_SEEDED_REJECTIONS); + for reason in DUPLEX_SEEDED_REJECTIONS { + let key = reason.tsv_key(); + let rows: Vec<_> = seeded.iter().filter(|m| m.key == key).collect(); + assert_eq!(rows.len(), 1, "{key} should be emitted exactly once even at zero"); + assert_eq!(rows[0].value, "0", "{key} should be seeded to zero"); + } + } + + #[test] + fn test_to_kv_metrics_seeded_reports_nonzero_counts() { + // A non-zero seeded reason is emitted with its count and not duplicated. + let mut metrics = ConsensusMetrics::new(); + metrics.rejected_non_paired_reads = 7; + + let seeded = metrics.to_kv_metrics_seeded(&DUPLEX_SEEDED_REJECTIONS); + let rows: Vec<_> = + seeded.iter().filter(|m| m.key == "raw_reads_rejected_for_non_paired_reads").collect(); + assert_eq!(rows.len(), 1, "non_paired_reads must not be duplicated"); + assert_eq!(rows[0].value, "7"); + } + #[test] fn test_consensus_kv_metric_new() { let kv = ConsensusKvMetric::new("test_key", 42u64.to_string(), "A test metric"); diff --git a/src/lib/commands/duplex.rs b/src/lib/commands/duplex.rs index 2206cd48e..2bf06e71b 100644 --- a/src/lib/commands/duplex.rs +++ b/src/lib/commands/duplex.rs @@ -822,7 +822,11 @@ impl Duplex { // Write statistics file if requested if let Some(stats_path) = &self.stats_opts.stats { - let kv_metrics = metrics.to_kv_metrics(); + // Duplex always-emits the fgbio `usedByDuplex` rejection rows + // (non_paired_reads, single_strand_only, potential_umi_collision), + // seeded to zero, for metric parity (R2-MET-05). + let kv_metrics = + metrics.to_kv_metrics_seeded(&fgumi_metrics::consensus::DUPLEX_SEEDED_REJECTIONS); DelimFile::default() .write_tsv(stats_path, kv_metrics) .with_context(|| format!("Failed to write statistics: {}", stats_path.display()))?; From 104d226c3304b04e051a83709620e4f6ec2c86bb Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Wed, 8 Jul 2026 23:07:07 -0700 Subject: [PATCH 4/4] fix(duplex): make metrics output thread-independent (single/multi parity) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The duplex statistics output depended on `--threads`: - Format flip: the single-threaded path wrote the wide `ConsensusMetrics` table (`write_tsv([metrics])`) while the threaded path wrote fgbio's KV vertical format (`to_kv_metrics_seeded`). Single-threaded duplex now writes the same seeded KV format, so the seeded `usedByDuplex` rejection rows (R2-MET-05) are emitted regardless of thread count. - `consensus_reads_emitted` mismatch: the threaded path reported this from the caller stat (which counts consensus events — one per molecule) instead of the number of reads actually written (R1+R2 per molecule), so a run emitting 6 consensus reads reported 3. The threaded path now sets `consensus_reads` from the pipeline's written-record count, matching the single-threaded path and the physical output. Adds single-vs-multi-threaded stats-equality assertions to the simplex, duplex, and codec threading-parity tests to lock the invariant that metrics output must not depend on thread count. simplex (KV) and codec (wide) were already self-consistent; duplex is fixed here. Note: codec still writes the wide (non-KV) format in both paths — self-consistent but not fgbio's KV format. Converting codec metrics to KV is a separate parity item, not addressed here. --- src/lib/commands/codec.rs | 11 ++++++++++ src/lib/commands/duplex.rs | 40 +++++++++++++++++++++++++++++++------ src/lib/commands/simplex.rs | 13 +++++++++++- 3 files changed, 57 insertions(+), 7 deletions(-) diff --git a/src/lib/commands/codec.rs b/src/lib/commands/codec.rs index 8887143b0..5c9a8f060 100644 --- a/src/lib/commands/codec.rs +++ b/src/lib/commands/codec.rs @@ -1306,20 +1306,31 @@ mod tests { // Run single-threaded let out_st = dir.path().join("out_st.bam"); + let stats_st = dir.path().join("stats_st.txt"); let mut cmd_st = create_codec_with_paths(input_path.clone(), out_st.clone()); cmd_st.outer_bases_length = 0; cmd_st.threading = ThreadingOptions::none(); + cmd_st.stats_opts.stats = Some(stats_st.clone()); cmd_st.execute("test")?; let records_st = read_bam_records(&out_st)?; // Run multi-threaded let out_mt = dir.path().join("out_mt.bam"); + let stats_mt = dir.path().join("stats_mt.txt"); let mut cmd_mt = create_codec_with_paths(input_path, out_mt.clone()); cmd_mt.outer_bases_length = 0; cmd_mt.threading = ThreadingOptions::new(2); + cmd_mt.stats_opts.stats = Some(stats_mt.clone()); cmd_mt.execute("test")?; let records_mt = read_bam_records(&out_mt)?; + // The metrics output must not depend on thread count. + assert_eq!( + std::fs::read_to_string(&stats_st)?, + std::fs::read_to_string(&stats_mt)?, + "single-threaded and multi-threaded codec must write identical stats" + ); + assert_eq!( records_st.len(), records_mt.len(), diff --git a/src/lib/commands/duplex.rs b/src/lib/commands/duplex.rs index 2bf06e71b..2251c84fa 100644 --- a/src/lib/commands/duplex.rs +++ b/src/lib/commands/duplex.rs @@ -530,11 +530,15 @@ impl Command for Duplex { metrics.consensus_reads = consensus_count as u64; log_consensus_summary(&metrics); - // Write statistics file if requested + // Write statistics file if requested. Emit the same fgbio KV format as + // the threaded path (execute_threads_mode), seeding the `usedByDuplex` + // rejection rows — the metrics output must not depend on thread count. if let Some(stats_path) = &self.stats_opts.stats { - DelimFile::default().write_tsv(stats_path, [metrics]).map_err(|e| { - anyhow::anyhow!("Failed to write statistics to {}: {}", stats_path.display(), e) - })?; + let kv_metrics = + metrics.to_kv_metrics_seeded(&fgumi_metrics::consensus::DUPLEX_SEEDED_REJECTIONS); + DelimFile::default() + .write_tsv(stats_path, kv_metrics) + .with_context(|| format!("Failed to write statistics: {}", stats_path.display()))?; info!("Statistics written to: {}", stats_path.display()); } @@ -816,7 +820,13 @@ impl Duplex { info!("Total MI groups processed: {total_groups}"); info!("Total consensus reads written by pipeline: {consensus_reads_written}"); - let metrics = merged_stats.to_metrics(); + // Report the actual number of consensus reads written by the pipeline + // (R1+R2 per molecule), matching the single-threaded path which sets + // `metrics.consensus_reads` from its streaming write count. The caller + // stat counts consensus events, not emitted reads, so the metric must + // not depend on thread count. + let mut metrics = merged_stats.to_metrics(); + metrics.consensus_reads = consensus_reads_written; let consensus_count = metrics.consensus_reads; log_consensus_summary(&metrics); @@ -1574,8 +1584,10 @@ mod tests { // Test with 1 thread let input1 = create_test_bam(records.clone())?; let paths = TestPaths::new()?; + let stats1_path = paths.output.with_extension("stats.txt"); - let cmd1 = create_duplex_with_paths(input1.path().to_path_buf(), paths.output.clone()); + let mut cmd1 = create_duplex_with_paths(input1.path().to_path_buf(), paths.output.clone()); + cmd1.stats_opts.stats = Some(stats1_path.clone()); cmd1.execute("test")?; let output1_records = read_bam_records(&paths.output)?; @@ -1583,13 +1595,29 @@ mod tests { // Test with 4 threads let input2 = create_test_bam(records)?; let output2_path = paths.output_n(2); + let stats2_path = output2_path.with_extension("stats.txt"); let mut cmd2 = create_duplex_with_paths(input2.path().to_path_buf(), output2_path.clone()); cmd2.threading = ThreadingOptions::new(4); + cmd2.stats_opts.stats = Some(stats2_path.clone()); cmd2.execute("test")?; let output2_records = read_bam_records(&output2_path)?; + // The metrics output must not depend on thread count: the single- and + // multi-threaded paths must write byte-identical stats (same fgbio KV + // format, same seeded rejection rows, same counts). + let stats1 = std::fs::read_to_string(&stats1_path)?; + let stats2 = std::fs::read_to_string(&stats2_path)?; + assert_eq!( + stats1, stats2, + "single-threaded and multi-threaded duplex must write identical stats" + ); + assert!( + stats1.contains("raw_reads_rejected_for_non_paired_reads"), + "duplex stats must always emit the seeded non_paired_reads row (KV format)" + ); + // Should produce the same number of consensus reads assert_eq!( output1_records.len(), diff --git a/src/lib/commands/simplex.rs b/src/lib/commands/simplex.rs index 30c1b59d2..5ef486a07 100644 --- a/src/lib/commands/simplex.rs +++ b/src/lib/commands/simplex.rs @@ -1165,16 +1165,20 @@ mod tests { let paths = TestPaths::new()?; let output_multi_path = paths.dir.path().join("output_multi.bam"); + let stats_single_path = paths.dir.path().join("stats_single.txt"); + let stats_multi_path = paths.dir.path().join("stats_multi.txt"); builder.write(&paths.input)?; // Run with single thread - let cmd_single = create_simplex_with_paths(paths.input.clone(), paths.output.clone()); + let mut cmd_single = create_simplex_with_paths(paths.input.clone(), paths.output.clone()); + cmd_single.stats_opts.stats = Some(stats_single_path.clone()); cmd_single.execute("test")?; // Run with multiple threads let mut cmd_multi = create_simplex_with_paths(paths.input.clone(), output_multi_path.clone()); cmd_multi.threading = ThreadingOptions::new(4); + cmd_multi.stats_opts.stats = Some(stats_multi_path.clone()); cmd_multi.execute("test")?; // Verify both outputs have same number of records @@ -1188,6 +1192,13 @@ mod tests { ); assert_eq!(records_single.len(), 100, "Should have 100 consensus reads"); + // The metrics output must not depend on thread count. + assert_eq!( + std::fs::read_to_string(&stats_single_path)?, + std::fs::read_to_string(&stats_multi_path)?, + "single-threaded and multi-threaded simplex must write identical stats" + ); + Ok(()) }