Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 8 additions & 6 deletions crates/fgumi-consensus/src/caller.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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,
}
}
}
Expand Down Expand Up @@ -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),
Expand All @@ -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),
];

Expand Down
20 changes: 19 additions & 1 deletion crates/fgumi-consensus/src/duplex_caller.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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<String> = 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| {
Expand All @@ -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.
Expand Down
107 changes: 103 additions & 4 deletions crates/fgumi-metrics/src/consensus.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -121,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.
Expand All @@ -132,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,
Expand All @@ -147,13 +157,24 @@ const OPTIONAL_REJECTIONS: [RejectionReason; 14] = [
RejectionReason::UmiTooShort,
RejectionReason::SameStrandOnly,
RejectionReason::DuplicateUmi,
RejectionReason::NonPairedReads,
RejectionReason::PotentialUmiCollision,
];

/// Returns an iterator over all rejection reasons (core + optional).
fn all_rejections() -> impl Iterator<Item = RejectionReason> {
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]
Expand Down Expand Up @@ -187,6 +208,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,
}
}

Expand All @@ -212,6 +235,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,
}
}

Expand Down Expand Up @@ -242,6 +267,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;
}
}
}

Expand Down Expand Up @@ -270,6 +299,17 @@ impl ConsensusMetrics {
/// metrics output format used by `CallMolecularConsensusReads`.
#[must_use]
pub fn to_kv_metrics(&self) -> Vec<ConsensusKvMetric> {
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<ConsensusKvMetric> {
let mut metrics = Vec::new();

let raw_reads_used = self.total_input_reads.saturating_sub(self.filtered_reads);
Expand Down Expand Up @@ -306,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(
Expand All @@ -320,8 +374,15 @@ impl ConsensusMetrics {
}

/// Appends fgumi-specific rejection metrics with non-zero counts.
fn push_optional_rejections(&self, metrics: &mut Vec<ConsensusKvMetric>) {
fn push_optional_rejections(
&self,
metrics: &mut Vec<ConsensusKvMetric>,
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(
Expand Down Expand Up @@ -563,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");
Expand Down
58 changes: 3 additions & 55 deletions crates/fgumi-metrics/src/correct.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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<S>(value: &f64, serializer: S) -> Result<S::Ok, S::Error>
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<f64, D::Error>
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
Expand Down Expand Up @@ -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,
}

Expand Down
Loading
Loading