diff --git a/src/lib/commands/duplex_metrics.rs b/src/lib/commands/duplex_metrics.rs index 657da6770..8ca6e9f90 100644 --- a/src/lib/commands/duplex_metrics.rs +++ b/src/lib/commands/duplex_metrics.rs @@ -351,10 +351,15 @@ impl DuplexMetrics { ds_groups.clear(); for m in &downsampled { let entry = ds_groups.entry(m.base_umi).or_default(); - if m.is_a_strand { - entry.0 += 1; - } else if m.is_b_strand { + if m.is_b_strand { entry.1 += 1; + } else { + // /A or an unsuffixed MI counts toward the AB strand. fgbio treats a + // single unsuffixed MI as Pair(ab=n, ba=0) — a valid single-strand DS + // family (CollectDuplexSeqMetrics). Without this, unsuffixed families + // get ds_size=0 and are silently dropped by the family-size histogram + // (which iterates 1..=max), undercounting ds_families (DXM3-03). + entry.0 += 1; } if is_full_fraction { entry.2.push(( @@ -803,6 +808,212 @@ mod tests { Ok(()) } + /// DXM3-07: the two strands of a duplex whose mates share an identical + /// (ref, unclipped-5') must co-group into ONE duplex family. Without a strand + /// tie-break in the canonical key they land in separate coordinate groups and are + /// reported as two single-strand (ab=1, ba=0) families instead of one (ab=1, ba=1). + #[test] + fn test_duplex_strands_cogroup_when_mates_share_five_prime() -> Result<()> { + // Palindromic fragment: both 5' ends map to reference position 100. + // AB strand: R1 fwd @100 (5'=100), R2 rev @1 (5'=end=100), MI .../A + // BA strand: R1 rev @1 (5'=100), R2 fwd @100 (5'=100), MI .../B + let mut records = Vec::new(); + let (r1, r2) = build_test_pair("ab", 0, 100, 1, "AAA-TTT", "AAA-TTT/A", true, false); + records.push(r1); + records.push(r2); + let (r1, r2) = build_test_pair("ba", 0, 1, 100, "TTT-AAA", "AAA-TTT/B", false, true); + records.push(r1); + records.push(r2); + + let input = create_test_bam(records)?; + let output_dir = TempDir::new()?; + let output = output_dir.path().join("output"); + + let cmd = DuplexMetrics { + input: input.path().to_path_buf(), + output: output.clone(), + min_ab_reads: 1, + min_ba_reads: 1, + duplex_umi_counts: false, + intervals: None, + description: None, + }; + cmd.execute("test")?; + + let duplex_family_path = format!("{}.duplex_family_sizes.txt", output.display()); + let metrics: Vec = + DelimFile::default().read_tsv(&duplex_family_path)?; + + // Exactly one duplex family, with both strands observed (ab=1, ba=1). + let total: usize = metrics.iter().map(|m| m.count).sum(); + assert_eq!(total, 1, "the two strands must form a single duplex family"); + let family = metrics.iter().find(|m| m.count > 0).expect("one duplex family"); + assert_eq!((family.ab_size, family.ba_size), (1, 1), "one read on each strand"); + Ok(()) + } + + /// DXM3-02: reads from different libraries (RG -> LB) at the same coordinate/strand + /// must form SEPARATE families, not merge into one. fgbio's `ReadInfo` carries the + /// library; fgumi's own group/dedup partition by it too. + #[test] + fn test_families_partitioned_by_library() -> Result<()> { + use bstr::BString; + use noodles::sam::alignment::record::data::field::Tag; + use noodles::sam::alignment::record_buf::data::field::Value as BufValue; + use noodles::sam::header::record::value::Map; + use noodles::sam::header::record::value::map::ReadGroup; + use noodles::sam::header::record::value::map::ReferenceSequence; + use noodles::sam::header::record::value::map::read_group::tag as rg_tag; + + // Header: chr1 + two read groups with distinct libraries. + let rg_a = Map::::builder() + .insert(rg_tag::LIBRARY, String::from("libA")) + .build() + .expect("read group A"); + let rg_b = Map::::builder() + .insert(rg_tag::LIBRARY, String::from("libB")) + .build() + .expect("read group B"); + let header = sam::Header::builder() + .add_reference_sequence( + BString::from("chr1"), + Map::::new( + NonZeroUsize::new(248_956_422).expect("non-zero length"), + ), + ) + .add_read_group(BString::from("rgA"), rg_a) + .add_read_group(BString::from("rgB"), rg_b) + .build(); + + // Two templates at the same coordinate/strand, tagged to different read groups. + let rg_tag_id = Tag::new(b'R', b'G'); + let mk = |name: &str, rg: &str| { + let (mut r1, mut r2) = + build_test_pair(name, 0, 100, 200, "AAA-TTT", "AAA-TTT/A", true, false); + r1.data_mut().insert(rg_tag_id, BufValue::String(rg.into())); + r2.data_mut().insert(rg_tag_id, BufValue::String(rg.into())); + [r1, r2] + }; + let mut records = Vec::new(); + records.extend(mk("q1", "rgA")); + records.extend(mk("q2", "rgB")); + + let temp_file = NamedTempFile::new()?; + let mut writer = bam::io::writer::Builder.build_from_path(temp_file.path())?; + writer.write_header(&header)?; + for record in &records { + writer.write_alignment_record(&header, record)?; + } + drop(writer); + + let output_dir = TempDir::new()?; + let output = output_dir.path().join("output"); + let cmd = DuplexMetrics { + input: temp_file.path().to_path_buf(), + output: output.clone(), + min_ab_reads: 1, + min_ba_reads: 0, + duplex_umi_counts: false, + intervals: None, + description: None, + }; + cmd.execute("test")?; + + let family_path = format!("{}.family_sizes.txt", output.display()); + let family_metrics: Vec = DelimFile::default().read_tsv(&family_path)?; + let cs_families: usize = family_metrics.iter().map(|m| m.cs_count).sum(); + assert_eq!(cs_families, 2, "different libraries must not merge into one CS family"); + let size1 = family_metrics.iter().find(|m| m.family_size == 1).expect("size-1 row present"); + assert_eq!(size1.cs_count, 2, "two CS families, each a single template"); + Ok(()) + } + + /// DXM3-02: reads with different cell barcodes (CB) at the same + /// coordinate/strand/library must form SEPARATE families, not merge into one. + /// Mirrors `test_families_partitioned_by_library` for the `cell_barcode` key field, + /// which is otherwise unexercised (a regression dropping CB from the key would pass). + #[test] + fn test_families_partitioned_by_cell_barcode() -> Result<()> { + use noodles::sam::alignment::record::data::field::Tag; + use noodles::sam::alignment::record_buf::data::field::Value as BufValue; + + // Two templates at the same coordinate/strand and (absent) library, tagged with + // distinct cell barcodes — they must not collapse into one family. + let cb_tag_id = Tag::new(b'C', b'B'); + let mk = |name: &str, cb: &str| { + let (mut r1, mut r2) = + build_test_pair(name, 0, 100, 200, "AAA-TTT", "AAA-TTT/A", true, false); + r1.data_mut().insert(cb_tag_id, BufValue::String(cb.into())); + r2.data_mut().insert(cb_tag_id, BufValue::String(cb.into())); + [r1, r2] + }; + let mut records = Vec::new(); + records.extend(mk("q1", "CELL_A")); + records.extend(mk("q2", "CELL_B")); + + let input = create_test_bam(records)?; + let output_dir = TempDir::new()?; + let output = output_dir.path().join("output"); + let cmd = DuplexMetrics { + input: input.path().to_path_buf(), + output: output.clone(), + min_ab_reads: 1, + min_ba_reads: 0, + duplex_umi_counts: false, + intervals: None, + description: None, + }; + cmd.execute("test")?; + + let family_path = format!("{}.family_sizes.txt", output.display()); + let family_metrics: Vec = DelimFile::default().read_tsv(&family_path)?; + let cs_families: usize = family_metrics.iter().map(|m| m.cs_count).sum(); + assert_eq!(cs_families, 2, "different cell barcodes must not merge into one CS family"); + let size1 = family_metrics.iter().find(|m| m.family_size == 1).expect("size-1 row present"); + assert_eq!(size1.cs_count, 2, "two CS families, each a single template"); + Ok(()) + } + + /// DXM3-03: non-duplex (unsuffixed MI) input must still be counted as a + /// single-strand DS family (`ab=n`, `ba=0`, `ds_size=n`), matching fgbio — not dropped + /// via a size-0 family-size bucket. + #[test] + fn test_unsuffixed_mi_counts_as_single_strand_ds_family() -> Result<()> { + // Two read pairs at the same coordinate with an unsuffixed MI ("1"). + let mut records = Vec::new(); + let (r1, r2) = build_test_pair("q1", 0, 100, 200, "AAA-TTT", "1", true, false); + records.push(r1); + records.push(r2); + let (r1, r2) = build_test_pair("q2", 0, 100, 200, "AAA-TTT", "1", true, false); + records.push(r1); + records.push(r2); + + let input = create_test_bam(records)?; + let output_dir = TempDir::new()?; + let output = output_dir.path().join("output"); + + let cmd = DuplexMetrics { + input: input.path().to_path_buf(), + output: output.clone(), + min_ab_reads: 1, + min_ba_reads: 0, + duplex_umi_counts: false, + intervals: None, + description: None, + }; + cmd.execute("test")?; + + let family_path = format!("{}.family_sizes.txt", output.display()); + let family_metrics: Vec = DelimFile::default().read_tsv(&family_path)?; + + // Exactly one DS family, of size 2 (both reads on the AB strand). + let ds_families: usize = family_metrics.iter().map(|m| m.ds_count).sum(); + assert_eq!(ds_families, 1, "unsuffixed-MI family must be counted, not dropped"); + let size2 = family_metrics.iter().find(|m| m.family_size == 2).expect("size-2 row present"); + assert_eq!(size2.ds_count, 1, "one DS family of size 2"); + Ok(()) + } + #[test] fn test_count_umis_once_per_read_pair() -> Result<()> { let mut records = Vec::new(); diff --git a/src/lib/commands/shared_metrics.rs b/src/lib/commands/shared_metrics.rs index 07db43096..329a8907c 100644 --- a/src/lib/commands/shared_metrics.rs +++ b/src/lib/commands/shared_metrics.rs @@ -5,6 +5,7 @@ //! deterministic downsampling. Both `duplex_metrics` and `simplex_metrics` commands //! build on these shared primitives. +use crate::read_info::LibraryIndex; use crate::sam::SamTag; use crate::template::TemplateIterator; use anyhow::{Context, Result}; @@ -61,7 +62,11 @@ pub struct TemplateInfo { /// Grouping key matching fgbio's `ReadInfo` structure. /// -/// Fields are ordered so the earlier-mapping read comes first. +/// The two mate positions are ordered so the earlier-mapping read comes first; +/// `library` and `cell_barcode` are template-level (identical for both mates) and are +/// part of the key so that reads from different libraries or cells at the same +/// coordinate/strand form separate families, matching fgbio's `ReadInfo` (which carries +/// `library` and `cellBarcode`) and fgumi's own `group`/`dedup` grouping. #[derive(Clone, PartialEq, Eq, Hash)] pub struct ReadInfoKey { /// Reference sequence index for read 1. @@ -76,6 +81,10 @@ pub struct ReadInfoKey { pub start2: i32, /// `true` if read 2 is reverse-complemented. pub strand2: bool, + /// Library index (from the read group's `LB`, via [`LibraryIndex`]); 0 = unknown. + pub library: u16, + /// Cell barcode (`CB` tag) value, or `None` when the tag is absent. + pub cell_barcode: Option>, } /// Pre-computed metadata for a template within a coordinate group. @@ -90,10 +99,14 @@ pub struct TemplateMetadata<'a> { pub is_b_strand: bool, } -/// Computes the unclipped 5' position for a read, matching fgbio's `positionOf`. +/// Computes the unclipped 5' position for a read. /// -/// Delegates to [`unclipped_5prime_from_raw_bam`]. Includes both soft- and hard-clip -/// bases on the 5' side (matching htsjdk / fgbio semantics). Returns `None` for +/// Delegates to [`unclipped_5prime_from_raw_bam`], which subtracts both soft- **and** +/// hard-clip bases on the 5' side. This deliberately matches samtools +/// (`unclipped_start` in `bam.c`, used by `markdup` and template-coordinate sort) and +/// htsjdk `getUnclippedStart`, keeping this grouping key consistent with how fgumi's +/// own `group`/`dedup` compute coordinates. It intentionally diverges from fgbio's +/// `positionOf`/`unSoftClippedStart`, which counts soft clips only. Returns `None` for /// unmapped records or records missing CIGAR ops. fn unclipped_five_prime_position_raw(record: &RawRecord) -> Option { let flags = record.flags(); @@ -468,6 +481,10 @@ where { let (reader, header) = create_raw_bam_reader(input, 1)?; + // Library index (RG -> LB) for partitioning families by library, matching fgbio's + // ReadInfo.library and fgumi's own group/dedup grouping. + let library_index = LibraryIndex::from_header(&header); + let template_iter = TemplateIterator::new(reader); // Streaming approach: process groups as they arrive (assumes consecutive ReadInfo grouping) @@ -549,25 +566,37 @@ where (r1_start, r1_end.unwrap_or(r1_start)) }; - // ReadInfoKey fields are ordered so the earlier-mapping read comes first. - let read_info_key = if (r1_ref, s1) <= (r2_ref, s2) { - ReadInfoKey { - ref_index1: r1_ref, - start1: s1, - strand1: r1_strand, - ref_index2: r2_ref, - start2: s2, - strand2: r2_strand, - } - } else { - ReadInfoKey { - ref_index1: r2_ref, - start1: s2, - strand1: r2_strand, - ref_index2: r1_ref, - start2: s1, - strand2: r1_strand, - } + // Template-level library (RG -> LB) and cell barcode (CB), taken from the primary + // R1, matching fgbio's ReadInfo(library, cellBarcode). Different libraries or cells + // at the same coordinate/strand form separate families (DXM3-02). The cell tag is + // hardcoded to CB, matching fgumi's opinionated group/dedup (no --cell-tag flag). + let library = find_string_tag_in_record(r1.as_ref(), SamTag::RG) + .map_or(0u16, |rg| library_index.get(LibraryIndex::hash_rg(rg))); + let cell_barcode: Option> = + find_string_tag_in_record(r1.as_ref(), SamTag::CB).map(Box::from); + + // Order the two mate positions so the earlier-mapping read comes first. The + // tie-break includes strand (positive sorts before negative, since Rust + // `false < true` and `strand == is_reverse`), matching fgumi's own group/dedup + // canonicalization (`unified_pipeline/base.rs`: `(ref, pos, strand) <= ...`) and + // fgbio's ReadInfo `r1Earlier` (GroupReadsByUmi.scala:105-111). Without the strand + // tie-break, the two strands of a duplex whose mates share an identical + // (ref, unclipped-5') canonicalize to different keys and fail to co-group. + let (ref_index1, start1, strand1, ref_index2, start2, strand2) = + if (r1_ref, s1, r1_strand) <= (r2_ref, s2, r2_strand) { + (r1_ref, s1, r1_strand, r2_ref, s2, r2_strand) + } else { + (r2_ref, s2, r2_strand, r1_ref, s1, r1_strand) + }; + let read_info_key = ReadInfoKey { + ref_index1, + start1, + strand1, + ref_index2, + start2, + strand2, + library, + cell_barcode, }; let hash_fraction = compute_hash_fraction(&read_name); diff --git a/src/lib/commands/simplex_metrics.rs b/src/lib/commands/simplex_metrics.rs index df5869578..82bd33934 100644 --- a/src/lib/commands/simplex_metrics.rs +++ b/src/lib/commands/simplex_metrics.rs @@ -139,7 +139,7 @@ impl Command for SimplexMetrics { &mut collectors, &mut umi_consensus_caller, fraction_counts, - ); + )?; // simplex UMIs are N-component by design, so there is no // wrong-segment-count guard (cf. DXM-03 for duplex). Ok(()) @@ -227,15 +227,34 @@ impl SimplexMetrics { collectors: &mut [SimplexMetricsCollector], umi_consensus_caller: &mut SimpleUmiConsensusCaller, fraction_template_counts: &mut [usize], - ) { + ) -> Result<()> { use std::collections::HashMap; if group.is_empty() { - return; + return Ok(()); } // Pre-compute metadata once for the entire group let metadata = compute_template_metadata(group); + + // SIMM3-01: simplex-metrics assumes single-strand (non-duplex) input. If a base + // UMI carries reads from BOTH the /A and /B strands, the input is duplex data: + // the per-base_umi RX consensus below would mix the two strands' swapped UMI + // orientations and produce garbage counts. Fail loud and point at duplex-metrics. + let mut base_umi_strands: HashMap<&str, (bool, bool)> = HashMap::new(); + for m in &metadata { + let seen = base_umi_strands.entry(m.base_umi).or_default(); + seen.0 |= m.is_a_strand; + seen.1 |= m.is_b_strand; + if seen.0 && seen.1 { + anyhow::bail!( + "simplex-metrics received duplex-UMI data: base UMI '{}' has reads on \ + both the /A and /B strands. Run duplex-metrics for duplex data.", + m.base_umi + ); + } + } + let last_fraction_idx = fractions.len() - 1; let mut ss_groups: HashMap<&str, usize> = HashMap::new(); @@ -301,6 +320,7 @@ impl SimplexMetrics { } } } + Ok(()) } /// Generates a yield metric from a collector at a specific downsampling fraction. @@ -436,6 +456,58 @@ mod tests { (r1, r2) } + /// Like [`build_test_pair`] but with an explicit R1 strand, for constructing realistic + /// duplex geometry (the BA strand has R1 reverse and its mate positions swapped). + fn build_test_pair_stranded( + name: &str, + ref_id: usize, + pos1: i32, + pos2: i32, + rx_umi: &str, + mi_tag: &str, + r1_reverse: bool, + ) -> (sam::alignment::RecordBuf, sam::alignment::RecordBuf) { + let seq = vec![b'A'; 100]; + let quals = vec![30u8; 100]; + let cigar = encode_op(0, 100); // 100M + let (r1_rev, r1_mate_rev) = + if r1_reverse { (flags::REVERSE, 0) } else { (0, flags::MATE_REVERSE) }; + let (r2_rev, r2_mate_rev) = + if r1_reverse { (0, flags::MATE_REVERSE) } else { (flags::REVERSE, 0) }; + + let mut b1 = RawSamBuilder::new(); + b1.read_name(name.as_bytes()) + .flags(flags::PAIRED | flags::FIRST_SEGMENT | r1_rev | r1_mate_rev) + .ref_id(ref_id as i32) + .pos(pos1 - 1) + .mapq(60) + .cigar_ops(&[cigar]) + .sequence(&seq) + .qualities(&quals) + .mate_ref_id(ref_id as i32) + .mate_pos(pos2 - 1); + b1.add_string_tag(SamTag::RX, rx_umi.as_bytes()); + b1.add_string_tag(SamTag::MI, mi_tag.as_bytes()); + let r1 = to_record_buf(b1.build()); + + let mut b2 = RawSamBuilder::new(); + b2.read_name(name.as_bytes()) + .flags(flags::PAIRED | flags::LAST_SEGMENT | r2_rev | r2_mate_rev) + .ref_id(ref_id as i32) + .pos(pos2 - 1) + .mapq(60) + .cigar_ops(&[cigar]) + .sequence(&seq) + .qualities(&quals) + .mate_ref_id(ref_id as i32) + .mate_pos(pos1 - 1); + b2.add_string_tag(SamTag::RX, rx_umi.as_bytes()); + b2.add_string_tag(SamTag::MI, mi_tag.as_bytes()); + let r2 = to_record_buf(b2.build()); + + (r1, r2) + } + fn create_test_bam(records: Vec) -> Result { let temp_file = NamedTempFile::new()?; let header = create_test_header(); @@ -467,6 +539,47 @@ mod tests { Ok(()) } + /// SIMM3-01: simplex-metrics must reject duplex-UMI input (a base UMI observed on + /// both the /A and /B strands at one coordinate), rather than silently producing + /// garbage UMI counts from the mixed-orientation RX consensus. + #[test] + fn test_simplex_metrics_rejects_duplex_input() { + // Real duplex geometry for base UMI "1": the AB strand (R1 forward @100, R2 + // reverse @200) and the BA strand (R1 reverse @200, R2 forward @100, UMI halves + // swapped) canonicalize to the same coordinate/strand key and so co-group. + let mut records = Vec::new(); + let (r1, r2) = build_test_pair_stranded("q1", 0, 100, 200, "AAA-TTT", "1/A", false); + records.push(r1); + records.push(r2); + let (r1, r2) = build_test_pair_stranded("q2", 0, 200, 100, "TTT-AAA", "1/B", true); + records.push(r1); + records.push(r2); + + let input = create_test_bam(records).expect("write test bam"); + let output_dir = TempDir::new().expect("tempdir"); + let output = output_dir.path().join("output"); + let cmd = SimplexMetrics { + input: input.path().to_path_buf(), + output, + min_reads: 1, + intervals: None, + description: None, + }; + let message = + cmd.execute("test").expect_err("duplex-UMI input must be rejected").to_string(); + // Pin the full diagnostic contract, not just the word "duplex": the message must both + // name the offending input ("duplex-UMI data") and give the actionable next step (run + // "duplex-metrics"), so a regression that drops either half is caught. + assert!( + message.contains("duplex-UMI data"), + "error should name duplex-UMI data: {message}" + ); + assert!( + message.contains("duplex-metrics"), + "error should point at duplex-metrics: {message}" + ); + } + #[test] fn test_generate_yield_metric_basic() { let mut collector = SimplexMetricsCollector::new();