From 30494b9e3e303e4e6ef3f3b8a3e60701c0b477af Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Sun, 4 Oct 2026 14:53:01 -0700 Subject: [PATCH 1/3] refactor(clip): share the end-clipping count and fix stale clipper comments clip_start_of_read_raw and clip_end_of_read_raw now use the clipping_at_end_raw helper that overlap clipping already uses. Comments that pointed at the removed typed SamRecordClipper now cite fgbio, and the clip_template_records comment no longer claims the per-read helpers never upgrade clipping (they upgrade the end they clip, as fgbio does). The strand-normalization overlap tests now also run in soft-with-mask mode. --- crates/fgumi-sam/src/clipper.rs | 23 +++++++---------------- src/lib/commands/clip.rs | 6 ++++-- 2 files changed, 11 insertions(+), 18 deletions(-) diff --git a/crates/fgumi-sam/src/clipper.rs b/crates/fgumi-sam/src/clipper.rs index f75c89647..16efa4cb9 100644 --- a/crates/fgumi-sam/src/clipper.rs +++ b/crates/fgumi-sam/src/clipper.rs @@ -755,12 +755,7 @@ impl RawRecordClipper { record: &mut fgumi_raw_bam::RawRecord, clip_length: usize, ) -> usize { - let ops = record.cigar_ops_vec(); - let existing_clipping: usize = ops - .iter() - .take_while(|&&op| matches!(op & 0xF, 4 | 5)) - .map(|&op| (op >> 4) as usize) - .sum(); + let existing_clipping = Self::clipping_at_end_raw(&record.cigar_ops_vec(), true); if clip_length > existing_clipping && !Self::lacks_seq_raw(record) { self.clip_start_of_alignment(record, clip_length - existing_clipping) @@ -776,13 +771,7 @@ impl RawRecordClipper { record: &mut fgumi_raw_bam::RawRecord, clip_length: usize, ) -> usize { - let ops = record.cigar_ops_vec(); - let existing_clipping: usize = ops - .iter() - .rev() - .take_while(|&&op| matches!(op & 0xF, 4 | 5)) - .map(|&op| (op >> 4) as usize) - .sum(); + let existing_clipping = Self::clipping_at_end_raw(&record.cigar_ops_vec(), false); if clip_length > existing_clipping && !Self::lacks_seq_raw(record) { self.clip_end_of_alignment(record, clip_length - existing_clipping) @@ -1002,7 +991,7 @@ impl RawRecordClipper { } // SoftWithMask: mask existing soft-clipped bases at both ends via upgrade_clipping_raw, - // leaving the CIGAR intact (see the typed upgrade_all_clipping for rationale). + // leaving the CIGAR intact, as fgbio `upgradeAllClipping` does via `clip{Start,End}OfRead`. if self.mode == ClippingMode::SoftWithMask { if leading_soft > 0 { self.upgrade_clipping_raw(record, leading_hard + leading_soft, true); @@ -3491,7 +3480,8 @@ mod tests { fn test_clip_overlapping_reads_normalizes_by_strand_typed( #[case] fwd_start: usize, #[case] rev_start: usize, - #[values(ClippingMode::Soft, ClippingMode::Hard)] mode: ClippingMode, + #[values(ClippingMode::Soft, ClippingMode::SoftWithMask, ClippingMode::Hard)] + mode: ClippingMode, ) { let clipper = RawClipperOnBuf::new(mode); let seq = "A".repeat(100); @@ -3530,7 +3520,8 @@ mod tests { fn test_clip_overlapping_reads_normalizes_by_strand_raw( #[case] fwd_start: usize, #[case] rev_start: usize, - #[values(ClippingMode::Soft, ClippingMode::Hard)] mode: ClippingMode, + #[values(ClippingMode::Soft, ClippingMode::SoftWithMask, ClippingMode::Hard)] + mode: ClippingMode, ) { use fgumi_raw_bam::encode_record_buf_to_raw; use noodles::sam::header::record::value::Map; diff --git a/src/lib/commands/clip.rs b/src/lib/commands/clip.rs index fa626892e..0872192c4 100644 --- a/src/lib/commands/clip.rs +++ b/src/lib/commands/clip.rs @@ -301,8 +301,10 @@ impl ClipParams { ) -> Result<(bool, bool)> { // Upgrade existing clipping on *every* read of the template first — including // secondary/supplementary alignments — matching fgbio ClipBam (ClipBam.scala:123) before - // clipping the primary pair. The per-read `clip_pair`/`clip_fragment` helpers deliberately - // do NOT upgrade, so this pre-pass is the sole upgrade site for both threading paths. + // clipping the primary pair. The per-read `clip_pair`/`clip_fragment` helpers do not run + // this whole-read upgrade, so this pre-pass is its sole site for both threading paths. + // (Those helpers still upgrade existing clipping at the specific end they clip, as fgbio's + // `clip{5,3}PrimeEndOfRead` and `clipOverlappingReads` do, regardless of this flag.) if self.upgrade_clipping { for record in records.iter_mut() { clipper.upgrade_all_clipping_raw(record)?; From b2574317ac34d138b780ccc4a94d8b12d623a1d2 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Sun, 4 Oct 2026 15:11:34 -0700 Subject: [PATCH 2/3] test(clip): port the fgbio SamRecordClipper and ClipBam tests fgumi was missing An audit of fgbio's SamRecordClipperTest and ClipBamTest (fgbio e51a661) against fgumi found 30 + 10 cases with no fgumi counterpart and 8 + 11 whose fgumi version dropped clipping modes or assertions (for example checking only the CIGAR, not bases, qualities, return values, mate info or metrics). This ports each of those cases with fgbio's inputs and expected values, in dedicated fgbio_sam_record_clipper_tests and fgbio_clip_bam_tests modules; every test cites its fgbio source line. fgumi already met every expectation, so no production code changes. Fixture-only deviations are documented on the affected tests: ClipBamTest L182's SEQ is resized to match its hard-clipped CIGAR (fgbio's fixture is malformed), and numBasesExtendingPastMate is exercised through the MC-based num_bases_extending_past_mate_raw. --- crates/fgumi-sam/src/clipper.rs | 551 ++++++++++++++++++++++ src/lib/commands/clip.rs | 803 ++++++++++++++++++++++++++++++++ 2 files changed, 1354 insertions(+) diff --git a/crates/fgumi-sam/src/clipper.rs b/crates/fgumi-sam/src/clipper.rs index 16efa4cb9..abdd0acad 100644 --- a/crates/fgumi-sam/src/clipper.rs +++ b/crates/fgumi-sam/src/clipper.rs @@ -1375,6 +1375,30 @@ mod clip_test_adapter { self.on_one(r, |c, raw| c.clip_end_of_read_raw(raw, n)) } + /// Round-trips [`RawRecordClipper::clip_5_prime_end_of_read_raw`] through `RecordBuf`. + pub(super) fn clip_5_prime_end_of_read(&self, r: &mut RecordBuf, n: usize) -> usize { + self.on_one(r, |c, raw| c.clip_5_prime_end_of_read_raw(raw, n)) + } + + /// Round-trips [`RawRecordClipper::clip_3_prime_end_of_read_raw`] through `RecordBuf`. + pub(super) fn clip_3_prime_end_of_read(&self, r: &mut RecordBuf, n: usize) -> usize { + self.on_one(r, |c, raw| c.clip_3_prime_end_of_read_raw(raw, n)) + } + + /// Round-trips [`RawRecordClipper::clip_extending_past_mate_ends`] through `RecordBuf`. + pub(super) fn clip_extending_past_mate_ends( + &self, + r1: &mut RecordBuf, + r2: &mut RecordBuf, + ) -> (usize, usize) { + let mut raw1 = to_raw(r1); + let mut raw2 = to_raw(r2); + let ret = self.inner.clip_extending_past_mate_ends(&mut raw1, &mut raw2); + *r1 = to_buf(&raw1); + *r2 = to_buf(&raw2); + ret + } + pub(super) fn clip_overlapping_reads( &self, r1: &mut RecordBuf, @@ -4797,6 +4821,533 @@ mod tests { nonzero — the past-mate properties risk passing vacuously" ); } + + /// Ports of fgbio `SamRecordClipperTest` (fgbio commit `e51a661`) that fgumi's own clipper + /// tests above either did not cover or covered more weakly (different inputs, missing + /// assertions). Inputs and expected values are copied verbatim from the fgbio tests; each + /// test cites its source as `SamRecordClipperTest.scala:`. + /// + /// fgbio's `r(start, cigar, strand, attrs)` fragment helper is [`fragment`] here, and its + /// `pair(...)` helper is the enclosing module's [`fr_pair`]. Clipping goes through the + /// [`RawClipperOnBuf`] façade, so every case exercises the live [`RawRecordClipper`]. + mod fgbio_sam_record_clipper_tests { + use super::*; + use ClippingMode::{Hard, Soft, SoftWithMask}; + use Strand::{Minus, Plus}; + use noodles::sam::alignment::record::data::field::Tag; + use noodles::sam::alignment::record_buf::data::field::value::Array; + + /// The 50-character `az` attribute fgbio attaches to a 50-base read in its + /// auto-clip-attribute tests. + const AZ_50: &str = "12345678901234567890123456789012345678901234567890"; + + /// fgbio `r(start, cigar, strand)`: a mapped fragment whose read length is the CIGAR's + /// query length. Bases cycle `ACGT` and qualities vary per base, so the tests that compare + /// bases/qualities before and after clipping can tell which bases were kept. + fn fragment(start: usize, cigar: &str, strand: Strand) -> RecordBuf { + let len = cigar_query_len(cigar); + let bases: String = (0..len).map(|i| char::from(b"ACGT"[i % 4])).collect(); + let quals: Vec = + (0..len).map(|i| u8::try_from(10 + i % 40).expect("quality fits in u8")).collect(); + RecordBuilder::mapped_read() + .sequence(&bases) + .qualities(&quals) + .cigar(cigar) + .alignment_start(start) + .reverse_complement(strand.is_reverse()) + .build() + } + + /// fgbio `r(start, cigar, attrs=Map("az" -> az))`: a forward-strand [`fragment`] that + /// carries a string `az` attribute. + fn fragment_with_az(start: usize, cigar: &str, az: &str) -> RecordBuf { + let mut rec = fragment(start, cigar, Plus); + rec.data_mut().insert(tag("az"), Value::from(az)); + rec + } + + /// 1-based alignment start, as fgbio's `rec.start`. + fn alignment_start(rec: &RecordBuf) -> Option { + rec.alignment_start().map(usize::from) + } + + /// 1-based inclusive alignment end, as fgbio's `rec.end`. + fn alignment_end(rec: &RecordBuf) -> Option { + rec.alignment_end().map(usize::from) + } + + /// The record's CIGAR rendered as a SAM string, as fgbio's `rec.cigar.toString`. + fn cigar_string(rec: &RecordBuf) -> String { + format_cigar(&rec.cigar()) + } + + /// The record's bases, as fgbio's `rec.bases`. + fn bases(rec: &RecordBuf) -> Vec { + rec.sequence().as_ref().to_vec() + } + + /// The record's base qualities, as fgbio's `rec.quals`. + fn quals(rec: &RecordBuf) -> Vec { + rec.quality_scores().as_ref().to_vec() + } + + /// The noodles [`Tag`] for a two-character tag name such as `"az"`. The fgbio fixtures + /// use opaque attribute names that have no `SamTag` constant. + fn tag(name: &str) -> Tag { + let &[first, second] = name.as_bytes() else { + panic!("tag name must be two characters: {name:?}") + }; + Tag::new(first, second) + } + + /// Reads a `Z` string tag, panicking if it is absent or of another type. + fn string_tag(rec: &RecordBuf, name: &str) -> String { + match rec.data().get(&tag(name)) { + Some(Value::String(s)) => String::from_utf8(s.to_vec()).expect("UTF-8 tag value"), + other => panic!("tag {name}: expected a string, got {other:?}"), + } + } + + /// Reads a `B:i` (Int32) array tag, panicking if it is absent or of another type. + fn int32_array_tag(rec: &RecordBuf, name: &str) -> Vec { + match rec.data().get(&tag(name)) { + Some(Value::Array(Array::Int32(values))) => values.clone(), + other => panic!("tag {name}: expected a B:i array, got {other:?}"), + } + } + + // ------------------------------------------------------------ clipStartOfAlignment + + /// `SamRecordClipperTest.scala:128` "mask bases and qualities when the clipping mode is + /// `SoftWithMask`". + #[test] + fn clip_start_of_alignment_soft_with_mask_masks_clipped_bases() { + let mut rec = fragment(10, "50M", Plus); + assert_eq!( + RawClipperOnBuf::new(SoftWithMask).clip_start_of_alignment(&mut rec, 10), + 10 + ); + assert_eq!(alignment_start(&rec), Some(20)); + assert_eq!(cigar_string(&rec), "10S40M"); + assert_eq!(bases(&rec)[..10], [NO_CALL_BASE; 10]); + assert_eq!(quals(&rec)[..10], [MIN_PHRED; 10]); + } + + /// `SamRecordClipperTest.scala:137` "mask bases and qualities when the clipping mode is + /// `SoftWithMask`, including existing soft-clips". + #[test] + fn clip_start_of_alignment_soft_with_mask_masks_existing_soft_clips() { + let mut rec = fragment(10, "10S40M", Plus); + assert_eq!( + RawClipperOnBuf::new(SoftWithMask).clip_start_of_alignment(&mut rec, 10), + 10 + ); + assert_eq!(alignment_start(&rec), Some(20)); + assert_eq!(cigar_string(&rec), "20S30M"); + assert_eq!(bases(&rec)[..20], [NO_CALL_BASE; 20]); + assert_eq!(quals(&rec)[..20], [MIN_PHRED; 20]); + } + + /// `SamRecordClipperTest.scala:146` "hard clip 10 bases". + #[test] + fn clip_start_of_alignment_hard_removes_clipped_bases() { + let mut rec = fragment(10, "50M", Plus); + let (prior_bases, prior_quals) = (bases(&rec), quals(&rec)); + assert_eq!(RawClipperOnBuf::new(Hard).clip_start_of_alignment(&mut rec, 10), 10); + assert_eq!(alignment_start(&rec), Some(20)); + assert_eq!(cigar_string(&rec), "10H40M"); + assert_eq!(bases(&rec), prior_bases[10..]); + assert_eq!(quals(&rec), prior_quals[10..]); + } + + /// A 20M fragment carrying fgbio's auto-trim fixture: `A1`/`A2` are per-base (length + /// 20) and so eligible for trimming; `B1`/`B2` (length 10) are not. + fn auto_trim_fragment() -> RecordBuf { + let mut rec = fragment(10, "20M", Plus); + let data = rec.data_mut(); + data.insert(tag("A1"), Value::from("AB".repeat(10))); + data.insert(tag("A2"), Value::from((1..=20).collect::>())); + data.insert(tag("B1"), Value::from("A".repeat(10))); + data.insert(tag("B2"), Value::from((1..=10).collect::>())); + rec + } + + /// `SamRecordClipperTest.scala:175` "correctly handle auto-trimming of attribute with + /// auto=$auto and mode=$mode" (`clipStartOfAlignment`). Per-base attributes are trimmed + /// only when hard clipping with auto-clip enabled; the others are never touched. + #[rstest] + fn clip_start_of_alignment_auto_trims_per_base_attributes( + #[values(Soft, SoftWithMask, Hard)] mode: ClippingMode, + #[values(true, false)] auto: bool, + ) { + let mut rec = auto_trim_fragment(); + let cut = mode == Hard && auto; + let clipper = RawClipperOnBuf::with_auto_clip(mode, auto); + assert_eq!(clipper.clip_start_of_alignment(&mut rec, 5), 5); + let expected_a1 = if cut { "BABABABABABABAB".to_string() } else { "AB".repeat(10) }; + assert_eq!(string_tag(&rec, "A1"), expected_a1); + let a2_from = if cut { 6 } else { 1 }; + assert_eq!(int32_array_tag(&rec, "A2"), (a2_from..=20).collect::>()); + assert_eq!(string_tag(&rec, "B1"), "A".repeat(10)); + assert_eq!(int32_array_tag(&rec, "B2"), (1..=10).collect::>()); + } + + // -------------------------------------------------------------- clipEndOfAlignment + + /// `SamRecordClipperTest.scala:319` "correctly handle auto-trimming of attribute with + /// auto=$auto and mode=$mode" (`clipEndOfAlignment`). + #[rstest] + fn clip_end_of_alignment_auto_trims_per_base_attributes( + #[values(Soft, SoftWithMask, Hard)] mode: ClippingMode, + #[values(true, false)] auto: bool, + ) { + let mut rec = auto_trim_fragment(); + let cut = mode == Hard && auto; + let clipper = RawClipperOnBuf::with_auto_clip(mode, auto); + assert_eq!(clipper.clip_end_of_alignment(&mut rec, 5), 5); + let expected_a1 = if cut { "ABABABABABABABA".to_string() } else { "AB".repeat(10) }; + assert_eq!(string_tag(&rec, "A1"), expected_a1); + let a2_to = if cut { 15 } else { 20 }; + assert_eq!(int32_array_tag(&rec, "A2"), (1..=a2_to).collect::>()); + assert_eq!(string_tag(&rec, "B1"), "A".repeat(10)); + assert_eq!(int32_array_tag(&rec, "B2"), (1..=10).collect::>()); + } + + // ------------------------------------------------- clip{Start,End}OfRead and 5'/3' + + /// `SamRecordClipperTest.scala:340` "expand existing clipping" (`clipStartOfRead`): + /// existing soft or hard clipping counts toward the requested 10 bases. + #[rstest] + #[case::unclipped("50M", 10, "10S40M")] + #[case::partly_soft_clipped("5S45M", 5, "10S40M")] + #[case::already_soft_clipped("20S30M", 0, "20S30M")] + #[case::already_hard_clipped("20H30M", 0, "20H30M")] + fn clip_start_of_read_expands_existing_clipping( + #[case] cigar: &str, + #[case] expected_clipped: usize, + #[case] expected_cigar: &str, + ) { + let mut rec = fragment(10, cigar, Plus); + assert_eq!( + RawClipperOnBuf::new(Soft).clip_start_of_read(&mut rec, 10), + expected_clipped + ); + assert_eq!(cigar_string(&rec), expected_cigar); + } + + /// `SamRecordClipperTest.scala:358` "convert soft-clipping to hard clipping or masked + /// clipping" (`clipStartOfRead`). + #[test] + fn clip_start_of_read_converts_soft_clipping_to_hard_or_masked() { + let mut hard = fragment_with_az(10, "2H8S40M", &AZ_50[2..]); + let mut mask = fragment_with_az(10, "10S40M", AZ_50); + RawClipperOnBuf::with_auto_clip(Hard, true).clip_start_of_read(&mut hard, 5); + RawClipperOnBuf::with_auto_clip(SoftWithMask, true).clip_start_of_read(&mut mask, 5); + + assert_eq!(cigar_string(&hard), "5H5S40M"); + assert_eq!(bases(&hard).len(), 45); + assert_eq!(string_tag(&hard, "az"), "678901234567890123456789012345678901234567890"); + + assert_eq!(cigar_string(&mask), "10S40M"); + assert_eq!(bases(&mask)[..5], *b"NNNNN"); + assert_eq!(string_tag(&mask, "az"), AZ_50); + } + + /// `SamRecordClipperTest.scala:373` "expand existing clipping" (`clipEndOfRead`). + #[rstest] + #[case::unclipped("50M", 10, "40M10S")] + #[case::partly_soft_clipped("45M5S", 5, "40M10S")] + #[case::already_soft_clipped("30M20S", 0, "30M20S")] + #[case::already_hard_clipped("30M20H", 0, "30M20H")] + fn clip_end_of_read_expands_existing_clipping( + #[case] cigar: &str, + #[case] expected_clipped: usize, + #[case] expected_cigar: &str, + ) { + let mut rec = fragment(10, cigar, Plus); + assert_eq!(RawClipperOnBuf::new(Soft).clip_end_of_read(&mut rec, 10), expected_clipped); + assert_eq!(cigar_string(&rec), expected_cigar); + } + + /// `SamRecordClipperTest.scala:391` "convert soft-clipping to hard clipping or masked + /// clipping" (`clipEndOfRead`). + #[test] + fn clip_end_of_read_converts_soft_clipping_to_hard_or_masked() { + let mut hard = fragment_with_az(10, "40M10S", AZ_50); + let mut mask = fragment_with_az(10, "40M10S", AZ_50); + RawClipperOnBuf::with_auto_clip(Hard, true).clip_end_of_read(&mut hard, 5); + RawClipperOnBuf::with_auto_clip(SoftWithMask, true).clip_end_of_read(&mut mask, 5); + + assert_eq!(cigar_string(&hard), "40M5S5H"); + assert_eq!(bases(&hard).len(), 45); + assert_eq!(string_tag(&hard, "az"), "123456789012345678901234567890123456789012345"); + + assert_eq!(cigar_string(&mask), "40M10S"); + assert_eq!(bases(&mask)[45..], *b"NNNNN"); + assert_eq!(string_tag(&mask, "az"), AZ_50); + } + + /// `SamRecordClipperTest.scala:442` "expand existing clipping at the 5' end" + /// (`clip5PrimeEndOfRead`): the 5' end is the leading end on `+` and trailing on `-`. + #[rstest] + #[case::plus_unclipped("50M", Plus, 10, "10S40M")] + #[case::plus_already_clipped("10S40M", Plus, 0, "10S40M")] + #[case::minus_unclipped("50M", Minus, 10, "40M10S")] + #[case::minus_already_clipped("40M10S", Minus, 0, "40M10S")] + fn clip_5_prime_end_of_read_expands_existing_clipping( + #[case] cigar: &str, + #[case] strand: Strand, + #[case] expected_clipped: usize, + #[case] expected_cigar: &str, + ) { + let mut rec = fragment(10, cigar, strand); + let clipped = RawClipperOnBuf::new(Soft).clip_5_prime_end_of_read(&mut rec, 10); + assert_eq!(clipped, expected_clipped); + assert_eq!(cigar_string(&rec), expected_cigar); + } + + /// `SamRecordClipperTest.scala:460` "convert soft-clipping to hard clipping or masked + /// clipping" (`clip5PrimeEndOfRead`). + #[test] + fn clip_5_prime_end_of_read_converts_soft_clipping_to_hard_or_masked() { + let mut hard = fragment_with_az(10, "10S40M", AZ_50); + let mut mask = fragment_with_az(10, "10S40M", AZ_50); + RawClipperOnBuf::with_auto_clip(Hard, true).clip_5_prime_end_of_read(&mut hard, 5); + RawClipperOnBuf::with_auto_clip(SoftWithMask, true) + .clip_5_prime_end_of_read(&mut mask, 5); + + assert_eq!(cigar_string(&hard), "5H5S40M"); + assert_eq!(bases(&hard).len(), 45); + assert_eq!(string_tag(&hard, "az"), "678901234567890123456789012345678901234567890"); + + assert_eq!(cigar_string(&mask), "10S40M"); + assert_eq!(bases(&mask)[..5], *b"NNNNN"); + assert_eq!(string_tag(&mask, "az"), AZ_50); + } + + /// `SamRecordClipperTest.scala:475` "expand existing clipping at the 3' end" + /// (`clip3PrimeEndOfRead`): the 3' end is the trailing end on `+` and leading on `-`. + #[rstest] + #[case::minus_unclipped("50M", Minus, 10, "10S40M")] + #[case::minus_already_clipped("10S40M", Minus, 0, "10S40M")] + #[case::plus_unclipped("50M", Plus, 10, "40M10S")] + #[case::plus_already_clipped("40M10S", Plus, 0, "40M10S")] + fn clip_3_prime_end_of_read_expands_existing_clipping( + #[case] cigar: &str, + #[case] strand: Strand, + #[case] expected_clipped: usize, + #[case] expected_cigar: &str, + ) { + let mut rec = fragment(10, cigar, strand); + let clipped = RawClipperOnBuf::new(Soft).clip_3_prime_end_of_read(&mut rec, 10); + assert_eq!(clipped, expected_clipped); + assert_eq!(cigar_string(&rec), expected_cigar); + } + + /// `SamRecordClipperTest.scala:493` "convert soft-clipping to hard clipping or masked + /// clipping" (`clip3PrimeEndOfRead`). + #[test] + fn clip_3_prime_end_of_read_converts_soft_clipping_to_hard_or_masked() { + let mut hard = fragment_with_az(10, "40M10S", AZ_50); + let mut mask = fragment_with_az(10, "40M10S", AZ_50); + RawClipperOnBuf::with_auto_clip(Hard, true).clip_3_prime_end_of_read(&mut hard, 5); + RawClipperOnBuf::with_auto_clip(SoftWithMask, true) + .clip_3_prime_end_of_read(&mut mask, 5); + + assert_eq!(cigar_string(&hard), "40M5S5H"); + assert_eq!(bases(&hard).len(), 45); + assert_eq!(string_tag(&hard, "az"), "123456789012345678901234567890123456789012345"); + + assert_eq!(cigar_string(&mask), "40M10S"); + assert_eq!(bases(&mask)[45..], *b"NNNNN"); + assert_eq!(string_tag(&mask, "az"), AZ_50); + } + + // -------------------------------------------------------------- upgradeAllClipping + + /// `SamRecordClipperTest.scala:538` "not convert reads that have no soft-clipping": + /// CIGAR, bases and the `az` attribute are all left untouched. + #[rstest] + #[case::no_clipping("55M")] + #[case::hard_clipping_only("5H55M10H")] + fn upgrade_all_clipping_leaves_reads_without_soft_clips_unchanged(#[case] cigar: &str) { + let mut rec = fragment_with_az(10, cigar, AZ_50); + let upgraded = RawClipperOnBuf::new(Hard) + .upgrade_all_clipping(&mut rec) + .expect("upgrade_all_clipping should succeed"); + assert_eq!(upgraded, (0, 0)); + assert_eq!(cigar_string(&rec), cigar); + assert_eq!(bases(&rec).len(), 55); + assert_eq!(string_tag(&rec, "az"), AZ_50); + } + + // ------------------------------------------------------------ clipOverlappingReads + + /// fgbio `clipOverlappingReads` cases, each a `pair(start1, cigar1, Plus, start2, + /// cigar2, Minus)`. Expected read/mate tuples are `(start, end, cigar)`; an end of + /// `None` means the fgbio test does not assert it. + #[rstest] + // SamRecordClipperTest.scala:561 "not clip if the reads are not overlapping" + #[case::abutting_not_overlapping(Soft, (1, "100M"), (101, "100M"), (0, 0), (1, None, "100M"), (101, None, "100M"))] + // SamRecordClipperTest.scala:601 "clip reads that overlap with soft-clipping on the forward read after the midpoint" + #[case::forward_soft_clip_after_midpoint(Hard, (1, "95M5S"), (50, "100M"), (20, 26), (1, Some(75), "75M25H"), (76, None, "26H74M"))] + // SamRecordClipperTest.scala:611 "clip reads that overlap with soft-clipping on the reverse read before the midpoint" + #[case::reverse_soft_clip_before_midpoint(Hard, (1, "100M"), (55, "5S95M"), (25, 21), (1, Some(75), "75M25H"), (76, None, "26H74M"))] + // SamRecordClipperTest.scala:621 "clip reads that overlap 1 bp directly in the middle of the pair" + #[case::one_bp_overlap_at_midpoint(Hard, (1, "99M1S"), (99, "1S99M"), (0, 1), (1, Some(99), "99M1H"), (100, None, "2H98M"))] + // SamRecordClipperTest.scala:633 "clip reads that overlap in the first half of the pair, not over the middle" + #[case::overlap_in_first_half(Hard, (1, "95M5S"), (90, "20S80M"), (6, 0), (1, Some(89), "89M11H"), (90, None, "20H80M"))] + // SamRecordClipperTest.scala:645 "clip reads that overlap in the second half of the pair, not over the middle" + #[case::overlap_in_second_half(Hard, (1, "80M20S"), (70, "5S95M"), (0, 11), (1, Some(80), "80M20H"), (81, None, "16H84M"))] + // SamRecordClipperTest.scala:685 "clip reads that overlap with one end having a deletion with mismatching cigars" + #[case::one_end_deletion_mismatching_cigars(Soft, (1, "100M"), (50, "10M10D80M10D10M"), (15, 26), (1, Some(85), "85M15S"), (86, None, "26S64M10D10M"))] + // SamRecordClipperTest.scala:695 "clip reads that fully overlap with both ends having deletions" + #[case::full_overlap_both_deletions(Soft, (1, "50M10D50M"), (1, "50M10D50M"), (50, 50), (1, Some(50), "50M50S"), (61, None, "50S50M"))] + // SamRecordClipperTest.scala:738 "clip reads that extend past each other with one read having deletions" + #[case::extend_past_one_read_deletions(Soft, (50, "100M"), (1, "10M10D80M10D10M"), (64, 75), (50, Some(85), "36M64S"), (86, Some(120), "75S15M10D10M"))] + // SamRecordClipperTest.scala:750 "clip reads that extend past each other with both read having deletions" + #[case::extend_past_both_reads_deletions(Soft, (50, "50M10D50M"), (1, "10M10D80M10D10M"), (64, 75), (50, Some(85), "36M64S"), (86, Some(120), "75S15M10D10M"))] + fn clip_overlapping_reads_matches_fgbio( + #[case] mode: ClippingMode, + #[case] read: (usize, &str), + #[case] mate: (usize, &str), + #[case] expected_clipped: (usize, usize), + #[case] expected_read: (usize, Option, &str), + #[case] expected_mate: (usize, Option, &str), + ) { + let (mut rec, mut mate_rec) = fr_pair((read.1, read.0, Plus), (mate.1, mate.0, Minus)); + let clipped = + RawClipperOnBuf::new(mode).clip_overlapping_reads(&mut rec, &mut mate_rec); + assert_eq!(clipped, expected_clipped); + for (label, actual, (start, end, cigar)) in + [("read", &rec, expected_read), ("mate", &mate_rec, expected_mate)] + { + assert_eq!(alignment_start(actual), Some(start), "{label} start"); + assert_eq!(cigar_string(actual), cigar, "{label} cigar"); + if let Some(end) = end { + assert_eq!(alignment_end(actual), Some(end), "{label} end"); + } + } + } + + // ------------------------------------------------------- clipExtendingPastMateEnds + + /// `SamRecordClipperTest.scala:771` "not clip reads that do not extend past each other", + /// run for both strand orders as fgbio does. + #[rstest] + #[case::plus_minus(Plus, Minus)] + #[case::minus_plus(Minus, Plus)] + fn clip_extending_past_mate_ends_does_not_clip_coincident_reads( + #[case] read_strand: Strand, + #[case] mate_strand: Strand, + ) { + let (mut rec, mut mate) = fr_pair(("100M", 1, read_strand), ("100M", 1, mate_strand)); + let clipped = + RawClipperOnBuf::new(Soft).clip_extending_past_mate_ends(&mut rec, &mut mate); + assert_eq!(clipped, (0, 0)); + assert_eq!((alignment_start(&rec), cigar_string(&rec)), (Some(1), "100M".to_string())); + assert_eq!( + (alignment_start(&mate), cigar_string(&mate)), + (Some(1), "100M".to_string()) + ); + } + + /// fgbio `clipExtendingPastMateEnds` cases, each a `pair(start1, cigar1, Plus, start2, + /// cigar2, Minus)`. Expected read/mate tuples are `(start, cigar)`. + #[rstest] + // SamRecordClipperTest.scala:782 "clip reads that extend one base past their mate's start" + #[case::extend_one_base(Soft, (2, "100M"), (1, "100M"), (1, 1), (2, "99M1S"), (2, "1S99M"))] + // SamRecordClipperTest.scala:791 "clip reads that extend two bases past their mate's start" + #[case::extend_two_bases(Soft, (3, "100M"), (1, "100M"), (2, 2), (3, "98M2S"), (3, "2S98M"))] + // SamRecordClipperTest.scala:800 "clip reads that where both ends extends their mate's start" + #[case::both_ends_extend(Soft, (51, "100M"), (1, "100M"), (50, 50), (51, "50M50S"), (51, "50S50M"))] + // SamRecordClipperTest.scala:809 "clip reads that where only one end extends their mate's start" + #[case::only_one_end_extends(Soft, (1, "100M"), (1, "50S50M"), (50, 0), (1, "50M50S"), (1, "50S50M"))] + // SamRecordClipperTest.scala:818 "clip reads where only one end extends their mate's start that has insertions" + #[case::only_one_end_extends_with_insertion(Soft, (1, "40M10I50M"), (1, "50S50M"), (40, 0), (1, "40M10I10M40S"), (1, "50S50M"))] + // SamRecordClipperTest.scala:827 "clip the forward read when it ends before the mate's start but soft-clipped bases extend past" + #[case::soft_clip_extends_past(Hard, (20, "30M20S"), (20, "10S40M"), (0, 0), (20, "30M10S10H"), (20, "10H40M"))] + // SamRecordClipperTest.scala:836 "... soft-clipped bases extend past while there is a deletion" + #[case::soft_clip_extends_past_with_deletion(Hard, (20, "15M1D15M20S"), (20, "10S15M1D25M"), (0, 0), (20, "15M1D15M10S10H"), (20, "10H15M1D25M"))] + // SamRecordClipperTest.scala:845 "... soft-clipped bases extend past while there is an insertion" + #[case::soft_clip_extends_past_with_insertion(Hard, (20, "15M1I15M20S"), (20, "10S15M1I25M"), (0, 0), (20, "15M1I15M10S10H"), (20, "10H15M1I25M"))] + // SamRecordClipperTest.scala:854 "clip the reverse read when it ends before the mate's start but soft-clipped bases extend past" + #[case::reverse_soft_clip_extends_past(Hard, (20, "40M10S"), (30, "20S30M"), (0, 0), (20, "40M10H"), (30, "10H10S30M"))] + // SamRecordClipperTest.scala:863 "not clip when the read pairs are mapped +/- with start(R1) > end(R2) but do not overlap" + #[case::disjoint_outward_facing(Soft, (1000, "100M"), (1, "100M"), (0, 0), (1000, "100M"), (1, "100M"))] + // SamRecordClipperTest.scala:872 "not clip when the reads do not extend past each other with insertions" + #[case::coincident_with_insertions(Soft, (1, "40M20I40M"), (1, "40M20I40M"), (0, 0), (1, "40M20I40M"), (1, "40M20I40M"))] + fn clip_extending_past_mate_ends_matches_fgbio( + #[case] mode: ClippingMode, + #[case] read: (usize, &str), + #[case] mate: (usize, &str), + #[case] expected_clipped: (usize, usize), + #[case] expected_read: (usize, &str), + #[case] expected_mate: (usize, &str), + ) { + let (mut rec, mut mate_rec) = fr_pair((read.1, read.0, Plus), (mate.1, mate.0, Minus)); + let clipped = + RawClipperOnBuf::new(mode).clip_extending_past_mate_ends(&mut rec, &mut mate_rec); + assert_eq!(clipped, expected_clipped); + assert_eq!( + (alignment_start(&rec), cigar_string(&rec)), + (Some(expected_read.0), expected_read.1.to_string()), + "read (start, cigar)" + ); + assert_eq!( + (alignment_start(&mate_rec), cigar_string(&mate_rec)), + (Some(expected_mate.0), expected_mate.1.to_string()), + "mate (start, cigar)" + ); + } + + // ------------------------------------------------------- numBasesExtendingPastMate + + /// fgbio's `numBasesExtendingPastMate` reads the mate's unclipped start/end from the + /// `MC` tag; fgumi's equivalent is the MC-based + /// [`fgumi_raw_bam::num_bases_extending_past_mate_raw`]. + fn num_bases_extending_past_mate(rec: &RecordBuf) -> usize { + fgumi_raw_bam::num_bases_extending_past_mate_raw(to_raw(rec).as_ref()) + } + + /// `SamRecordClipperTest.scala:881` "return zero when reads do not extend past the end + /// or are not FR pairs". Each case is `(start, cigar, strand)` for read one and two. + #[rstest] + #[case::fr_r1_more_soft_clipped((100, "20S80M", Plus), (200, "10S90M", Minus))] + #[case::fr_r2_more_soft_clipped((100, "10S90M", Plus), (200, "20S80M", Minus))] + #[case::forward_forward((100, "10S90M", Plus), (100, "20S80M", Plus))] + #[case::reverse_reverse((100, "10S90M", Minus), (100, "20S80M", Minus))] + fn num_bases_extending_past_mate_is_zero_when_not_extending_or_not_fr( + #[case] read_one: (usize, &str, Strand), + #[case] read_two: (usize, &str, Strand), + ) { + let (r1, r2) = + fr_pair((read_one.1, read_one.0, read_one.2), (read_two.1, read_two.0, read_two.2)); + assert_eq!(num_bases_extending_past_mate(&r1), 0, "read one"); + assert_eq!(num_bases_extending_past_mate(&r2), 0, "read two"); + } + + /// `SamRecordClipperTest.scala:904` "return a return a positive value when reads extend + /// past its mate". Each case is an FR pair `(start, cigar)` for read one (`+`) and read + /// two (`-`), with the expected count for each read. + #[rstest] + #[case::both_extend_fifty((100, "100M"), (50, "100M"), (50, 50))] + #[case::same_unclipped_ends((100, "50S50M"), (100, "50S50M"), (0, 0))] + #[case::same_aligned_ends((100, "50M50S"), (100, "50M50S"), (0, 0))] + #[case::read_one_inside_read_two((100, "50S50M"), (100, "30S70M"), (0, 0))] + #[case::read_one_past_read_two((100, "30S70M"), (100, "50S50M"), (20, 20))] + fn num_bases_extending_past_mate_counts_bases_past_the_mate( + #[case] read_one: (usize, &str), + #[case] read_two: (usize, &str), + #[case] expected: (usize, usize), + ) { + let (r1, r2) = fr_pair((read_one.1, read_one.0, Plus), (read_two.1, read_two.0, Minus)); + assert_eq!( + (num_bases_extending_past_mate(&r1), num_bases_extending_past_mate(&r2)), + expected + ); + } + } } // ============================================================================ diff --git a/src/lib/commands/clip.rs b/src/lib/commands/clip.rs index 0872192c4..82a66444f 100644 --- a/src/lib/commands/clip.rs +++ b/src/lib/commands/clip.rs @@ -3505,4 +3505,807 @@ mod tests { "mapped r1 mate-reverse must reflect unmapped r2's actual REVERSE flag" ); } + + /// Ports of fgbio `ClipBamTest` (fgbio commit `e51a661`) that fgumi's own `clip` tests above + /// either did not cover or covered more weakly (different inputs, missing assertions). + /// Inputs and expected values are copied verbatim from the fgbio tests; each test cites its + /// source as `ClipBamTest.scala:`, and any fixture deviation is explained there. + /// + /// fgbio's `clipper.clipPair(r1, r2)` cases run [`ClipParams::clip_pair`] with a `Hard` + /// clipper (fgbio `ClipBam`'s default mode); its `.execute()` cases run [`Clip::execute`] + /// end to end. + mod fgbio_clip_bam_tests { + use super::*; + use crate::metrics::clip::{ClippingMetrics, ReadType}; + use crate::sam::builder::PairBuilder; + use fgumi_raw_bam::{encode_record_buf_to_raw, raw_record_to_record_buf}; + use noodles::sam::alignment::RecordBuf; + use noodles::sam::alignment::record::cigar::Op; + use noodles::sam::alignment::record::cigar::op::Kind; + use noodles::sam::alignment::record::data::field::Tag; + use noodles::sam::alignment::record_buf::data::field::Value; + + /// Length of fgbio's `chr1` test reference, which is all `A`. + const REFERENCE_LENGTH: usize = 5000; + + /// Writes fgbio's test reference (`chr1`: 5000 `A`s) plus its `.fai`; `NM`/`UQ`/`MD` + /// expectations in these tests are computed against it. + fn write_all_a_reference(dir: &TempDir) -> PathBuf { + let path = dir.path().join("ref.fa"); + std::fs::write(&path, format!(">chr1\n{}\n", "A".repeat(REFERENCE_LENGTH))) + .expect("write reference"); + std::fs::write( + dir.path().join("ref.fa.fai"), + format!( + "chr1\t{REFERENCE_LENGTH}\t6\t{REFERENCE_LENGTH}\t{}\n", + REFERENCE_LENGTH + 1 + ), + ) + .expect("write reference index"); + path + } + + /// fgbio `new SamBuilder(readLength).addPair(...)`: builds an FR pair (by default) of + /// all-`A` reads of `read_length` bases, letting `configure` set positions, strands and + /// CIGARs, and encodes both reads as [`RawRecord`]s. + fn fgbio_pair( + read_length: usize, + configure: impl for<'a> FnOnce(PairBuilder<'a>) -> PairBuilder<'a>, + ) -> (RawRecord, RawRecord) { + let mut builder = SamBuilder::new(); + let bases = "A".repeat(read_length); + let pair = builder.add_pair().name("q").bases1(&bases).bases2(&bases); + let (r1, r2) = configure(pair).build(); + let header = builder.header.clone(); + ( + encode_record_buf_to_raw(&r1, &header).expect("encode r1"), + encode_record_buf_to_raw(&r2, &header).expect("encode r2"), + ) + } + + /// fgbio `clipper.clipPair(r1, r2)`: runs the per-template clipping that `clip`'s + /// options configure on one pair, in fgbio `ClipBam`'s default `Hard` mode. + fn clip_pair_hard(clip: &Clip, r1: &mut RawRecord, r2: &mut RawRecord) { + ClipParams::from_clip(clip) + .clip_pair(&RawRecordClipper::new(ClippingMode::Hard), r1, r2, None) + .expect("clip_pair should succeed"); + } + + /// A `Clip` with only `--clip-overlapping-reads` (and no fixed clipping) enabled. + fn overlap_only() -> Clip { + let mut clip = make_clip(0, 0, 0, 0); + clip.clip_overlapping_reads = true; + clip + } + + /// 1-based alignment start of a mapped raw record, as fgbio's `rec.start`. + fn start(rec: &RawRecord) -> usize { + rec.alignment_start_1based().expect("record should be mapped") + } + + /// 1-based inclusive alignment end of a mapped raw record, as fgbio's `rec.end`. + fn end(rec: &RawRecord) -> usize { + rec.alignment_end_1based().expect("record should be mapped") + } + + /// fgbio `StartAndEnd.checkClipping`: asserts `rec` lost exactly `five_prime` aligned + /// bases at its 5' end and `three_prime` at its 3' end relative to `prior` `(start, + /// end)`, accounting for strand. + fn assert_clipped_by( + prior: (usize, usize), + rec: &RawRecord, + five_prime: usize, + three_prime: usize, + label: &str, + ) { + let (prior_start, prior_end) = prior; + let expected = if rec.is_reverse() { + (prior_start + three_prime, prior_end - five_prime) + } else { + (prior_start + five_prime, prior_end - three_prime) + }; + assert_eq!((start(rec), end(rec)), expected, "{label} (start, end)"); + } + + /// `ClipBamTest.scala:86` "not clip reads where either read is unaligned". + #[test] + fn clip_pair_does_not_clip_when_a_read_is_unaligned() { + let (mut r1, mut r2) = fgbio_pair(50, |p| p.start1(100).unmapped2()); + let expected = r1.cigar_to_string(); + clip_pair_hard(&overlap_only(), &mut r1, &mut r2); + assert_eq!(r1.cigar_to_string(), expected); + } + + /// `ClipBamTest.scala:95` "not clip reads that are on different chromosomes". + #[test] + fn clip_pair_does_not_clip_reads_on_different_chromosomes() { + let (mut r1, mut r2) = fgbio_pair(50, |p| p.start1(100).start2(100).contig2(1)); + let expected = (r1.cigar_to_string(), r2.cigar_to_string()); + clip_pair_hard(&overlap_only(), &mut r1, &mut r2); + assert_eq!((r1.cigar_to_string(), r2.cigar_to_string()), expected); + } + + /// `ClipBamTest.scala:108` "not clip reads that are abutting but not overlapped". + #[test] + fn clip_pair_does_not_clip_abutting_reads() { + let (mut r1, mut r2) = fgbio_pair(50, |p| p.start1(100).start2(150)); + let expected = (r1.cigar_to_string(), r2.cigar_to_string()); + clip_pair_hard(&overlap_only(), &mut r1, &mut r2); + assert_eq!((r1.cigar_to_string(), r2.cigar_to_string()), expected); + } + + /// `ClipBamTest.scala:119` "not clip non-FR reads". + #[test] + fn clip_pair_does_not_clip_non_fr_reads() { + let (mut r1, mut r2) = + fgbio_pair(50, |p| p.start1(100).start2(100).strand2(Strand::Plus)); + let expected = (r1.cigar_to_string(), r2.cigar_to_string()); + clip_pair_hard(&overlap_only(), &mut r1, &mut r2); + assert_eq!((r1.cigar_to_string(), r2.cigar_to_string()), expected); + } + + /// `ClipBamTest.scala:130` "clip reads that are fully overlapped". + #[test] + fn clip_pair_hard_clips_fully_overlapped_reads() { + let (mut r1, mut r2) = fgbio_pair(50, |p| p.start1(100).start2(100)); + clip_pair_hard(&overlap_only(), &mut r1, &mut r2); + assert_eq!(r1.cigar_to_string(), "25M25H"); + assert_eq!(r2.cigar_to_string(), "25H25M"); + } + + /// `ClipBamTest.scala:160` "handle reads that contain deletions". + #[test] + fn clip_pair_removes_overlap_of_reads_with_deletions() { + let (mut r1, mut r2) = + fgbio_pair(50, |p| p.start1(100).start2(130).cigar1("40M2D10M").cigar2("10M2D40M")); + assert!(end(&r1) >= start(&r2), "fixture should overlap"); + clip_pair_hard(&overlap_only(), &mut r1, &mut r2); + assert!(end(&r1) < start(&r2), "r1.end={} r2.start={}", end(&r1), start(&r2)); + } + + /// `ClipBamTest.scala:170` "clip a fixed amount on the ends of the reads with reads that + /// do not overlap". + #[test] + fn clip_pair_fixed_clipping_without_overlap() { + let (mut r1, mut r2) = fgbio_pair(50, |p| p.start1(100).start2(150)); + let (prior1, prior2) = ((start(&r1), end(&r1)), (start(&r2), end(&r2))); + assert_eq!(end(&r1), start(&r2) - 1); + clip_pair_hard(&make_clip(1, 2, 3, 4), &mut r1, &mut r2); + assert_clipped_by(prior1, &r1, 1, 2, "r1"); + assert_clipped_by(prior2, &r2, 3, 4, "r2"); + } + + /// `ClipBamTest.scala:182` "clip a fixed amount on the ends of the reads with reads with + /// clipping present": existing hard clipping counts toward the fixed amounts. + /// + /// Deviation: fgbio sets `4H46M` / `44M6H` on 50-base reads, so its SEQ is longer than + /// the CIGAR's query length. fgumi encodes records to raw BAM, which rejects that, so + /// the reads here carry 46 / 44 bases to match their CIGARs. The expected clipping is + /// unchanged. + #[test] + fn clip_pair_fixed_clipping_counts_existing_hard_clips() { + let (mut r1, mut r2) = fgbio_pair(50, |p| { + p.start1(104) + .start2(150) + .cigar1("4H46M") + .bases1(&"A".repeat(46)) + .cigar2("44M6H") + .bases2(&"A".repeat(44)) + }); + let (prior1, prior2) = ((start(&r1), end(&r1)), (start(&r2), end(&r2))); + assert_eq!(end(&r1), start(&r2) - 1); + clip_pair_hard(&make_clip(5, 2, 3, 4), &mut r1, &mut r2); + // R1: one more 5' base (4H already counts toward 5), two 3' bases. + assert_clipped_by(prior1, &r1, 1, 2, "r1"); + // R2: no more 5' bases (6H already exceeds 3), four 3' bases. + assert_clipped_by(prior2, &r2, 0, 4, "r2"); + } + + /// `ClipBamTest.scala:201` "clip a fixed amount on the ends of the reads then clip + /// overlapping reads": fixed clipping is applied first, then the remaining overlap. + #[test] + fn clip_pair_fixed_clipping_then_overlap_clipping() { + let (mut r1, mut r2) = fgbio_pair(50, |p| p.start1(100).start2(146)); + let (prior1, prior2) = ((start(&r1), end(&r1)), (start(&r2), end(&r2))); + assert_eq!(end(&r1), start(&r2) + 3, "four bases overlap"); + let mut clip = make_clip(0, 1, 0, 1); + clip.clip_overlapping_reads = true; + clip_pair_hard(&clip, &mut r1, &mut r2); + assert_eq!(end(&r1), start(&r2) - 1); + assert_clipped_by(prior1, &r1, 0, 2, "r1"); + assert_clipped_by(prior2, &r2, 0, 2, "r2"); + } + + /// `ClipBamTest.scala:216` "clip a fixed amount on the ends of the reads in + /// $strand1/$strand2", for all four strand combinations. + #[rstest] + #[case::plus_plus(Strand::Plus, Strand::Plus)] + #[case::plus_minus(Strand::Plus, Strand::Minus)] + #[case::minus_plus(Strand::Minus, Strand::Plus)] + #[case::minus_minus(Strand::Minus, Strand::Minus)] + fn clip_pair_fixed_clipping_is_strand_aware( + #[case] strand1: Strand, + #[case] strand2: Strand, + ) { + let (mut r1, mut r2) = + fgbio_pair(50, |p| p.start1(100).start2(150).strand1(strand1).strand2(strand2)); + let (prior1, prior2) = ((start(&r1), end(&r1)), (start(&r2), end(&r2))); + assert_eq!(end(&r1), start(&r2) - 1); + clip_pair_hard(&make_clip(1, 2, 3, 4), &mut r1, &mut r2); + assert_clipped_by(prior1, &r1, 1, 2, "r1"); + assert_clipped_by(prior2, &r2, 3, 4, "r2"); + } + + /// All 15 counters of a [`ClippingMetrics`], in fgbio's column order. + fn metric_values(m: &ClippingMetrics) -> [usize; 15] { + [ + m.reads, + m.reads_unmapped, + m.reads_clipped_pre, + m.reads_clipped_post, + m.reads_clipped_five_prime, + m.reads_clipped_three_prime, + m.reads_clipped_overlapping, + m.reads_clipped_extending, + m.bases, + m.bases_clipped_pre, + m.bases_clipped_post, + m.bases_clipped_five_prime, + m.bases_clipped_three_prime, + m.bases_clipped_overlapping, + m.bases_clipped_extending, + ] + } + + /// A [`ClippingMetrics`] whose 15 counters are `first, first + 1, ..., first + 14`, as + /// in fgbio's `ClippingMetrics(readType, 1, 2, 3, ...)` fixtures. + fn metrics_counting_up_from(read_type: ReadType, first: usize) -> ClippingMetrics { + let mut m = ClippingMetrics::new(read_type); + let fields = [ + &mut m.reads, + &mut m.reads_unmapped, + &mut m.reads_clipped_pre, + &mut m.reads_clipped_post, + &mut m.reads_clipped_five_prime, + &mut m.reads_clipped_three_prime, + &mut m.reads_clipped_overlapping, + &mut m.reads_clipped_extending, + &mut m.bases, + &mut m.bases_clipped_pre, + &mut m.bases_clipped_post, + &mut m.bases_clipped_five_prime, + &mut m.bases_clipped_three_prime, + &mut m.bases_clipped_overlapping, + &mut m.bases_clipped_extending, + ]; + for (offset, field) in fields.into_iter().enumerate() { + *field = first + offset; + } + m + } + + /// `ClipBamTest.scala:230` "add two metrics" (`ClippingMetrics.add`): every counter is + /// summed, and the sum differs from each input in every field. fgbio builds `readOne` + /// with `ReadType.ReadTwo`; that is kept. + #[test] + fn clipping_metrics_add_sums_every_field() { + let fragment = metrics_counting_up_from(ReadType::Fragment, 1); + let read_one = metrics_counting_up_from(ReadType::ReadTwo, 2); + let read_two = metrics_counting_up_from(ReadType::ReadTwo, 3); + let mut added = ClippingMetrics::new(ReadType::All); + added.add(&fragment); + added.add(&read_one); + added.add(&read_two); + + assert_eq!(added.read_type, ReadType::All); + assert_eq!( + metric_values(&added), + [6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48] + ); + assert_ne!(added.read_type, fragment.read_type); + for (field, (sum, input)) in + metric_values(&added).into_iter().zip(metric_values(&fragment)).enumerate() + { + assert_ne!(sum, input, "field {field} should differ from the fragment input"); + } + } + + /// The `(NM, UQ, MD)` htsjdk's `calculateMdAndNmTags` / `sumQualitiesOfMismatches` give + /// `rec` against fgbio's all-`A` reference. Supports only the CIGAR operators these + /// fixtures use. + fn expected_nm_uq_md(rec: &RecordBuf) -> (i64, i64, String) { + use std::fmt::Write as _; + let bases: &[u8] = rec.sequence().as_ref(); + let quals: &[u8] = rec.quality_scores().as_ref(); + let (mut offset, mut nm, mut uq, mut matches) = (0usize, 0i64, 0i64, 0usize); + let mut md = String::new(); + for op in rec.cigar().as_ref() { + match op.kind() { + Kind::Match | Kind::SequenceMatch | Kind::SequenceMismatch => { + for _ in 0..op.len() { + if bases[offset].eq_ignore_ascii_case(&b'A') { + matches += 1; + } else { + nm += 1; + uq += i64::from(quals[offset]); + write!(md, "{matches}A").expect("write to String"); + matches = 0; + } + offset += 1; + } + } + Kind::Insertion | Kind::SoftClip => offset += op.len(), + Kind::HardClip => {} + kind => panic!("unsupported CIGAR operator in fixture: {kind:?}"), + } + } + write!(md, "{matches}").expect("write to String"); + (nm, uq, md) + } + + /// Sets `NM`, `UQ` and `MD` on `rec` to their correct values, as fgbio's fixtures do + /// before running `ClipBam`. + fn set_nm_uq_md(rec: &mut RecordBuf) { + let (nm, uq, md) = expected_nm_uq_md(rec); + let data = rec.data_mut(); + data.insert( + Tag::from(SamTag::NM), + Value::from(i32::try_from(nm).expect("NM fits i32")), + ); + data.insert( + Tag::from(SamTag::UQ), + Value::from(i32::try_from(uq).expect("UQ fits i32")), + ); + data.insert(Tag::from(SamTag::MD), Value::from(md.as_str())); + } + + /// Asserts `rec` carries `NM`, `UQ` and `MD` equal to [`expected_nm_uq_md`]. + fn assert_nm_uq_md_recomputed(rec: &RecordBuf) { + let (nm, uq, md) = expected_nm_uq_md(rec); + let int_tag = |tag: SamTag| { + rec.data() + .get(&Tag::from(tag)) + .and_then(Value::as_int) + .unwrap_or_else(|| panic!("{tag:?} missing or not an integer")) + }; + let md_tag = match rec.data().get(&Tag::from(SamTag::MD)) { + Some(Value::String(s)) => s.to_string(), + other => panic!("MD missing or not a string: {other:?}"), + }; + let label = cigar_string(rec); + assert_eq!((int_tag(SamTag::NM), int_tag(SamTag::UQ), md_tag), (nm, uq, md), "{label}"); + } + + /// The 50-base all-`A` read fgbio's NM/UQ/MD fixtures use, with a single `C` at + /// `mismatch_at`. + fn all_a_with_mismatch(mismatch_at: usize) -> String { + let mut bases = vec![b'A'; 50]; + bases[mismatch_at] = b'C'; + String::from_utf8(bases).expect("ASCII bases") + } + + /// The first 12 values of Java's `new Random(1).nextInt(50)`: the mismatch positions + /// fgbio's NM/UQ/MD fixtures draw, in the order its `SamBuilder` yields the reads. + const FGBIO_MISMATCH_POSITIONS: [usize; 12] = [35, 38, 47, 13, 4, 4, 34, 6, 28, 48, 19, 23]; + + /// A `Clip` reading `input` and writing `output` (with `metrics`) against `reference`, + /// in `mode`, with all clipping disabled; callers enable what their case needs. + fn end_to_end_clip( + input: PathBuf, + output: PathBuf, + reference: PathBuf, + metrics: PathBuf, + mode: ClippingMode, + ) -> Clip { + let mut clip = make_clip(0, 0, 0, 0); + clip.io.input = input; + clip.io.output = output; + clip.reference = reference; + clip.metrics = Some(metrics); + clip.clipping_mode = mode; + clip + } + + /// Reads `clip`'s metrics file back. + fn read_clipping_metrics(path: &std::path::Path) -> Vec { + crate::metrics::read_metrics(path, "clipping").expect("read clipping metrics") + } + + /// The first and last CIGAR operations of a mapped record. + fn first_and_last_ops(rec: &RecordBuf) -> (Op, Op) { + let ops = rec.cigar().as_ref(); + (*ops.first().expect("non-empty CIGAR"), *ops.last().expect("non-empty CIGAR")) + } + + /// `ClipBamTest.scala:245` "clip overlapping reads, update mate info, and reset NM, UQ & + /// MD": six pairs overlapping by 10, 8, 6, 4, 2 and 0 bases, with fixed 5' clipping (R1 + /// 2, R2 3) and overlap clipping in `Hard` mode. Checks clipping, coordinate order, mate + /// info (`MC`), recomputed `NM`/`UQ`/`MD`, and every metrics row. + #[test] + fn clip_execute_pairs_updates_mate_info_tags_and_metrics() { + let dir = TempDir::new().expect("temp dir"); + let reference = write_all_a_reference(&dir); + let mut source = SamBuilder::with_single_ref("chr1", REFERENCE_LENGTH); + let mut input = SamBuilder::with_single_ref("chr1", REFERENCE_LENGTH); + input.set_queryname_sort_order(); + let starts = [(100, 140), (200, 242), (300, 344), (400, 446), (500, 548), (600, 650)]; + for (i, (start1, start2)) in starts.into_iter().enumerate() { + let (mut r1, mut r2) = source + .add_pair() + .name(&format!("q{}", i + 1)) + .bases1(&all_a_with_mismatch(FGBIO_MISMATCH_POSITIONS[2 * i])) + .bases2(&all_a_with_mismatch(FGBIO_MISMATCH_POSITIONS[2 * i + 1])) + .start1(start1) + .start2(start2) + .build(); + set_nm_uq_md(&mut r1); + set_nm_uq_md(&mut r2); + input.push_record(r1); + input.push_record(r2); + } + let (input_path, output, metrics) = + (dir.path().join("in.bam"), dir.path().join("out.bam"), dir.path().join("m.txt")); + input.write(&input_path).expect("write input"); + let mut clip = end_to_end_clip( + input_path, + output.clone(), + reference, + metrics.clone(), + ClippingMode::Hard, + ); + clip.read_one_five_prime = 2; + clip.read_two_five_prime = 3; + clip.clip_overlapping_reads = true; + clip.execute("test").expect("clip should succeed"); + + let clipped = read_bam_records(&output).expect("read output"); + assert_eq!(clipped.len(), 12); + let is_hard = |op: Op| op.kind() == Kind::HardClip; + let is_hard_of = |op: Op, len: usize| is_hard(op) && op.len() == len; + let (minus, plus): (Vec<_>, Vec<_>) = + clipped.iter().partition(|r| r.flags().is_reverse_complemented()); + let count = |recs: &[&RecordBuf], pred: &dyn Fn((Op, Op)) -> bool| { + recs.iter().filter(|r| pred(first_and_last_ops(r))).count() + }; + // Overlap clipping hit every pair but q6. + assert_eq!(count(&minus, &|(first, _)| is_hard(first)), 5); + assert_eq!(count(&plus, &|(_, last)| is_hard(last)), 5); + // Fixed 5' clipping hit every read. + assert_eq!(count(&minus, &|(_, last)| is_hard_of(last, 3)), 6); + assert_eq!(count(&plus, &|(first, _)| is_hard_of(first, 2)), 6); + + for pair in clipped.windows(2) { + assert!(pair[0].alignment_start() <= pair[1].alignment_start(), "coordinate order"); + } + let mut templates: std::collections::BTreeMap<_, Vec<&RecordBuf>> = + std::collections::BTreeMap::new(); + for rec in &clipped { + templates.entry(rec.name()).or_default().push(rec); + } + assert_eq!(templates.len(), 6); + for template in templates.values() { + let [lhs, rhs] = template.as_slice() else { + panic!("expected a pair, got {template:?}") + }; + for (a, b) in [(lhs, rhs), (rhs, lhs)] { + assert_eq!(a.mate_alignment_start(), b.alignment_start(), "mate start"); + assert_eq!( + a.flags().is_mate_reverse_complemented(), + b.flags().is_reverse_complemented(), + "mate strand" + ); + match a.data().get(&Tag::from(SamTag::MC)) { + Some(Value::String(mc)) => assert_eq!(mc.to_string(), cigar_string(b)), + other => panic!("MC missing or not a string: {other:?}"), + } + } + } + for rec in &clipped { + assert_nm_uq_md_recomputed(rec); + } + + let rows = read_clipping_metrics(&metrics); + let row = |read_type: ReadType| { + let row = rows.iter().find(|m| m.read_type == read_type); + metric_values(row.unwrap_or_else(|| panic!("no {read_type:?} row"))) + }; + // Columns: reads, unmapped, clipped pre/post/5'/3'/overlapping/extending, bases, + // bases clipped pre/post/5'/3'/overlapping/extending. + let read_one = [6, 0, 0, 6, 6, 0, 5, 0, 273, 0, 27, 12, 0, 15, 0]; + let read_two = [6, 0, 0, 6, 6, 0, 5, 0, 267, 0, 33, 18, 0, 15, 0]; + let pair: [usize; 15] = std::array::from_fn(|i| read_one[i] + read_two[i]); + assert_eq!(rows.len(), 5); + assert_eq!(row(ReadType::Fragment), [0; 15], "Fragment"); + assert_eq!(row(ReadType::ReadOne), read_one, "ReadOne"); + assert_eq!(row(ReadType::ReadTwo), read_two, "ReadTwo"); + assert_eq!(row(ReadType::Pair), pair, "Pair"); + assert_eq!(row(ReadType::All), pair, "All"); + } + + /// `ClipBamTest.scala:351` "clip fragment reads, and reset NM, UQ & MD": three fragments + /// with R1 fixed clipping (5' 2, 3' 10); R2 clipping and overlap clipping do not apply to + /// fragments, and the `40S10M` fragment is clipped away entirely and unmapped. + #[test] + fn clip_execute_fragments_resets_tags_and_metrics() { + let dir = TempDir::new().expect("temp dir"); + let reference = write_all_a_reference(&dir); + let mut source = SamBuilder::with_single_ref("chr1", REFERENCE_LENGTH); + let mut input = SamBuilder::with_single_ref("chr1", REFERENCE_LENGTH); + input.set_queryname_sort_order(); + let fragments = [ + (100, Strand::Plus, "50M"), + (200, Strand::Minus, "50M"), + (300, Strand::Plus, "40S10M"), + ]; + for (i, (start, strand, cigar)) in fragments.into_iter().enumerate() { + let mut rec = source + .add_frag() + .name(&format!("f{}", i + 1)) + .bases(&all_a_with_mismatch(FGBIO_MISMATCH_POSITIONS[i])) + .start(start) + .strand(strand) + .cigar(cigar) + .build(); + set_nm_uq_md(&mut rec); + input.push_record(rec); + } + let (input_path, output, metrics) = + (dir.path().join("in.bam"), dir.path().join("out.bam"), dir.path().join("m.txt")); + input.write(&input_path).expect("write input"); + let mut clip = end_to_end_clip( + input_path, + output.clone(), + reference, + metrics.clone(), + ClippingMode::Hard, + ); + (clip.read_one_five_prime, clip.read_one_three_prime) = (2, 10); + (clip.read_two_five_prime, clip.read_two_three_prime) = (5, 5); + clip.clip_overlapping_reads = true; + clip.execute("test").expect("clip should succeed"); + + let clipped = read_bam_records(&output).expect("read output"); + assert_eq!(clipped.len(), 3); + let mapped: Vec<_> = clipped.iter().filter(|r| !r.flags().is_unmapped()).collect(); + let is_2h = |op: Op| op.kind() == Kind::HardClip && op.len() == 2; + let minus_5p_clipped = mapped + .iter() + .filter(|r| r.flags().is_reverse_complemented() && is_2h(first_and_last_ops(r).1)) + .count(); + let plus_5p_clipped = mapped + .iter() + .filter(|r| !r.flags().is_reverse_complemented() && is_2h(first_and_last_ops(r).0)) + .count(); + assert_eq!((minus_5p_clipped, plus_5p_clipped), (1, 1)); + for pair in clipped.windows(2) { + let (lhs, rhs) = (&pair[0], &pair[1]); + if !lhs.flags().is_unmapped() && !rhs.flags().is_unmapped() { + assert!(lhs.alignment_start() <= rhs.alignment_start(), "coordinate order"); + } else if lhs.flags().is_unmapped() { + assert!(rhs.flags().is_unmapped(), "unmapped reads sort last"); + } + } + for rec in &mapped { + assert_nm_uq_md_recomputed(rec); + } + + let rows = read_clipping_metrics(&metrics); + // Columns as in `clip_execute_pairs_updates_mate_info_tags_and_metrics`. + let fragment = [3, 1, 1, 3, 2, 3, 0, 0, 76, 40, 74, 4, 30, 0, 0]; + assert_eq!(rows.len(), 5); + for row in &rows { + let expected = match row.read_type { + ReadType::Fragment | ReadType::All => fragment, + ReadType::ReadOne | ReadType::ReadTwo | ReadType::Pair => [0; 15], + }; + assert_eq!(metric_values(row), expected, "{:?}", row.read_type); + } + } + + /// fgbio's `--upgrade-clipping` fixture: two 50 bp fragments, `q1` (`+`, start 100) and + /// `q2` (`-`, start 200), each clipped 10 bases at the 5' end and 4 at the 3' end in + /// `prior` mode (with auto-clip attributes). With `reclip_in_prior_mode`, each is then + /// clipped again (5' 5, 3' 2), which must be a no-op, as in fgbio's upgrade cases. + /// Returns the fragments as written to `clip`'s input, and `clip`'s output in `mode`. + fn run_upgrade_clipping( + prior: ClippingMode, + mode: ClippingMode, + reclip_in_prior_mode: bool, + ) -> (Vec, Vec) { + let dir = TempDir::new().expect("temp dir"); + let reference = write_all_a_reference(&dir); + let mut source = SamBuilder::with_single_ref("chr1", REFERENCE_LENGTH); + let mut input = SamBuilder::with_single_ref("chr1", REFERENCE_LENGTH); + input.set_queryname_sort_order(); + let quals: Vec = (0..50u8).map(|i| 10 + i % 30).collect(); + let fragments = [ + ("q1", "ACGTTGCAAC".repeat(5), 100, Strand::Plus), + ("q2", "TTGACCAGTA".repeat(5), 200, Strand::Minus), + ]; + let header = source.header.clone(); + let clipper = RawRecordClipper::with_auto_clip(prior, true); + for (name, bases, start, strand) in fragments { + let frag = source + .add_frag() + .name(name) + .bases(&bases) + .quals(&quals) + .start(start) + .strand(strand) + .build(); + let mut raw = encode_record_buf_to_raw(&frag, &header).expect("encode"); + assert_eq!(clipper.clip_5_prime_end_of_read_raw(&mut raw, 10), 10); + assert_eq!(clipper.clip_3_prime_end_of_read_raw(&mut raw, 4), 4); + if reclip_in_prior_mode { + assert_eq!(clipper.clip_5_prime_end_of_read_raw(&mut raw, 5), 0); + assert_eq!(clipper.clip_3_prime_end_of_read_raw(&mut raw, 2), 0); + } + input.push_record(raw_record_to_record_buf(&raw, &header).expect("decode")); + } + let (input_path, output) = (dir.path().join("in.bam"), dir.path().join("out.bam")); + input.write(&input_path).expect("write input"); + let mut clip = end_to_end_clip( + input_path, + output.clone(), + reference, + dir.path().join("m.txt"), + mode, + ); + clip.upgrade_clipping = true; + clip.execute("test").expect("clip should succeed"); + let clipped = read_bam_records(&output).expect("read output"); + assert_eq!(clipped.len(), 2); + (input.records().to_vec(), clipped) + } + + /// fgbio's `maskBases`/`maskQuals`: `values` with its first `leading` and last + /// `trailing` entries replaced by `fill`. + fn masked(values: &[u8], leading: usize, trailing: usize, fill: u8) -> Vec { + let mut out = values.to_vec(); + let len = out.len(); + out[..leading].fill(fill); + out[len - trailing..].fill(fill); + out + } + + /// Asserts `clipped` is `prior` with its clipping in `expected` mode, per fgbio's + /// upgrade expectations: `q1` (`+`) carries 10 5' / 4 3' clipped bases (`10?36M4?`), + /// `q2` (`-`) the mirror (`4?36M10?`). `prior_mode` says whether `prior`'s bases were + /// already hard-clipped away. + fn assert_clipping_in_mode( + prior_mode: ClippingMode, + expected: ClippingMode, + prior: &[RecordBuf], + clipped: &[RecordBuf], + ) { + let seq = |r: &RecordBuf| r.sequence().as_ref().to_vec(); + let qual = |r: &RecordBuf| r.quality_scores().as_ref().to_vec(); + for (prior, clipped, (leading, trailing)) in + [(&prior[0], &clipped[0], (10, 4)), (&prior[1], &clipped[1], (4, 10))] + { + let op = if expected == ClippingMode::Hard { 'H' } else { 'S' }; + assert_eq!(cigar_string(clipped), format!("{leading}{op}36M{trailing}{op}")); + let (want_seq, want_qual) = match expected { + ClippingMode::Soft => (seq(prior), qual(prior)), + ClippingMode::SoftWithMask => ( + masked(&seq(prior), leading, trailing, b'N'), + masked(&qual(prior), leading, trailing, fgumi_dna::MIN_PHRED), + ), + ClippingMode::Hard if prior_mode == ClippingMode::Hard => { + (seq(prior), qual(prior)) + } + ClippingMode::Hard => ( + seq(prior)[leading..50 - trailing].to_vec(), + qual(prior)[leading..50 - trailing].to_vec(), + ), + }; + assert_eq!( + seq(clipped), + want_seq, + "bases of {cigar}", + cigar = cigar_string(clipped) + ); + assert_eq!( + qual(clipped), + want_qual, + "quals of {cigar}", + cigar = cigar_string(clipped) + ); + } + } + + /// `ClipBamTest.scala:428` "upgrade existing clipping from $prior to $mode with + /// --upgrade-clipping". + #[rstest] + #[case::soft_to_soft_with_mask(ClippingMode::Soft, ClippingMode::SoftWithMask)] + #[case::soft_to_hard(ClippingMode::Soft, ClippingMode::Hard)] + #[case::soft_with_mask_to_hard(ClippingMode::SoftWithMask, ClippingMode::Hard)] + fn clip_execute_upgrades_existing_clipping( + #[case] prior: ClippingMode, + #[case] mode: ClippingMode, + ) { + let (input, clipped) = run_upgrade_clipping(prior, mode, true); + assert_clipping_in_mode(prior, mode, &input, &clipped); + } + + /// `ClipBamTest.scala:473` "not upgrade existing clipping from $prior to $mode with + /// --upgrade-clipping": clipping already at or above `mode` is left as it was. + #[rstest] + #[case::soft_to_soft(ClippingMode::Soft, ClippingMode::Soft)] + #[case::soft_with_mask_to_soft(ClippingMode::SoftWithMask, ClippingMode::Soft)] + #[case::soft_with_mask_to_soft_with_mask( + ClippingMode::SoftWithMask, + ClippingMode::SoftWithMask + )] + #[case::hard_to_hard(ClippingMode::Hard, ClippingMode::Hard)] + #[case::hard_to_soft_with_mask(ClippingMode::Hard, ClippingMode::SoftWithMask)] + #[case::hard_to_soft(ClippingMode::Hard, ClippingMode::Soft)] + fn clip_execute_does_not_downgrade_existing_clipping( + #[case] prior: ClippingMode, + #[case] mode: ClippingMode, + ) { + let (input, clipped) = run_upgrade_clipping(prior, mode, false); + assert_clipping_in_mode(prior, prior, &input, &clipped); + } + + /// `ClipBamTest.scala:518` "clip FR reads that extend past the mate". + #[test] + fn clip_pair_clips_reads_extending_past_mate() { + let (mut r1, mut r2) = fgbio_pair(50, |p| p.start1(100).start2(90)); + assert_eq!((end(&r1), end(&r2)), (149, 139)); + let mut clip = make_clip(0, 0, 0, 0); + clip.clip_extending_past_mate = true; + clip_pair_hard(&clip, &mut r1, &mut r2); + assert_eq!(start(&r1), start(&r2)); + assert_eq!(end(&r1), end(&r2)); + } + + /// A `Clip` with fixed clipping `(r1 5', r1 3', r2 5', r2 3')`, past-mate clipping, and + /// optionally overlap clipping. + fn past_mate_clip(fixed: [usize; 4], clip_overlapping: bool) -> Clip { + let mut clip = make_clip(fixed[0], fixed[1], fixed[2], fixed[3]); + clip.clip_extending_past_mate = true; + clip.clip_overlapping_reads = clip_overlapping; + clip + } + + /// fgbio's past-mate cases: an FR pair of `read_length` reads at `starts`, whose prior + /// ends must be `prior_ends`, clipped by `clip`, giving `(r1 start, r1 end, r2 start, r2 + /// end)`. + #[rstest] + // ClipBamTest.scala:530 "clip FR reads that extend past their mate and remove overlap" + #[case::past_mate_and_overlap(100, (100, 90), (199, 189), past_mate_clip([0, 0, 0, 0], true), (100, 144, 145, 189))] + // ClipBamTest.scala:551 "clip FR reads that extend past their mate with asymmetrical five prime hard clipping" + #[case::asymmetric_five_prime(200, (100, 90), (299, 289), past_mate_clip([10, 0, 50, 0], false), (110, 239, 110, 239))] + // ClipBamTest.scala:574 "clip FR reads that extend past their mate with some irrelevant three prime clipping and removal of overlap" + #[case::three_prime_and_overlap(200, (100, 90), (299, 289), past_mate_clip([0, 0, 0, 50], true), (100, 194, 195, 289))] + // ClipBamTest.scala:597 "clip FR reads that extend past their mate, overlap, and have clipping on the 3-prime side of one and the 5-prime side of another" + #[case::overlap_with_mixed_fixed(200, (140, 90), (339, 289), past_mate_clip([25, 0, 0, 175], true), (165, 264, 265, 289))] + fn clip_pair_past_mate_with_fixed_and_overlap_clipping( + #[case] read_length: usize, + #[case] starts: (usize, usize), + #[case] prior_ends: (usize, usize), + #[case] clip: Clip, + #[case] expected: (usize, usize, usize, usize), + ) { + let (mut r1, mut r2) = fgbio_pair(read_length, |p| p.start1(starts.0).start2(starts.1)); + assert_eq!((end(&r1), end(&r2)), prior_ends); + clip_pair_hard(&clip, &mut r1, &mut r2); + assert_eq!((start(&r1), end(&r1), start(&r2), end(&r2)), expected); + } + + /// `ClipBamTest.scala:621` "unmap reads when the hard clipping length requested is + /// greater than the length of the reads". fgbio's `UnmappedStart` (0) corresponds to BAM + /// `POS` -1 with no alignment start or end. + #[test] + fn clip_pair_unmaps_reads_when_clipping_exceeds_read_length() { + let (mut r1, mut r2) = fgbio_pair(100, |p| p.start1(100).start2(300)); + assert_eq!((end(&r1), end(&r2)), (199, 399)); + clip_pair_hard(&past_mate_clip([101, 0, 101, 0], true), &mut r1, &mut r2); + assert_eq!((r1.is_unmapped(), r2.is_unmapped()), (true, true)); + assert_eq!((r1.pos(), r2.pos()), (-1, -1)); + assert_eq!((r1.alignment_start_1based(), r2.alignment_start_1based()), (None, None)); + assert_eq!((r1.alignment_end_1based(), r2.alignment_end_1based()), (None, None)); + } + } } From 5c14575357f208889f98a17ebc36b46082a15566 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Sun, 4 Oct 2026 15:11:51 -0700 Subject: [PATCH 3/3] test(clip): fail instead of skipping when an asserted tag is missing 28 clipper tests checked per-base tag values inside `if let Some(Value::...)`, so a missing tag or a changed value type skipped the assertion and the test passed. They now use let-else and panic. --- crates/fgumi-sam/src/clipper.rs | 256 ++++++++++++++++++-------------- 1 file changed, 142 insertions(+), 114 deletions(-) diff --git a/crates/fgumi-sam/src/clipper.rs b/crates/fgumi-sam/src/clipper.rs index abdd0acad..ab26ea081 100644 --- a/crates/fgumi-sam/src/clipper.rs +++ b/crates/fgumi-sam/src/clipper.rs @@ -2602,10 +2602,11 @@ mod tests { clipper.clip_start_of_alignment(&mut record, 3); // Attribute should remain unchanged in Soft mode - if let Some(Value::String(s)) = record.data().get(&tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, b"0123456789"); - } + let Some(Value::String(s)) = record.data().get(&tag) else { + panic!("tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, b"0123456789"); } #[test] @@ -2623,10 +2624,11 @@ mod tests { clipper.clip_start_of_alignment(&mut record, 3); // Attribute should remain unchanged when auto-clip is disabled - if let Some(Value::String(s)) = record.data().get(&tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, b"0123456789"); - } + let Some(Value::String(s)) = record.data().get(&tag) else { + panic!("tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, b"0123456789"); } #[test] @@ -2648,16 +2650,18 @@ mod tests { clipper.clip_start_of_alignment(&mut record, 3); // Check tag1 was clipped - if let Some(Value::String(s)) = record.data().get(&tag1) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, b"3456789"); - } + let Some(Value::String(s)) = record.data().get(&tag1) else { + panic!("tag1 missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, b"3456789"); // Check tag2 was NOT clipped - if let Some(Value::String(s)) = record.data().get(&tag2) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, b"01234"); - } + let Some(Value::String(s)) = record.data().get(&tag2) else { + panic!("tag2 missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, b"01234"); } // =================================================================== @@ -2935,10 +2939,11 @@ mod tests { assert_eq!(clipped, 5); // In Soft mode with auto=false, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); } #[test] @@ -2956,10 +2961,11 @@ mod tests { assert_eq!(clipped, 5); // In Soft mode, even with auto=true, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); } #[test] @@ -2977,10 +2983,11 @@ mod tests { assert_eq!(clipped, 5); // In SoftWithMask mode with auto=false, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); } #[test] @@ -2998,10 +3005,11 @@ mod tests { assert_eq!(clipped, 5); // In SoftWithMask mode, even with auto=true, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); } #[test] @@ -3023,14 +3031,16 @@ mod tests { assert_eq!(clipped, 5); // In Hard mode with auto=false, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } - if let Some(Value::Array(Array::Int32(arr))) = record.data().get(&a2_tag) { - let vec: Vec = arr.clone(); - assert_eq!(vec, (1..=20).collect::>()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); + let Some(Value::Array(Array::Int32(arr))) = record.data().get(&a2_tag) else { + panic!("a2_tag missing or of unexpected type"); + }; + let vec: Vec = arr.clone(); + assert_eq!(vec, (1..=20).collect::>()); } #[test] @@ -3056,24 +3066,28 @@ mod tests { assert_eq!(clipped, 5); // In Hard mode with auto=true, attributes matching read length should be clipped - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - // "ABABABABABABABABABAB" -> remove first 5 -> "BABABABABABABAB" - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "BABABABABABABAB".as_bytes()); - } - if let Some(Value::Array(Array::Int32(arr))) = record.data().get(&a2_tag) { - let vec: Vec = arr.clone(); - assert_eq!(vec, (6..=20).collect::>()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + // "ABABABABABABABABABAB" -> remove first 5 -> "BABABABABABABAB" + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "BABABABABABABAB".as_bytes()); + let Some(Value::Array(Array::Int32(arr))) = record.data().get(&a2_tag) else { + panic!("a2_tag missing or of unexpected type"); + }; + let vec: Vec = arr.clone(); + assert_eq!(vec, (6..=20).collect::>()); // B1 and B2 should NOT be modified (length doesn't match) - if let Some(Value::String(s)) = record.data().get(&b1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "A".repeat(10).as_bytes()); - } - if let Some(Value::Array(Array::Int32(arr))) = record.data().get(&b2_tag) { - let vec: Vec = arr.clone(); - assert_eq!(vec, (1..=10).collect::>()); - } + let Some(Value::String(s)) = record.data().get(&b1_tag) else { + panic!("b1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "A".repeat(10).as_bytes()); + let Some(Value::Array(Array::Int32(arr))) = record.data().get(&b2_tag) else { + panic!("b2_tag missing or of unexpected type"); + }; + let vec: Vec = arr.clone(); + assert_eq!(vec, (1..=10).collect::>()); } #[test] @@ -3091,10 +3105,11 @@ mod tests { assert_eq!(clipped, 5); // In Soft mode with auto=false, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); } #[test] @@ -3112,10 +3127,11 @@ mod tests { assert_eq!(clipped, 5); // In Soft mode, even with auto=true, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); } #[test] @@ -3133,10 +3149,11 @@ mod tests { assert_eq!(clipped, 5); // In SoftWithMask mode with auto=false, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); } #[test] @@ -3154,10 +3171,11 @@ mod tests { assert_eq!(clipped, 5); // In SoftWithMask mode, even with auto=true, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); } #[test] @@ -3179,14 +3197,16 @@ mod tests { assert_eq!(clipped, 5); // In Hard mode with auto=false, attributes should NOT be modified - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "AB".repeat(10).as_bytes()); - } - if let Some(Value::Array(Array::Int32(arr))) = record.data().get(&a2_tag) { - let vec: Vec = arr.clone(); - assert_eq!(vec, (1..=20).collect::>()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "AB".repeat(10).as_bytes()); + let Some(Value::Array(Array::Int32(arr))) = record.data().get(&a2_tag) else { + panic!("a2_tag missing or of unexpected type"); + }; + let vec: Vec = arr.clone(); + assert_eq!(vec, (1..=20).collect::>()); } #[test] @@ -3212,24 +3232,28 @@ mod tests { assert_eq!(clipped, 5); // In Hard mode with auto=true, attributes matching read length should be clipped - if let Some(Value::String(s)) = record.data().get(&a1_tag) { - // "ABABABABABABABABABAB" -> remove last 5 -> "ABABABABABABABA" - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "ABABABABABABABA".as_bytes()); - } - if let Some(Value::Array(Array::Int32(arr))) = record.data().get(&a2_tag) { - let vec: Vec = arr.clone(); - assert_eq!(vec, (1..=15).collect::>()); - } + let Some(Value::String(s)) = record.data().get(&a1_tag) else { + panic!("a1_tag missing or of unexpected type"); + }; + // "ABABABABABABABABABAB" -> remove last 5 -> "ABABABABABABABA" + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "ABABABABABABABA".as_bytes()); + let Some(Value::Array(Array::Int32(arr))) = record.data().get(&a2_tag) else { + panic!("a2_tag missing or of unexpected type"); + }; + let vec: Vec = arr.clone(); + assert_eq!(vec, (1..=15).collect::>()); // B1 and B2 should NOT be modified (length doesn't match) - if let Some(Value::String(s)) = record.data().get(&b1_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "A".repeat(10).as_bytes()); - } - if let Some(Value::Array(Array::Int32(arr))) = record.data().get(&b2_tag) { - let vec: Vec = arr.clone(); - assert_eq!(vec, (1..=10).collect::>()); - } + let Some(Value::String(s)) = record.data().get(&b1_tag) else { + panic!("b1_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "A".repeat(10).as_bytes()); + let Some(Value::Array(Array::Int32(arr))) = record.data().get(&b2_tag) else { + panic!("b2_tag missing or of unexpected type"); + }; + let vec: Vec = arr.clone(); + assert_eq!(vec, (1..=10).collect::>()); } // =================================================================== @@ -3259,10 +3283,11 @@ mod tests { assert_eq!(no_auto.sequence().len(), 35); // Attributes should NOT be modified without auto-clip - if let Some(Value::String(s)) = no_auto.data().get(&az_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "12345678901234567890123456789012345678901234567890".as_bytes()); - } + let Some(Value::String(s)) = no_auto.data().get(&az_tag) else { + panic!("az_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "12345678901234567890123456789012345678901234567890".as_bytes()); // Test with auto-clip let clipper_auto = RawClipperOnBuf::with_auto_clip(ClippingMode::Hard, true); @@ -3281,10 +3306,11 @@ mod tests { assert_eq!(with_auto.sequence().len(), 35); // Attributes SHOULD be modified with auto-clip (remove first 5 and last 10) - if let Some(Value::String(s)) = with_auto.data().get(&az_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "67890123456789012345678901234567890".as_bytes()); - } + let Some(Value::String(s)) = with_auto.data().get(&az_tag) else { + panic!("az_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "67890123456789012345678901234567890".as_bytes()); } /// The soft→hard upgrade path (`upgrade_all_clipping`) must also clip a @@ -3341,10 +3367,11 @@ mod tests { assert_eq!(no_auto.sequence().len(), 35); // Attributes should NOT be modified without auto-clip - if let Some(Value::String(s)) = no_auto.data().get(&az_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "12345678901234567890123456789012345678901234567890".as_bytes()); - } + let Some(Value::String(s)) = no_auto.data().get(&az_tag) else { + panic!("az_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "12345678901234567890123456789012345678901234567890".as_bytes()); // Test with auto-clip let clipper_auto = RawClipperOnBuf::with_auto_clip(ClippingMode::Hard, true); @@ -3363,10 +3390,11 @@ mod tests { assert_eq!(with_auto.sequence().len(), 35); // Attributes SHOULD be modified with auto-clip - if let Some(Value::String(s)) = with_auto.data().get(&az_tag) { - let bytes: &[u8] = s.as_ref(); - assert_eq!(bytes, "67890123456789012345678901234567890".as_bytes()); - } + let Some(Value::String(s)) = with_auto.data().get(&az_tag) else { + panic!("az_tag missing or of unexpected type"); + }; + let bytes: &[u8] = s.as_ref(); + assert_eq!(bytes, "67890123456789012345678901234567890".as_bytes()); } #[test]