diff --git a/src/lib/commands/clip.rs b/src/lib/commands/clip.rs index a267a5d17..caf4edbe7 100644 --- a/src/lib/commands/clip.rs +++ b/src/lib/commands/clip.rs @@ -13,7 +13,9 @@ use crate::metrics::writer::write_metrics as write_metrics_tsv; use crate::per_thread_accumulator::PerThreadAccumulator; use crate::reference::ReferenceReader; use crate::sam::SamTag; -use crate::template::{TemplateBatch, TemplateIterator}; +use crate::template::{ + InsertSizeEnd, TemplateBatch, TemplateIterator, compute_insert_size_from_ends, +}; use crate::unified_pipeline::{ GroupKeyConfig, Grouper, MemoryEstimate, run_bam_pipeline_from_reader, }; @@ -853,6 +855,23 @@ struct MateSnap { } impl MateSnap { + /// Snapshots only the fields `compute_insert_size_raw` reads, skipping the `cigar_to_string` + /// allocation. `cigar` and `mapq` are left empty/zero, so the result must not be used to write + /// a mate's MC or MQ tag. + fn coords_of(rec: &RawRecord) -> Self { + use fgumi_raw_bam::flags as rflags; + Self { + ref_id: rec.ref_id(), + pos: rec.pos(), + neg: rec.flags() & rflags::REVERSE != 0, + unmapped: rec.flags() & rflags::UNMAPPED != 0, + mapq: 0, + cigar: String::new(), + aln_start: rec.alignment_start_1based(), + aln_end: rec.alignment_end_1based(), + } + } + fn of(rec: &RawRecord) -> Self { use fgumi_raw_bam::flags as rflags; Self { @@ -902,17 +921,21 @@ fn clear_mate_mq_mc_raw(rec: &mut RawRecord) { /// Computes the TLEN (inferred insert size) for the first read of a pair, mirroring /// htsjdk `SamPairUtil.computeInsertSize`. Returns 0 unless both reads are mapped to /// the same reference; the second read's TLEN is the negation of this value. +/// +/// Only the 5'-position extraction lives here — `MateSnap` caches the alignment ends, and a +/// read whose CIGAR yields none has no insert size. The arithmetic itself is +/// [`compute_insert_size_from_ends`], shared with `template`'s supplementary-TLEN fix-up so +/// the two commands cannot drift apart on ties, overflow, or the htsjdk adjustment. fn compute_insert_size_raw(s1: &MateSnap, s2: &MateSnap) -> i32 { - if s1.unmapped || s2.unmapped || s1.ref_id != s2.ref_id { - return 0; - } let five_prime = |s: &MateSnap| if s.neg { s.aln_end } else { s.aln_start }; let (Some(p1), Some(p2)) = (five_prime(s1), five_prime(s2)) else { return 0; }; - let (p1, p2) = (p1 as i64, p2 as i64); - let adjustment = if p2 >= p1 { 1 } else { -1 }; - i32::try_from(p2 - p1 + adjustment).unwrap_or(0) + let to_i64 = |p: usize| i64::try_from(p).unwrap_or(i64::MAX); + compute_insert_size_from_ends( + InsertSizeEnd::new(s1.ref_id, s1.unmapped, to_i64(p1)), + InsertSizeEnd::new(s2.ref_id, s2.unmapped, to_i64(p2)), + ) } /// Sets full mate-pair information on a read pair, mirroring htsjdk @@ -989,19 +1012,28 @@ fn set_mate_info_raw(r1: &mut RawRecord, r2: &mut RawRecord) { } /// Ports htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment(supp, matePrimary, -/// setMateCigar=true)`. +/// setMateCigar=true)` **except for TLEN**, where it intentionally diverges. +/// +/// htsjdk sets the supplementary's TLEN to the negation of the mate primary's. That assumes the +/// supplementary sits where its own primary sits, which is false by construction for split +/// alignments: the primary's TLEN describes coordinates this record does not occupy. Copying it +/// yields a non-zero TLEN across references, and the wrong sign and magnitude when the +/// supplementary lies beyond its mate. This computes TLEN from the supplementary's own alignment +/// against the mate primary instead, matching bwa-mem and minibwa. See issue #673 and +/// samtools/htsjdk#1795. /// -/// fgbio `ClipBam` calls this on every supplementary alignment after clipping the primary pair, -/// so a supplemental's mate fields point at its mate *primary* read. It sets mate ref/pos/strand, -/// the mate-unmapped flag, TLEN (the negation of the mate primary's already-updated TLEN), the -/// mate CIGAR (MC, only when the mate is mapped) and — as of htsjdk 5.0.0 — the mate mapping -/// quality (MQ, unconditionally). `mate` and `mate_tlen` are snapshotted from the post-clip -/// primary so the caller can fix several supplementals without re-borrowing the primaries. -fn set_supplemental_mate_info_raw(supp: &mut RawRecord, mate: &MateSnap, mate_tlen: i32) { +/// Everything else follows htsjdk: fgbio `ClipBam` calls this on every supplementary alignment +/// after clipping the primary pair, so a supplemental's mate fields point at its mate *primary* +/// read. It sets mate ref/pos/strand, the mate-unmapped flag, the mate CIGAR (MC, only when the +/// mate is mapped) and — as of htsjdk 5.0.0 — the mate mapping quality (MQ, unconditionally). +/// `mate` is snapshotted from the post-clip primary so the caller can fix several supplementals +/// without re-borrowing the primaries. +fn set_supplemental_mate_info_raw(supp: &mut RawRecord, mate: &MateSnap) { + let tlen = compute_insert_size_raw(&MateSnap::coords_of(supp), mate); supp.set_mate_ref_id(mate.ref_id); supp.set_mate_pos(mate.pos); set_mate_flags_raw(supp, mate.neg, mate.unmapped); - supp.set_template_length(-mate_tlen); + supp.set_template_length(tlen); let mut editor = supp.tags_editor(); if mate.unmapped { editor.remove(SamTag::MC); @@ -1054,9 +1086,7 @@ fn find_primary_pair_indices(records: &[RawRecord]) -> Result<(Option, Op fn fix_supplemental_mate_info(records: &mut [RawRecord], r1_idx: usize, r2_idx: usize) { // Snapshot the post-clip primaries so the per-supplemental updates don't re-borrow them. let r1_snap = MateSnap::of(&records[r1_idx]); - let r1_tlen = records[r1_idx].template_length(); let r2_snap = MateSnap::of(&records[r2_idx]); - let r2_tlen = records[r2_idx].template_length(); for rec in records.iter_mut() { if !rec.is_supplementary() { @@ -1064,9 +1094,9 @@ fn fix_supplemental_mate_info(records: &mut [RawRecord], r1_idx: usize, r2_idx: } // R1 supplementals (unpaired or first-of-pair) take R2 as their mate; R2 supplementals R1. if !rec.is_paired() || rec.is_first_segment() { - set_supplemental_mate_info_raw(rec, &r2_snap, r2_tlen); + set_supplemental_mate_info_raw(rec, &r2_snap); } else if rec.is_last_segment() { - set_supplemental_mate_info_raw(rec, &r1_snap, r1_tlen); + set_supplemental_mate_info_raw(rec, &r1_snap); } } } @@ -1234,7 +1264,8 @@ mod tests { } // set_supplemental_mate_info_raw copies a mapped mate's coordinate/strand/MAPQ onto a - // supplementary read, writes MC from the mate CIGAR, sets MQ, and negates the mate TLEN. + // supplementary read, writes MC from the mate CIGAR, sets MQ, and computes TLEN from the + // supplementary's own alignment against the mate. #[test] fn test_set_supplemental_mate_info_raw_mapped_mate() { use fgumi_raw_bam::flags as rflags; @@ -1243,11 +1274,13 @@ mod tests { let mate = raw_read(false, false, true, 301, 40); let mut supp = raw_read(true, true, false, 700, 30); - set_supplemental_mate_info_raw(&mut supp, &MateSnap::of(&mate), 120); + set_supplemental_mate_info_raw(&mut supp, &MateSnap::of(&mate)); assert_eq!(supp.mate_ref_id(), 0); assert_eq!(supp.mate_pos(), 300); - assert_eq!(supp.template_length(), -120); + // The supplementary (700..749, forward, 5' = 700) is the rightmost segment; the mate + // (301..350, reverse, 5' = 350) is leftmost. TLEN = 350 - 700 - 1. + assert_eq!(supp.template_length(), -351); assert_ne!(supp.flags() & rflags::MATE_REVERSE, 0, "mate is reverse"); assert_eq!(supp.flags() & rflags::MATE_UNMAPPED, 0, "mate is mapped"); assert_eq!(supp.tags().find_mc(), Some("50M")); @@ -1281,8 +1314,9 @@ mod tests { supp.tags_editor().update_string(SamTag::MC, b"10M"); assert!(supp.tags().contains(SamTag::MC), "MC present before"); - set_supplemental_mate_info_raw(&mut supp, &MateSnap::of(&mate), 0); + set_supplemental_mate_info_raw(&mut supp, &MateSnap::of(&mate)); + assert_eq!(supp.template_length(), 0, "TLEN is 0 when the mate is unmapped"); assert_ne!(supp.flags() & rflags::MATE_UNMAPPED, 0, "mate is unmapped"); assert!(!supp.tags().contains(SamTag::MC), "MC dropped when mate unmapped"); // MQ is set unconditionally to the mate's mapping quality, even for an unmapped mate. @@ -1290,7 +1324,8 @@ mod tests { } // fix_supplemental_mate_info points R1 supplementals at the primary R2 and R2 supplementals - // at the primary R1, inheriting each primary's coordinate/strand/MAPQ and negated TLEN. + // at the primary R1, inheriting each primary's coordinate/strand/MAPQ and computing TLEN from + // the supplementary's own alignment. #[test] fn test_fix_supplemental_mate_info() { use fgumi_raw_bam::flags as rflags; @@ -1311,14 +1346,62 @@ mod tests { assert_eq!(recs[2].mate_pos(), 300); assert_ne!(recs[2].flags() & rflags::MATE_REVERSE, 0, "primary R2 is reverse"); assert_eq!(recs[2].tags().find_int(SamTag::MQ), Some(40)); - assert_eq!(recs[2].template_length(), 200); // -(-200) + // Supp R1 (701..750, forward, 5' = 701) lies right of primary R2 (5' = 350): 350-701-1. + assert_eq!(recs[2].template_length(), -352); // Supp R2 (idx 3) takes primary R1 (idx 0) as its mate. assert_eq!(recs[3].mate_ref_id(), 0); assert_eq!(recs[3].mate_pos(), 100); assert_eq!(recs[3].flags() & rflags::MATE_REVERSE, 0, "primary R1 is forward"); assert_eq!(recs[3].tags().find_int(SamTag::MQ), Some(60)); - assert_eq!(recs[3].template_length(), -200); + // Supp R2 (901..950, forward, 5' = 901) lies right of primary R1 (5' = 101): 101-901-1. + assert_eq!(recs[3].template_length(), -801); + } + + /// Builds a mapped record on an explicit reference for the supplementary-TLEN cases below. + /// + /// The builder validates `reference_sequence_id` against its single-contig test header, so the + /// record is built on reference 0 and the id is stamped onto the encoded record afterwards. + fn raw_read_on_ref( + first: bool, + supplementary: bool, + reverse: bool, + ref_id: i32, + start: usize, + ) -> RawRecord { + let mut rec = raw_read(first, supplementary, reverse, start, 60); + rec.set_ref_id(ref_id); + rec + } + + /// A supplementary's TLEN is computed from its own alignment against the mate primary, not + /// copied from the mate primary's TLEN. See issue #673 and samtools/htsjdk#1795. + /// + /// The mate primary is always a reverse-strand read at 1-based 301 (50M, so 5' = 350). + #[rstest] + // Supplementary right of the mate: it is the rightmost segment, so TLEN is negative. + #[case::beyond_mate(0, 700, -351)] + // Supplementary left of the mate: it is the leftmost segment, so TLEN is positive. + #[case::before_mate(0, 100, 251)] + // Coincident 5' ends: htsjdk's convention gives the two ends differing signs via the +1/-1 + // adjustment, so the leftmost-by-tie-break gets +1. + #[case::coincident_five_prime(0, 350, 1)] + // Different reference: the information is unavailable, so TLEN is 0. + #[case::cross_reference(1, 700, 0)] + fn test_supplemental_tlen_is_computed_not_copied( + #[case] supp_ref_id: i32, + #[case] supp_start: usize, + #[case] expected_tlen: i32, + ) { + let mate = raw_read_on_ref(false, false, true, 0, 301); + let mut supp = raw_read_on_ref(true, true, false, supp_ref_id, supp_start); + // Seed a value that the old copy-the-mate's-TLEN behaviour would have propagated, so a + // regression cannot pass by coincidence. + supp.set_template_length(-9999); + + set_supplemental_mate_info_raw(&mut supp, &MateSnap::of(&mate)); + + assert_eq!(supp.template_length(), expected_tlen); } // The RecordBuf unsoftclipped_start/end helpers (used by the RecordBuf clip path) subtract or diff --git a/src/lib/commands/zipper.rs b/src/lib/commands/zipper.rs index 4b5a14a20..b7ea06dd8 100644 --- a/src/lib/commands/zipper.rs +++ b/src/lib/commands/zipper.rs @@ -3245,16 +3245,19 @@ mod tests { "R1 supp ms should be R2 primary's AS (55), got {supp_ms:?}" ); - // TLEN on supplementary should be negative of R2 primary's TLEN - let r1_primary_tlen = r1_primary.template_length(); - let r1_supp_tlen = r1_supp.template_length(); - // Supplementary TLEN = -(R2 primary TLEN) = -(-R1 primary TLEN) = R1 primary TLEN... no - // Actually: supplementary's TLEN = -(mate primary's TLEN) = -(R2's TLEN) = R1's TLEN - // Wait, the code says: *self.records[i].template_length_mut() = -r2_tlen; - // R2's TLEN is the negative of R1's TLEN, so -r2_tlen = R1's TLEN + // TLEN on the supplementary is computed from its own alignment against the mate primary, + // not copied from that primary's TLEN. See issue #673. + // + // R1 primary 100..149 forward -> 5' = 100 + // R2 primary 300..374 reverse -> 5' = 374 => primaries are +275 / -275 + // R1 supp 500..559 reverse -> 5' = 559 + // + // The supplementary is the rightmost segment, so its TLEN is negative: 374 - 559 - 1. + assert_eq!(r1_primary.template_length(), 275, "R1 primary TLEN"); assert_eq!( - r1_supp_tlen, r1_primary_tlen, - "R1 supp TLEN should equal R1 primary TLEN (both = -R2_TLEN)" + r1_supp.template_length(), + -186, + "R1 supp TLEN is computed against the mate primary, not copied from it" ); Ok(()) diff --git a/src/lib/template.rs b/src/lib/template.rs index f2524fb03..921de6bc4 100644 --- a/src/lib/template.rs +++ b/src/lib/template.rs @@ -518,7 +518,7 @@ impl Template { let r2_flags = RawRecordView::new(&rr[r2_i]).flags(); let r2_is_reverse = (r2_flags & fgumi_raw_bam::flags::REVERSE) != 0; let r2_is_unmapped = (r2_flags & fgumi_raw_bam::flags::UNMAPPED) != 0; - let r2_tlen = fgumi_raw_bam::template_length(&rr[r2_i]); + let r2_end = insert_size_end(&rr[r2_i]); let r2_mapq = fgumi_raw_bam::mapq(&rr[r2_i]); let r2_cigar_str = fgumi_raw_bam::cigar_to_string_from_raw(&rr[r2_i]); let r2_as = @@ -530,7 +530,12 @@ impl Template { fgumi_raw_bam::set_mate_ref_id(rec, r2_ref_id); fgumi_raw_bam::set_mate_pos(rec, r2_pos); set_mate_flags(rec, r2_is_reverse, r2_is_unmapped); - fgumi_raw_bam::set_template_length(rec, -r2_tlen); + // TLEN is computed from the supplementary's own alignment against the mate + // primary, not copied from the mate primary's TLEN: a supplementary sits at a + // different locus than its own primary, so the primary's value describes + // coordinates this record does not occupy. See issue #673. + let tlen = compute_insert_size_from_ends(insert_size_end(rec), r2_end); + fgumi_raw_bam::set_template_length(rec, tlen); let mq_val = if r2_mapq == 255 { 255 } else { i32::from(r2_mapq) }; fgumi_raw_bam::update_int_tag(rec.as_mut_vec(), SamTag::MQ, mq_val); @@ -563,7 +568,7 @@ impl Template { let r1_flags = RawRecordView::new(&rr[r1_i]).flags(); let r1_is_reverse = (r1_flags & fgumi_raw_bam::flags::REVERSE) != 0; let r1_is_unmapped = (r1_flags & fgumi_raw_bam::flags::UNMAPPED) != 0; - let r1_tlen = fgumi_raw_bam::template_length(&rr[r1_i]); + let r1_end = insert_size_end(&rr[r1_i]); let r1_mapq = fgumi_raw_bam::mapq(&rr[r1_i]); let r1_cigar_str = fgumi_raw_bam::cigar_to_string_from_raw(&rr[r1_i]); let r1_as = @@ -575,7 +580,9 @@ impl Template { fgumi_raw_bam::set_mate_ref_id(rec, r1_ref_id); fgumi_raw_bam::set_mate_pos(rec, r1_pos); set_mate_flags(rec, r1_is_reverse, r1_is_unmapped); - fgumi_raw_bam::set_template_length(rec, -r1_tlen); + // See the R1 supplemental loop above: computed, not copied. Issue #673. + let tlen = compute_insert_size_from_ends(insert_size_end(rec), r1_end); + fgumi_raw_bam::set_template_length(rec, tlen); let mq_val = if r1_mapq == 255 { 255 } else { i32::from(r1_mapq) }; fgumi_raw_bam::update_int_tag(rec.as_mut_vec(), SamTag::MQ, mq_val); @@ -813,41 +820,82 @@ fn set_mate_flags(record: &mut [u8], mate_is_reverse: bool, mate_is_unmapped: bo fgumi_raw_bam::set_flags(record, f); } -/// Computes insert size (TLEN) from two raw BAM records. +/// The fields of one record needed to compute an insert size against another. /// -/// Uses 0-based pos from BAM; converts to 1-based for the calculation. -fn compute_insert_size_raw(rec1: &[u8], rec2: &[u8]) -> i32 { - use fgumi_raw_bam; - - let f1 = RawRecordView::new(rec1).flags(); - let f2 = RawRecordView::new(rec2).flags(); +/// Snapshotting these lets the supplementary fix-up loops compute a TLEN against the +/// mate primary without holding a second borrow of `self.records` inside the loop. +/// +/// This is also the shared input to [`compute_insert_size_from_ends`], the one copy of the +/// htsjdk 5'-difference-plus-adjustment formula. `zipper` (through [`insert_size_end`]) and +/// `clip` (through its `MateSnap`) reach it from different record representations — raw BAM +/// bytes versus a pre-mutation snapshot — but must not each carry the arithmetic: a future +/// tie-break or overflow fix applied to only one copy would silently diverge that command's +/// output from fgbio/htsjdk. +#[derive(Clone, Copy)] +pub(crate) struct InsertSizeEnd { + ref_id: i32, + is_unmapped: bool, + /// 1-based 5' position: alignment start for forward reads, alignment end for reverse. + five_prime: i32, +} - // If either read is unmapped, return 0 - if (f1 & fgumi_raw_bam::flags::UNMAPPED) != 0 || (f2 & fgumi_raw_bam::flags::UNMAPPED) != 0 { - return 0; +impl InsertSizeEnd { + /// Builds an end from a 1-based 5' position, saturating one that does not fit in `i32`. + /// + /// A position past `i32::MAX` is reachable only from a malformed record; saturating keeps + /// [`compute_insert_size_from_ends`] from overflowing rather than wrapping. + pub(crate) fn new(ref_id: i32, is_unmapped: bool, five_prime: i64) -> Self { + Self { ref_id, is_unmapped, five_prime: i32::try_from(five_prime).unwrap_or(i32::MAX) } } +} - // If reads are on different references, return 0 - if fgumi_raw_bam::ref_id(rec1) != fgumi_raw_bam::ref_id(rec2) { - return 0; - } +/// Snapshots the insert-size-relevant fields of a raw BAM record. +/// +/// Uses 0-based pos from BAM; converts to 1-based for the calculation. +fn insert_size_end(rec: &[u8]) -> InsertSizeEnd { + use fgumi_raw_bam; - // pos is 0-based in BAM; convert to 1-based for the calculation - let pos1 = fgumi_raw_bam::pos(rec1) + 1; - let pos2 = fgumi_raw_bam::pos(rec2) + 1; + let flags = RawRecordView::new(rec).flags(); + // pos is 0-based in BAM; convert to 1-based for the calculation. Widened to i64 so a malformed + // record with an extreme POS cannot overflow before the unmapped/cross-reference guards run. + let start = i64::from(fgumi_raw_bam::pos(rec)) + 1; // alignment end (1-based inclusive) = pos_1based + ref_len - 1 - let ref_len1 = fgumi_raw_bam::reference_length_from_raw_bam(rec1); - let ref_len2 = fgumi_raw_bam::reference_length_from_raw_bam(rec2); - let end1 = pos1 + ref_len1 - 1; - let end2 = pos2 + ref_len2 - 1; + let end = start + i64::from(fgumi_raw_bam::reference_length_from_raw_bam(rec)) - 1; // 5' position: forward=start, reverse=end - let first_5prime = if (f1 & fgumi_raw_bam::flags::REVERSE) != 0 { end1 } else { pos1 }; - let second_5prime = if (f2 & fgumi_raw_bam::flags::REVERSE) != 0 { end2 } else { pos2 }; + let five_prime = if (flags & fgumi_raw_bam::flags::REVERSE) != 0 { end } else { start }; + + InsertSizeEnd::new( + fgumi_raw_bam::ref_id(rec), + (flags & fgumi_raw_bam::flags::UNMAPPED) != 0, + five_prime, + ) +} + +/// Computes the insert size (TLEN) `first` should carry relative to `second`. +/// +/// Returns 0 when either end is unmapped or the two are on different references, matching +/// htsjdk `SamPairUtil.computeInsertSize`. For a primary pair the other end carries the negation +/// of this value; supplementary alignments instead take their own value from a second call with +/// the supplementary as `first`. +/// +/// This is the single copy of the formula: `commands::clip::compute_insert_size_raw` delegates +/// here rather than repeating it. See [`InsertSizeEnd`]. +pub(crate) fn compute_insert_size_from_ends(first: InsertSizeEnd, second: InsertSizeEnd) -> i32 { + if first.is_unmapped || second.is_unmapped || first.ref_id != second.ref_id { + return 0; + } + + // Widen to i64 so a malformed record with an extreme POS cannot overflow the subtraction. + let (first_5p, second_5p) = (i64::from(first.five_prime), i64::from(second.five_prime)); + let adjustment = if second_5p >= first_5p { 1 } else { -1 }; + i32::try_from(second_5p - first_5p + adjustment).unwrap_or(0) +} - let adjustment = if second_5prime >= first_5prime { 1 } else { -1 }; - second_5prime - first_5prime + adjustment +/// Computes insert size (TLEN) from two raw BAM records. +fn compute_insert_size_raw(rec1: &[u8], rec2: &[u8]) -> i32 { + compute_insert_size_from_ends(insert_size_end(rec1), insert_size_end(rec2)) } /// Determines the pair orientation for a paired read using raw BAM bytes. @@ -1006,6 +1054,7 @@ impl MemoryEstimate for Template { mod tests { use super::*; use fgumi_raw_bam::{RawRecord, SamBuilder as RawSamBuilder, flags as raw_flags}; + use rstest::rstest; // SAM flag constants const FLAG_PAIRED: u16 = 0x1; @@ -1547,7 +1596,7 @@ mod tests { } /// Tests that `fix_mate_info` sets TLEN on R1 supplementary alignments. - /// TLEN should be set to negative of mate primary's (R2) TLEN. + /// TLEN is computed from the supplementary's own alignment against the mate primary. #[test] fn test_fix_mate_info_sets_tlen_on_r1_supplementals() -> Result<()> { let r1 = create_mapped_record_with_tlen(b"read1", FLAG_PAIRED | FLAG_READ1, 100, 30, 200); @@ -1563,19 +1612,19 @@ mod tests { let mut template = Template::from_records(vec![r1, r2, r1_supp])?; template.fix_mate_info()?; - // R1(pos=100,forward) and R2(pos=200,forward) → insert_size = 101 - // R1 supplementary TLEN = -(-101) = 101 + // The supplementary (pos=300, forward) lies right of its mate primary R2 (pos=200), so it + // is the rightmost segment: TLEN = 200 - 300 - 1 = -101. assert_eq!( template.records()[2].template_length(), - 101, - "R1 supplementary TLEN should be negative of R2's TLEN" + -101, + "R1 supplementary TLEN is computed against the mate primary, not copied from it" ); Ok(()) } /// Tests that `fix_mate_info` sets TLEN on R2 supplementary alignments. - /// TLEN should be set to negative of mate primary's (R1) TLEN. + /// TLEN is computed from the supplementary's own alignment against the mate primary. #[test] fn test_fix_mate_info_sets_tlen_on_r2_supplementals() -> Result<()> { let r1 = create_mapped_record_with_tlen(b"read1", FLAG_PAIRED | FLAG_READ1, 100, 30, 300); @@ -1591,12 +1640,12 @@ mod tests { let mut template = Template::from_records(vec![r1, r2, r2_supp])?; template.fix_mate_info()?; - // R1(pos=100,forward) and R2(pos=200,forward) → insert_size = 101 - // R2 supplementary TLEN = -(101) = -101 + // The supplementary (pos=400, forward) lies right of its mate primary R1 (pos=100), so it + // is the rightmost segment: TLEN = 100 - 400 - 1 = -301. assert_eq!( template.records()[2].template_length(), - -101, - "R2 supplementary TLEN should be negative of R1's TLEN" + -301, + "R2 supplementary TLEN is computed against the mate primary, not copied from it" ); Ok(()) @@ -1641,18 +1690,109 @@ mod tests { Template::from_records(vec![r1, r2, r1_supp1, r1_supp2, r2_supp1, r2_supp2])?; template.fix_mate_info()?; - // R1(pos=100,forward) and R2(pos=200,forward) → insert_size = 101 - // R1 supplementaries should have TLEN = -(-101) = 101 - assert_eq!(template.records()[2].template_length(), 101, "R1 supp1 TLEN"); - assert_eq!(template.records()[3].template_length(), 101, "R1 supp2 TLEN"); + // Each supplementary gets its own TLEN, computed against its mate primary — so unlike the + // old copy-the-mate's-TLEN behaviour, the two R1 supplementaries no longer share a value. + // R1 supplementaries are measured against primary R2 (pos=200); order is reversed. + assert_eq!(template.records()[2].template_length(), -201, "R1 supp@400 TLEN"); + assert_eq!(template.records()[3].template_length(), -101, "R1 supp@300 TLEN"); + + // R2 supplementaries are measured against primary R1 (pos=100); order is reversed. + assert_eq!(template.records()[4].template_length(), -501, "R2 supp@600 TLEN"); + assert_eq!(template.records()[5].template_length(), -401, "R2 supp@500 TLEN"); + + Ok(()) + } + + /// Builds a mapped 100M record on an explicit reference at a 1-based position. + fn create_mapped_record_on_ref(name: &[u8], flags: u16, ref_id: i32, pos: usize) -> RawRecord { + let mut rec = create_mapped_record_with_tlen(name, flags, pos, 60, 0); + fgumi_raw_bam::set_ref_id(&mut rec, ref_id); + rec + } + + /// A supplementary's TLEN is computed from its own alignment against the mate primary rather + /// than copied from the mate primary's TLEN. See issue #673 and samtools/htsjdk#1795. + /// + /// R1 and R2 primaries are forward 100M at 1-based 1000 and 1350; the R1 supplementary is + /// forward 100M and varies by reference and position. + #[rstest] + // Supplementary right of the mate primary: it is the rightmost segment, so TLEN is negative. + #[case::beyond_mate(0, 5_000, -3_651)] + // Supplementary left of the mate primary: it is the leftmost segment, so TLEN is positive. + #[case::before_mate(0, 500, 851)] + // Coincident 5' ends: the two ends must differ in sign, so the +1 adjustment applies. + #[case::coincident_five_prime(0, 1_350, 1)] + // Different reference: the information is unavailable, so TLEN is 0 rather than the mate's. + #[case::cross_reference(1, 5_000, 0)] + fn test_supplemental_tlen_is_computed_not_copied( + #[case] supp_ref_id: i32, + #[case] supp_pos: usize, + #[case] expected_tlen: i32, + ) -> Result<()> { + let r1 = create_mapped_record_on_ref(b"read1", FLAG_PAIRED | FLAG_READ1, 0, 1_000); + let r2 = create_mapped_record_on_ref(b"read1", FLAG_PAIRED | FLAG_READ2, 0, 1_350); + let r1_supp = create_mapped_record_on_ref( + b"read1", + FLAG_PAIRED | FLAG_READ1 | FLAG_SUPPLEMENTARY, + supp_ref_id, + supp_pos, + ); + + let mut template = Template::from_records(vec![r1, r2, r1_supp])?; + template.fix_mate_info()?; + + assert_eq!(template.records()[2].template_length(), expected_tlen); - // R2 supplementaries should have TLEN = -(101) = -101 - assert_eq!(template.records()[4].template_length(), -101, "R2 supp1 TLEN"); - assert_eq!(template.records()[5].template_length(), -101, "R2 supp2 TLEN"); + // The primaries are unaffected by where the supplementary landed. + assert_eq!(template.records()[0].template_length(), 351, "R1 primary TLEN"); + assert_eq!(template.records()[1].template_length(), -351, "R2 primary TLEN"); Ok(()) } + /// `compute_insert_size_from_ends` is the one copy of the htsjdk formula — `zipper`'s + /// supplementary fix-up and `commands::clip::compute_insert_size_raw` both route through + /// it — so its contract is pinned directly against htsjdk + /// `SamPairUtil.computeInsertSize`: 0 if either end is unmapped or the ends are on + /// different references, otherwise `second5p - first5p + (second5p >= first5p ? 1 : -1)`. + #[rstest] + // Second end to the right: positive, with the +1 adjustment. + #[case::second_end_right(0, false, 1_000, 0, false, 1_350, 351)] + // Second end to the left: negative, with the -1 adjustment. + #[case::second_end_left(0, false, 1_350, 0, false, 1_000, -351)] + // Coincident 5' ends: `>=` selects +1, so the two ends of a pair still differ in sign. + #[case::coincident_five_prime(0, false, 1_000, 0, false, 1_000, 1)] + // One base apart in each direction, pinning the sign flip at the boundary. + #[case::one_base_right(0, false, 1_000, 0, false, 1_001, 2)] + #[case::one_base_left(0, false, 1_000, 0, false, 999, -2)] + // Either end unmapped, or the ends on different references: no insert size. + #[case::first_unmapped(0, true, 1_000, 0, false, 1_350, 0)] + #[case::second_unmapped(0, false, 1_000, 0, true, 1_350, 0)] + #[case::cross_reference(0, false, 1_000, 1, false, 1_350, 0)] + fn test_compute_insert_size_from_ends_matches_htsjdk( + #[case] first_ref_id: i32, + #[case] first_unmapped: bool, + #[case] first_five_prime: i64, + #[case] second_ref_id: i32, + #[case] second_unmapped: bool, + #[case] second_five_prime: i64, + #[case] expected_tlen: i32, + ) { + let first = InsertSizeEnd::new(first_ref_id, first_unmapped, first_five_prime); + let second = InsertSizeEnd::new(second_ref_id, second_unmapped, second_five_prime); + assert_eq!(compute_insert_size_from_ends(first, second), expected_tlen); + } + + /// A 5' position past `i32::MAX` is only reachable from a malformed record. + /// `InsertSizeEnd::new` saturates it so the subtraction cannot overflow; the pair of + /// saturated ends then reads as coincident rather than wrapping to a bogus TLEN. + #[test] + fn test_insert_size_end_saturates_an_out_of_range_five_prime() { + let huge = InsertSizeEnd::new(0, false, i64::from(i32::MAX) + 10); + let at_max = InsertSizeEnd::new(0, false, i64::from(i32::MAX)); + assert_eq!(compute_insert_size_from_ends(huge, at_max), 1); + } + /// Tests that `fix_mate_info` recalculates TLEN for supplementaries based on primary positions #[test] fn test_fix_mate_info_tlen_recalculated() -> Result<()> { @@ -1669,10 +1809,11 @@ mod tests { let mut template = Template::from_records(vec![r1, r2, r1_supp])?; template.fix_mate_info()?; - // R1(pos=100,forward) and R2(pos=200,forward) → insert_size = 101 + // The stale TLEN of 999 is discarded and recomputed: the supplementary (pos=300) lies + // right of its mate primary R2 (pos=200), giving 200 - 300 - 1 = -101. assert_eq!( template.records()[2].template_length(), - 101, + -101, "R1 supplementary TLEN should be recalculated from positions" ); @@ -1868,8 +2009,9 @@ mod tests { assert!(supp.is_mate_reverse(), "R1 supp should have mate_reverse since R2 is reverse"); // Mate unmapped should NOT be set assert!(!supp.is_mate_unmapped(), "R1 supp should NOT have mate_unmapped"); - // TLEN should be -(-200) = 200 - assert_eq!(supp.template_length(), 200, "R1 supp TLEN should be -(-200) = 200"); + // TLEN is computed against R2 (200..299 reverse, 5' = 299), not copied from it: the + // supplementary (500..599 forward, 5' = 500) is rightmost, so 299 - 500 - 1 = -202. + assert_eq!(supp.template_length(), -202, "R1 supp TLEN computed against R2"); // MQ tag should be R2's mapq (40) assert_eq!(supp.tags().find_int(SamTag::MQ), Some(40), "R1 supp MQ should be R2's mapq"); // MC tag should be present @@ -1936,9 +2078,9 @@ mod tests { !supp.is_mate_reverse(), "R2 supp should NOT have mate_reverse since R1 is forward" ); - // TLEN = -(101) = -101 - // R1(pos=100,forward) and R2(pos=200,forward) → insert_size = 101 - assert_eq!(supp.template_length(), -101, "R2 supp TLEN should be -101"); + // TLEN is computed against R1 (forward, 5' = 100), not copied from it: the supplementary + // (600..699 forward, 5' = 600) is rightmost, so 100 - 600 - 1 = -501. + assert_eq!(supp.template_length(), -501, "R2 supp TLEN computed against R1"); // MQ tag should be R1's mapq (35) assert_eq!(supp.tags().find_int(SamTag::MQ), Some(35), "R2 supp MQ should be R1's mapq"); // MC tag should be present @@ -2106,14 +2248,18 @@ mod tests { let mut template = Template::from_records(vec![r1, r2, r1_supp1, r1_supp2])?; template.fix_mate_info()?; - // R1(pos=100,forward) and R2(pos=300,100M,reverse) → 5' are 100 and 399 → insert_size = 300 + // Mate info is shared by both supplementaries; TLEN is not, since each is computed from + // its own alignment against R2 (300..399 reverse, 5' = 399). for i in 2..=3 { let supp = &template.records()[i]; assert_eq!(supp.mate_pos() + 1, 300, "Supp {i} should have mate pos 300"); assert!(supp.is_mate_reverse(), "Supp {i} should have mate_reverse"); - assert_eq!(supp.template_length(), 300, "Supp {i} TLEN should be 300"); } + // Supplementaries are stored in reverse order: [2] is the one at 700, [3] at 500. + assert_eq!(template.records()[2].template_length(), -302, "supp@700: 399 - 700 - 1"); + assert_eq!(template.records()[3].template_length(), -102, "supp@500: 399 - 500 - 1"); + Ok(()) } diff --git a/tests/integration/test_clip_command.rs b/tests/integration/test_clip_command.rs index c92eb0b73..f85df7aac 100644 --- a/tests/integration/test_clip_command.rs +++ b/tests/integration/test_clip_command.rs @@ -866,10 +866,17 @@ fn test_clip_command_threads_mode_supplementary_mate_repair() { ); assert_eq!(mate_cigar(supp).as_deref(), Some("98M2H"), "MC = primary R2 post-clip CIGAR"); assert_eq!(mate_mapq(supp), Some(40), "MQ = primary R2 MAPQ"); + // TLEN is computed from the supplementary's own alignment against the post-clip primary R2, + // not copied from R2's TLEN (see issue #673). Hand-derived from the inputs: + // R1 primary 100..199 forward -> 5' = 100 + // R2 primary 300..397 reverse -> 5' = 397 => primaries are +298 / -298 + // R1 supp 5000..5049 forward -> 5' = 5000 + // The supplementary is the rightmost segment, so its TLEN is negative: 397 - 5000 - 1. + assert_eq!(r2.template_length(), -298, "primary R2 TLEN"); assert_eq!( supp.template_length(), - -r2.template_length(), - "supplementary TLEN = negation of primary R2 TLEN" + -4604, + "supplementary TLEN is computed against primary R2, not copied from it" ); // Cross-check: the repaired mate fields agree with the actual primary R2 record, proving the