Skip to content
Merged
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
44 changes: 26 additions & 18 deletions crates/fgumi-metrics/src/dedup.rs
Original file line number Diff line number Diff line change
Expand Up @@ -30,24 +30,28 @@ use crate::library_size::estimate_library_size;
/// dropped at write time, after they have been counted here.
/// 3. **Diagnostics** — `secondary_reads`, `supplementary_reads`,
/// `missing_tc_tag`.
/// 4. **Pair/orphan breakdown** — the Picard `DuplicationMetrics` split. Because
/// `dedup` dedups each mate in its own position group, reads are classified by
/// their SAM flags: `mapped_pairs`/`duplicate_pairs` count mapped read *pairs*
/// 4. **Pair/orphan breakdown** — the Picard `DuplicationMetrics` split. A read
/// pair is one *template* but two *reads*; `dedup` classifies each primary read
/// by its SAM flags: `mapped_pairs`/`duplicate_pairs` count mapped read *pairs*
/// (both mates mapped; Picard `READ_PAIRS_EXAMINED`/`READ_PAIR_DUPLICATES`,
/// already halved from the two mates to whole pairs), `mapped_orphans`/
/// halved from the two mates to whole pairs), `mapped_orphans`/
/// `duplicate_orphans` count mapped reads with no mapped mate (Picard
/// `UNPAIRED_READS_EXAMINED`/`UNPAIRED_READ_DUPLICATES`), and `unmapped_pairs`,
/// `unmapped_orphans`, `unmated_templates` are read-unit diagnostics. Every
/// counted primary read falls in exactly one mapping bucket, so (ignoring the
/// at-most-one truncated half per library from integer halving)
/// `2*mapped_pairs + mapped_orphans + unmapped_pairs + unmapped_orphans ==
/// total_templates` and `2*duplicate_pairs + duplicate_orphans ==
/// duplicate_templates` — but only for templates holding a single primary
/// read, which is every deduplicated template (each mate is dedup'd in its
/// own position group). The one exception is a `--include-unmapped`
/// pass-through: a fully unmapped pair is emitted verbatim as one template
/// holding both mates, so it adds 1 to `total_templates` but 2 to
/// `unmapped_pairs`, and the left-hand side then exceeds `total_templates`.
/// counted primary read falls in exactly one mapping bucket, so in *read*
/// units (ignoring the at-most-one truncated half per library from integer
/// halving) `2*mapped_pairs + mapped_orphans + unmapped_pairs +
/// unmapped_orphans == total_reads` and `2*duplicate_pairs + duplicate_orphans
/// == duplicate_reads`. This read-unit reconciliation is the robust one; there
/// is no general *template*-unit identity, because a single template can spread
/// its two primary reads across two buckets. A mixed-mapping pair (one mate
/// mapped, one unmapped) adds 1 to `mapped_orphans` and 1 to `unmapped_pairs`
/// yet 1 to `total_templates`, and a `--include-unmapped` pass-through adds 2 to
/// `unmapped_pairs` yet 1 to `total_templates`. The template-unit sum
/// `mapped_pairs + mapped_orphans + unmapped_pairs + unmapped_orphans` therefore
/// equals `total_templates` only for an all-mapped fixture — every surviving
/// template a fully-mapped pair or a mapped single-end read — and exceeds it
/// otherwise.
/// 5. **Complexity** — `percent_duplication` (Picard-parity pair/orphan units) and
/// the Lander-Waterman `estimated_library_size`.
///
Expand Down Expand Up @@ -576,6 +580,8 @@ mod tests {
library: "lib1".to_string(),
total_templates: 100,
duplicate_templates: 30,
total_reads: 100,
duplicate_reads: 30,
mapped_pairs: 35,
duplicate_pairs: 12,
mapped_orphans: 20,
Expand All @@ -592,18 +598,20 @@ mod tests {
assert_eq!(metrics.unmapped_pairs, 6);
assert_eq!(metrics.unmapped_orphans, 4);
assert_eq!(metrics.unmated_templates, 24);
// Documented read-unit invariants: each pair is two counted primary
// reads, so pairs count double against the template total.
// Documented read-unit reconciliation: each pair is two counted primary
// reads, so pairs count double against the read total. This is the robust
// identity (it holds regardless of mixed-mapping templates), so it is
// asserted against the read-unit fields, not the template counts.
assert_eq!(
2 * metrics.mapped_pairs
+ metrics.mapped_orphans
+ metrics.unmapped_pairs
+ metrics.unmapped_orphans,
metrics.total_templates
metrics.total_reads
);
assert_eq!(
2 * metrics.duplicate_pairs + metrics.duplicate_orphans,
metrics.duplicate_templates
metrics.duplicate_reads
);
}

Expand Down
75 changes: 66 additions & 9 deletions src/lib/commands/dedup.rs
Original file line number Diff line number Diff line change
Expand Up @@ -101,8 +101,8 @@ pub struct DedupCounts {
/// Templates emitted untouched by `--include-unmapped` (bypass the filter)
pub passthrough_templates: u64,
/// Mapped primary reads whose mate is also mapped, i.e. one half of a mapped
/// pair. Accumulated in read (pair-half) units because a pair's two mates are
/// dedup'd in separate position groups; halved to `mapped_pairs`
/// pair. Accumulated in read (pair-half) units — a mapped pair's two primary
/// reads are each counted here — and halved to `mapped_pairs`
/// (Picard `READ_PAIRS_EXAMINED`) only when the final row is built.
pub mapped_pair_reads: u64,
/// Duplicate-marked half of a mapped pair; halved to `duplicate_pairs`
Expand Down Expand Up @@ -855,11 +855,11 @@ fn mark_template_as_duplicate(template: &mut Template, dedup_counts: &mut DedupC
/// Classify each primary read of a counted template into the Picard-style
/// pair/orphan breakdown (#804).
///
/// fgumi `dedup` dedups each mate in its own position group, so a template here
/// almost always holds a single primary read; classifying **per primary read**
/// (rather than per template) is what lets a normal read-pair — split across two
/// templates in two groups — be recognised as a pair at all. Each read is bucketed
/// by its own SAM flags, mirroring Picard `DuplicationMetrics`:
/// A template holds a read pair's primary reads (a read pair is one template but
/// two reads). Classifying **per primary read** (rather than per template) is
/// what mirrors Picard `DuplicationMetrics`: a mapped pair contributes two
/// pair-halves, one per mate, later halved to whole pairs. Each read is bucketed
/// by its own SAM flags:
///
/// - mapped with a mapped mate -> `mapped_pair_reads` (a pair half; Picard
/// `READ_PAIRS_EXAMINED` after halving), and `duplicate_pair_reads` if marked.
Expand Down Expand Up @@ -2321,8 +2321,8 @@ mod tests {
// ========================================================================

/// Build a single-primary-read template whose read carries exactly `flags`.
/// `dedup` sees single-mate templates, so classifying one read is the real
/// unit of the pair/orphan pass.
/// The pair/orphan pass classifies per primary read, so a one-read template
/// exercises exactly one bucket in isolation.
fn template_with_primary_flags(flags: u16) -> Template {
let mut b = RawSamBuilder::new();
b.read_name(b"q1").sequence(b"ACGT").qualities(&[30, 30, 30, 30]).flags(flags);
Expand Down Expand Up @@ -2361,6 +2361,63 @@ mod tests {
assert_eq!(breakdown(&counts), expected);
}

/// A mixed-mapping pair (one mate mapped, the other unmapped) is one template
/// but two primary reads that land in *different* buckets: the mapped mate is a
/// `mapped_orphan` (mapped, no mapped mate) and the unmapped mate is an
/// `unmapped_pair`. This is why the pair/orphan breakdown reconciles in **read**
/// units and NOT in template units — the read-unit sum for this template is 2
/// (== its reads), while the template-unit sum is also 2 yet the template counts
/// once. Pins that the documented reconciliation is read-unit, not template-unit
/// (which the `--include-unmapped` pass-through is not the only exception to).
#[test]
fn test_count_pair_orphan_mixed_mapping_pair_reconciles_in_read_units() {
// R1 mapped with its mate unmapped -> mapped_orphan; R2 unmapped, paired ->
// unmapped_pair. Both primary reads of ONE template.
let mut r1 = RawSamBuilder::new();
r1.read_name(b"q1")
.sequence(b"ACGT")
.qualities(&[30, 30, 30, 30])
.flags(flags::PAIRED | flags::FIRST_SEGMENT | flags::MATE_UNMAPPED);
let mut r2 = RawSamBuilder::new();
r2.read_name(b"q1")
.sequence(b"TGCA")
.qualities(&[30, 30, 30, 30])
.flags(flags::PAIRED | flags::LAST_SEGMENT | flags::UNMAPPED);
let template = Template::from_records(vec![r1.build(), r2.build()])
.expect("test template construction should not fail");

let mut counts = DedupCounts::default();
count_template_pair_orphan(&template, &mut counts);
// mapped_pair_reads, mapped_orphans, unmapped_pairs, unmapped_orphans, unmated.
assert_eq!(breakdown(&counts), (0, 1, 1, 0, 0));

// One mixed pair: one surviving template holding two primary reads.
counts.total_templates = 1;
counts.total_reads = 2;
let metrics = to_deduplication_metrics("smpl".to_string(), "lib1".to_string(), &counts);

// Read-unit reconciliation holds: each primary read is counted in exactly
// one bucket, so the read-unit sum equals total_reads.
assert_eq!(
2 * metrics.mapped_pairs
+ metrics.mapped_orphans
+ metrics.unmapped_pairs
+ metrics.unmapped_orphans,
metrics.total_reads,
"the pair/orphan breakdown must reconcile against total_reads",
);
// The template-unit identity does NOT hold: a mixed pair adds 2 to the
// breakdown sum but only 1 to total_templates.
assert_ne!(
metrics.mapped_pairs
+ metrics.mapped_orphans
+ metrics.unmapped_pairs
+ metrics.unmapped_orphans,
metrics.total_templates,
"a mixed-mapping pair breaks the template-unit reconciliation",
);
}

/// A duplicate-marked read increments the duplicate sub-counter of its bucket:
/// pair-half reads feed `duplicate_pair_reads`, orphans feed `duplicate_orphans`.
#[rstest]
Expand Down
Loading
Loading