From fce7ec2176d4c97db28725ca26fcbae327a573cf Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 01:47:18 -0700 Subject: [PATCH 1/3] Preserve independent MM group progression --- modkit-core/src/pileup/base_mods_adapter.rs | 40 +++++++++++++++++++-- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/modkit-core/src/pileup/base_mods_adapter.rs b/modkit-core/src/pileup/base_mods_adapter.rs index 3dbf3d5..c790b06 100644 --- a/modkit-core/src/pileup/base_mods_adapter.rs +++ b/modkit-core/src/pileup/base_mods_adapter.rs @@ -221,6 +221,9 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { let mut pos = self.left_to_right_seq_pos; let mut done = false; let mut inferred = true; + // A different MM group can count down to zero while this position is + // scanned. Only code slots already at zero have an explicit call here. + let mut explicit_codes_at_position = [false; SIZE]; while !done && pos < self.seq.len() { let base = if self.reverse { @@ -236,6 +239,7 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { done = true; mod_pos = Some(pos); inferred = false; + explicit_codes_at_position[i] = true; } } Some(skip_count) => { @@ -273,6 +277,7 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { let mut mod_code = ModCodeRepr::Code(base as char); for i in 0..self.n_codes { if self.canonical_bases[i] == base + && explicit_codes_at_position[i] && self.mm_next[i].map(|x| x == 0).unwrap_or(false) { let q = self.ml.get(self.ml_pos[i]).unwrap(); @@ -319,7 +324,7 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { )) }; - self.move_forward(mod_pos, base); + self.move_forward(mod_pos, base, &explicit_codes_at_position); mod_state } else { @@ -344,13 +349,18 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { } #[inline] - fn move_forward(&mut self, last_pos: usize, base: u8) { + fn move_forward( + &mut self, + last_pos: usize, + base: u8, + explicit_codes_at_position: &[bool; SIZE], + ) { self.left_to_right_seq_pos = last_pos.saturating_add(1); for i in (0..self.n_codes).filter(|i| self.canonical_bases[*i] == base) { assert!(i < SIZE); match &mut self.mm_next[i] { - Some(x) if *x == 0 => { + Some(x) if *x == 0 && explicit_codes_at_position[i] => { let num_explicit_positions = self.num_explicit_positions[i].saturating_sub(1); if num_explicit_positions == 0 { @@ -694,6 +704,30 @@ mod base_mods_adapter_tests { assert!(mod_state.inferred); } + #[test] + fn test_independent_same_base_groups_keep_their_delta_progression() { + let mut record = + make_record("C+m?,0;C+h?,1;", &[200, 250], "CC", None, false); + record.push_aux(b"MN", Aux::I32(2)).unwrap(); + + let mut scanner = BaseModsAdapter::<2>::new(&record).unwrap(); + let mod_states = std::iter::from_fn(|| { + scanner.next_modified_position_no_thresh().unwrap() + }) + .map(|state| { + (state.mod_position, state.mod_code, state.mod_qual, state.modified) + }) + .collect::>(); + + assert_eq!( + mod_states, + vec![ + (0, METHYL_CYTOSINE, 200, true), + (1, HYDROXY_METHYL_CYTOSINE, 250, true), + ] + ); + } + #[test] fn test_implicit_calls() { let seq = "ATCATCATTCCTACCGCTATAGCCT"; From 573a28aa02e757d67305bef3ed0d50ee82bd2185 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 02:11:29 -0700 Subject: [PATCH 2/3] Support MM tags with more than 16 modification codes --- modkit-core/Cargo.toml | 1 + modkit-core/src/pileup/base_mods_adapter.rs | 453 +++++++++++++++----- 2 files changed, 340 insertions(+), 114 deletions(-) diff --git a/modkit-core/Cargo.toml b/modkit-core/Cargo.toml index 881c7d8..d28076e 100644 --- a/modkit-core/Cargo.toml +++ b/modkit-core/Cargo.toml @@ -47,6 +47,7 @@ rust-htslib = "0.46.0" rust-lapper = "1.1.0" rustc-hash = "1.1.0" rv = "=0.16.0" +smallvec = { version = "1.15.2", features = ["const_generics"] } statrs = "0.16.0" substring = "1.4.5" thiserror = "2.0.11" diff --git a/modkit-core/src/pileup/base_mods_adapter.rs b/modkit-core/src/pileup/base_mods_adapter.rs index c790b06..2911359 100644 --- a/modkit-core/src/pileup/base_mods_adapter.rs +++ b/modkit-core/src/pileup/base_mods_adapter.rs @@ -7,6 +7,7 @@ use rust_htslib::bam::{ ext::{BamRecordExtensions, IterAlignedPairs}, record::AuxArray, }; +use smallvec::SmallVec; use crate::{ errs::MkResult, @@ -25,21 +26,25 @@ pub(crate) struct ModState { pub mod_qual: u8, } +#[derive(Debug, Copy, Clone)] +struct ModCodeState { + mm_pos: usize, + ml_pos: usize, + ml_stride: usize, + mod_code: ModCodeRepr, + mm_next: Option, + num_explicit_positions: u32, + canonical_base: u8, + implicit: bool, + explicit_at_position: bool, +} + #[derive(Debug)] pub(crate) struct BaseModsAdapter<'a, const SIZE: usize = 16> { mm: &'a [u8], seq: bam::record::Seq<'a>, ml: AuxArray<'a, u8>, - mod_codes: [ModCodeRepr; SIZE], - canonical_bases: [u8; SIZE], - mm_pos: [usize; SIZE], - ml_pos: [usize; SIZE], - // strands: [u8; SIZE], // could be bitvec - implicits: [bool; SIZE], - ml_strides: [usize; SIZE], - n_codes: usize, - mm_next: [Option; SIZE], - num_explicit_positions: [u32; SIZE], + code_states: SmallVec<[ModCodeState; SIZE]>, reverse: bool, left_to_right_seq_pos: usize, } @@ -88,18 +93,8 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { let mut i = 0; let mut ml_start = 0usize; - let mut mod_codes = [ModCodeRepr::Code('N'); SIZE]; - let mut canonical_bases = [0u8; SIZE]; - // let mut strands = [0u8; SIZE]; - let mut ml_strides = [0usize; SIZE]; - let mut implicits = [false; SIZE]; - let mut mm_pos = [0usize; SIZE]; - let mut ml_pos = [0usize; SIZE]; - let mut num_explicit_positions = [0u32; SIZE]; - let mut n_codes = 0; - let mut mm_next = [None; SIZE]; + let mut code_states = SmallVec::<[ModCodeState; SIZE]>::new(); while i < mm.len() { - assert!(n_codes < SIZE); let base = mm[i]; i += 1; let strand = mm[i]; @@ -107,8 +102,8 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { bail!("duplex data not currently supported") } i += 1; - let (mods_in_rec, offset) = - parse_mod_code(&mm[i..], &mut mod_codes, n_codes); + let (parsed_mod_codes, offset) = parse_mod_code(&mm[i..]); + let mods_in_rec = parsed_mod_codes.len(); i += offset; assert!(i < mm.len()); let implicit_mode = match mm[i] { @@ -173,37 +168,23 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { }; for j in 0..mods_in_rec { - mm_pos[j + n_codes] = mm_idx; - canonical_bases[j + n_codes] = base; - // strands[j + n_codes] = strand; - ml_strides[j + n_codes] = mods_in_rec; - implicits[j + n_codes] = implicit_mode; - mm_next[j + n_codes] = delta; - num_explicit_positions[j + n_codes] = n_deltas; - ml_pos[j + n_codes] = ml_idx + j; + code_states.push(ModCodeState { + mm_pos: mm_idx, + ml_pos: ml_idx + j, + ml_stride: mods_in_rec, + mod_code: parsed_mod_codes.get(j), + mm_next: delta, + num_explicit_positions: n_deltas, + canonical_base: base, + implicit: implicit_mode, + explicit_at_position: false, + }); } - n_codes += mods_in_rec; ml_start += n_deltas as usize * mods_in_rec; i += record_end + 1; } - Ok(Self { - mm, - ml, - seq, - mod_codes, - canonical_bases, - n_codes, - implicits, - mm_pos, - ml_pos, - // strands, - ml_strides, - mm_next, - num_explicit_positions, - reverse, - left_to_right_seq_pos: 0, - }) + Ok(Self { mm, ml, seq, code_states, reverse, left_to_right_seq_pos: 0 }) } pub fn next_modified_position_no_thresh( @@ -222,8 +203,12 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { let mut done = false; let mut inferred = true; // A different MM group can count down to zero while this position is - // scanned. Only code slots already at zero have an explicit call here. - let mut explicit_codes_at_position = [false; SIZE]; + // scanned. Only code states already at zero have an explicit call here. + // Store readiness with each dynamically sized state so this scan does + // not allocate, including when the inline capacity is exceeded. + self.code_states + .iter_mut() + .for_each(|state| state.explicit_at_position = false); while !done && pos < self.seq.len() { let base = if self.reverse { @@ -231,32 +216,34 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { } else { self.seq[pos] }; - for i in 0..self.n_codes { - match &mut self.mm_next[i] { - Some(skip_count) if *skip_count == 0u32 => { + for code_state in &mut self.code_states { + match code_state.mm_next { + Some(0) => { // there is at least one explicit modification here - if self.canonical_bases[i] == base { + if code_state.canonical_base == base { done = true; mod_pos = Some(pos); inferred = false; - explicit_codes_at_position[i] = true; + code_state.explicit_at_position = true; } } Some(skip_count) => { - if self.canonical_bases[i] == base { - if self.implicits[i] { + if code_state.canonical_base == base { + if code_state.implicit { done = true; mod_pos = Some(pos); } else { - *skip_count = skip_count - .checked_sub(1u32) - .expect("should not go off the end"); + code_state.mm_next = Some( + skip_count + .checked_sub(1u32) + .expect("should not go off the end"), + ); } } } None => { - if self.canonical_bases[i] == base { - if self.implicits[i] { + if code_state.canonical_base == base { + if code_state.implicit { done = true; mod_pos = Some(pos); } @@ -275,14 +262,14 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { self.seq[mod_pos] }; let mut mod_code = ModCodeRepr::Code(base as char); - for i in 0..self.n_codes { - if self.canonical_bases[i] == base - && explicit_codes_at_position[i] - && self.mm_next[i].map(|x| x == 0).unwrap_or(false) + for code_state in &self.code_states { + if code_state.canonical_base == base + && code_state.explicit_at_position + && code_state.mm_next.map(|x| x == 0).unwrap_or(false) { - let q = self.ml.get(self.ml_pos[i]).unwrap(); + let q = self.ml.get(code_state.ml_pos).unwrap(); if q > mod_qual { - mod_code = self.mod_codes[i]; + mod_code = code_state.mod_code; mod_qual = q; } total_mod_qual = total_mod_qual.saturating_add(q); @@ -324,7 +311,7 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { )) }; - self.move_forward(mod_pos, base, &explicit_codes_at_position); + self.move_forward(mod_pos, base); mod_state } else { @@ -336,8 +323,8 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { pub fn primary_bases_in_record(&self) -> u8 { let mut bs = 0u8; - for &raw_can_base in self.canonical_bases.iter().take(self.n_codes) { - match raw_can_base { + for code_state in &self.code_states { + match code_state.canonical_base { b'A' => bs.view_bits_mut::().set(0, true), b'C' => bs.view_bits_mut::().set(1, true), b'G' => bs.view_bits_mut::().set(2, true), @@ -349,25 +336,22 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { } #[inline] - fn move_forward( - &mut self, - last_pos: usize, - base: u8, - explicit_codes_at_position: &[bool; SIZE], - ) { + fn move_forward(&mut self, last_pos: usize, base: u8) { self.left_to_right_seq_pos = last_pos.saturating_add(1); - for i in (0..self.n_codes).filter(|i| self.canonical_bases[*i] == base) + for code_state in self + .code_states + .iter_mut() + .filter(|state| state.canonical_base == base) { - assert!(i < SIZE); - match &mut self.mm_next[i] { - Some(x) if *x == 0 && explicit_codes_at_position[i] => { + match code_state.mm_next { + Some(0) if code_state.explicit_at_position => { let num_explicit_positions = - self.num_explicit_positions[i].saturating_sub(1); + code_state.num_explicit_positions.saturating_sub(1); if num_explicit_positions == 0 { // done - self.mm_next[i] = None; + code_state.mm_next = None; } else { - let mm_pos = self.mm_pos[i]; + let mm_pos = code_state.mm_pos; if self.reverse { let mut p = mm_pos.saturating_sub(1); while self.mm[p] != b',' { @@ -378,27 +362,30 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { let c = self.mm[idx]; val = val * 10 + (c - b'0') as u32; } - self.mm_next[i] = Some(val); - self.mm_pos[i] = p; - self.ml_pos[i] -= self.ml_strides[i]; + code_state.mm_next = Some(val); + code_state.mm_pos = p; + code_state.ml_pos -= code_state.ml_stride; // if num_explicit_positions > 1 { // } } else { let (delta, offset) = parse_int::(&self.mm[mm_pos..]); - self.mm_next[i] = Some(delta); - self.mm_pos[i] = offset + mm_pos; - self.ml_pos[i] += self.ml_strides[i]; + code_state.mm_next = Some(delta); + code_state.mm_pos = offset + mm_pos; + code_state.ml_pos += code_state.ml_stride; // if num_explicit_positions > 1 { // } } - self.num_explicit_positions[i] = num_explicit_positions; + code_state.num_explicit_positions = + num_explicit_positions; } } - Some(x) if self.implicits[i] => { - *x = x - .checked_sub(1u32) - .expect("should not go less than zero"); + Some(skip_count) if code_state.implicit => { + code_state.mm_next = Some( + skip_count + .checked_sub(1u32) + .expect("should not go less than zero"), + ); } _ => {} } @@ -406,12 +393,32 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { } } +enum ParsedModCodes<'a> { + Alphabetic(&'a [u8]), + ChEbi(u32), +} + +impl ParsedModCodes<'_> { + fn len(&self) -> usize { + match self { + Self::Alphabetic(codes) => codes.len(), + Self::ChEbi(_) => 1, + } + } + + fn get(&self, index: usize) -> ModCodeRepr { + match self { + Self::Alphabetic(codes) => ModCodeRepr::Code(codes[index] as char), + Self::ChEbi(code) => { + debug_assert_eq!(index, 0); + ModCodeRepr::ChEbi(*code) + } + } + } +} + #[inline(always)] -fn parse_mod_code( - bytes: &[u8], - mod_codes: &mut [ModCodeRepr; SIZE], - mut n_mods_so_far: usize, -) -> (usize, usize) { +fn parse_mod_code(bytes: &[u8]) -> (ParsedModCodes<'_>, usize) { let mut idx = 0; let mut val: u32 = 0; @@ -434,25 +441,18 @@ fn parse_mod_code( } if !all_digits { - let mut mods_parsed = 0usize; + let start = idx; while bytes[idx] != b'.' && bytes[idx] != b'?' && bytes[idx] != b',' && bytes[idx] != b';' { - let c = bytes[idx]; - assert!(n_mods_so_far < SIZE); - mod_codes[n_mods_so_far] = ModCodeRepr::Code(c as char); - n_mods_so_far += 1; idx += 1; - mods_parsed += 1; } - return (mods_parsed, idx); + (ParsedModCodes::Alphabetic(&bytes[start..idx]), idx) } else { - assert!(n_mods_so_far < SIZE); - mod_codes[n_mods_so_far] = ModCodeRepr::ChEbi(val); - (1usize, idx) + (ParsedModCodes::ChEbi(val), idx) } } @@ -585,6 +585,38 @@ mod base_mods_adapter_tests { record } + fn make_forward_record_with_mn( + mm: &str, + ml: &[u8], + seq: &str, + ) -> bam::Record { + let mut record = make_record(mm, ml, seq, None, false); + record.push_aux(b"MN", Aux::I32(seq.len() as i32)).unwrap(); + record + } + + fn collect_states( + scanner: &mut BaseModsAdapter, + ) -> Vec<(usize, ModCodeRepr, u8, bool, bool)> { + std::iter::from_fn(|| { + scanner.next_modified_position_no_thresh().unwrap() + }) + .map(|state| { + ( + state.mod_position, + state.mod_code, + state.mod_qual, + state.modified, + state.inferred, + ) + }) + .collect() + } + + fn alphabetic_codes(count: usize) -> String { + (b'a'..).take(count).map(char::from).collect() + } + #[test] fn test_reverse() { let mm = "C+m.,0,0;C+h.,0,0;"; @@ -803,7 +835,14 @@ mod base_mods_adapter_tests { let record = make_record(mm, &ml, seq, None, true); let thresholds = [0f32; 4]; let mut scanner = BaseModsAdapter::<2>::new(&record).unwrap(); - assert_eq!(scanner.ml_pos, [2, 3]); + assert_eq!( + scanner + .code_states + .iter() + .map(|state| state.ml_pos) + .collect::>(), + vec![2, 3] + ); let mod_state = scanner.next_modified_position(thresholds, &[]).unwrap().unwrap(); assert_eq!(mod_state.mod_position, 0); @@ -966,6 +1005,188 @@ mod base_mods_adapter_tests { assert_eq!(mod_state.mod_code, ANY_CYTOSINE); } + #[test] + fn test_single_group_grows_beyond_inline_code_capacity() { + for count in [16, 17] { + let codes = alphabetic_codes(count); + let mm = format!("C+{codes}?,0;"); + let mut ml = vec![0; count]; + ml[count - 1] = 200; + let record = make_forward_record_with_mn(&mm, &ml, "C"); + + let mut scanner = BaseModsAdapter::<16>::new(&record).unwrap(); + assert_eq!( + scanner + .code_states + .iter() + .map(|state| state.mod_code) + .collect::>(), + codes + .bytes() + .map(|code| ModCodeRepr::Code(code as char)) + .collect::>() + ); + assert!(scanner + .code_states + .iter() + .all(|state| state.ml_stride == count)); + assert_eq!(scanner.code_states.spilled(), count > 16); + let expected_code = + ModCodeRepr::Code(codes.as_bytes()[count - 1] as char); + assert_eq!( + collect_states(&mut scanner), + vec![(0, expected_code, 200, true, false)] + ); + } + } + + #[test] + fn test_spilled_single_group_reverse_cursor_and_ml_stride() { + let codes = alphabetic_codes(17); + let mm = format!("C+{codes}?,0,0;"); + let mut ml = vec![0; 34]; + ml[0] = 200; + ml[33] = 220; + let mut record = make_record(&mm, &ml, "GG", None, true); + record.push_aux(b"MN", Aux::I32(2)).unwrap(); + + let mut scanner = BaseModsAdapter::<16>::new(&record).unwrap(); + assert!(scanner.code_states.spilled()); + assert_eq!(scanner.code_states.len(), 17); + assert!(scanner.code_states.iter().all(|state| state.ml_stride == 17)); + assert_eq!( + scanner + .code_states + .iter() + .map(|state| state.ml_pos) + .collect::>(), + (17..34).collect::>() + ); + let initial_mm_pos = scanner.code_states[0].mm_pos; + + let state = + scanner.next_modified_position_no_thresh().unwrap().unwrap(); + assert_eq!( + ( + state.mod_position, + state.mod_code, + state.mod_qual, + state.modified, + state.inferred, + ), + (0, ModCodeRepr::Code('q'), 220, true, false) + ); + assert!(scanner.code_states.iter().all(|state| { + state.mm_pos == initial_mm_pos - 2 + && state.mm_next == Some(0) + && state.num_explicit_positions == 1 + })); + assert_eq!( + scanner + .code_states + .iter() + .map(|state| state.ml_pos) + .collect::>(), + (0..17).collect::>() + ); + + let state = + scanner.next_modified_position_no_thresh().unwrap().unwrap(); + assert_eq!( + ( + state.mod_position, + state.mod_code, + state.mod_qual, + state.modified, + state.inferred, + ), + (1, ModCodeRepr::Code('a'), 200, true, false) + ); + assert!(scanner.next_modified_position_no_thresh().unwrap().is_none()); + } + + #[test] + fn test_one_code_groups_grow_beyond_inline_capacity() { + for count in [16, 17] { + // Separate MM groups for the same canonical base must target + // different positions to avoid conflicting calls. + let codes = + (0..count).map(|idx| 10_000 + idx as u32).collect::>(); + let mm = codes + .iter() + .enumerate() + .map(|(idx, code)| format!("C+{code}?,{idx};")) + .collect::(); + let ml = vec![200; count]; + let seq = "C".repeat(count); + let record = make_forward_record_with_mn(&mm, &ml, &seq); + + let mut scanner = BaseModsAdapter::<16>::new(&record).unwrap(); + assert_eq!( + scanner + .code_states + .iter() + .map(|state| state.mod_code) + .collect::>(), + codes + .iter() + .copied() + .map(ModCodeRepr::ChEbi) + .collect::>() + ); + assert_eq!( + scanner + .code_states + .iter() + .map(|state| state.mm_next) + .collect::>(), + (0..count).map(|idx| Some(idx as u32)).collect::>() + ); + assert_eq!( + scanner + .code_states + .iter() + .map(|state| state.ml_pos) + .collect::>(), + (0..count).collect::>() + ); + assert!(scanner + .code_states + .iter() + .all(|state| state.ml_stride == 1)); + assert_eq!(scanner.code_states.spilled(), count > 16); + let expected = codes + .iter() + .enumerate() + .map(|(position, code)| { + (position, ModCodeRepr::ChEbi(*code), 200, true, false) + }) + .collect::>(); + assert_eq!(collect_states(&mut scanner), expected); + } + } + + #[test] + fn test_explicit_position_count_is_not_code_capacity() { + let seq = "C".repeat(17); + let mm = format!("C+m?,{};", vec!["0"; 17].join(",")); + let ml = (200..217).collect::>(); + let record = make_forward_record_with_mn(&mm, &ml, &seq); + + let mut scanner = BaseModsAdapter::<16>::new(&record).unwrap(); + assert_eq!(scanner.code_states.len(), 1); + assert_eq!(scanner.code_states[0].num_explicit_positions, 17); + assert!(!scanner.code_states.spilled()); + let expected = ml + .iter() + .enumerate() + .map(|(position, qual)| { + (position, METHYL_CYTOSINE, *qual, true, false) + }) + .collect::>(); + assert_eq!(collect_states(&mut scanner), expected); + } + #[test] fn test_chebi() { let mm = "C+m.,0,0;C+76792.,0,0;"; @@ -974,8 +1195,12 @@ mod base_mods_adapter_tests { let record = make_record(mm, &ml, seq, None, true); let scanner = BaseModsAdapter::<2>::new(&record).unwrap(); assert_eq!( - scanner.mod_codes, - [METHYL_CYTOSINE, ModCodeRepr::ChEbi(76792)] + scanner + .code_states + .iter() + .map(|state| state.mod_code) + .collect::>(), + vec![METHYL_CYTOSINE, ModCodeRepr::ChEbi(76792)] ) } From 9fb9aea763aa1b78ac1172fa6d6734e7955ccc70 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Mon, 10 Aug 2026 17:22:50 -0700 Subject: [PATCH 3/3] Test independent MM groups through pileup --- modkit/tests/test_pileup.rs | 137 +++++++++++++++++++++++++++++++++++- 1 file changed, 134 insertions(+), 3 deletions(-) diff --git a/modkit/tests/test_pileup.rs b/modkit/tests/test_pileup.rs index 68f2cac..8f42c4a 100644 --- a/modkit/tests/test_pileup.rs +++ b/modkit/tests/test_pileup.rs @@ -1,11 +1,16 @@ use anyhow::Context; use itertools::Itertools; -use rust_htslib::bam; +use rust_htslib::bam::{ + self, + header::HeaderRecord, + record::{Aux, Cigar, CigarString}, + Format, Header, Record, Writer as BamWriter, +}; use std::cmp::Ordering; use std::collections::HashMap; use std::fs::File; -use std::io::{BufRead, BufReader}; -use std::path::PathBuf; +use std::io::{BufRead, BufReader, Write}; +use std::path::{Path, PathBuf}; use common::{check_against_expected_text_file, run_modkit}; use mod_kit::dmr::bedmethyl::BedMethylLine; @@ -13,6 +18,132 @@ use mod_kit::mod_base_code::{ModCodeRepr, METHYL_CYTOSINE}; mod common; +fn write_independent_mm_group_fixture(root: &Path) -> (PathBuf, PathBuf) { + let bam_path = root.join("independent-mm-groups.bam"); + let fasta_path = root.join("reference.fa"); + + let mut header = Header::new(); + let mut hd = HeaderRecord::new(b"HD"); + hd.push_tag(b"VN", "1.6").push_tag(b"SO", "coordinate"); + header.push_record(&hd); + let mut sq = HeaderRecord::new(b"SQ"); + sq.push_tag(b"SN", "chr1").push_tag(b"LN", 2); + header.push_record(&sq); + + let mut record = Record::new(); + record.set( + b"independent-groups", + Some(&CigarString(vec![Cigar::Match(2)])), + b"CC", + &[30; 2], + ); + record.set_tid(0); + record.set_pos(0); + record.set_mapq(60); + record.push_aux(b"MM", Aux::String("C+m?,0;C+h?,1;")).unwrap(); + record.push_aux(b"ML", Aux::ArrayU8((&[200, 250][..]).into())).unwrap(); + record.push_aux(b"MN", Aux::U32(2)).unwrap(); + record.push_aux(b"NM", Aux::U32(0)).unwrap(); + + let mut writer = + BamWriter::from_path(&bam_path, &header, Format::Bam).unwrap(); + writer.write(&record).unwrap(); + drop(writer); + bam::index::build(&bam_path, None, bam::index::Type::Bai, 1).unwrap(); + + File::create(&fasta_path).unwrap().write_all(b">chr1\nCC\n").unwrap(); + File::create(root.join("reference.fa.fai")) + .unwrap() + .write_all(b"chr1\t2\t6\t2\t3\n") + .unwrap(); + + (bam_path, fasta_path) +} + +fn collect_nonzero_modified_events( + output_path: &Path, +) -> Vec<(u32, u32, char, u16, String, u16, u16, u16)> { + BufReader::new(File::open(output_path).unwrap()) + .lines() + .map(|line| line.unwrap()) + .filter_map(|line| { + let fields = line.split('\t').collect::>(); + let n_modified = fields[11].parse::().unwrap(); + (n_modified > 0).then(|| { + ( + fields[1].parse::().unwrap(), + fields[2].parse::().unwrap(), + fields[3] + .split(',') + .next() + .unwrap() + .parse::() + .unwrap(), + fields[9].parse::().unwrap(), + fields[10].to_string(), + n_modified, + fields[12].parse::().unwrap(), + fields[13].parse::().unwrap(), + ) + }) + }) + .collect() +} + +#[test] +fn test_pileup_independent_same_base_mm_groups_match_generic() { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + let (bam_path, fasta_path) = write_independent_mm_group_fixture(root); + let optimized_path = root.join("optimized.bed"); + let generic_path = root.join("generic.bed"); + + run_modkit(&[ + "pileup", + bam_path.to_str().unwrap(), + optimized_path.to_str().unwrap(), + "--ref", + fasta_path.to_str().unwrap(), + "--modified-bases", + "C:m", + "C:h", + "--no-filtering", + "--threads", + "1", + "--io-threads", + "1", + "--suppress-progress", + ]) + .unwrap(); + run_modkit(&[ + "pileup", + bam_path.to_str().unwrap(), + generic_path.to_str().unwrap(), + "--ref", + fasta_path.to_str().unwrap(), + "--motif", + "CC", + "0", + "--motif", + "CC", + "1", + "--no-filtering", + "--threads", + "1", + "--io-threads", + "1", + "--suppress-progress", + ]) + .unwrap(); + + let expected = vec![ + (0, 1, 'm', 1, "100.00".to_string(), 1, 0, 0), + (1, 2, 'h', 1, "100.00".to_string(), 1, 0, 0), + ]; + assert_eq!(collect_nonzero_modified_events(&optimized_path), expected); + assert_eq!(collect_nonzero_modified_events(&generic_path), expected); +} + #[test] fn test_pileup_help() { let pileup_help_args = ["pileup", "--help"];