diff --git a/.coderabbit.yaml b/.coderabbit.yaml index fe94227c9..de5fd8616 100644 --- a/.coderabbit.yaml +++ b/.coderabbit.yaml @@ -329,7 +329,7 @@ reviews: require it to be called out. Do NOT file wording-only nits on help text or doc comments unless the text states a behavior the code does not implement. - - path: "src/lib/{grouper,template,mi_group,read_info}.rs" + - path: "src/lib/{grouper,template,template_filter,mi_group,read_info}.rs" instructions: >- Core read-grouping and template-assembly logic that sits outside the crates but feeds the same fgbio-validated output as `fgumi-umi` and diff --git a/CHANGELOG.md b/CHANGELOG.md index ed245d0c5..affd5e53b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,16 @@ All notable changes to this project will be documented in this file. ### Changed +- **Breaking:** `dedup --metrics` gains eleven columns: one `filtered_*` column per template + filter reason, `filtered_templates` for the total, and `passthrough_templates` for the + `--include-unmapped` templates that bypass the filter. Consumers that parse the file + positionally must be updated; consumers keyed on column names are unaffected. Note that + `fgumi compare metrics` treats a column-set difference as `DIFFER`, so comparing a + pre-change dedup metrics file against a post-change one will fail by design. + + These columns count *templates*, whereas the `discarded_*` columns of `group`'s + fgbio-shaped `.grouping_metrics.txt` count *primary records*. `group`'s output is unchanged. + - **Breaking:** consensus downsampling now selects which reads to retain by a hash of the read name rather than by shuffling a seeded random number generator. This affects `simplex --max-reads`, `duplex --max-reads-per-strand`, and `codec --max-reads`. @@ -28,6 +38,14 @@ All notable changes to this project will be documented in this file. ### Bug Fixes +- `dedup --metrics` now reports the templates dropped by its template filter, broken out by + reason. `dedup` is a read filter as well as a duplicate marker, but the drop counts were + collected, merged across workers, and then discarded at the serialization boundary: a run + that filtered out every input template wrote a row of zeros and logged only "Deduplication + complete: 0 templates", with nothing to indicate records had been dropped or why. A + per-reason summary is now also logged, so the drops are visible without `--metrics` + ([#739](https://github.com/fulcrumgenomics/fgumi/issues/739)). + - `duplex --max-reads-per-strand 0` is now rejected at startup. It previously exited 0 having written an empty BAM, because a cap of zero empties every strand. `simplex` and `codec` already validated their equivalents. diff --git a/crates/fgumi-metrics/src/dedup.rs b/crates/fgumi-metrics/src/dedup.rs new file mode 100644 index 000000000..351d23a46 --- /dev/null +++ b/crates/fgumi-metrics/src/dedup.rs @@ -0,0 +1,164 @@ +//! Metrics for the `dedup` command. + +use serde::{Deserialize, Serialize}; + +use crate::Metric; + +/// Metrics written by `fgumi dedup --metrics`. +/// +/// Three sections, in column order: +/// +/// 1. **Filtering** — `filtered_templates` through `passthrough_templates`. +/// `dedup` is a read *filter* as well as a duplicate marker: templates failing +/// the pre-grouping filter are dropped and never reach the output. The +/// `filtered_*` columns give the count per reason in **template** units — note +/// that `group`'s fgbio-shaped `.grouping_metrics.txt` counts primary +/// *records*, so the two files' discard counts are not comparable. +/// 2. **Deduplication** — `total_templates` through `duplicate_reads`, covering +/// only what survived filtering. These count what *passed the filter*, not what +/// reaches the output file: under `--remove-duplicates` the duplicates are +/// dropped at write time, after they have been counted here. +/// 3. **Diagnostics** — `secondary_reads`, `supplementary_reads`, +/// `missing_tc_tag`. +/// +/// Templates read from the input are `filtered_templates + total_templates`; that +/// sum is not emitted as its own column because it is derivable, matching +/// `UmiGroupingMetrics`' treatment of its derivable fields. +// NB: plain code spans above, never intra-doc links — this doc is rendered verbatim +// into docs/src/metrics/deduplication-metrics.md, where a rustdoc link would come +// out as literal markdown. Keep field docs to a single short line: they become the +// Description column of that page's table. +#[derive(Debug, Clone, Default, Serialize, Deserialize)] +pub struct DeduplicationMetrics { + /// Templates dropped by the filter (any reason) + pub filtered_templates: u64, + /// Templates dropped (record shorter than the minimum BAM record length) + pub filtered_malformed_record: u64, + /// Templates dropped (no primary R1 or R2) + pub filtered_no_primary_reads: u64, + /// Templates dropped (no mapped primary read) + pub filtered_unmapped: u64, + /// Templates dropped (QC-fail flag) + pub filtered_not_passing_filter: u64, + /// Templates dropped (mapping quality below `--min-map-q`) + pub filtered_low_mapping_quality: u64, + /// Templates dropped (mate `MQ` below `--min-map-q`) + pub filtered_low_mate_mapping_quality: u64, + /// Templates dropped (no UMI tag) + pub filtered_missing_umi: u64, + /// Templates dropped (N base in the UMI) + pub filtered_ns_in_umi: u64, + /// Templates dropped (UMI shorter than `--min-umi-length`) + pub filtered_umi_too_short: u64, + /// Templates emitted untouched by `--include-unmapped` + pub passthrough_templates: u64, + /// Templates that passed the filter + pub total_templates: u64, + /// Templates kept as their family's representative + pub unique_templates: u64, + /// Templates marked as duplicates + pub duplicate_templates: u64, + /// `duplicate_templates` / `total_templates` + #[serde(with = "crate::float")] + pub duplicate_rate: f64, + /// Reads in templates that passed the filter + pub total_reads: u64, + /// Reads not marked as duplicates + pub unique_reads: u64, + /// Reads marked as duplicates + pub duplicate_reads: u64, + /// Secondary alignments that passed the filter + pub secondary_reads: u64, + /// Supplementary alignments that passed the filter + pub supplementary_reads: u64, + /// Secondary/supplementary reads lacking the `tc` tag + pub missing_tc_tag: u64, +} + +impl DeduplicationMetrics { + /// Creates a metrics struct with all counts initialized to zero. + #[must_use] + pub fn new() -> Self { + Self::default() + } +} + +impl Metric for DeduplicationMetrics { + fn metric_name() -> &'static str { + "deduplication" + } +} + +impl crate::ProcessingMetrics for DeduplicationMetrics { + fn total_input(&self) -> u64 { + self.filtered_templates + self.total_templates + } + + fn total_output(&self) -> u64 { + self.total_templates + } + + fn total_filtered(&self) -> u64 { + self.filtered_templates + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::template_filter::TemplateFilterReason; + use fgoxide::io::DelimFile; + use tempfile::NamedTempFile; + + fn header_of(metrics: DeduplicationMetrics) -> String { + let file = NamedTempFile::new().expect("temp file"); + DelimFile::default().write_tsv(file.path(), [metrics]).expect("write"); + std::fs::read_to_string(file.path()) + .expect("read") + .lines() + .next() + .expect("header") + .to_string() + } + + /// Every filter reason must have a column. This is the test that would have + /// caught fgumi#739: the counts existed and were merged correctly, then were + /// dropped at the serialization boundary because no column held them. + #[test] + fn every_filter_reason_has_a_column() { + let header = header_of(DeduplicationMetrics::default()); + let columns: Vec<&str> = header.split('\t').collect(); + for reason in TemplateFilterReason::ALL { + let column = format!("filtered_{}", reason.column_suffix()); + assert!(columns.contains(&column.as_str()), "missing {column} for {reason:?}"); + } + } + + /// Column order is the file's contract; pin it whole so a field reorder is a + /// deliberate, reviewed change rather than a silent schema break. + #[test] + fn serializes_expected_columns_in_order() { + assert_eq!( + header_of(DeduplicationMetrics::default()), + "filtered_templates\tfiltered_malformed_record\tfiltered_no_primary_reads\t\ + filtered_unmapped\tfiltered_not_passing_filter\tfiltered_low_mapping_quality\t\ + filtered_low_mate_mapping_quality\tfiltered_missing_umi\tfiltered_ns_in_umi\t\ + filtered_umi_too_short\tpassthrough_templates\ttotal_templates\tunique_templates\t\ + duplicate_templates\tduplicate_rate\ttotal_reads\tunique_reads\tduplicate_reads\t\ + secondary_reads\tsupplementary_reads\tmissing_tc_tag" + ); + } + + #[test] + fn processing_metrics_reconciles_input_from_the_columns() { + use crate::ProcessingMetrics; + let metrics = DeduplicationMetrics { + filtered_templates: 10, + total_templates: 90, + ..DeduplicationMetrics::default() + }; + assert_eq!(metrics.total_input(), 100); + assert_eq!(metrics.total_output(), 90); + assert_eq!(metrics.total_filtered(), 10); + } +} diff --git a/crates/fgumi-metrics/src/group.rs b/crates/fgumi-metrics/src/group.rs index 8179977d9..88bb42921 100644 --- a/crates/fgumi-metrics/src/group.rs +++ b/crates/fgumi-metrics/src/group.rs @@ -6,6 +6,7 @@ use serde::{Deserialize, Serialize}; use crate::Metric; +use crate::template_filter::{TemplateFilterCounts, TemplateFilterReason}; /// Build a size distribution from (size, count) pairs. /// @@ -107,6 +108,45 @@ impl UmiGroupingMetrics { pub fn new() -> Self { Self::default() } + + /// Builds the fgbio filter columns from [`TemplateFilterCounts`]. + /// + /// fgbio's `UmiGroupingMetric` counts **primary records**, not templates, so + /// this reads the primary-read side of `counts`. The exhaustive match makes + /// adding a new [`TemplateFilterReason`] a compile error until it is assigned + /// a column, so a rejection can never be silently dropped from the output. + /// + /// Six reasons share `discarded_poor_alignment` because fgbio's schema has no + /// finer bucket for them. That collapse is lossy and is retained only for + /// fgbio parity; `dedup` reports each reason separately. + /// + /// Leaves the family-size and summary fields at their defaults; the caller + /// fills those. + #[must_use] + pub fn from_filter_counts(counts: &TemplateFilterCounts) -> Self { + let mut metrics = Self { + total_records: counts.total_primary_reads(), + accepted_records: counts.accepted_primary_reads(), + ..Self::default() + }; + + for reason in TemplateFilterReason::ALL { + let reads = counts.rejected_primary_reads(reason); + let column = match reason { + TemplateFilterReason::NotPassingFilter => &mut metrics.discarded_non_pf, + TemplateFilterReason::NsInUmi => &mut metrics.discarded_ns_in_umi, + TemplateFilterReason::UmiTooShort => &mut metrics.discarded_umi_too_short, + TemplateFilterReason::MalformedRecord + | TemplateFilterReason::NoPrimaryReads + | TemplateFilterReason::Unmapped + | TemplateFilterReason::LowMappingQuality + | TemplateFilterReason::LowMateMappingQuality + | TemplateFilterReason::MissingUmi => &mut metrics.discarded_poor_alignment, + }; + *column += reads; + } + metrics + } } impl Metric for UmiGroupingMetrics { @@ -258,6 +298,60 @@ mod tests { assert_eq!(metrics.unique_molecule_ids, 0); } + /// fgbio columns are in **primary-read** units, and the six reasons fgbio has + /// no bucket for must all land in `discarded_poor_alignment`. This mapping is + /// what keeps group's output byte-identical across the `TemplateFilterCounts` + /// refactor. + #[rstest::rstest] + #[case::malformed(TemplateFilterReason::MalformedRecord, 0, 2, 0, 0)] + #[case::no_primary(TemplateFilterReason::NoPrimaryReads, 0, 2, 0, 0)] + #[case::unmapped(TemplateFilterReason::Unmapped, 0, 2, 0, 0)] + #[case::non_pf(TemplateFilterReason::NotPassingFilter, 2, 0, 0, 0)] + #[case::low_mapq(TemplateFilterReason::LowMappingQuality, 0, 2, 0, 0)] + #[case::low_mate_mapq(TemplateFilterReason::LowMateMappingQuality, 0, 2, 0, 0)] + #[case::missing_umi(TemplateFilterReason::MissingUmi, 0, 2, 0, 0)] + #[case::ns_in_umi(TemplateFilterReason::NsInUmi, 0, 0, 2, 0)] + #[case::umi_too_short(TemplateFilterReason::UmiTooShort, 0, 0, 0, 2)] + fn set_filter_counts_maps_reasons_to_fgbio_columns( + #[case] reason: TemplateFilterReason, + #[case] expected_non_pf: u64, + #[case] expected_poor_alignment: u64, + #[case] expected_ns_in_umi: u64, + #[case] expected_umi_too_short: u64, + ) { + let mut counts = TemplateFilterCounts::new(); + counts.record_accepted(2); + counts.record_rejected(reason, 2); + + let metrics = UmiGroupingMetrics::from_filter_counts(&counts); + + assert_eq!(metrics.accepted_records, 2, "accepted is in primary-read units"); + assert_eq!(metrics.total_records, 4); + assert_eq!(metrics.discarded_non_pf, expected_non_pf); + assert_eq!(metrics.discarded_poor_alignment, expected_poor_alignment); + assert_eq!(metrics.discarded_ns_in_umi, expected_ns_in_umi); + assert_eq!(metrics.discarded_umi_too_short, expected_umi_too_short); + } + + /// No rejection may be dropped: the four fgbio columns must sum to the total. + #[test] + fn set_filter_counts_conserves_every_rejection() { + let mut counts = TemplateFilterCounts::new(); + counts.record_accepted(2); + for reason in TemplateFilterReason::ALL { + counts.record_rejected(reason, 2); + } + + let metrics = UmiGroupingMetrics::from_filter_counts(&counts); + + let discarded = metrics.discarded_non_pf + + metrics.discarded_poor_alignment + + metrics.discarded_ns_in_umi + + metrics.discarded_umi_too_short; + assert_eq!(discarded, counts.total_rejected_primary_reads()); + assert_eq!(metrics.total_records, metrics.accepted_records + discarded); + } + #[test] fn test_umi_grouping_metrics_serializes_fgbio_five_columns() { use fgoxide::io::DelimFile; diff --git a/crates/fgumi-metrics/src/lib.rs b/crates/fgumi-metrics/src/lib.rs index d26520e8d..aaeace398 100644 --- a/crates/fgumi-metrics/src/lib.rs +++ b/crates/fgumi-metrics/src/lib.rs @@ -12,12 +12,14 @@ pub mod clip; pub mod consensus; pub mod correct; +pub mod dedup; pub mod duplex; pub mod float; pub mod group; pub mod rejection; pub mod shared; pub mod simplex; +pub mod template_filter; pub mod writer; use serde::{Deserialize, Serialize}; @@ -98,6 +100,7 @@ pub trait ProcessingMetrics { pub use clip::{ClipCounts, ClippingMetrics, ClippingMetricsCollection, ReadType}; pub use consensus::{ConsensusCallerKind, ConsensusKvMetric, ConsensusMetrics}; pub use correct::UmiCorrectionMetrics; +pub use dedup::DeduplicationMetrics; pub use duplex::{ DuplexFamilySizeMetric, DuplexMetricsCollector, DuplexUmiMetric, DuplexYieldMetric, FamilySizeMetric, @@ -106,6 +109,7 @@ pub use group::{FamilySizeMetrics, PositionGroupSizeMetrics, UmiGroupingMetrics} pub use rejection::{RejectionReason, format_count}; pub use shared::UmiMetric; pub use simplex::{SimplexFamilySizeMetric, SimplexMetricsCollector, SimplexYieldMetric}; +pub use template_filter::{TemplateFilterCounts, TemplateFilterReason}; pub use writer::{read_metrics, read_metrics_auto, write_metrics}; #[cfg(test)] diff --git a/crates/fgumi-metrics/src/template_filter.rs b/crates/fgumi-metrics/src/template_filter.rs new file mode 100644 index 000000000..18ddde326 --- /dev/null +++ b/crates/fgumi-metrics/src/template_filter.rs @@ -0,0 +1,298 @@ +//! Template-level filter accounting shared by the `group` and `dedup` commands. +//! +//! Both commands run the same pre-grouping template filter but report it +//! differently: `group` emits fgbio's five-column `UmiGroupingMetric` in +//! **primary-record** units, `dedup` emits fgumi-native `filtered_*` columns in +//! **template** units. [`TemplateFilterCounts`] is the single measurement those +//! two views render from, so it records every decision in both units and keys +//! rejections by [`TemplateFilterReason`] — one variant per rejection site in the +//! filter, independent of any output file's column vocabulary. + +/// Why the pre-grouping template filter rejected a template. +/// +/// One variant per rejection site, not per output column: fgbio's +/// `UmiGroupingMetric` collapses six of these into a single +/// `discarded_poor_alignment` column, which is a property of that file format +/// rather than of the filter. Adding a rejection site requires adding a variant +/// here, which breaks the exhaustive matches in every renderer until each decides +/// what to do with it. +/// +/// `#[repr(usize)]` is load-bearing: [`TemplateFilterCounts`] indexes its arrays +/// by `reason as usize`, and [`Self::ALL`] is ordered to match. +// `EnumIter` (test-only) auto-generates `iter()` so coverage assertions enumerate +// every variant without a hand-maintained list that can drift. +#[cfg_attr(test, derive(strum::EnumIter))] +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +#[repr(usize)] +pub enum TemplateFilterReason { + /// A record was shorter than the minimum BAM record length (corrupt input). + MalformedRecord, + /// The template had neither an R1 nor an R2 primary record. + NoPrimaryReads, + /// Every primary read in the template was unmapped. + Unmapped, + /// A primary read carried the QC-fail (`0x200`) flag. + NotPassingFilter, + /// A mapped primary read's mapping quality was below the threshold. + LowMappingQuality, + /// A primary read's mate mapping quality (`MQ` tag) was below the threshold. + LowMateMappingQuality, + /// A primary read had no UMI tag. + MissingUmi, + /// The UMI contained at least one `N` base. + NsInUmi, + /// The UMI was shorter than the configured minimum length. + UmiTooShort, +} + +impl TemplateFilterReason { + /// Number of variants. + pub const COUNT: usize = 9; + + /// Every variant, in discriminant order. + pub const ALL: [Self; Self::COUNT] = [ + Self::MalformedRecord, + Self::NoPrimaryReads, + Self::Unmapped, + Self::NotPassingFilter, + Self::LowMappingQuality, + Self::LowMateMappingQuality, + Self::MissingUmi, + Self::NsInUmi, + Self::UmiTooShort, + ]; + + /// Snake-case identifier for this reason, used as the suffix of the + /// `filtered_*` column in `dedup`'s metrics file and in log summaries. + #[must_use] + pub fn column_suffix(self) -> &'static str { + match self { + Self::MalformedRecord => "malformed_record", + Self::NoPrimaryReads => "no_primary_reads", + Self::Unmapped => "unmapped", + Self::NotPassingFilter => "not_passing_filter", + Self::LowMappingQuality => "low_mapping_quality", + Self::LowMateMappingQuality => "low_mate_mapping_quality", + Self::MissingUmi => "missing_umi", + Self::NsInUmi => "ns_in_umi", + Self::UmiTooShort => "umi_too_short", + } + } +} + +/// Per-reason accounting for the pre-grouping template filter. +/// +/// The filter is the only writer: each template that enters is recorded exactly +/// once, via [`Self::record_accepted`] or [`Self::record_rejected`]. That makes +/// `total == accepted + sum(rejected)` true by construction in both units. +/// +/// Counts are tracked in two units simultaneously: +/// - **templates**: one per template, regardless of how many records it holds +/// - **primary reads**: the number of primary (non-secondary, non-supplementary) +/// records in the template, which is fgbio's unit for `UmiGroupingMetric` +/// +/// A template with neither an R1 nor an R2 contributes `1` template and `0` +/// primary reads, so the units are not proportional and neither is derivable +/// from the other. +#[derive(Debug, Clone, Default, PartialEq, Eq)] +pub struct TemplateFilterCounts { + accepted_templates: u64, + accepted_primary_reads: u64, + rejected_templates: [u64; TemplateFilterReason::COUNT], + rejected_primary_reads: [u64; TemplateFilterReason::COUNT], +} + +impl TemplateFilterCounts { + /// Creates an empty set of counts. + #[must_use] + pub fn new() -> Self { + Self::default() + } + + /// Records a template that passed the filter, carrying `primary_reads` + /// primary records. + pub fn record_accepted(&mut self, primary_reads: u64) { + self.accepted_templates += 1; + self.accepted_primary_reads += primary_reads; + } + + /// Records a template rejected for `reason`, carrying `primary_reads` primary + /// records. + pub fn record_rejected(&mut self, reason: TemplateFilterReason, primary_reads: u64) { + self.rejected_templates[reason as usize] += 1; + self.rejected_primary_reads[reason as usize] += primary_reads; + } + + /// Adds every count in `other` into `self`. + pub fn merge(&mut self, other: &Self) { + self.accepted_templates += other.accepted_templates; + self.accepted_primary_reads += other.accepted_primary_reads; + for index in 0..TemplateFilterReason::COUNT { + self.rejected_templates[index] += other.rejected_templates[index]; + self.rejected_primary_reads[index] += other.rejected_primary_reads[index]; + } + } + + /// Templates that passed the filter. + #[must_use] + pub fn accepted_templates(&self) -> u64 { + self.accepted_templates + } + + /// Primary reads in templates that passed the filter. + #[must_use] + pub fn accepted_primary_reads(&self) -> u64 { + self.accepted_primary_reads + } + + /// Templates rejected for `reason`. + #[must_use] + pub fn rejected_templates(&self, reason: TemplateFilterReason) -> u64 { + self.rejected_templates[reason as usize] + } + + /// Primary reads in templates rejected for `reason`. + #[must_use] + pub fn rejected_primary_reads(&self, reason: TemplateFilterReason) -> u64 { + self.rejected_primary_reads[reason as usize] + } + + /// Templates rejected for any reason. + #[must_use] + pub fn total_rejected_templates(&self) -> u64 { + self.rejected_templates.iter().sum() + } + + /// Primary reads in templates rejected for any reason. + #[must_use] + pub fn total_rejected_primary_reads(&self) -> u64 { + self.rejected_primary_reads.iter().sum() + } + + /// Templates that entered the filter (accepted plus rejected). + #[must_use] + pub fn total_templates(&self) -> u64 { + self.accepted_templates + self.total_rejected_templates() + } + + /// Primary reads that entered the filter (accepted plus rejected). + #[must_use] + pub fn total_primary_reads(&self) -> u64 { + self.accepted_primary_reads + self.total_rejected_primary_reads() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use strum::IntoEnumIterator; + + /// `TemplateFilterCounts` indexes its arrays by `reason as usize`, so `ALL` + /// must list every variant in discriminant order or counts land in the wrong + /// bucket. + #[test] + fn all_covers_every_variant_in_discriminant_order() { + assert_eq!(TemplateFilterReason::ALL.len(), TemplateFilterReason::COUNT); + assert_eq!(TemplateFilterReason::iter().count(), TemplateFilterReason::COUNT); + for (position, &reason) in TemplateFilterReason::ALL.iter().enumerate() { + assert_eq!(reason as usize, position, "{reason:?} is out of order in ALL"); + } + } + + /// Column suffixes become metric column names: unique and `snake_case`. + #[test] + fn column_suffixes_are_unique_snake_case() { + let mut seen = std::collections::HashSet::new(); + for reason in TemplateFilterReason::ALL { + let suffix = reason.column_suffix(); + assert!(seen.insert(suffix), "duplicate column suffix {suffix}"); + assert!( + suffix.bytes().all(|b| b.is_ascii_lowercase() || b == b'_'), + "{suffix} is not snake_case" + ); + } + } + + /// The central invariant: every template entering the filter is recorded + /// exactly once, so totals reconcile in both units. + #[rstest::rstest] + #[case::only_accepts(&[], 3, 2)] + #[case::only_rejects(&[(TemplateFilterReason::NsInUmi, 2)], 0, 0)] + #[case::mixed(&[(TemplateFilterReason::NotPassingFilter, 2), (TemplateFilterReason::MalformedRecord, 1)], 4, 2)] + #[case::zero_primary_reads(&[(TemplateFilterReason::NoPrimaryReads, 0)], 1, 2)] + fn totals_reconcile( + #[case] rejects: &[(TemplateFilterReason, u64)], + #[case] accepts: u64, + #[case] primary_reads_per_accept: u64, + ) { + let mut counts = TemplateFilterCounts::new(); + for _ in 0..accepts { + counts.record_accepted(primary_reads_per_accept); + } + for &(reason, primary_reads) in rejects { + counts.record_rejected(reason, primary_reads); + } + + assert_eq!( + counts.total_templates(), + counts.accepted_templates() + counts.total_rejected_templates() + ); + assert_eq!( + counts.total_primary_reads(), + counts.accepted_primary_reads() + counts.total_rejected_primary_reads() + ); + assert_eq!(counts.accepted_templates(), accepts); + assert_eq!(counts.accepted_primary_reads(), accepts * primary_reads_per_accept); + assert_eq!(counts.total_rejected_templates(), rejects.len() as u64); + assert_eq!( + counts.total_rejected_primary_reads(), + rejects.iter().map(|&(_, n)| n).sum::() + ); + } + + /// A rejection lands in its own bucket and nowhere else. + #[rstest::rstest] + #[case::malformed(TemplateFilterReason::MalformedRecord)] + #[case::no_primary(TemplateFilterReason::NoPrimaryReads)] + #[case::unmapped(TemplateFilterReason::Unmapped)] + #[case::non_pf(TemplateFilterReason::NotPassingFilter)] + #[case::low_mapq(TemplateFilterReason::LowMappingQuality)] + #[case::low_mate_mapq(TemplateFilterReason::LowMateMappingQuality)] + #[case::missing_umi(TemplateFilterReason::MissingUmi)] + #[case::ns_in_umi(TemplateFilterReason::NsInUmi)] + #[case::umi_too_short(TemplateFilterReason::UmiTooShort)] + fn rejection_lands_in_exactly_one_bucket(#[case] recorded: TemplateFilterReason) { + let mut counts = TemplateFilterCounts::new(); + counts.record_rejected(recorded, 2); + + for reason in TemplateFilterReason::ALL { + let expected_templates = u64::from(reason == recorded); + let expected_reads = if reason == recorded { 2 } else { 0 }; + assert_eq!(counts.rejected_templates(reason), expected_templates, "{reason:?}"); + assert_eq!(counts.rejected_primary_reads(reason), expected_reads, "{reason:?}"); + } + } + + /// Merge aggregates per-worker counters; it must be additive on every field. + #[test] + fn merge_is_additive() { + let mut left = TemplateFilterCounts::new(); + left.record_accepted(2); + left.record_rejected(TemplateFilterReason::NsInUmi, 2); + + let mut right = TemplateFilterCounts::new(); + right.record_accepted(1); + right.record_rejected(TemplateFilterReason::NsInUmi, 2); + right.record_rejected(TemplateFilterReason::Unmapped, 2); + + left.merge(&right); + + assert_eq!(left.accepted_templates(), 2); + assert_eq!(left.accepted_primary_reads(), 3); + assert_eq!(left.rejected_templates(TemplateFilterReason::NsInUmi), 2); + assert_eq!(left.rejected_primary_reads(TemplateFilterReason::NsInUmi), 4); + assert_eq!(left.rejected_templates(TemplateFilterReason::Unmapped), 1); + assert_eq!(left.total_templates(), 5); + assert_eq!(left.total_primary_reads(), 9); + } +} diff --git a/crates/fgumi-raw-bam/src/lib.rs b/crates/fgumi-raw-bam/src/lib.rs index f4e650312..186a9f7c3 100644 --- a/crates/fgumi-raw-bam/src/lib.rs +++ b/crates/fgumi-raw-bam/src/lib.rs @@ -123,8 +123,8 @@ pub use fgumi_tag::{AsTagBytes, SamTag}; pub use tags::{ ArrayTagRef, AuxStringTags, AuxTagsIter, RawTagsEditor, RawTagsMut, RawTagsView, TagEntry, TagValue, TemplateAuxTags, append_i32_array_tag, append_signed_int_tag, array_tag_element_u16, - array_tag_to_vec_u16, extract_aux_string_tags, extract_template_aux_tags, find_array_tag, - find_float_tag, find_int_tag, find_mc_tag_in_record, find_string_tag, + array_tag_to_vec_u16, extract_aux_string_tags, extract_int_value, extract_template_aux_tags, + find_array_tag, find_float_tag, find_int_tag, find_mc_tag_in_record, find_string_tag, find_string_tag_in_record, find_string_tag_position, find_tag_type, find_uint8_tag, normalize_int_tag_to_smallest_signed, read_tc_template_coordinate, remove_tag, reverse_array_tag_in_place, reverse_complement_string_tag_in_place, diff --git a/crates/fgumi-raw-bam/src/tags.rs b/crates/fgumi-raw-bam/src/tags.rs index e68d73955..3fa6f1a1d 100644 --- a/crates/fgumi-raw-bam/src/tags.rs +++ b/crates/fgumi-raw-bam/src/tags.rs @@ -146,31 +146,50 @@ pub fn find_int_tag(aux_data: &[u8], tag: impl AsTagBytes) -> Option { extract_int_value(aux_data, p, val_type) } +/// Borrow the `width` value bytes of the tag entry that starts at `p`. +/// +/// The value begins at `p + 3`, after the two tag bytes and the type byte. Returns +/// `None` if that range overflows `usize` or runs past the end of `aux_data`. +/// +/// Every offset is computed with checked arithmetic because `p` reaches +/// `extract_int_value` from the caller rather than from a `find_tag_position` scan: +/// a position near `usize::MAX` would otherwise wrap and decode unrelated bytes. +#[inline] +fn tag_value_bytes(aux_data: &[u8], p: usize, width: usize) -> Option<&[u8]> { + let start = p.checked_add(3)?; + let end = start.checked_add(width)?; + aux_data.get(start..end) +} + /// Extract an integer value at position `p` with the given type byte. /// -/// Shared by [`find_int_tag`] and [`find_mi_tag`]. -fn extract_int_value(aux_data: &[u8], p: usize, val_type: u8) -> Option { +/// Shared by [`find_int_tag`] and `find_mi_tag`, and by callers that walk the +/// aux block themselves and already hold the tag's position and type byte — they +/// get the decode ladder without paying for a second `find_tag_position` scan. +/// +/// Returns `None` for a non-integer type byte, and for any `p` whose value bytes +/// do not lie wholly within `aux_data`. +#[must_use] +pub fn extract_int_value(aux_data: &[u8], p: usize, val_type: u8) -> Option { match val_type { - b'c' if p + 4 <= aux_data.len() => Some(i64::from(aux_data[p + 3].cast_signed())), - b'C' if p + 4 <= aux_data.len() => Some(i64::from(aux_data[p + 3])), - b's' if p + 5 <= aux_data.len() => { - Some(i64::from(i16::from_le_bytes([aux_data[p + 3], aux_data[p + 4]]))) + b'c' => Some(i64::from(tag_value_bytes(aux_data, p, 1)?[0].cast_signed())), + b'C' => Some(i64::from(tag_value_bytes(aux_data, p, 1)?[0])), + b's' => { + let bytes = tag_value_bytes(aux_data, p, 2)?; + Some(i64::from(i16::from_le_bytes([bytes[0], bytes[1]]))) } - b'S' if p + 5 <= aux_data.len() => { - Some(i64::from(u16::from_le_bytes([aux_data[p + 3], aux_data[p + 4]]))) + b'S' => { + let bytes = tag_value_bytes(aux_data, p, 2)?; + Some(i64::from(u16::from_le_bytes([bytes[0], bytes[1]]))) + } + b'i' => { + let bytes = tag_value_bytes(aux_data, p, 4)?; + Some(i64::from(i32::from_le_bytes([bytes[0], bytes[1], bytes[2], bytes[3]]))) + } + b'I' => { + let bytes = tag_value_bytes(aux_data, p, 4)?; + Some(i64::from(u32::from_le_bytes([bytes[0], bytes[1], bytes[2], bytes[3]]))) } - b'i' if p + 7 <= aux_data.len() => Some(i64::from(i32::from_le_bytes([ - aux_data[p + 3], - aux_data[p + 4], - aux_data[p + 5], - aux_data[p + 6], - ]))), - b'I' if p + 7 <= aux_data.len() => Some(i64::from(u32::from_le_bytes([ - aux_data[p + 3], - aux_data[p + 4], - aux_data[p + 5], - aux_data[p + 6], - ]))), _ => None, } } @@ -2413,6 +2432,55 @@ mod tests { assert_eq!(find_int_tag(&aux, b"ZZ"), None); } + // ======================================================================== + // extract_int_value bounds tests + // ======================================================================== + + /// `extract_int_value` is public and takes `p` from the caller rather than from a + /// `find_tag_position` scan, so a position near `usize::MAX` reaches the decode + /// ladder. With unchecked `p + 3` / `p + 4` / `p + 7` the guard wraps to a small + /// number, passes the length check, and the ladder then decodes bytes from a + /// wrapped offset — the wrong value, silently, in a release build. + #[rstest] + #[case::signed_byte(b'c')] + #[case::unsigned_byte(b'C')] + #[case::signed_short(b's')] + #[case::unsigned_short(b'S')] + #[case::signed_int(b'i')] + #[case::unsigned_int(b'I')] + fn test_extract_int_value_rejects_positions_that_overflow(#[case] type_byte: u8) { + let aux = [b'X', b'Y', type_byte, 1, 2, 3, 4, 5]; + // Every `p` within a value width of the top of the address space: `p + 3`, + // `p + 4`, and `p + 7` all wrap for at least one of these. + for delta in 0..8usize { + let p = usize::MAX - delta; + assert_eq!( + extract_int_value(&aux, p, type_byte), + None, + "type {} at p = usize::MAX - {delta} must not decode", + type_byte as char + ); + } + } + + /// The in-bounds ladder is unchanged: a value that ends exactly at the end of the + /// aux block still decodes, and one byte short of that still does not. + #[rstest] + #[case::signed_byte(b'c', 1)] + #[case::unsigned_byte(b'C', 1)] + #[case::signed_short(b's', 2)] + #[case::unsigned_short(b'S', 2)] + #[case::signed_int(b'i', 4)] + #[case::unsigned_int(b'I', 4)] + fn test_extract_int_value_boundary(#[case] type_byte: u8, #[case] width: usize) { + let mut exact = vec![b'X', b'Y', type_byte]; + exact.extend(std::iter::repeat_n(0u8, width)); + assert_eq!(extract_int_value(&exact, 0, type_byte), Some(0)); + + let truncated = &exact[..exact.len() - 1]; + assert_eq!(extract_int_value(truncated, 0, type_byte), None); + } + // --- Per-type find_int_tag tests --- #[test] diff --git a/docs/src/guide/working-with-metrics.md b/docs/src/guide/working-with-metrics.md index 2f814979d..92b6da320 100644 --- a/docs/src/guide/working-with-metrics.md +++ b/docs/src/guide/working-with-metrics.md @@ -27,10 +27,18 @@ A header row followed by a single data row. Used by `dedup`, `codec`, `duplex-me `simplex-metrics`, and `group`. ```text -total_templates unique_templates duplicate_templates duplicate_rate -25000 18750 6250 0.25 +filtered_templates filtered_low_mapping_quality … total_templates unique_templates duplicate_templates duplicate_rate +1200 1200 … 25000 18750 6250 0.25 ``` +The example above is abridged; see +[DeduplicationMetrics](../metrics/deduplication-metrics.md) for the full column list. + +Note that `dedup` and `group` count their discards in **different units**. `dedup`'s +`filtered_*` columns count *templates*, while `group`'s `discarded_*` columns count *primary +records* — the latter match fgbio's `UmiGroupingMetric` schema exactly, so they are readable +by fgbio's `Metric.read`. Do not compare the two files' discard counts directly. + ### Vertical Key-Value (Simplex/Duplex) The `simplex` and `duplex` commands use a three-column format with one metric per row: diff --git a/src/lib/commands/dedup.rs b/src/lib/commands/dedup.rs index 24266b763..e1b1f1ccb 100644 --- a/src/lib/commands/dedup.rs +++ b/src/lib/commands/dedup.rs @@ -23,16 +23,19 @@ use std::sync::Arc; use std::sync::atomic::{AtomicU64, Ordering}; use crate::assigner::{PairedUmiAssigner, Strategy, UmiAssigner}; -use crate::grouper::{ - FilterMetrics, RawPositionGroup, RecordPositionGrouper, build_templates_from_records, -}; +use crate::grouper::{RawPositionGroup, RecordPositionGrouper, build_templates_from_records}; use crate::logging::OperationTimer; +use crate::metrics::DeduplicationMetrics; use crate::metrics::group::FamilySizeMetrics; +use crate::metrics::{TemplateFilterCounts, TemplateFilterReason}; use crate::read_info::LibraryIndex; use crate::sam::SamTag; use crate::sam::is_template_coordinate_sorted; use crate::template::Template; -use crate::umi::{UmiValidation, validate_umi}; +use crate::template_filter::{ + TemplateFilterConfig, filter_template, template_has_malformed_record, + template_is_fully_unmapped, +}; use crate::unified_pipeline::{ BatchWeight, GroupKeyConfig, Grouper, MemoryEstimate, run_bam_pipeline_from_reader_with_mi_assign, @@ -48,7 +51,6 @@ use log::info; use noodles::sam::Header; use noodles::sam::alignment::record::data::field::Tag; use parking_lot::Mutex; -use serde::{Deserialize, Serialize}; use crate::commands::command::Command; use crate::commands::common::{ @@ -57,7 +59,7 @@ use crate::commands::common::{ }; use crate::sam::TC_TAG; use fgumi_raw_bam; -use fgumi_raw_bam::RawRecordView; +use fgumi_raw_bam::{RawRecord, RawRecordView}; /// Duplicate flag bit in SAM flags (0x400) const DUPLICATE_FLAG: u16 = 0x400; @@ -66,16 +68,23 @@ const DUPLICATE_FLAG: u16 = 0x400; // Metrics ////////////////////////////////////////////////////////////////////////////// -/// Metrics collected during deduplication. +/// Per-worker deduplication counters, merged into a single total and then +/// rendered into the serializable metrics schema for output. +/// +/// Named `Counts` rather than `Metrics` to keep it distinct from the file +/// schema: this is the accumulator, that is the file. #[derive(Debug, Default, Clone)] -pub struct DedupMetrics { - /// Total templates processed +pub struct DedupCounts { + /// Templates that passed the filter (excludes filtered) + /// + /// Not "written": under `--remove-duplicates` the duplicates are dropped at + /// serialization, after this counter has already taken them. pub total_templates: u64, /// Templates marked as duplicates pub duplicate_templates: u64, /// Templates kept (not duplicates) pub unique_templates: u64, - /// Total reads processed + /// Reads in templates that passed the filter (excludes filtered) pub total_reads: u64, /// Reads marked as duplicates pub duplicate_reads: u64, @@ -87,13 +96,15 @@ pub struct DedupMetrics { pub supplementary_reads: u64, /// Secondary/supplementary without `tc` tag pub missing_tc_tag: u64, - /// Filter metrics from position grouping - pub filter_metrics: FilterMetrics, + /// Per-reason counts from the template filter + pub filter_counts: TemplateFilterCounts, + /// Templates emitted untouched by `--include-unmapped` (bypass the filter) + pub passthrough_templates: u64, } -impl DedupMetrics { - /// Merge another `DedupMetrics` into this one. - pub fn merge(&mut self, other: &DedupMetrics) { +impl DedupCounts { + /// Merge another `DedupCounts` into this one. + pub fn merge(&mut self, other: &DedupCounts) { self.total_templates += other.total_templates; self.duplicate_templates += other.duplicate_templates; self.unique_templates += other.unique_templates; @@ -103,7 +114,8 @@ impl DedupMetrics { self.secondary_reads += other.secondary_reads; self.supplementary_reads += other.supplementary_reads; self.missing_tc_tag += other.missing_tc_tag; - self.filter_metrics.merge(&other.filter_metrics); + self.filter_counts.merge(&other.filter_counts); + self.passthrough_templates += other.passthrough_templates; } /// Calculate duplicate rate. @@ -117,37 +129,33 @@ impl DedupMetrics { } } -/// Serializable version of `DedupMetrics` for file output. -#[derive(Debug, Clone, Serialize, Deserialize)] -struct DedupMetricsOutput { - total_templates: u64, - unique_templates: u64, - duplicate_templates: u64, - // DXM3-04: format the one fraction column like fgbio's DecimalFormat, matching - // the float serialization #504 wired onto the fgumi-metrics-crate structs. - #[serde(with = "fgumi_metrics::float")] - duplicate_rate: f64, - total_reads: u64, - unique_reads: u64, - duplicate_reads: u64, - secondary_reads: u64, - supplementary_reads: u64, - missing_tc_tag: u64, -} - -impl From<&DedupMetrics> for DedupMetricsOutput { - fn from(m: &DedupMetrics) -> Self { +impl From<&DedupCounts> for DeduplicationMetrics { + fn from(counts: &DedupCounts) -> Self { + use TemplateFilterReason as Reason; + let filter = &counts.filter_counts; Self { - total_templates: m.total_templates, - unique_templates: m.unique_templates, - duplicate_templates: m.duplicate_templates, - duplicate_rate: m.duplicate_rate(), - total_reads: m.total_reads, - unique_reads: m.unique_reads, - duplicate_reads: m.duplicate_reads, - secondary_reads: m.secondary_reads, - supplementary_reads: m.supplementary_reads, - missing_tc_tag: m.missing_tc_tag, + filtered_templates: filter.total_rejected_templates(), + filtered_malformed_record: filter.rejected_templates(Reason::MalformedRecord), + filtered_no_primary_reads: filter.rejected_templates(Reason::NoPrimaryReads), + filtered_unmapped: filter.rejected_templates(Reason::Unmapped), + filtered_not_passing_filter: filter.rejected_templates(Reason::NotPassingFilter), + filtered_low_mapping_quality: filter.rejected_templates(Reason::LowMappingQuality), + filtered_low_mate_mapping_quality: filter + .rejected_templates(Reason::LowMateMappingQuality), + filtered_missing_umi: filter.rejected_templates(Reason::MissingUmi), + filtered_ns_in_umi: filter.rejected_templates(Reason::NsInUmi), + filtered_umi_too_short: filter.rejected_templates(Reason::UmiTooShort), + passthrough_templates: counts.passthrough_templates, + total_templates: counts.total_templates, + unique_templates: counts.unique_templates, + duplicate_templates: counts.duplicate_templates, + duplicate_rate: counts.duplicate_rate(), + total_reads: counts.total_reads, + unique_reads: counts.unique_reads, + duplicate_reads: counts.duplicate_reads, + secondary_reads: counts.secondary_reads, + supplementary_reads: counts.supplementary_reads, + missing_tc_tag: counts.missing_tc_tag, } } } @@ -158,9 +166,9 @@ impl From<&DedupMetrics> for DedupMetricsOutput { /// Metrics collected per position group, aggregated after pipeline completion. #[derive(Default, Debug)] -struct CollectedDedupMetrics { +struct CollectedDedupCounts { /// Dedup-specific metrics - dedup_metrics: DedupMetrics, + dedup_counts: DedupCounts, /// Family size counts family_sizes: AHashMap, } @@ -176,7 +184,7 @@ pub struct ProcessedDedupGroup { /// Family size counts for this group. pub family_sizes: AHashMap, /// Dedup metrics for this group. - pub dedup_metrics: DedupMetrics, + pub dedup_counts: DedupCounts, /// Total input records processed (for progress tracking). pub input_record_count: u64, /// Number of distinct numeric `MoleculeId`s assigned in this group. @@ -197,7 +205,6 @@ impl MemoryEstimate for ProcessedDedupGroup { self.templates.iter().map(MemoryEstimate::estimate_heap_size).sum::() + self.templates.capacity() * std::mem::size_of::