diff --git a/crates/fgumi-sam/src/clipper.rs b/crates/fgumi-sam/src/clipper.rs index f75c89647..ab26ea081 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); @@ -1386,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, @@ -2589,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] @@ -2610,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] @@ -2635,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"); } // =================================================================== @@ -2922,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] @@ -2943,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] @@ -2964,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] @@ -2985,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] @@ -3010,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] @@ -3043,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] @@ -3078,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] @@ -3099,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] @@ -3120,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] @@ -3141,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] @@ -3166,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] @@ -3199,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::>()); } // =================================================================== @@ -3246,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); @@ -3268,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 @@ -3328,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); @@ -3350,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] @@ -3491,7 +3532,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 +3572,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; @@ -4806,6 +4849,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 fa626892e..82a66444f 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)?; @@ -3503,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)); + } + } }