diff --git a/modkit-core/src/pileup/pileup_processor.rs b/modkit-core/src/pileup/pileup_processor.rs index 262371a..c7c7204 100644 --- a/modkit-core/src/pileup/pileup_processor.rs +++ b/modkit-core/src/pileup/pileup_processor.rs @@ -18,8 +18,7 @@ use crate::{ MmTagInfo, }, mod_base_code::{ - DnaBase, ModCodeRepr, ANY_CYTOSINE, HYDROXY_METHYL_CYTOSINE, - METHYL_CYTOSINE, SIX_METHYL_ADENINE, + DnaBase, ModCodeRepr, ANY_CYTOSINE, METHYL_CYTOSINE, SIX_METHYL_ADENINE, }, motifs::motif_bed::MotifInfo, pileup::{ @@ -257,10 +256,6 @@ impl< "bs:{bs:?} st:{st} end:{end} \ num_motifs:{num_motifs}" ); - #[cfg(debug_assertions)] - { - assert!(bs.count_ones() <= 1); - } if bs[0] { Some((qpos, rpos, self.motif_bases[0])) } else if *num_motifs >= 2u8 && bs[1] { @@ -1012,11 +1007,63 @@ const DYN_CAN_C: usize = 5usize; const DYN_CAN_G: usize = 6usize; const DYN_CAN_T: usize = 7usize; const DYN_OTHER_MOD_A: usize = 8usize; -// const DYN_OTHER_MOD_C: usize = 9usize; -// const DYN_OTHER_MOD_G: usize = 10usize; -// const DYN_OTHER_MOD_T: usize = 11usize; +const DYN_OTHER_MOD_C: usize = 9usize; +const DYN_OTHER_MOD_G: usize = 10usize; +const DYN_OTHER_MOD_T: usize = 11usize; const DYN_N_CONSTANT_COUNTS: usize = 12usize; +fn dynamic_canonical_offset(canonical_base: DnaBase) -> usize { + match canonical_base { + DnaBase::A => DYN_CAN_A, + DnaBase::C => DYN_CAN_C, + DnaBase::G => DYN_CAN_G, + DnaBase::T => DYN_CAN_T, + } +} + +fn dynamic_other_mod_offset(canonical_base: DnaBase) -> usize { + match canonical_base { + DnaBase::A => DYN_OTHER_MOD_A, + DnaBase::C => DYN_OTHER_MOD_C, + DnaBase::G => DYN_OTHER_MOD_G, + DnaBase::T => DYN_OTHER_MOD_T, + } +} + +#[inline] +fn dynamic_anchor_rpos(rpos: u32, reverse: bool, motif_offset: u32) -> u32 { + if reverse { + rpos.saturating_sub(motif_offset) + } else { + rpos + } +} + +fn dynamic_mod_offset( + mod_codes: &[(DnaBase, ModCodeRepr)], + canonical_base: DnaBase, + mod_code: ModCodeRepr, + combine_mods: bool, +) -> usize { + let compact_slot = if combine_mods { + mod_codes.iter().position(|(base, _)| *base == canonical_base) + } else { + mod_codes.iter().position(|(base, code)| { + *base == canonical_base && *code == mod_code + }) + }; + let mod_offset = compact_slot + .map(|slot| DYN_N_CONSTANT_COUNTS + slot) + .unwrap_or_else(|| dynamic_other_mod_offset(canonical_base)); + let row_stride = DYN_N_CONSTANT_COUNTS + mod_codes.len(); + assert!( + mod_offset < row_stride, + "dynamic modification offset {mod_offset} exceeds row stride \ + {row_stride}" + ); + mod_offset +} + pub(super) trait ACountsMatrix { fn new( @@ -1460,7 +1507,7 @@ impl ACountsMatrix for CountsMatrix { &mut self, rpos: u32, canonical_base: DnaBase, - mod_code: ModCodeRepr, + _mod_code: ModCodeRepr, reverse: bool, haplotype: u8, _combine_mods: bool, @@ -1473,12 +1520,8 @@ impl ACountsMatrix for CountsMatrix { ); match canonical_base { DnaBase::C => { - if mod_code == METHYL_CYTOSINE - || mod_code == HYDROXY_METHYL_CYTOSINE - { - self.inner[offset + METH_C_OFFSET] = - self.inner[offset + METH_C_OFFSET].saturating_add(1); - } + self.inner[offset + METH_C_OFFSET] = + self.inner[offset + METH_C_OFFSET].saturating_add(1); } _ => {} } @@ -1501,7 +1544,7 @@ impl ACountsMatrix for CountsMatrix { match canonical_base { DnaBase::C => { self.inner[offset + FILT_C_OFFSET] = - self.inner[offset + FILT_C_OFFSET]; + self.inner[offset + FILT_C_OFFSET].saturating_add(1); } _ => {} } @@ -1821,6 +1864,7 @@ impl return; } + let rpos = dynamic_anchor_rpos(rpos, reverse, self.motif_offset); let offset = calc_offset_dyn::( rpos, self.strand_width, @@ -1844,8 +1888,7 @@ impl if base != reference_base { return; } - let rpos = - if reverse { rpos.saturating_sub(self.motif_offset) } else { rpos }; + let rpos = dynamic_anchor_rpos(rpos, reverse, self.motif_offset); let offset = calc_offset_dyn::( rpos, self.strand_width, @@ -1853,12 +1896,7 @@ impl reverse, haplotype, ); - let base_offset = match base { - DnaBase::A => DYN_CAN_A, - DnaBase::C => DYN_CAN_C, - DnaBase::G => DYN_CAN_G, - DnaBase::T => DYN_CAN_T, - }; + let base_offset = dynamic_canonical_offset(base); assert!(offset + base_offset < self.inner.len()); self.inner[offset + base_offset] = self.inner[offset + base_offset].saturating_add(1); @@ -1873,8 +1911,7 @@ impl haplotype: u8, combine_mods: bool, ) { - let rpos = - if reverse { rpos.saturating_sub(self.motif_offset) } else { rpos }; + let rpos = dynamic_anchor_rpos(rpos, reverse, self.motif_offset); let offset = calc_offset_dyn::( rpos, self.strand_width, @@ -1882,15 +1919,12 @@ impl reverse, haplotype, ); - let mod_offset = if combine_mods { - canonical_base as usize + DYN_N_CONSTANT_COUNTS - } else { - self.mod_codes - .iter() - .position(|(_b, c)| c == &mod_code) - .map(|p| p + DYN_N_CONSTANT_COUNTS) - .unwrap_or(DYN_OTHER_MOD_A + canonical_base as usize) - }; + let mod_offset = dynamic_mod_offset( + &self.mod_codes, + canonical_base, + mod_code, + combine_mods, + ); assert!(offset + mod_offset < self.inner.len(), "STRANDS {STRANDS}"); self.inner[offset + mod_offset] = self.inner[offset + mod_offset].saturating_add(1); @@ -1903,6 +1937,7 @@ impl reverse: bool, haplotype: u8, ) { + let rpos = dynamic_anchor_rpos(rpos, reverse, self.motif_offset); let offset = calc_offset_dyn::( rpos, self.strand_width, @@ -1926,6 +1961,7 @@ impl if >::reached_max_depth(self, rpos, reverse, haplotype) { return; } + let rpos = dynamic_anchor_rpos(rpos, reverse, self.motif_offset); let offset = calc_offset_dyn::( rpos, self.strand_width, @@ -1951,8 +1987,7 @@ impl reverse: bool, haplotype: u8, ) -> bool { - let rpos = - if reverse { rpos.saturating_sub(self.motif_offset) } else { rpos }; + let rpos = dynamic_anchor_rpos(rpos, reverse, self.motif_offset); let offset = calc_offset_dyn::( rpos, self.strand_width, @@ -2018,9 +2053,10 @@ fn slice_to_counts_dynamic<'a>( mod_codes.iter().enumerate().map( move |(j, (primary_base, code))| { let n_modified = chunk[DYN_N_CONSTANT_COUNTS + j]; - let n_canonical = chunk[DYN_CAN_A + *primary_base as usize]; + let n_canonical = + chunk[dynamic_canonical_offset(*primary_base)]; let n_anon_modified = - chunk[DYN_OTHER_MOD_A + *primary_base as usize]; + chunk[dynamic_other_mod_offset(*primary_base)]; let n_other_modified = mod_codes .iter() .enumerate() @@ -2507,3 +2543,388 @@ fn indices_to_byte(motif_idxs: &BitSlice) -> u8 { } agg } + +#[cfg(test)] +mod tests { + use crate::{ + mod_base_code::{ + DnaBase, ModCodeRepr, ANY_CYTOSINE, FORMYL_CYTOSINE, + HYDROXY_METHYL_CYTOSINE, METHYL_CYTOSINE, TWO_OME_CYTOSINE, + }, + pileup::base_mods_adapter::ModState, + }; + + use super::{ + ACountsMatrix, CountsMatrix, DnaCytosineCombine, DnaModOption, Dynamic, + }; + + type DynamicCounts = + (u32, char, u16, ModCodeRepr, u16, u16, u16, u16, u16, u16, u16, u8); + + fn increment_dynamic_call( + matrix: &mut CountsMatrix, + rpos: u32, + primary_base: DnaBase, + reference_base: DnaBase, + mod_code: ModCodeRepr, + modified: bool, + filtered: bool, + reverse: bool, + haplotype: u8, + ) { + >::incr_call( + matrix, + ModState { + mod_position: rpos as usize, + modified, + filtered, + mod_code, + primary_base, + inferred: false, + mod_qual: 255, + }, + rpos, + reference_base, + reverse, + haplotype, + false, + ); + } + + fn valid_dynamic_counts( + counts: &[crate::pileup::PileupFeatureCounts2], + ) -> Vec { + counts + .iter() + .filter(|counts| counts.is_valid()) + .map(|counts| { + ( + counts.position, + counts.raw_strand, + counts.filtered_coverage, + counts.mod_code, + counts.n_canonical, + counts.n_modified, + counts.n_other_modified, + counts.n_delete, + counts.n_filtered, + counts.n_diff, + counts.n_nocall, + counts.motif_idxs, + ) + }) + .collect() + } + + fn new_cytosine_combine_matrix() -> CountsMatrix { + >::new( + 1, + false, + Vec::new(), + 0, + u16::MAX, + ) + } + + fn increment_cytosine_call( + matrix: &mut CountsMatrix, + modified: bool, + filtered: bool, + mod_code: ModCodeRepr, + ) { + >::incr_call( + matrix, + ModState { + mod_position: 0, + modified, + filtered, + mod_code, + primary_base: DnaBase::C, + inferred: false, + mod_qual: 255, + }, + 0, + DnaBase::C, + false, + 0, + true, + ); + } + + fn decode_cytosine_combine( + matrix: &CountsMatrix, + ) -> Vec { + >::decode( + matrix, + 100, + DnaModOption::Combine, + false, + ) + } + + #[test] + fn dynamic_explicit_mod_offsets_are_keyed_by_base_and_code() { + let mod_codes = + [(DnaBase::A, METHYL_CYTOSINE), (DnaBase::C, METHYL_CYTOSINE)]; + + assert_eq!( + super::dynamic_mod_offset( + &mod_codes, + DnaBase::A, + METHYL_CYTOSINE, + false, + ), + super::DYN_N_CONSTANT_COUNTS + ); + assert_eq!( + super::dynamic_mod_offset( + &mod_codes, + DnaBase::C, + METHYL_CYTOSINE, + false, + ), + super::DYN_N_CONSTANT_COUNTS + 1 + ); + assert_eq!( + super::dynamic_mod_offset( + &mod_codes, + DnaBase::T, + METHYL_CYTOSINE, + false, + ), + super::DYN_OTHER_MOD_T + ); + } + + #[test] + fn dynamic_explicit_slots_preserve_statuses_across_strand_and_phase() { + let mod_codes = + vec![(DnaBase::A, METHYL_CYTOSINE), (DnaBase::C, METHYL_CYTOSINE)]; + let mut matrix = >::new( + 2, + true, + mod_codes, + 0, + u16::MAX, + ); + + // The unphased partition proves that identical modification codes on + // different primary bases retain independent requested slots. + increment_dynamic_call( + &mut matrix, + 0, + DnaBase::A, + DnaBase::A, + METHYL_CYTOSINE, + true, + false, + false, + 0, + ); + increment_dynamic_call( + &mut matrix, + 1, + DnaBase::C, + DnaBase::C, + METHYL_CYTOSINE, + true, + false, + false, + 0, + ); + + // Exercise every dynamic status in the reverse HP1 partition. A zero + // motif offset intentionally isolates slot mapping from motif shifts. + for (mod_code, modified, filtered) in [ + (METHYL_CYTOSINE, true, false), + (METHYL_CYTOSINE, false, false), + (HYDROXY_METHYL_CYTOSINE, true, false), + (METHYL_CYTOSINE, false, true), + ] { + increment_dynamic_call( + &mut matrix, + 1, + DnaBase::C, + DnaBase::C, + mod_code, + modified, + filtered, + true, + 1, + ); + } + >::incr_diff_call( + &mut matrix, + 1, + DnaBase::A, + DnaBase::C, + true, + 1, + ); + >::incr_diff_call( + &mut matrix, + 1, + DnaBase::C, + DnaBase::C, + true, + 1, + ); + >::incr_delete( + &mut matrix, + 1, + true, + 1, + ); + + // Repeat the same status oracle for adenine in the forward HP2 + // partition, including an unrequested modification code. + for (mod_code, modified, filtered) in [ + (METHYL_CYTOSINE, true, false), + (METHYL_CYTOSINE, false, false), + (ModCodeRepr::Code('x'), true, false), + (METHYL_CYTOSINE, false, true), + ] { + increment_dynamic_call( + &mut matrix, + 0, + DnaBase::A, + DnaBase::A, + mod_code, + modified, + filtered, + false, + 2, + ); + } + >::incr_diff_call( + &mut matrix, + 0, + DnaBase::C, + DnaBase::A, + false, + 2, + ); + >::incr_diff_call( + &mut matrix, + 0, + DnaBase::A, + DnaBase::A, + false, + 2, + ); + >::incr_delete( + &mut matrix, + 0, + false, + 2, + ); + + let decoded = >::decode( + &matrix, + 100, + DnaModOption::Other(ANY_CYTOSINE), + true, + ); + let partition_width = 2 * 2 * 2; + assert_eq!(decoded.len(), partition_width * 3); + assert_eq!( + valid_dynamic_counts(&decoded[..partition_width]), + vec![ + (100, '+', 1, METHYL_CYTOSINE, 0, 1, 0, 0, 0, 0, 0, 0), + (101, '+', 1, METHYL_CYTOSINE, 0, 1, 0, 0, 0, 0, 0, 0), + ] + ); + assert_eq!( + valid_dynamic_counts( + &decoded[partition_width..partition_width * 2] + ), + vec![(101, '-', 3, METHYL_CYTOSINE, 1, 1, 1, 1, 1, 1, 1, 0,)] + ); + assert_eq!( + valid_dynamic_counts(&decoded[partition_width * 2..]), + vec![(100, '+', 3, METHYL_CYTOSINE, 1, 1, 1, 1, 1, 1, 1, 0,)] + ); + } + + #[test] + fn dna_cytosine_combine_counts_every_c_modification() { + let mut matrix = new_cytosine_combine_matrix(); + increment_cytosine_call(&mut matrix, false, false, ANY_CYTOSINE); + for mod_code in [ + METHYL_CYTOSINE, + HYDROXY_METHYL_CYTOSINE, + FORMYL_CYTOSINE, + TWO_OME_CYTOSINE, + ] { + increment_cytosine_call(&mut matrix, true, false, mod_code); + } + + let decoded = decode_cytosine_combine(&matrix); + let counts = &decoded[0]; + assert_eq!( + ( + counts.position, + counts.raw_strand, + counts.filtered_coverage, + counts.mod_code, + counts.n_canonical, + counts.n_modified, + counts.n_other_modified, + counts.n_delete, + counts.n_filtered, + counts.n_diff, + counts.n_nocall, + counts.motif_idxs, + ), + (100, '+', 5, ANY_CYTOSINE, 1, 4, 0, 0, 0, 0, 0, 0) + ); + assert_eq!( + counts.filtered_coverage, + counts + .n_canonical + .saturating_add(counts.n_modified) + .saturating_add(counts.n_other_modified) + ); + } + + #[test] + fn dna_cytosine_combine_counts_filtered_c_calls_saturating() { + let mut matrix = new_cytosine_combine_matrix(); + increment_cytosine_call(&mut matrix, false, false, ANY_CYTOSINE); + increment_cytosine_call(&mut matrix, false, true, ANY_CYTOSINE); + let decoded = decode_cytosine_combine(&matrix); + let counts = &decoded[0]; + assert_eq!( + ( + counts.position, + counts.raw_strand, + counts.filtered_coverage, + counts.mod_code, + counts.n_canonical, + counts.n_modified, + counts.n_other_modified, + counts.n_delete, + counts.n_filtered, + counts.n_diff, + counts.n_nocall, + counts.motif_idxs, + ), + (100, '+', 1, ANY_CYTOSINE, 1, 0, 0, 0, 1, 0, 0, 0) + ); + assert_eq!( + counts.filtered_coverage, + counts + .n_canonical + .saturating_add(counts.n_modified) + .saturating_add(counts.n_other_modified) + ); + + matrix.inner[super::FILT_C_OFFSET] = u16::MAX; + increment_cytosine_call(&mut matrix, false, true, ANY_CYTOSINE); + assert_eq!(decode_cytosine_combine(&matrix)[0].n_filtered, u16::MAX); + } +} diff --git a/modkit-core/src/pileup/subcommand.rs b/modkit-core/src/pileup/subcommand.rs index 4503f99..bf10577 100644 --- a/modkit-core/src/pileup/subcommand.rs +++ b/modkit-core/src/pileup/subcommand.rs @@ -424,6 +424,31 @@ pub struct ModBamPileup { } impl ModBamPileup { + fn normalized_modified_base_selections( + &self, + ) -> Option> { + self.modified_bases.as_ref().map(|modified_bases| { + let mut seen = HashSet::with_capacity(modified_bases.len()); + modified_bases + .iter() + .filter_map(|option| { + let selection = (option.primary_base, option.mod_code); + if seen.insert(selection) { + Some(option.clone()) + } else { + warn!( + "ignoring duplicate --modified-bases selection \ + '{}:{}'", + option.primary_base.char(), + option.mod_code + ); + None + } + }) + .collect() + }) + } + fn parse_user_motifs(&self) -> Option>> { if let Some(raw_motif_parts) = &self.motif { Some(parse_raw_motifs(raw_motif_parts, self.cpg, false)) @@ -540,6 +565,7 @@ impl ModBamPileup { fn determine_preset( &self, + modified_base_selections: Option<&[ModifiedBasesOptions]>, ) -> anyhow::Result<(Option, Option>)> { let mut regex_motifs = self.parse_user_motifs().transpose()?.unwrap_or_else(Vec::new); @@ -557,7 +583,7 @@ impl ModBamPileup { .iter() .map(|mot| mot.motif_info.primary_base) .collect::>(); - if let Some(modified_bases) = self.modified_bases.as_ref() { + if let Some(modified_bases) = modified_base_selections { for primary_base in modified_bases.iter().map(|x| x.primary_base) { @@ -578,7 +604,7 @@ impl ModBamPileup { } return Ok((None, Some(regex_motifs))); } - if let Some(modified_bases) = self.modified_bases.as_ref() { + if let Some(modified_bases) = modified_base_selections { let modified_bases = modified_bases .iter() .map(|x| (x.primary_base, x.mod_code)) @@ -930,6 +956,9 @@ impl ModBamPileup { .set_draw_target(indicatif::ProgressDrawTarget::stderr()); } + let modified_base_selections = + self.normalized_modified_base_selections(); + let io_thread_pool = rust_htslib::tpool::ThreadPool::new(self.io_threads)?; @@ -980,7 +1009,7 @@ impl ModBamPileup { .transpose()?; let reference_records = get_targets(&header, region.as_ref()); if self.include_bed.is_some() { - if !(self.modified_bases.is_some() || self.motif.is_some()) { + if !(modified_base_selections.is_some() || self.motif.is_some()) { bail!( "currently, --include-bed requires a --motif or \ --modified-bases. This limitation may be removed in the \ @@ -1022,7 +1051,8 @@ impl ModBamPileup { reference_records }; - let (preset, regex_motifs) = self.determine_preset()?; + let (preset, regex_motifs) = + self.determine_preset(modified_base_selections.as_deref())?; let (pileup_options, combine_strands) = match &preset { Some(preset) => match preset { @@ -1102,7 +1132,7 @@ impl ModBamPileup { self.with_header, &self.bedrmodargs, &header, - self.modified_bases.as_ref(), + modified_base_selections.as_ref(), empties_tx.clone(), master_progress.clone(), )?) @@ -1112,7 +1142,7 @@ impl ModBamPileup { self.with_header, &self.bedrmodargs, &header, - self.modified_bases.as_ref(), + modified_base_selections.as_ref(), master_progress.clone(), empties_tx.clone(), )?) @@ -1133,7 +1163,7 @@ impl ModBamPileup { self.with_header, &self.bedrmodargs, &header, - self.modified_bases.as_ref(), + modified_base_selections.as_ref(), master_progress.clone(), empties_tx.clone(), self.bgzf_threads, @@ -1144,7 +1174,7 @@ impl ModBamPileup { self.with_header, &self.bedrmodargs, &header, - self.modified_bases.as_ref(), + modified_base_selections.as_ref(), master_progress.clone(), empties_tx.clone(), )?) @@ -1161,7 +1191,7 @@ impl ModBamPileup { self.with_header, &self.bedrmodargs, &header, - self.modified_bases.as_ref(), + modified_base_selections.as_ref(), master_progress.clone(), empties_tx.clone(), self.bgzf_threads, @@ -1173,7 +1203,7 @@ impl ModBamPileup { self.with_header, &self.bedrmodargs, &header, - self.modified_bases.as_ref(), + modified_base_selections.as_ref(), master_progress.clone(), empties_tx.clone(), )?) diff --git a/modkit/tests/test_pileup.rs b/modkit/tests/test_pileup.rs index 68f2cac..c845116 100644 --- a/modkit/tests/test_pileup.rs +++ b/modkit/tests/test_pileup.rs @@ -1,11 +1,14 @@ use anyhow::Context; use itertools::Itertools; use rust_htslib::bam; +use rust_htslib::bam::header::HeaderRecord; +use rust_htslib::bam::record::{Aux, Cigar, CigarString}; +use rust_htslib::bam::{Format, Header, Record, Writer as BamWriter}; use std::cmp::Ordering; -use std::collections::HashMap; +use std::collections::{HashMap, HashSet}; 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 +16,332 @@ use mod_kit::mod_base_code::{ModCodeRepr, METHYL_CYTOSINE}; mod common; +fn read_bed_sites(path: &str) -> HashSet<(String, u32, char)> { + BufReader::new(File::open(path).unwrap()) + .lines() + .map(|line| { + let line = line.unwrap(); + let fields = line.split('\t').collect::>(); + ( + fields[0].to_string(), + fields[1].parse::().unwrap(), + fields[5].parse::().unwrap(), + ) + }) + .collect() +} + +fn write_dynamic_slot_fixture_with_probs( + root: &Path, + name: &str, + probabilities: &[[u8; 3]], +) -> (PathBuf, PathBuf) { + let bam_path = root.join(format!("{name}.bam")); + let fasta_path = root.join("reference.fa"); + + let mut header = Header::new(); + let mut sq = HeaderRecord::new(b"SQ"); + sq.push_tag(b"SN", "chr1"); + sq.push_tag(b"LN", 3); + header.push_record(&sq); + + let mut writer = + BamWriter::from_path(&bam_path, &header, Format::Bam).unwrap(); + for (i, probabilities) in probabilities.iter().enumerate() { + let cigar = CigarString(vec![Cigar::Match(3)]); + let mut record = Record::new(); + record.set( + format!("read-{i}").as_bytes(), + Some(&cigar), + b"ACT", + &[30; 3], + ); + record.set_tid(0); + record.set_pos(0); + record.set_mapq(60); + record.push_aux(b"MM", Aux::String("A+m?,0;C+m?,0;T+g?,0;")).unwrap(); + record + .push_aux(b"ML", Aux::ArrayU8((&probabilities[..]).into())) + .unwrap(); + record.push_aux(b"MN", Aux::U32(3)).unwrap(); + record.push_aux(b"NM", Aux::U32(0)).unwrap(); + writer.write(&record).unwrap(); + } + drop(writer); + bam::index::build(bam_path.clone(), None, bam::index::Type::Bai, 1) + .unwrap(); + + File::create(&fasta_path).unwrap().write_all(b">chr1\nACT\n").unwrap(); + File::create(root.join("reference.fa.fai")) + .unwrap() + .write_all(b"chr1\t3\t6\t3\t4\n") + .unwrap(); + + (bam_path, fasta_path) +} + +fn write_dynamic_slot_fixture(root: &Path) -> (PathBuf, PathBuf) { + write_dynamic_slot_fixture_with_probs( + root, + "dynamic-slots", + &[[255, 255, 255]], + ) +} + +fn write_reverse_motif_record( + writer: &mut BamWriter, + name: &str, + sequence: &[u8], + cigar: CigarString, + mm_tag: &str, + probability: u8, + edit_distance: u32, +) { + let qualities = vec![30; sequence.len()]; + let probabilities = [probability]; + let mut record = Record::new(); + record.set(name.as_bytes(), Some(&cigar), sequence, &qualities); + record.set_tid(0); + record.set_pos(0); + record.set_mapq(60); + record.set_flags(16); + record.push_aux(b"MM", Aux::String(mm_tag)).unwrap(); + record.push_aux(b"ML", Aux::ArrayU8((&probabilities[..]).into())).unwrap(); + record.push_aux(b"MN", Aux::U32(sequence.len() as u32)).unwrap(); + record.push_aux(b"NM", Aux::U32(edit_distance)).unwrap(); + writer.write(&record).unwrap(); +} + +fn write_combined_non_cpg_anchor_fixture(root: &Path) -> (PathBuf, PathBuf) { + let bam_path = root.join("combined-non-cpg-anchor.bam"); + let fasta_path = root.join("combined-non-cpg-anchor.fa"); + + let mut header = Header::new(); + let mut sq = HeaderRecord::new(b"SQ"); + sq.push_tag(b"SN", "chr1"); + sq.push_tag(b"LN", 4); + header.push_record(&sq); + + let mut writer = + BamWriter::from_path(&bam_path, &header, Format::Bam).unwrap(); + let full_match = || CigarString(vec![Cigar::Match(4)]); + + // GATC is reverse-complement palindromic. Focusing offset 1 makes the + // forward A anchor position 1 and the reverse T anchor position 2. + for (name, probability) in + [("modified", 255), ("canonical", 0), ("filtered", 128)] + { + write_reverse_motif_record( + &mut writer, + name, + b"GATC", + full_match(), + "A+m?,0;", + probability, + 0, + ); + } + write_reverse_motif_record( + &mut writer, + "no-call", + b"GATC", + full_match(), + "C+m?,0;", + 255, + 0, + ); + write_reverse_motif_record( + &mut writer, + "mismatch", + b"GACC", + full_match(), + "C+m?,0;", + 255, + 1, + ); + write_reverse_motif_record( + &mut writer, + "deletion", + b"GAC", + CigarString(vec![Cigar::Match(2), Cigar::Del(1), Cigar::Match(1)]), + "C+m?,0;", + 255, + 1, + ); + drop(writer); + bam::index::build(bam_path.clone(), None, bam::index::Type::Bai, 1) + .unwrap(); + + File::create(&fasta_path).unwrap().write_all(b">chr1\nGATC\n").unwrap(); + File::create(root.join("combined-non-cpg-anchor.fa.fai")) + .unwrap() + .write_all(b"chr1\t4\t6\t4\t5\n") + .unwrap(); + + (bam_path, fasta_path) +} + +fn run_combined_non_cpg_anchor_pileup( + root: &Path, + bam_path: &Path, + fasta_path: &Path, + high_depth: bool, + threads: usize, +) -> Vec { + let mode = if high_depth { "high" } else { "normal" }; + let output_path = + root.join(format!("combined-non-cpg-{mode}-{threads}.bed")); + let threads = threads.to_string(); + let mut args = vec![ + "pileup", + bam_path.to_str().unwrap(), + output_path.to_str().unwrap(), + "--ref", + fasta_path.to_str().unwrap(), + "--motif", + "GATC", + "1", + "--combine-strands", + "--modified-bases", + "A:m", + "--filter-threshold", + "0.9", + "--interval-size", + "4", + "--threads", + &threads, + "--suppress-progress", + ]; + if high_depth { + args.extend_from_slice(&["--high-depth", "--max-depth", "10"]); + } + + run_modkit(&args) + .with_context(|| { + format!( + "combined non-CpG pileup failed for mode {mode}, threads \ + {threads}" + ) + }) + .unwrap(); + std::fs::read(output_path).unwrap() +} + +fn run_combined_slot_pileup( + root: &Path, + bam_path: &Path, + fasta_path: &Path, + bases: &[&str], + high_depth: bool, + threads: usize, +) -> Vec { + let mode = if high_depth { "high" } else { "normal" }; + let output_path = + root.join(format!("combine-{}-{mode}-{threads}.bed", bases.join(""))); + let threads = threads.to_string(); + let mut args = vec![ + "pileup", + bam_path.to_str().unwrap(), + output_path.to_str().unwrap(), + "--ref", + fasta_path.to_str().unwrap(), + "--modified-bases", + ]; + args.extend_from_slice(bases); + args.extend_from_slice(&[ + "--combine-mods", + "--no-filtering", + "--interval-size", + "1", + "--threads", + &threads, + "--suppress-progress", + ]); + if high_depth { + args.extend_from_slice(&["--high-depth", "--max-depth", "10"]); + } + + run_modkit(&args) + .with_context(|| { + format!( + "combine pileup failed for bases {}, mode {mode}, threads {}", + bases.join(","), + threads + ) + }) + .unwrap(); + std::fs::read(output_path).unwrap() +} + +fn assert_combined_slot_parity(bases: &[&str], expected: &str) { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + let (bam_path, fasta_path) = write_dynamic_slot_fixture(root); + let baseline = + run_combined_slot_pileup(root, &bam_path, &fasta_path, bases, false, 1); + assert_eq!(baseline, expected.as_bytes()); + + for threads in [1, 2, 3, 8] { + for high_depth in [false, true] { + let observed = run_combined_slot_pileup( + root, + &bam_path, + &fasta_path, + bases, + high_depth, + threads, + ); + assert_eq!( + observed, + baseline, + "bases {}, high_depth {high_depth}, threads {threads}", + bases.join(",") + ); + } + } +} + +fn run_explicit_pair_slot_pileup( + root: &Path, + bam_path: &Path, + fasta_path: &Path, + high_depth: bool, + threads: usize, +) -> Vec { + let mode = if high_depth { "high" } else { "normal" }; + let output_path = root.join(format!("explicit-{mode}-{threads}.bed")); + let threads = threads.to_string(); + let mut args = vec![ + "pileup", + bam_path.to_str().unwrap(), + output_path.to_str().unwrap(), + "--ref", + fasta_path.to_str().unwrap(), + "--modified-bases", + "A:m", + "C:m", + "--no-filtering", + "--interval-size", + "1", + "--threads", + &threads, + "--suppress-progress", + ]; + if high_depth { + args.extend_from_slice(&["--high-depth", "--max-depth", "10"]); + } + + run_modkit(&args) + .with_context(|| { + format!( + "explicit pair pileup failed for mode {mode}, threads \ + {threads}" + ) + }) + .unwrap(); + std::fs::read(output_path).unwrap() +} + #[test] fn test_pileup_help() { let pileup_help_args = ["pileup", "--help"]; @@ -251,6 +580,166 @@ fn test_pileup_cpg_motif_filtering() { ); } +#[test] +fn test_pileup_cpg_combined_cytosine_debug_regression() { + let temp_file = std::env::temp_dir() + .join("test_pileup_cpg_combined_cytosine_debug_regression.bed"); + run_modkit(&[ + "pileup", + "../tests/resources/bc_anchored_10_reads.sorted.bam", + temp_file.to_str().unwrap(), + "--cpg", + "--modified-bases", + "C", + "--combine-mods", + "--no-filtering", + "--ref", + "../tests/resources/CGI_ladder_3.6kb_ref.fa", + "--threads", + "1", + ]) + .expect("valid overlapping CpG and internal C masks should not panic"); + + let observed = read_bed_sites(temp_file.to_str().unwrap()); + let expected = read_bed_sites( + "../tests/resources/bc_anchored_10_reads_nofilt_cg_motif.bed", + ); + assert!(!observed.is_empty()); + assert_eq!(observed, expected); +} + +#[test] +fn test_pileup_combined_c_compact_slot_matches_high_depth() { + assert_combined_slot_parity( + &["C"], + "chr1\t1\t2\tC\t1\t+\t1\t2\t255,0,0\t1\t100.00\t1\t0\t0\t0\t0\t0\t0\n", + ); +} + +#[test] +fn test_pileup_combined_act_compact_slots_match_high_depth() { + assert_combined_slot_parity( + &["A", "C", "T"], + concat!( + "chr1\t0\t1\tA\t1\t+\t0\t1\t255,0,0\t1\t100.00\t1\t0\t0\t0\t0\t0\t0\n", + "chr1\t1\t2\tC\t1\t+\t1\t2\t255,0,0\t1\t100.00\t1\t0\t0\t0\t0\t0\t0\n", + "chr1\t2\t3\tT\t1\t+\t2\t3\t255,0,0\t1\t100.00\t1\t0\t0\t0\t0\t0\t0\n", + ), + ); +} + +#[test] +fn test_pileup_explicit_same_code_slots_are_keyed_by_base() { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + let (bam_path, fasta_path) = write_dynamic_slot_fixture_with_probs( + root, + "explicit-pairs", + &[[255, 255, 255], [0, 0, 0]], + ); + let expected = concat!( + "chr1\t0\t1\tm\t2\t+\t0\t1\t255,0,0\t2\t50.00\t1\t1\t0\t0\t0\t0\t0\n", + "chr1\t1\t2\tm\t2\t+\t1\t2\t255,0,0\t2\t50.00\t1\t1\t0\t0\t0\t0\t0\n", + ); + + for threads in [1, 2, 3, 8] { + for high_depth in [false, true] { + let observed = run_explicit_pair_slot_pileup( + root, + &bam_path, + &fasta_path, + high_depth, + threads, + ); + assert_eq!( + observed, + expected.as_bytes(), + "high_depth {high_depth}, threads {threads}" + ); + } + } +} + +#[test] +fn test_pileup_duplicate_modified_base_selections_are_idempotent() { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + let (bam_path, fasta_path) = write_dynamic_slot_fixture(root); + for (mode, mode_args) in [ + ("optimized", &[][..]), + ("dynamic", &["--use-dynamic"][..]), + ("high-depth", &["--high-depth", "--max-depth", "100"][..]), + ] { + let unique_output = root.join(format!("{mode}-unique.bed")); + let duplicate_output = root.join(format!("{mode}-duplicate.bed")); + let mut unique_args = vec![ + "pileup", + bam_path.to_str().unwrap(), + unique_output.to_str().unwrap(), + "--ref", + fasta_path.to_str().unwrap(), + "--modified-bases", + "C:m", + "--no-filtering", + "--threads", + "1", + "--suppress-progress", + ]; + unique_args.extend_from_slice(mode_args); + run_modkit(&unique_args).unwrap(); + + let mut duplicate_args = vec![ + "pileup", + bam_path.to_str().unwrap(), + duplicate_output.to_str().unwrap(), + "--ref", + fasta_path.to_str().unwrap(), + "--modified-bases", + "C:m", + "C:m", + "--no-filtering", + "--threads", + "1", + "--suppress-progress", + ]; + duplicate_args.extend_from_slice(mode_args); + run_modkit(&duplicate_args).unwrap(); + + let expected = std::fs::read(unique_output).unwrap(); + assert!(!expected.is_empty(), "{mode} control emitted no rows"); + assert_eq!( + std::fs::read(duplicate_output).unwrap(), + expected, + "duplicate selection changed {mode} output" + ); + } +} + +#[test] +fn test_pileup_combined_non_cpg_reverse_statuses_share_anchor() { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + let (bam_path, fasta_path) = write_combined_non_cpg_anchor_fixture(root); + let expected = + b"chr1\t1\t2\tm\t2\t.\t1\t2\t255,0,0\t2\t50.00\t1\t1\t0\t1\t1\t1\t1\n"; + + for threads in [1, 2] { + for high_depth in [false, true] { + let observed = run_combined_non_cpg_anchor_pileup( + root, + &bam_path, + &fasta_path, + high_depth, + threads, + ); + assert_eq!( + observed, expected, + "high_depth {high_depth}, threads {threads}" + ); + } + } +} + #[test] fn test_pileup_cpg_motif_filtering_compressed_ref() { let temp_file = std::env::temp_dir()