diff --git a/modkit-core/src/repair_tags.rs b/modkit-core/src/repair_tags.rs index 3aa8f2a..cc78006 100644 --- a/modkit-core/src/repair_tags.rs +++ b/modkit-core/src/repair_tags.rs @@ -6,6 +6,7 @@ use clap::Args; use derive_new::new; use indicatif::{MultiProgress, ProgressBar}; use log::{debug, error, info, warn}; +use memchr::memmem::Finder; use rayon::prelude::*; use rust_htslib::bam::record::{Aux, AuxArray}; use rust_htslib::bam::{self, Read}; @@ -324,17 +325,14 @@ fn repair_record_pair(record_pair: RecordPair) -> anyhow::Result { if donor_seq.len() < acceptor_seq.len() { bail!("donor sequence for {read_name} is longer than acceptor sequence") } - let matches = donor_seq.match_indices(&acceptor_seq); - - let starts = - matches.into_iter().map(|(start, _)| start).collect::>(); - if starts.len() > 1 { - bail!("multiple potential corrections found for {read_name}") - } else if starts.is_empty() { + let finder = Finder::new(acceptor_seq.as_bytes()); + let Some(start) = finder.find(donor_seq.as_bytes()) else { bail!("acceptor sequence is not a substring of the donor sequence") + }; + if finder.find(&donor_seq.as_bytes()[start + 1..]).is_some() { + bail!("multiple potential corrections found for {read_name}") } else { let acceptor_seq_len = acceptor_seq.len(); - let start = *starts.get(0).unwrap(); let end = start + acceptor_seq_len; let mm_style = modbase_info.mm_style; @@ -392,3 +390,102 @@ fn repair_record_pair(record_pair: RecordPair) -> anyhow::Result { Ok(repaired_record) } } + +#[cfg(test)] +mod tests { + use std::sync::Arc; + + use rust_htslib::bam::{self, record::Aux}; + + use super::{repair_record_pair, RecordPair}; + + fn make_record(name: &[u8], sequence: &str) -> bam::Record { + let mut record = bam::Record::new(); + record.set(name, None, sequence.as_bytes(), &vec![255; sequence.len()]); + record + } + + fn make_donor(name: &[u8], sequence: &str) -> bam::Record { + let mut record = make_record(name, sequence); + record.push_aux(b"MM", Aux::String("A+a.,0;")).unwrap(); + record.push_aux(b"ML", Aux::ArrayU8((&[200u8][..]).into())).unwrap(); + record + } + + fn repair(donor: &str, acceptor: &str) -> anyhow::Result { + repair_record_pair(RecordPair::new( + Arc::new(make_donor(b"read", donor)), + make_record(b"read", acceptor), + )) + } + + fn repair_reverse( + donor_stored_sequence: &str, + acceptor_stored_sequence: &str, + ) -> anyhow::Result { + let mut donor = make_donor(b"reverse", donor_stored_sequence); + donor.set_reverse(); + let mut acceptor = make_record(b"reverse", acceptor_stored_sequence); + acceptor.set_reverse(); + repair_record_pair(RecordPair::new(Arc::new(donor), acceptor)) + } + + #[test] + fn repair_rejects_overlapping_ambiguous_placements() { + let error = repair("ACACAC", "ACAC").unwrap_err(); + assert!( + error.to_string().contains("multiple potential corrections"), + "unexpected error: {error:#}" + ); + + assert!(repair("ACACGGACAC", "ACAC") + .unwrap_err() + .to_string() + .contains("multiple potential corrections")); + } + + #[test] + fn repair_preserves_unique_and_absent_placement_behavior() { + let repaired = repair("TACACG", "ACAC").unwrap(); + assert!(matches!(repaired.aux(b"MM").unwrap(), Aux::String("A+a.,0;"))); + match repaired.aux(b"ML").unwrap() { + Aux::ArrayU8(qualities) => { + assert_eq!(qualities.iter().collect::>(), vec![200]); + } + other => panic!("unexpected ML tag: {other:?}"), + } + assert!(matches!(repaired.aux(b"MN").unwrap(), Aux::U32(4))); + + assert!(repair("ACACAC", "CGCG") + .unwrap_err() + .to_string() + .contains("not a substring")); + } + + #[test] + fn repair_searches_forward_sequences_for_reverse_records() { + // Stored placement starts at 3, while the forward-sequence placement + // TACGAAA -> ACG starts at 1. + let repaired = repair_reverse("TTTCGTA", "CGT").unwrap(); + assert!(repaired.is_reverse()); + assert!(matches!(repaired.aux(b"MM").unwrap(), Aux::String("A+a.,0;"))); + match repaired.aux(b"ML").unwrap() { + Aux::ArrayU8(qualities) => { + assert_eq!(qualities.iter().collect::>(), vec![200]); + } + other => panic!("unexpected ML tag: {other:?}"), + } + assert!(matches!(repaired.aux(b"MN").unwrap(), Aux::U32(3))); + + assert!(repair_reverse("GTGTGT", "GTGT") + .unwrap_err() + .to_string() + .contains("multiple potential corrections")); + } + + #[test] + fn repair_returns_errors_for_empty_sequences() { + assert!(repair("ACAC", "").is_err()); + assert!(repair("", "ACAC").is_err()); + } +}