diff --git a/crates/fgumi-sam/src/clipper.rs b/crates/fgumi-sam/src/clipper.rs index 4e2ba6aa3..c7dcce7b4 100644 --- a/crates/fgumi-sam/src/clipper.rs +++ b/crates/fgumi-sam/src/clipper.rs @@ -871,11 +871,13 @@ impl SamRecordClipper { return (0, 0); } - // Get unclipped positions for both reads - let r1_unclipped_start = Self::unclipped_start(r1); - let r1_unclipped_end = Self::unclipped_end(r1); - let r2_unclipped_start = Self::unclipped_start(r2); - let r2_unclipped_end = Self::unclipped_end(r2); + // Bound each read's window by the mate's *un-soft-clipped* span (soft clips only; + // hard-clipped bases are physically absent), matching fgbio's use of + // `mate.unSoftClippedStart` / `unSoftClippedEnd` in `clipExtendingPastMateEnds`. + let r1_unclipped_start = record_utils::unsoftclipped_start(r1); + let r1_unclipped_end = record_utils::unsoftclipped_end(r1); + let r2_unclipped_start = record_utils::unsoftclipped_start(r2); + let r2_unclipped_end = record_utils::unsoftclipped_end(r2); let (Some(r1_start), Some(r1_end), Some(r2_start), Some(r2_end)) = (r1_unclipped_start, r1_unclipped_end, r2_unclipped_start, r2_unclipped_end) @@ -923,19 +925,6 @@ impl SamRecordClipper { } } - /// Get unclipped start position (delegates to `record_utils`) - // RecordBuf kept: thin wrapper over record_utils::unclipped_start which uses noodles - // typed flags and Position; serves RecordBuf callers within this module. - fn unclipped_start(rec: &RecordBuf) -> Option { - record_utils::unclipped_start(rec) - } - - /// Get unclipped end position (delegates to `record_utils`) - // RecordBuf kept: same as unclipped_start; delegates to record_utils::unclipped_end. - fn unclipped_end(rec: &RecordBuf) -> Option { - record_utils::unclipped_end(rec) - } - /// Count leading soft clips // RecordBuf kept: accepts noodles CigarBuf (typed CIGAR) — converts ops to (Kind, usize) // pairs and delegates to record_utils::leading_soft_clipping. @@ -2108,10 +2097,12 @@ impl RawRecordClipper { return (0, 0); } - let r1_unclipped_start = crate::record_utils::unclipped_start_raw(r1.as_ref()); - let r1_unclipped_end = crate::record_utils::unclipped_end_raw(r1.as_ref()); - let r2_unclipped_start = crate::record_utils::unclipped_start_raw(r2.as_ref()); - let r2_unclipped_end = crate::record_utils::unclipped_end_raw(r2.as_ref()); + // Soft-only mate window (see the RecordBuf sibling above and fgbio + // `SamRecordClipper.clipExtendingPastMateEnds`); hard-clipped bases are absent. + let r1_unclipped_start = crate::record_utils::unsoftclipped_start_raw(r1.as_ref()); + let r1_unclipped_end = crate::record_utils::unsoftclipped_end_raw(r1.as_ref()); + let r2_unclipped_start = crate::record_utils::unsoftclipped_start_raw(r2.as_ref()); + let r2_unclipped_end = crate::record_utils::unsoftclipped_end_raw(r2.as_ref()); let (Some(r1_start), Some(r1_end), Some(r2_start), Some(r2_end)) = (r1_unclipped_start, r1_unclipped_end, r2_unclipped_start, r2_unclipped_end) diff --git a/crates/fgumi-sam/src/record_utils.rs b/crates/fgumi-sam/src/record_utils.rs index c01610c4f..578a7dcd0 100644 --- a/crates/fgumi-sam/src/record_utils.rs +++ b/crates/fgumi-sam/src/record_utils.rs @@ -494,6 +494,49 @@ pub fn unclipped_end(record: &RecordBuf) -> Option { Some(start + ref_len.saturating_sub(1) + trailing) } +/// Gets the un-*soft*-clipped start position of a read (alignment start minus leading +/// **soft** clips only; hard clips are ignored because the hard-clipped bases are +/// physically absent from the read). +/// +/// This matches fgbio's `SamRecord.unSoftClippedStart`, used when bounding a mate's +/// window for `clipExtendingPastMate` (see fgbio `SamRecordClipper`). Contrast +/// [`unclipped_start`], which includes hard clips (htsjdk `getUnclippedStart`). +/// +/// Returns `None` for unmapped reads. +#[must_use] +pub fn unsoftclipped_start(record: &RecordBuf) -> Option { + if record.flags().is_unmapped() { + return None; + } + let start = usize::from(record.alignment_start()?); + let leading = leading_soft_clipping(&cigar_to_ops(record)); + Some(start.saturating_sub(leading)) +} + +/// Gets the un-*soft*-clipped end position of a read (alignment end plus trailing +/// **soft** clips only; hard clips are ignored because the hard-clipped bases are +/// physically absent from the read). +/// +/// This matches fgbio's `SamRecord.unSoftClippedEnd`. Contrast [`unclipped_end`], +/// which includes hard clips (htsjdk `getUnclippedEnd`). +/// +/// Returns `None` for unmapped reads or records with no CIGAR ops (matching the raw +/// sibling [`unsoftclipped_end_raw`]). +#[must_use] +pub fn unsoftclipped_end(record: &RecordBuf) -> Option { + if record.flags().is_unmapped() { + return None; + } + let start = usize::from(record.alignment_start()?); + let ops = cigar_to_ops(record); + if ops.is_empty() { + return None; + } + let ref_len = reference_length(&record.cigar()); + let trailing = trailing_soft_clipping(&ops); + Some(start + ref_len.saturating_sub(1) + trailing) +} + /// Gets the unclipped 5' position of a read. /// /// For forward strand reads, returns the unclipped start position. @@ -771,19 +814,46 @@ pub fn unclipped_end_raw(bam: &[u8]) -> Option { if cigar_ops.is_empty() { return None; } - let ref_len: usize = cigar_ops - .iter() - .map(|&op| { - let t = op & 0xF; - let len = (op >> 4) as usize; - if consumes_ref(t) { len } else { 0 } - }) - .sum(); + let ref_len = cigar_reference_length_raw(&cigar_ops); let trailing = trailing_clipping_raw(&cigar_ops); // alignment_end = start + ref_len - 1; unclipped_end = alignment_end + trailing Some(start + ref_len.saturating_sub(1) + trailing) } +/// Raw-byte sibling of [`unsoftclipped_start`]: alignment start minus leading **soft** +/// clips only (hard clips ignored). Matches fgbio `SamRecord.unSoftClippedStart`. +/// +/// Returns `None` for unmapped reads. +#[must_use] +pub fn unsoftclipped_start_raw(bam: &[u8]) -> Option { + if RawRecordView::new(bam).flags() & flags::UNMAPPED != 0 { + return None; + } + let start = alignment_start_from_raw(bam)?; + let cigar_ops = get_cigar_ops(bam); + let leading = leading_soft_clipping_raw(&cigar_ops); + Some(start.saturating_sub(leading)) +} + +/// Raw-byte sibling of [`unsoftclipped_end`]: alignment end plus trailing **soft** +/// clips only (hard clips ignored). Matches fgbio `SamRecord.unSoftClippedEnd`. +/// +/// Returns `None` for unmapped reads or records with no CIGAR ops. +#[must_use] +pub fn unsoftclipped_end_raw(bam: &[u8]) -> Option { + if RawRecordView::new(bam).flags() & flags::UNMAPPED != 0 { + return None; + } + let start = alignment_start_from_raw(bam)?; + let cigar_ops = get_cigar_ops(bam); + if cigar_ops.is_empty() { + return None; + } + let ref_len = cigar_reference_length_raw(&cigar_ops); + let trailing = trailing_soft_clipping_raw(&cigar_ops); + Some(start + ref_len.saturating_sub(1) + trailing) +} + /// Returns the reference length consumed by raw u32 CIGAR operations. /// /// This is the raw-byte sibling of [`cigar_reference_length`]. @@ -1637,6 +1707,51 @@ mod tests { assert_eq!(unclipped_end(&read), Some(149)); } + // R2-CLIP-03: un-*soft*-clipped positions ignore hard clips (fgbio `unSoftClipped*`), + // unlike `unclipped_*` (htsjdk `getUnclipped*`, which counts hard clips too). + #[test] + fn test_unsoftclipped_start_ignores_leading_hard_clip() { + // 5H10S35M at position 100: unclipped_start = 100 - 5 - 10 = 85 (counts hard); + // unsoftclipped_start = 100 - 10 = 90 (soft only). + let read = create_cigar_test_read("hs", 100, "5H10S35M"); + assert_eq!(unclipped_start(&read), Some(85)); + assert_eq!(unsoftclipped_start(&read), Some(90)); + } + + #[test] + fn test_unsoftclipped_end_ignores_trailing_hard_clip() { + // 35M10S5H at position 100: alignment_end = 134; + // unclipped_end = 134 + 10 + 5 = 149 (counts hard); + // unsoftclipped_end = 134 + 10 = 144 (soft only). + let read = create_cigar_test_read("sh", 100, "35M10S5H"); + assert_eq!(unclipped_end(&read), Some(149)); + assert_eq!(unsoftclipped_end(&read), Some(144)); + } + + #[test] + fn test_unsoftclipped_no_hard_matches_unclipped() { + // With no hard clips, soft-only equals the full unclipped positions. + let read = create_cigar_test_read("soft", 100, "10S30M10S"); + assert_eq!(unsoftclipped_start(&read), unclipped_start(&read)); + assert_eq!(unsoftclipped_end(&read), unclipped_end(&read)); + } + + #[test] + fn test_unsoftclipped_unmapped_returns_none() { + let read = RecordBuilder::new().name("u").sequence("ACGT").flags(Flags::UNMAPPED).build(); + assert_eq!(unsoftclipped_start(&read), None); + assert_eq!(unsoftclipped_end(&read), None); + } + + #[test] + fn test_unsoftclipped_end_empty_cigar_returns_none() { + // A mapped record with no CIGAR ops has no end position: the typed function must + // return None, matching the raw sibling `unsoftclipped_end_raw` (not `Some(start)`). + let read = create_cigar_test_read("empty", 100, ""); + assert!(read.cigar().as_ref().is_empty(), "test setup: CIGAR must be empty"); + assert_eq!(unsoftclipped_end(&read), None); + } + #[test] fn test_unclipped_five_prime_forward_strand() { // Forward strand: 5' is at unclipped_start @@ -1885,6 +2000,33 @@ mod tests { assert_eq!(unclipped_end_raw(&raw), None); } + // R2-CLIP-03 (raw): soft-only unclipped positions ignore hard clips. + #[test] + fn test_unsoftclipped_start_raw_ignores_hard_clip() { + let record = create_cigar_test_read("hs", 100, "5H10S35M"); + let raw = to_raw(&record); + assert_eq!(unsoftclipped_start_raw(&raw), unsoftclipped_start(&record)); + assert_eq!(unsoftclipped_start_raw(&raw), Some(90)); // 100 - 10S (hard ignored) + assert_eq!(unclipped_start_raw(&raw), Some(85)); // 100 - 5H - 10S + } + + #[test] + fn test_unsoftclipped_end_raw_ignores_hard_clip() { + let record = create_cigar_test_read("sh", 100, "35M10S5H"); + let raw = to_raw(&record); + assert_eq!(unsoftclipped_end_raw(&raw), unsoftclipped_end(&record)); + assert_eq!(unsoftclipped_end_raw(&raw), Some(144)); // 134 + 10S (hard ignored) + assert_eq!(unclipped_end_raw(&raw), Some(149)); // 134 + 10S + 5H + } + + #[test] + fn test_unsoftclipped_raw_unmapped_returns_none() { + let record = RecordBuilder::new().name("u").sequence("ACGT").flags(Flags::UNMAPPED).build(); + let raw = to_raw(&record); + assert_eq!(unsoftclipped_start_raw(&raw), None); + assert_eq!(unsoftclipped_end_raw(&raw), None); + } + // ========================================================================= // Tests for mate_unclipped_start_raw / mate_unclipped_end_raw // ========================================================================= diff --git a/src/lib/commands/clip.rs b/src/lib/commands/clip.rs index 85d655701..c9b337ab4 100644 --- a/src/lib/commands/clip.rs +++ b/src/lib/commands/clip.rs @@ -58,8 +58,8 @@ done in streaming fashion with, for example: fgumi sort -i in.bam --order queryname | fgumi clip -i /dev/stdin ... -The output sort order may be specified with --sort-order. If not given, then the output will be in the same -order as input. +The output is written in the same order as the input. To produce coordinate-sorted output, pipe the +result through a separate `fgumi sort`. Any existing NM, UQ and MD tags are repaired, and mate-pair information is updated. @@ -87,14 +87,6 @@ pub struct Clip { #[arg(short = 'c', long = "clipping-mode", default_value_t = ClippingMode::Hard)] pub clipping_mode: ClippingMode, - /// Output sort order (if not specified, output is in same order as input) - #[arg( - short = 'S', - long = "sort-order", - value_parser = ["unknown", "unsorted", "queryname", "coordinate"] - )] - pub sort_order: Option, - /// Clip overlapping read pairs #[arg(long = "clip-overlapping-reads", default_value = "false", num_args = 0..=1, default_missing_value = "true", action = clap::ArgAction::Set, value_parser = parse_bool)] pub clip_overlapping_reads: bool, @@ -248,9 +240,6 @@ impl Command for Clip { // Load reference (always required for clip) let reference = Arc::new(ReferenceReader::new(&self.reference)?); - // Update header sort order if specified - let header = self.update_header_sort_order(header)?; - // Add @PG record with PP chaining to input's last program let header = crate::commands::common::add_pg_record(header, command_line)?; @@ -286,9 +275,6 @@ impl Clip { let reader_threads = self.threading.num_threads(); let (reader, header) = create_raw_bam_reader(&self.io.input, reader_threads)?; - // Update header sort order if specified - let header = self.update_header_sort_order(header)?; - // Add @PG record with PP chaining to input's last program let header = crate::commands::common::add_pg_record(header, command_line)?; @@ -316,31 +302,34 @@ impl Clip { let template = template?; let mut records: Vec = template.into_records(); - // Process based on template type - #[allow(clippy::len_zero)] // We specifically want len() == 1, not !is_empty() - if records.len() == 1 { - // Fragment - let record = &mut records[0]; - self.clip_fragment(&clipper, record, metrics_collection.as_mut())?; - } else if records.len() == 2 { - // Paired reads - let (r1, r2) = records.split_at_mut(1); - let r1 = &mut r1[0]; - let r2 = &mut r2[0]; - - let (overlap_clip, extend_clip) = - self.clip_pair(&clipper, r1, r2, metrics_collection.as_mut())?; - - if overlap_clip { - total_clipped_overlap += 1; + // Clip the primary R1/R2 (or a lone fragment) — never positional records[0]/[1] — + // so templates with secondary/supplementary reads (len > 2) are clipped too, matching + // fgbio ClipBam's Template.r1/r2 handling. + match find_primary_pair_indices(&records)? { + (Some(i1), Some(i2)) => { + let [r1, r2] = + records.get_disjoint_mut([i1, i2]).expect("distinct primary indices"); + + let (overlap_clip, extend_clip) = + self.clip_pair(&clipper, r1, r2, metrics_collection.as_mut())?; + + if overlap_clip { + total_clipped_overlap += 1; + } + if extend_clip { + total_clipped_mate_extension += 1; + } + + // Fix full mate-pair info (mate coords/strand, mate-unmapped flag, MQ/MC, + // TLEN), matching fgbio ClipBam's SamPairUtil.setMateInfo call, then repair + // mate info on any supplementary alignments in the template. + set_mate_info_raw(r1, r2); + fix_supplemental_mate_info(&mut records, i1, i2); } - if extend_clip { - total_clipped_mate_extension += 1; + (Some(i1), None) => { + self.clip_fragment(&clipper, &mut records[i1], metrics_collection.as_mut())?; } - - // Fix full mate-pair info (mate coords/strand, mate-unmapped flag, MQ/MC, - // TLEN), matching fgbio ClipBam's SamPairUtil.setMateInfo call. - set_mate_info_raw(r1, r2); + _ => {} } // Regenerate alignment tags (always done to match Scala fgbio behavior) @@ -575,49 +564,6 @@ impl Clip { Ok((overlap_clipped, extend_clipped)) } - /// Updates the header sort order if specified via command line option. - fn update_header_sort_order(&self, header: Header) -> Result
{ - if let Some(ref sort_order) = self.sort_order { - use bstr::BString; - use noodles::sam::header::record::value::Map; - use noodles::sam::header::record::value::map::header::tag::SORT_ORDER; - - // Get or create the header map - let mut header_map = if let Some(hd) = header.header() { - hd.clone() - } else { - Map::::default() - }; - - // Update sort order - *header_map.other_fields_mut().entry(SORT_ORDER).or_insert(BString::from("")) = - BString::from(sort_order.as_str()); - - // Rebuild header with new header map - let mut builder = noodles::sam::Header::builder(); - - // Copy existing components - for (name, rg) in header.read_groups() { - builder = builder.add_read_group(name.clone(), rg.clone()); - } - for (name, reference) in header.reference_sequences() { - builder = builder.add_reference_sequence(name.clone(), reference.clone()); - } - for (id, pg) in header.programs().as_ref() { - builder = builder.add_program(id.clone(), pg.clone()); - } - for comment in header.comments() { - builder = builder.add_comment(comment.clone()); - } - - // Set the modified header - builder = builder.set_header(header_map); - Ok(builder.build()) - } else { - Ok(header) - } - } - /// Execute using the 7-step unified pipeline with multi-threading. /// /// Uses `TemplateGrouper` to batch records by template (QNAME) for parallel processing. @@ -686,77 +632,80 @@ impl Clip { templates_count += 1; let mut records: Vec = template.into_records(); - #[allow(clippy::len_zero)] - if records.len() == 1 { - // Fragment - apply fixed-position clipping - let record = &mut records[0]; - if upgrade_clipping { - clipper.upgrade_all_clipping_raw(record).map_err(io::Error::other)?; - } - if read_one_five_prime > 0 { - clipper.clip_5_prime_end_of_alignment(record, read_one_five_prime); - } - if read_one_three_prime > 0 { - clipper.clip_3_prime_end_of_alignment(record, read_one_three_prime); - } - } else if records.len() == 2 { - // Paired reads - let (r1_slice, r2_slice) = records.split_at_mut(1); - let r1 = &mut r1_slice[0]; - let r2 = &mut r2_slice[0]; - - if upgrade_clipping { - clipper.upgrade_all_clipping_raw(r1).map_err(io::Error::other)?; - clipper.upgrade_all_clipping_raw(r2).map_err(io::Error::other)?; - } + // Clip the primary R1/R2 (or a lone fragment) by flags, never positional + // records[0]/[1], so templates with secondary/supplementary reads are clipped + // too (matches fgbio ClipBam's Template.r1/r2 handling). + match find_primary_pair_indices(&records).map_err(io::Error::other)? { + (Some(i1), Some(i2)) => { + let [r1, r2] = + records.get_disjoint_mut([i1, i2]).expect("distinct primary indices"); + + if upgrade_clipping { + clipper.upgrade_all_clipping_raw(r1).map_err(io::Error::other)?; + clipper.upgrade_all_clipping_raw(r2).map_err(io::Error::other)?; + } - // Determine read types (raw flags) - let is_r1_first = r1.is_first_segment(); - let is_r2_last = r2.is_last_segment(); + // Determine read types (raw flags) + let is_r1_first = r1.is_first_segment(); - // Apply fixed-position clipping for R1 - if is_r1_first && read_one_five_prime > 0 { - clipper.clip_5_prime_end_of_alignment(r1, read_one_five_prime); - } else if !is_r1_first && read_two_five_prime > 0 { - clipper.clip_5_prime_end_of_alignment(r1, read_two_five_prime); - } - if is_r1_first && read_one_three_prime > 0 { - clipper.clip_3_prime_end_of_alignment(r1, read_one_three_prime); - } else if !is_r1_first && read_two_three_prime > 0 { - clipper.clip_3_prime_end_of_alignment(r1, read_two_three_prime); - } + // Apply fixed-position clipping for R1 + if is_r1_first && read_one_five_prime > 0 { + clipper.clip_5_prime_end_of_alignment(r1, read_one_five_prime); + } else if !is_r1_first && read_two_five_prime > 0 { + clipper.clip_5_prime_end_of_alignment(r1, read_two_five_prime); + } + if is_r1_first && read_one_three_prime > 0 { + clipper.clip_3_prime_end_of_alignment(r1, read_one_three_prime); + } else if !is_r1_first && read_two_three_prime > 0 { + clipper.clip_3_prime_end_of_alignment(r1, read_two_three_prime); + } - // Apply fixed-position clipping for R2 - if is_r2_last && read_two_five_prime > 0 { - clipper.clip_5_prime_end_of_alignment(r2, read_two_five_prime); - } else if !is_r2_last && read_one_five_prime > 0 { - clipper.clip_5_prime_end_of_alignment(r2, read_one_five_prime); - } - if is_r2_last && read_two_three_prime > 0 { - clipper.clip_3_prime_end_of_alignment(r2, read_two_three_prime); - } else if !is_r2_last && read_one_three_prime > 0 { - clipper.clip_3_prime_end_of_alignment(r2, read_one_three_prime); - } + // Apply fixed-position clipping for R2. The primary R2 slot is always a + // last-segment read (find_primary_pair_indices only fills it from an + // is_last_segment read), so it always uses the read-two thresholds. + if read_two_five_prime > 0 { + clipper.clip_5_prime_end_of_alignment(r2, read_two_five_prime); + } + if read_two_three_prime > 0 { + clipper.clip_3_prime_end_of_alignment(r2, read_two_three_prime); + } - // Clip overlapping reads - if clip_overlapping_reads { - let (num_r1, num_r2) = clipper.clip_overlapping_reads(r1, r2); - if num_r1 > 0 || num_r2 > 0 { - overlap_clipped_count += 1; + // Clip overlapping reads + if clip_overlapping_reads { + let (num_r1, num_r2) = clipper.clip_overlapping_reads(r1, r2); + if num_r1 > 0 || num_r2 > 0 { + overlap_clipped_count += 1; + } } - } - // Clip reads extending past mate - if clip_extending_past_mate { - let (num_r1, num_r2) = clipper.clip_extending_past_mate_ends(r1, r2); - if num_r1 > 0 || num_r2 > 0 { - extend_clipped_count += 1; + // Clip reads extending past mate + if clip_extending_past_mate { + let (num_r1, num_r2) = clipper.clip_extending_past_mate_ends(r1, r2); + if num_r1 > 0 || num_r2 > 0 { + extend_clipped_count += 1; + } } - } - // Fix full mate-pair info (mate coords/strand, mate-unmapped flag, - // MQ/MC, TLEN), matching fgbio ClipBam's SamPairUtil.setMateInfo call. - set_mate_info_raw(r1, r2); + // Fix full mate-pair info (mate coords/strand, mate-unmapped flag, + // MQ/MC, TLEN), matching fgbio ClipBam's SamPairUtil.setMateInfo call, + // then repair mate info on any supplementary alignments. + set_mate_info_raw(r1, r2); + fix_supplemental_mate_info(&mut records, i1, i2); + } + (Some(i1), None) => { + // Fragment - apply fixed-position clipping (R1 thresholds) + let record = &mut records[i1]; + if upgrade_clipping { + clipper.upgrade_all_clipping_raw(record).map_err(io::Error::other)?; + } + if read_one_five_prime > 0 { + clipper.clip_5_prime_end_of_alignment(record, read_one_five_prime); + } + if read_one_three_prime > 0 { + clipper.clip_3_prime_end_of_alignment(record, read_one_three_prime); + } + } + _ => {} } // Regenerate alignment tags (always done to match fgbio behavior) @@ -984,6 +933,89 @@ fn set_mate_info_raw(r1: &mut RawRecord, r2: &mut RawRecord) { } } +/// Ports htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment(supp, matePrimary, +/// setMateCigar=true)`. +/// +/// 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) { + 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); + let mut editor = supp.tags_editor(); + if mate.unmapped { + editor.remove(SamTag::MC); + } else { + editor.update_string(SamTag::MC, mate.cigar.as_bytes()); + } + // htsjdk 5.0.0 sets MQ unconditionally from the mate primary's mapping quality. + editor.update_int(SamTag::MQ, i32::from(mate.mapq)); +} + +/// Finds the primary R1 and R2 record indices in a template, following fgbio `Template.r1`/`r2`: +/// the first non-secondary, non-supplementary read that is unpaired or first-of-pair, and the +/// first that is paired and second-of-pair, respectively. Either may be `None`. +/// +/// Like fgbio's `Bams.Template` (`Bams.scala:161,167`), a template carrying more than one primary +/// (non-secondary, non-supplementary) R1 — or more than one primary R2 — is malformed and rejected +/// with an error rather than silently keeping the first and passing the extra through unclipped and +/// with unrepaired mate info. The error message mirrors fgbio verbatim. +fn find_primary_pair_indices(records: &[RawRecord]) -> Result<(Option, Option)> { + let mut r1_idx = None; + let mut r2_idx = None; + for (i, rec) in records.iter().enumerate() { + if rec.is_secondary() || rec.is_supplementary() { + continue; + } + if !rec.is_paired() || rec.is_first_segment() { + if r1_idx.is_some() { + anyhow::bail!( + "Multiple non-secondary, non-supplemental R1s for {}", + String::from_utf8_lossy(rec.read_name()).trim_end_matches('\0') + ); + } + r1_idx = Some(i); + } else if rec.is_last_segment() { + if r2_idx.is_some() { + anyhow::bail!( + "Multiple non-secondary, non-supplemental R2s for {}", + String::from_utf8_lossy(rec.read_name()).trim_end_matches('\0') + ); + } + r2_idx = Some(i); + } + } + Ok((r1_idx, r2_idx)) +} + +/// Repairs mate information on the template's supplementary alignments after the primary pair +/// (at `r1_idx`/`r2_idx`) has been clipped and had its own mate info set. Mirrors fgbio +/// `ClipBam`: R1 supplementals point at the primary R2, R2 supplementals at the primary R1. +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() { + continue; + } + // 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); + } else if rec.is_last_segment() { + set_supplemental_mate_info_raw(rec, &r1_snap, r1_tlen); + } + } +} + /// Returns the number of clipped bases (soft + hard) in a raw record's CIGAR. fn clipped_bases_raw(record: &RawRecord) -> usize { record @@ -1004,6 +1036,261 @@ mod tests { use rstest::rstest; use std::path::PathBuf; + // R2-CLIP-02: the primary R1/R2 are located by SAM flags, not by position, so templates + // that carry secondary/supplementary alignments (len > 2) still resolve their primary pair. + #[test] + fn test_find_primary_pair_indices_ignores_secondary_and_supplementary() { + use crate::sam::RecordBuilder; + use fgumi_raw_bam::encode_record_buf_to_raw; + use noodles::sam::header::record::value::Map; + use noodles::sam::header::record::value::map::ReferenceSequence; + use std::num::NonZeroUsize; + + let ref_seq = Map::::new( + NonZeroUsize::new(100_000).expect("ref length must be nonzero"), + ); + let header = + noodles::sam::Header::builder().add_reference_sequence(b"chr1", ref_seq).build(); + let enc = |b: &noodles::sam::alignment::RecordBuf| { + encode_record_buf_to_raw(b, &header).expect("encode") + }; + let mapped = |first: bool, secondary: bool, supplementary: bool, start: usize| { + RecordBuilder::mapped_read() + .name("q") + .paired(true) + .first_segment(first) + .secondary(secondary) + .supplementary(supplementary) + .reference_sequence_id(0) + .alignment_start(start) + .cigar("50M") + .sequence(&"A".repeat(50)) + .build() + }; + + let r1 = mapped(true, false, false, 100); + let r2 = mapped(false, false, false, 300); + let supp_r1 = mapped(true, false, true, 700); + let sec_r2 = mapped(false, true, false, 900); + + // Primaries at positions 0/1 with trailing secondary + supplementary reads. + let recs = vec![enc(&r1), enc(&r2), enc(&supp_r1), enc(&sec_r2)]; + assert_eq!(find_primary_pair_indices(&recs).unwrap(), (Some(0), Some(1))); + + // Order-independent: a supplementary read first must not be mistaken for a primary. + let recs2 = vec![enc(&supp_r1), enc(&r2), enc(&r1)]; + assert_eq!(find_primary_pair_indices(&recs2).unwrap(), (Some(2), Some(1))); + + // A lone fragment (unpaired primary) resolves R1 only. + let frag = RecordBuilder::mapped_read() + .name("f") + .paired(false) + .reference_sequence_id(0) + .alignment_start(100) + .cigar("50M") + .sequence(&"A".repeat(50)) + .build(); + assert_eq!(find_primary_pair_indices(&[enc(&frag)]).unwrap(), (Some(0), None)); + } + + // A malformed template with two primary (non-secondary, non-supplementary) R1s — or two + // primary R2s — is rejected loudly, matching fgbio `Bams.Template` (`Bams.scala:161,167`) + // rather than silently keeping the first and passing the extra through unclipped. + #[test] + fn test_find_primary_pair_indices_rejects_duplicate_primaries() { + use crate::sam::RecordBuilder; + use fgumi_raw_bam::encode_record_buf_to_raw; + use noodles::sam::header::record::value::Map; + use noodles::sam::header::record::value::map::ReferenceSequence; + use std::num::NonZeroUsize; + + let ref_seq = Map::::new( + NonZeroUsize::new(100_000).expect("ref length must be nonzero"), + ); + let header = + noodles::sam::Header::builder().add_reference_sequence(b"chr1", ref_seq).build(); + let enc = |b: &noodles::sam::alignment::RecordBuf| { + encode_record_buf_to_raw(b, &header).expect("encode") + }; + let mapped = |first: bool, start: usize| { + RecordBuilder::mapped_read() + .name("dup") + .paired(true) + .first_segment(first) + .reference_sequence_id(0) + .alignment_start(start) + .cigar("50M") + .sequence(&"A".repeat(50)) + .build() + }; + + // Two primary R1s (both first-of-pair, neither secondary/supplementary). + let two_r1 = + vec![enc(&mapped(true, 100)), enc(&mapped(false, 300)), enc(&mapped(true, 500))]; + let err = find_primary_pair_indices(&two_r1).unwrap_err().to_string(); + assert_eq!(err, "Multiple non-secondary, non-supplemental R1s for dup"); + + // Two primary R2s (both last-of-pair). + let two_r2 = + vec![enc(&mapped(true, 100)), enc(&mapped(false, 300)), enc(&mapped(false, 500))]; + let err = find_primary_pair_indices(&two_r2).unwrap_err().to_string(); + assert_eq!(err, "Multiple non-secondary, non-supplemental R2s for dup"); + } + + /// Encodes a `RecordBuf` to a `RawRecord` with a shared single-contig header. Used by the + /// raw mate-info tests below. + fn encode_raw(rec: &noodles::sam::alignment::RecordBuf) -> RawRecord { + use fgumi_raw_bam::encode_record_buf_to_raw; + use noodles::sam::header::record::value::Map; + use noodles::sam::header::record::value::map::ReferenceSequence; + use std::num::NonZeroUsize; + + let ref_seq = Map::::new( + NonZeroUsize::new(100_000).expect("ref length must be nonzero"), + ); + let header = + noodles::sam::Header::builder().add_reference_sequence(b"chr1", ref_seq).build(); + encode_record_buf_to_raw(rec, &header).expect("encode") + } + + /// Builds a mapped raw record from the common fields the mate-info tests need. + fn raw_read( + first: bool, + supplementary: bool, + reverse: bool, + start: usize, + mapq: u8, + ) -> RawRecord { + use crate::sam::RecordBuilder; + encode_raw( + &RecordBuilder::mapped_read() + .name("q") + .paired(true) + .first_segment(first) + .supplementary(supplementary) + .reverse_complement(reverse) + .reference_sequence_id(0) + .alignment_start(start) + .mapping_quality(mapq) + .cigar("50M") + .sequence(&"A".repeat(50)) + .build(), + ) + } + + // 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. + #[test] + fn test_set_supplemental_mate_info_raw_mapped_mate() { + use fgumi_raw_bam::flags as rflags; + + // A reverse-strand primary mate at 1-based 301 (0-based 300), MAPQ 40, CIGAR 50M. + 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); + + assert_eq!(supp.mate_ref_id(), 0); + assert_eq!(supp.mate_pos(), 300); + assert_eq!(supp.template_length(), -120); + 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")); + assert_eq!(supp.tags().find_int(SamTag::MQ), Some(40)); + } + + // With an unmapped mate, set_supplemental_mate_info_raw flags the mate unmapped and drops MC. + #[test] + fn test_set_supplemental_mate_info_raw_unmapped_mate() { + use crate::sam::RecordBuilder; + use fgumi_raw_bam::flags as rflags; + + let mate = encode_raw( + &RecordBuilder::mapped_read() + .name("q") + .paired(true) + .first_segment(false) + .unmapped(true) + .reference_sequence_id(0) + .alignment_start(301) + .cigar("50M") + .sequence(&"A".repeat(50)) + .build(), + ); + // Seed an MC tag so we can confirm it is removed when the mate is unmapped. + let mut supp = raw_read(true, true, false, 700, 30); + 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); + + assert_ne!(supp.flags() & rflags::MATE_UNMAPPED, 0, "mate is unmapped"); + assert!(!supp.tags().contains(SamTag::MC), "MC dropped when mate unmapped"); + } + + // 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. + #[test] + fn test_fix_supplemental_mate_info() { + use fgumi_raw_bam::flags as rflags; + + let mut recs = vec![ + raw_read(true, false, false, 101, 60), // 0: primary R1 (forward, MAPQ 60) + raw_read(false, false, true, 301, 40), // 1: primary R2 (reverse, MAPQ 40) + raw_read(true, true, false, 701, 30), // 2: supplementary R1 + raw_read(false, true, false, 901, 20), // 3: supplementary R2 + ]; + recs[0].set_template_length(200); + recs[1].set_template_length(-200); + + fix_supplemental_mate_info(&mut recs, 0, 1); + + // Supp R1 (idx 2) takes primary R2 (idx 1) as its mate. + assert_eq!(recs[2].mate_ref_id(), 0); + 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 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); + } + + // The RecordBuf unsoftclipped_start/end helpers (used by the RecordBuf clip path) subtract or + // add only *soft* clips, ignore hard clips, and return None for unmapped reads. + #[test] + fn test_unsoftclipped_recordbuf_helpers() { + use crate::sam::RecordBuilder; + use crate::sam::record_utils::{unsoftclipped_end, unsoftclipped_start}; + + // 5H10S30M10S at 1-based 100: start = 100 - 10 (leading soft) = 90; hard clips ignored. + // end = 100 + 30 (ref span) - 1 + 10 (trailing soft) = 139. + let mapped = RecordBuilder::mapped_read() + .name("q") + .reference_sequence_id(0) + .alignment_start(100) + .cigar("5H10S30M10S") + .sequence(&"A".repeat(50)) + .build(); + assert_eq!(unsoftclipped_start(&mapped), Some(90)); + assert_eq!(unsoftclipped_end(&mapped), Some(139)); + + let unmapped = RecordBuilder::mapped_read() + .name("q") + .unmapped(true) + .reference_sequence_id(0) + .alignment_start(100) + .cigar("50M") + .sequence(&"A".repeat(50)) + .build(); + assert_eq!(unsoftclipped_start(&unmapped), None); + assert_eq!(unsoftclipped_end(&unmapped), None); + } + #[test] fn test_default_clip_parameters() { let clip = Clip { @@ -1024,7 +1311,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1056,7 +1342,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1089,7 +1374,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1121,7 +1405,6 @@ mod tests { upgrade_clipping: true, auto_clip_attributes: false, metrics: Some(PathBuf::from("metrics.txt")), - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1153,7 +1436,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: true, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1184,7 +1466,6 @@ mod tests { upgrade_clipping: true, auto_clip_attributes: true, metrics: Some(PathBuf::from("metrics.txt")), - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1220,7 +1501,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1230,36 +1510,6 @@ mod tests { assert_eq!(soft.clipping_mode, ClippingMode::Soft); } - #[test] - fn test_clip_with_sort_order_specification() { - let clip = Clip { - io: BamIoOptions { - input: PathBuf::from("input.bam"), - output: PathBuf::from("output.bam"), - async_reader: false, - }, - reference: PathBuf::from("reference.fa"), - clipping_mode: ClippingMode::Hard, - clip_overlapping_reads: false, - clip_extending_past_mate: false, - - read_one_five_prime: 0, - read_one_three_prime: 0, - read_two_five_prime: 0, - read_two_three_prime: 0, - upgrade_clipping: false, - auto_clip_attributes: false, - metrics: None, - sort_order: Some("coordinate".to_string()), - threading: ThreadingOptions::none(), - compression: CompressionOptions { compression_level: 1 }, - scheduler_opts: SchedulerOptions::default(), - queue_memory: QueueMemoryOptions::default(), - }; - - assert_eq!(clip.sort_order, Some("coordinate".to_string())); - } - #[test] fn test_clip_asymmetric_fixed_positions() { let clip = Clip { @@ -1280,7 +1530,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1314,7 +1563,6 @@ mod tests { upgrade_clipping: true, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1346,7 +1594,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1378,7 +1625,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1410,7 +1656,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: true, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1442,7 +1687,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1476,7 +1720,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1508,7 +1751,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1542,7 +1784,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1575,7 +1816,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1600,7 +1840,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1625,7 +1864,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1638,66 +1876,6 @@ mod tests { assert_eq!(hard.clipping_mode, ClippingMode::Hard); } - #[test] - fn test_clip_with_queryname_sort_order() { - let clip = Clip { - io: BamIoOptions { - input: PathBuf::from("input.bam"), - output: PathBuf::from("output.bam"), - async_reader: false, - }, - reference: PathBuf::from("reference.fa"), - clipping_mode: ClippingMode::Hard, - clip_overlapping_reads: true, - clip_extending_past_mate: true, - - read_one_five_prime: 0, - read_one_three_prime: 0, - read_two_five_prime: 0, - read_two_three_prime: 0, - upgrade_clipping: false, - auto_clip_attributes: false, - metrics: None, - sort_order: Some("queryname".to_string()), - threading: ThreadingOptions::none(), - compression: CompressionOptions { compression_level: 1 }, - scheduler_opts: SchedulerOptions::default(), - queue_memory: QueueMemoryOptions::default(), - }; - - assert_eq!(clip.sort_order, Some("queryname".to_string())); - } - - #[test] - fn test_clip_with_unsorted_sort_order() { - let clip = Clip { - io: BamIoOptions { - input: PathBuf::from("input.bam"), - output: PathBuf::from("output.bam"), - async_reader: false, - }, - reference: PathBuf::from("reference.fa"), - clipping_mode: ClippingMode::Hard, - clip_overlapping_reads: true, - clip_extending_past_mate: true, - - read_one_five_prime: 0, - read_one_three_prime: 0, - read_two_five_prime: 0, - read_two_three_prime: 0, - upgrade_clipping: false, - auto_clip_attributes: false, - metrics: None, - sort_order: Some("unsorted".to_string()), - threading: ThreadingOptions::none(), - compression: CompressionOptions { compression_level: 1 }, - scheduler_opts: SchedulerOptions::default(), - queue_memory: QueueMemoryOptions::default(), - }; - - assert_eq!(clip.sort_order, Some("unsorted".to_string())); - } - #[test] fn test_clip_single_read_end_clipping() { let clip = Clip { @@ -1718,7 +1896,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1752,7 +1929,6 @@ mod tests { upgrade_clipping: true, auto_clip_attributes: false, metrics: Some(PathBuf::from("metrics.txt")), - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1784,7 +1960,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1857,7 +2032,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1908,7 +2082,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -1958,7 +2131,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2008,7 +2180,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2058,7 +2229,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2108,7 +2278,6 @@ mod tests { upgrade_clipping: true, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2159,7 +2328,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: Some(metrics_path.clone()), - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2173,56 +2341,6 @@ mod tests { Ok(()) } - #[test] - fn test_clip_execute_with_sort_order() -> Result<()> { - let dir = TempDir::new()?; - let ref_path = create_test_reference(&dir); - let input_path = dir.path().join("input.bam"); - let output_path = dir.path().join("output.bam"); - - let mut builder = SamBuilder::with_single_ref("chr1", 200); - let _ = builder - .add_pair() - .name("read1") - .bases1("ACGTACGTACGTACGTACGT") - .contig(0) - .start1(10) - .start2(30) - .build(); - - builder.write(&input_path)?; - - let clip = Clip { - io: BamIoOptions { - input: input_path, - output: output_path.clone(), - async_reader: false, - }, - reference: ref_path, - clipping_mode: ClippingMode::Hard, - clip_overlapping_reads: true, - clip_extending_past_mate: false, - - read_one_five_prime: 0, - read_one_three_prime: 0, - read_two_five_prime: 0, - read_two_three_prime: 0, - upgrade_clipping: false, - auto_clip_attributes: false, - metrics: None, - sort_order: Some("queryname".to_string()), - threading: ThreadingOptions::none(), - compression: CompressionOptions { compression_level: 1 }, - scheduler_opts: SchedulerOptions::default(), - queue_memory: QueueMemoryOptions::default(), - }; - clip.execute("test")?; - - assert!(output_path.exists()); - - Ok(()) - } - #[test] fn test_clip_execute_with_fragment() -> Result<()> { let dir = TempDir::new()?; @@ -2260,7 +2378,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2311,7 +2428,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: true, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2362,7 +2478,6 @@ mod tests { upgrade_clipping: true, auto_clip_attributes: true, metrics: Some(metrics_path.clone()), - sort_order: Some("queryname".to_string()), threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2409,7 +2524,6 @@ mod tests { upgrade_clipping: false, // No clipping option auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2466,7 +2580,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading, compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), @@ -2529,7 +2642,6 @@ mod tests { upgrade_clipping: false, auto_clip_attributes: false, metrics: None, - sort_order: None, threading: ThreadingOptions::none(), compression: CompressionOptions { compression_level: 1 }, scheduler_opts: SchedulerOptions::default(), diff --git a/tests/integration/test_clip_command.rs b/tests/integration/test_clip_command.rs index 78f79dd57..f0f5968f0 100644 --- a/tests/integration/test_clip_command.rs +++ b/tests/integration/test_clip_command.rs @@ -10,7 +10,9 @@ use fgumi_lib::commands::clip::Clip; use fgumi_lib::commands::command::Command; use fgumi_raw_bam::{RawRecord, SamBuilder, flags}; use noodles::bam; +use noodles::sam::alignment::RecordBuf; use noodles::sam::alignment::io::Write as AlignmentWrite; +use noodles::sam::alignment::record::cigar::op::Kind as CigarKind; use std::fs; use std::path::Path; use tempfile::TempDir; @@ -162,3 +164,375 @@ fn test_clip_command_with_metrics() { cmd.execute("fgumi clip").expect("Clip command with metrics failed"); assert!(metrics_path.exists(), "Metrics file not created"); } + +/// Write a BAM from an explicit list of records (any flags), using the shared minimal header. +fn create_bam_from_records(path: &Path, records: &[RawRecord]) { + let header = create_minimal_header("chr1", 10000); + let mut writer = + bam::io::Writer::new(fs::File::create(path).expect("Failed to create BAM file")); + writer.write_header(&header).expect("Failed to write header"); + for r in records { + writer.write_alignment_record(&header, &to_record_buf(r)).expect("Failed to write record"); + } + writer.try_finish().expect("Failed to finish BAM"); +} + +/// One decoded output record: `(name, cigar-ops, raw-flags)`. +type OutputRecord = (String, Vec<(CigarKind, usize)>, u16); + +/// Read output records so tests can assert the *actual* clip state (CIGAR changed / +/// unchanged), not merely the record count. +fn read_output_records(path: &Path) -> Vec { + let mut reader = bam::io::Reader::new(fs::File::open(path).unwrap()); + let header = reader.read_header().unwrap(); + reader + .record_bufs(&header) + .map(|r| { + let rec = r.expect("Failed to read output record"); + let name = rec + .name() + .map(|n| String::from_utf8_lossy(n.as_ref()).into_owned()) + .unwrap_or_default(); + let ops: Vec<(CigarKind, usize)> = + rec.cigar().as_ref().iter().map(|op| (op.kind(), op.len())).collect(); + (name, ops, u16::from(rec.flags())) + }) + .collect() +} + +/// Read output records as full noodles `RecordBuf`s so tests can assert mate metadata +/// (mate ref/pos/strand, MC, MQ, TLEN) — not just CIGAR — after clipping. +fn read_output_record_bufs(path: &Path) -> Vec { + let mut reader = bam::io::Reader::new(fs::File::open(path).unwrap()); + let header = reader.read_header().unwrap(); + reader.record_bufs(&header).map(|r| r.expect("Failed to read output record")).collect() +} + +/// The record's CIGAR as `(kind, len)` ops (e.g. `[(Match, 98), (HardClip, 2)]`). +fn cigar_ops(rec: &RecordBuf) -> Vec<(CigarKind, usize)> { + rec.cigar().as_ref().iter().map(|op| (op.kind(), op.len())).collect() +} + +/// The record's mate-CIGAR (`MC`) tag as a string, if present. +fn mate_cigar(rec: &RecordBuf) -> Option { + use noodles::sam::alignment::record::data::field::Tag; + use noodles::sam::alignment::record_buf::data::field::Value; + match rec.data().get(&Tag::MATE_CIGAR) { + Some(Value::String(s)) => Some(s.to_string()), + _ => None, + } +} + +/// The record's mate-mapping-quality (`MQ`) tag as an integer, if present. +fn mate_mapq(rec: &RecordBuf) -> Option { + use noodles::sam::alignment::record::data::field::Tag; + use noodles::sam::alignment::record_buf::data::field::Value; + rec.data().get(&Tag::MATE_MAPPING_QUALITY).and_then(Value::as_int) +} + +/// Build a mapped read with the given flags/position/CIGAR-length (all `M`). +fn mapped_read(name: &[u8], flags: u16, pos: i32, match_len: u32, mapq: u8) -> RawRecord { + let len = match_len as usize; + let mut b = SamBuilder::new(); + b.read_name(name) + .sequence(&vec![b'A'; len]) + .qualities(&vec![30; len]) + .flags(flags) + .ref_id(0) + .pos(pos) + .mapq(mapq) + .cigar_ops(&[match_len << 4]) // M + .mate_ref_id(0) + .mate_pos(pos); + b.build() +} + +/// `--threads` mode routes through `execute_threads_mode`; exercise the `(Some, Some)` primary-pair +/// branch with fixed R1/R2 clipping, overlap clipping, mate-extension clipping and clip upgrading +/// all enabled. R1 (fwd, 150M) and R2 (rev, 100M) fully overlap, so overlap clipping engages. +#[test] +fn test_clip_command_threads_mode_primary_pair_all_options() { + let temp_dir = TempDir::new().unwrap(); + let input_bam = temp_dir.path().join("input.bam"); + let output_bam = temp_dir.path().join("output.bam"); + let ref_path = create_test_reference(temp_dir.path()); + + let mut r1 = + mapped_read(b"p", flags::PAIRED | flags::FIRST_SEGMENT | flags::MATE_REVERSE, 99, 150, 60); + let mut r2 = + mapped_read(b"p", flags::PAIRED | flags::LAST_SEGMENT | flags::REVERSE, 99, 100, 60); + // A non-zero TLEN is required for `is_fr_pair_raw` to recognize the pair (htsjdk keys FR + // detection off the insert size); without it overlap and mate-extension clipping would + // silently no-op and only the fixed 5'/3' clips would apply. + fgumi_raw_bam::set_template_length(r1.as_mut_vec(), 150); + fgumi_raw_bam::set_template_length(r2.as_mut_vec(), -150); + create_bam_from_records(&input_bam, &[r1, r2]); + + let cmd = Clip::try_parse_from([ + "clip", + "--input", + input_bam.to_str().unwrap(), + "--output", + output_bam.to_str().unwrap(), + "--reference", + ref_path.to_str().unwrap(), + "--threads", + "2", + "--read-one-five-prime", + "1", + "--read-one-three-prime", + "1", + "--read-two-five-prime", + "1", + "--read-two-three-prime", + "1", + "--clip-overlapping-reads", + "true", + "--clip-bases-past-mate", + "true", + "--upgrade-clipping", + "true", + "--compression-level", + "1", + ]) + .expect("failed to parse clip args"); + cmd.execute("fgumi clip").expect("Clip command failed"); + + // Exact post-clip CIGARs, derived independently from the clip pipeline (Hard mode), + // not just "some hard clip exists". Inputs are fixed: R1 fwd 150M@100 (1-based), + // R2 rev 100M@100. + // + // R1: fixed 5'/3' = 1 each -> 1H148M1H, ref 101..248 + // overlap clip trims R1's 3' end down to the pair midpoint + // (midpoint(r1_start=101, r2_end=198) = 149), clipping ref 150..248 (99 M bases) + // -> 1H49M100H, ref 101..149 + // R2: fixed 5'/3' = 1 each -> 1H98M1H, ref 101..198 + // overlap clip trims R2's 5' (low-coord) end up to midpoint+1 = 150, + // clipping ref 101..149 (49 M bases) -> 50H49M1H, ref 150..198 + // + // Mate-extension clipping then finds nothing past the mate (the reads now meet at the + // midpoint), so it is a no-op here — but its code path still runs. Totals check out: + // R1 = 1+49+100 = 150 read bases, R2 = 50+49+1 = 100 read bases. + let recs = read_output_records(&output_bam); + assert_eq!(recs.len(), 2, "both reads retained"); + let expected_r1 = + vec![(CigarKind::HardClip, 1), (CigarKind::Match, 49), (CigarKind::HardClip, 100)]; + let expected_r2 = + vec![(CigarKind::HardClip, 50), (CigarKind::Match, 49), (CigarKind::HardClip, 1)]; + for (name, ops, flag_bits) in &recs { + assert_eq!(name, "p", "unexpected read name"); + let is_r1 = flag_bits & flags::FIRST_SEGMENT != 0; + let (label, expected) = if is_r1 { ("R1", &expected_r1) } else { ("R2", &expected_r2) }; + assert_eq!(ops, expected, "{label} exact post-clip CIGAR mismatch; cigar={ops:?}"); + } +} + +/// `--threads` mode with a lone unpaired fragment exercises the `(Some, None)` fragment branch. +#[test] +fn test_clip_command_threads_mode_fragment() { + let temp_dir = TempDir::new().unwrap(); + let input_bam = temp_dir.path().join("input.bam"); + let output_bam = temp_dir.path().join("output.bam"); + let ref_path = create_test_reference(temp_dir.path()); + + // No PAIRED flag: find_primary_pair_indices resolves R1 only. + let frag = mapped_read(b"f", 0, 99, 100, 60); + create_bam_from_records(&input_bam, &[frag]); + + let cmd = Clip::try_parse_from([ + "clip", + "--input", + input_bam.to_str().unwrap(), + "--output", + output_bam.to_str().unwrap(), + "--reference", + ref_path.to_str().unwrap(), + "--threads", + "2", + "--read-one-five-prime", + "2", + "--read-one-three-prime", + "2", + "--upgrade-clipping", + "true", + "--compression-level", + "1", + ]) + .expect("failed to parse clip args"); + cmd.execute("fgumi clip").expect("Clip command failed"); + + // The lone forward-strand fragment is hard-clipped 2bp at each end (read-one 5'/3' = 2, + // Hard mode) => 2H96M2H, proving the (Some, None) branch actually clipped rather than + // just passing the record through. + let recs = read_output_records(&output_bam); + assert_eq!(recs.len(), 1, "fragment retained"); + assert_eq!( + recs[0].1, + vec![(CigarKind::HardClip, 2), (CigarKind::Match, 96), (CigarKind::HardClip, 2)], + "fragment must be hard-clipped 2bp at each end (2H96M2H); cigar={:?}", + recs[0].1 + ); +} + +/// `--threads` mode with a template that has no primary read (secondary only) exercises the +/// `(None, None)` no-op branch — the record passes through untouched. +#[test] +fn test_clip_command_threads_mode_secondary_only() { + let temp_dir = TempDir::new().unwrap(); + let input_bam = temp_dir.path().join("input.bam"); + let output_bam = temp_dir.path().join("output.bam"); + let ref_path = create_test_reference(temp_dir.path()); + + let secondary = + mapped_read(b"s", flags::PAIRED | flags::FIRST_SEGMENT | flags::SECONDARY, 99, 100, 60); + create_bam_from_records(&input_bam, &[secondary]); + + let cmd = Clip::try_parse_from([ + "clip", + "--input", + input_bam.to_str().unwrap(), + "--output", + output_bam.to_str().unwrap(), + "--reference", + ref_path.to_str().unwrap(), + "--threads", + "2", + "--read-one-five-prime", + "1", + "--compression-level", + "1", + ]) + .expect("failed to parse clip args"); + cmd.execute("fgumi clip").expect("Clip command failed"); + + // The (None, None) branch is a true no-op: the secondary read passes through with its + // CIGAR unchanged (still 100M, not clipped) and its SECONDARY flag preserved. + let recs = read_output_records(&output_bam); + assert_eq!(recs.len(), 1, "secondary read passed through"); + let (_name, ops, flag_bits) = &recs[0]; + assert_eq!( + ops, + &vec![(CigarKind::Match, 100)], + "secondary read must be unchanged (100M no-op); cigar={ops:?}" + ); + assert_ne!(flag_bits & flags::SECONDARY, 0, "SECONDARY flag must be preserved"); +} + +/// `--threads` mode with a chimeric/split template — primary R1 + primary R2 + a supplementary +/// R1 — exercises the PR's headline fix end-to-end: after the primary pair is clipped, +/// `fix_supplemental_mate_info` must repair the supplementary alignment's mate metadata to point +/// at the *post-clip* primary R2. The `clip.rs` unit test covers the helper in isolation; this +/// verifies the `execute_threads_mode` wiring (correct r1/r2 indices, and that the supplemental +/// mate snapshot is taken *after* the primary is clipped, not before). A wiring bug (wrong +/// index, or records mutated out of order) would be invisible to every other test in this file. +#[test] +fn test_clip_command_threads_mode_supplementary_mate_repair() { + let temp_dir = TempDir::new().unwrap(); + let input_bam = temp_dir.path().join("input.bam"); + let output_bam = temp_dir.path().join("output.bam"); + let ref_path = create_test_reference(temp_dir.path()); + + // Chimeric template: a primary FR pair (R1 fwd 100M @100, R2 rev 100M @300) plus a + // supplementary R1 mapped far away (50M @5000). `mapped_read` seeds each record's mate fields + // to its *own* position with no MC/MQ, so the supplementary's input mate info is deliberately + // stale (points at 5000, not the primary R2) — only a correct repair yields the expected + // primary-R2-derived values. + let r1 = mapped_read( + b"chim", + flags::PAIRED | flags::FIRST_SEGMENT | flags::MATE_REVERSE, + 99, + 100, + 60, + ); + let r2 = + mapped_read(b"chim", flags::PAIRED | flags::LAST_SEGMENT | flags::REVERSE, 299, 100, 40); + let supp = mapped_read( + b"chim", + flags::PAIRED | flags::FIRST_SEGMENT | flags::SUPPLEMENTARY, + 4999, + 50, + 30, + ); + create_bam_from_records(&input_bam, &[r1, r2, supp]); + + // Clip only read-two's 5' end by 2bp. R2 is reverse, so its 5' end is the high-coordinate + // (right) end: 100M @300 -> 98M2H, alignment start unchanged at 300 (1-based). The + // supplementary R1 is never touched by the clip loop (only the primary pair is), so it keeps + // 50M — but its mate fields must be rewritten to the post-clip primary R2. + let cmd = Clip::try_parse_from([ + "clip", + "--input", + input_bam.to_str().unwrap(), + "--output", + output_bam.to_str().unwrap(), + "--reference", + ref_path.to_str().unwrap(), + "--threads", + "2", + "--read-two-five-prime", + "2", + "--compression-level", + "1", + ]) + .expect("failed to parse clip args"); + cmd.execute("fgumi clip").expect("Clip command failed"); + + let recs = read_output_record_bufs(&output_bam); + assert_eq!(recs.len(), 3, "all three records retained"); + + let flag_bits = |r: &RecordBuf| u16::from(r.flags()); + let supp = recs + .iter() + .find(|r| flag_bits(r) & flags::SUPPLEMENTARY != 0) + .expect("supplementary record present"); + let r2 = recs + .iter() + .find(|r| { + let f = flag_bits(r); + f & flags::LAST_SEGMENT != 0 + && f & flags::SUPPLEMENTARY == 0 + && f & flags::SECONDARY == 0 + }) + .expect("primary R2 present"); + + // Independent oracle for the clip itself: primary R2 (reverse, 5' clip of 2) -> 98M2H @300, + // start unchanged, MAPQ 40. Derived by hand from the fixed inputs, not from the code under test. + assert_eq!( + cigar_ops(r2), + vec![(CigarKind::Match, 98), (CigarKind::HardClip, 2)], + "primary R2 must be 98M2H after a 5' clip of 2" + ); + assert_eq!(usize::from(r2.alignment_start().unwrap()), 300, "R2 alignment start unchanged"); + assert_eq!(r2.mapping_quality().map(u8::from), Some(40), "R2 MAPQ unchanged"); + assert_ne!(flag_bits(r2) & flags::REVERSE, 0, "primary R2 is reverse"); + + // The supplementary alignment itself is not clipped (only the primary pair is). + assert_eq!(cigar_ops(supp), vec![(CigarKind::Match, 50)], "supplementary unchanged (50M)"); + + // Headline fix: the supplementary's mate metadata now points at the post-clip primary R2. + // Absolute values are hand-derived from the inputs (mate ref 0, mate pos 300, MC 98M2H, MQ 40, + // mate-reverse set) — so this is an independent oracle, not a self-consistency check. + assert_eq!(supp.mate_reference_sequence_id(), Some(0), "mate ref = primary R2 ref"); + assert_eq!( + usize::from(supp.mate_alignment_start().unwrap()), + 300, + "mate pos = primary R2 start" + ); + assert_ne!( + flag_bits(supp) & flags::MATE_REVERSE, + 0, + "mate-reverse must be set (primary R2 is reverse)" + ); + 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"); + assert_eq!( + supp.template_length(), + -r2.template_length(), + "supplementary TLEN = negation of primary R2 TLEN" + ); + + // Cross-check: the repaired mate fields agree with the actual primary R2 record, proving the + // snapshot was taken from the correct (post-clip) primary and not a stale/wrong index. + assert_eq!(supp.mate_reference_sequence_id(), r2.reference_sequence_id()); + assert_eq!(supp.mate_alignment_start(), r2.alignment_start()); +}