diff --git a/book/src/advanced_usage.md b/book/src/advanced_usage.md index 6745d3a..6c380eb 100644 --- a/book/src/advanced_usage.md +++ b/book/src/advanced_usage.md @@ -269,9 +269,12 @@ Modified Base Options: is used, the resulting output BED file will indicate the motif in the "name" field as ,,. For example, given `--motif CGCG 2 --motif CG 0` there will be output lines with name - fields such as "m,CG,0" and "m,CGCG,2". To use this option with - `--combine-strands`, all motifs must be reverse-complement palindromic - or an error will be raised. + fields such as "m,CG,0" and "m,CGCG,2". At most eight motifs are + supported, including motifs added by `--cpg` or `--modified-bases`. + To use `--combine-strands`, exactly one motif must be supplied. It + must be reverse-complement palindromic, and the reverse-strand anchor + must not precede the forward-strand anchor, or an error will be + raised. --cpg Only output counts at CpG motifs diff --git a/book/src/intro_pileup.md b/book/src/intro_pileup.md index f813900..0f488a0 100644 --- a/book/src/intro_pileup.md +++ b/book/src/intro_pileup.md @@ -141,8 +141,11 @@ oligo_741_adapters 39 40 m,CG,0 4 - 39 40 255,0,0 4 100.00 4 0 0 0 0 0 0 oligo_741_adapters 39 40 m,CGCG,2 4 - 39 40 255,0,0 4 100.00 4 0 0 0 0 0 0 ``` -The `--combine-strands` flag can be combined with `--motif` however all motifs must be reverse-complement palindromic (`CG` _is_ a palindrome but `CHH` is not). -Only one motif at a time is supported with `--combine-strands` is used (see [limitations](./limitations.md) for details). +Pileup supports at most eight motifs, including motifs added by `--cpg` or `--modified-bases`. + +The `--combine-strands` flag requires exactly one motif, supplied with `--motif` or `--cpg`, and that motif must be reverse-complement palindromic (`CG` _is_ a palindrome but `CHH` is not). +The reverse-strand anchor must not precede the forward-strand anchor; for example, `--motif CGCG 0 --combine-strands` is supported, while `--motif CGCG 2 --combine-strands` is not currently supported. +See [limitations](./limitations.md) for more details. ## Partitioning reads based on phasing information with `--phased` diff --git a/modkit-core/src/fasta.rs b/modkit-core/src/fasta.rs index 4cac943..58acbdf 100644 --- a/modkit-core/src/fasta.rs +++ b/modkit-core/src/fasta.rs @@ -5,7 +5,6 @@ use bitvec::bitvec; use bitvec::vec::BitVec; use itertools::Itertools; use log::debug; -use rayon::prelude::*; use rust_htslib::faidx; use rustc_hash::FxHashMap; use substring::Substring; @@ -22,6 +21,10 @@ struct HtsFastaHandle { contigs: FxHashMap, preloaded: bool, sequences: FxHashMap, + #[cfg(test)] + sequence_fetches: std::sync::Arc, + #[cfg(test)] + sequence_bases_fetched: std::sync::Arc, } impl HtsFastaHandle { @@ -66,6 +69,10 @@ impl HtsFastaHandle { contigs, sequences, preloaded: preload, + #[cfg(test)] + sequence_fetches: Default::default(), + #[cfg(test)] + sequence_bases_fetched: Default::default(), }) } @@ -76,9 +83,20 @@ impl HtsFastaHandle { end: u64, ) -> MkResult { if let Some(length) = self.contigs.get(contig) { - if end > *length { + if start > end || end > *length { Err(MkError::InvalidReferenceCoordinates) + } else if start == end { + Ok(String::new()) } else { + #[cfg(test)] + { + self.sequence_fetches + .fetch_add(1, std::sync::atomic::Ordering::Relaxed); + self.sequence_bases_fetched.fetch_add( + end - start, + std::sync::atomic::Ordering::Relaxed, + ); + } let seq = if self.preloaded && self.sequences.get(contig).is_some() { @@ -88,7 +106,11 @@ impl HtsFastaHandle { let tmp_reader = faidx::Reader::from_path(&self.fasta_fp) .map_err(|e| MkError::HtsLibError(e))?; tmp_reader - .fetch_seq_string(contig, start as usize, end as usize) + .fetch_seq_string( + contig, + start as usize, + end.saturating_sub(1) as usize, + ) .map_err(|e| MkError::HtsLibError(e))? }; Ok(seq) @@ -101,6 +123,15 @@ impl HtsFastaHandle { fn has_sequence(&self, seq_name: &str) -> bool { self.contigs.contains_key(seq_name) } + + #[cfg(test)] + fn sequence_fetch_stats(&self) -> (usize, u64) { + ( + self.sequence_fetches.load(std::sync::atomic::Ordering::Relaxed), + self.sequence_bases_fetched + .load(std::sync::atomic::Ordering::Relaxed), + ) + } } pub struct HtsLibFastaRecords { @@ -175,111 +206,244 @@ impl MotifLocationsLookup { } #[inline] - fn get_motifs_on_seq( + fn get_motifs_for_owner( &self, - seq: &str, - start: u64, + contig: &str, + owner: std::ops::Range, tid: u32, stranded_position_filter: Option<&StrandedPositionFilter<()>>, - ) -> BitVec { + ) -> MkResult { debug_assert!(!self.motifs.is_empty()); let num_motifs = self.motifs.len(); assert!(num_motifs < u8::MAX as usize); let bits_per_pos = self.motifs.len() * 2; // two strands - let mut mask = bitvec![0; seq.len() * bits_per_pos]; + let ref_end = *self.reader.contigs.get(contig).ok_or_else(|| { + MkError::ContigMissing(format!("{contig} not in FASTA index")) + })?; + if owner.start > owner.end || owner.end > ref_end { + return Err(MkError::InvalidReferenceCoordinates); + } + let owner_len = (owner.end - owner.start) as usize; + let mut mask = bitvec![0; owner_len * bits_per_pos]; + if owner.is_empty() { + return Ok(mask); + } + + // Motif anchors belong to exactly one half-open owner interval, but + // discovering an anchor may require bases on either side of it. Fetch + // enough clipped context for the longest motif, then translate every + // local hit back to its global reference coordinate before filtering + // and writing it into the owner-local mask. + let context = self.longest_motif_length.saturating_sub(1); + let fetch_start = owner.start.saturating_sub(context); + let fetch_end = owner.end.saturating_add(context).min(ref_end); + let seq = self.reader.get_sequence(contig, fetch_start, fetch_end)?; + let seq = if self.mask { seq } else { seq.to_ascii_uppercase() }; for (offset, motif) in self.motifs.iter().enumerate() { - let positions = if let Some(spf) = stranded_position_filter { - find_motif_hits(seq, motif) - .into_par_iter() - .filter(|(pos, strand)| { - spf.contains( - tid as i32, - (*pos as u64).saturating_add(start), - *strand, - ) + for (local_pos, strand) in find_motif_hits(&seq, motif) { + let global_pos = fetch_start.saturating_add(local_pos as u64); + if !owner.contains(&global_pos) + || stranded_position_filter.is_some_and(|spf| { + !spf.contains(tid as i32, global_pos, strand) }) - .collect::>() - } else { - find_motif_hits(seq, motif) - }; - for (pos, strand) in positions { + { + continue; + } + let owner_pos = (global_pos - owner.start) as usize; match strand { Strand::Positive => { - mask.set((pos * bits_per_pos) + offset, true) + mask.set((owner_pos * bits_per_pos) + offset, true) } Strand::Negative => mask.set( - ((pos * bits_per_pos) + num_motifs) + offset, + ((owner_pos * bits_per_pos) + num_motifs) + offset, true, ), } } } - mask + Ok(mask) + } + + /// Find complete palindromic-motif anchor pairs whose positive anchor is + /// in `positive_scan`. A combined-strand output row is written at that + /// positive coordinate, so any stranded position filter must select it; + /// the paired negative anchor is then admitted atomically even when a + /// strand-specific filter does not select it. + fn get_combined_motif_pairs( + &self, + contig: &str, + positive_scan: std::ops::Range, + owner_start: u64, + owner_limit: u64, + tid: u32, + stranded_position_filter: Option<&StrandedPositionFilter<()>>, + ) -> MkResult> { + let ref_end = *self.reader.contigs.get(contig).ok_or_else(|| { + MkError::ContigMissing(format!("{contig} not in FASTA index")) + })?; + if owner_start > positive_scan.start + || positive_scan.start > positive_scan.end + || positive_scan.end > owner_limit + || owner_limit > ref_end + { + return Err(MkError::InvalidReferenceCoordinates); + } + if positive_scan.is_empty() { + return Ok(Vec::new()); + } + + let context = self.longest_motif_length.saturating_sub(1); + let fetch_start = positive_scan.start.saturating_sub(context); + let fetch_end = positive_scan.end.saturating_add(context).min(ref_end); + let seq = self.reader.get_sequence(contig, fetch_start, fetch_end)?; + let seq = if self.mask { seq } else { seq.to_ascii_uppercase() }; + + let mut pairs = Vec::new(); + for (motif_idx, motif) in self.motifs.iter().enumerate() { + debug_assert!(motif.is_palendrome()); + let forward_offset = motif.forward_offset() as u64; + let reverse_offset = motif.reverse_offset() as u64; + for (local_pos, strand) in find_motif_hits(&seq, motif) { + if strand != Strand::Positive { + continue; + } + let Some(positive_position) = + fetch_start.checked_add(local_pos as u64) + else { + continue; + }; + if !positive_scan.contains(&positive_position) + || stranded_position_filter.is_some_and(|spf| { + !spf.contains( + tid as i32, + positive_position, + Strand::Positive, + ) + }) + { + continue; + } + let negative_position = if reverse_offset >= forward_offset { + positive_position + .checked_add(reverse_offset - forward_offset) + } else { + positive_position + .checked_sub(forward_offset - reverse_offset) + }; + let Some(negative_position) = negative_position else { + continue; + }; + + // Both anchors must be in the same selected reference region, + // even when this scan window starts inside that region. + if negative_position < owner_start + || negative_position >= owner_limit + { + continue; + } + pairs.push((positive_position, negative_position, motif_idx)); + } + } + Ok(pairs) } fn get_motif_positions_combine_strands( &mut self, contig: &str, tid: u32, - _ref_end: u64, + ref_end: u64, range: std::ops::Range, stranded_position_filter: Option<&StrandedPositionFilter<()>>, ) -> MkResult<(FocusPositions2, u32)> { - let ref_end = *self.reader.contigs.get(contig).ok_or_else(|| { + let contig_end = *self.reader.contigs.get(contig).ok_or_else(|| { MkError::ContigMissing(format!("{contig} not in FASTA index")) })?; - let buffer_size = self.longest_motif_length * 5; + let owner_limit = ref_end.min(contig_end); + if range.start > range.end || range.end > owner_limit { + return Err(MkError::InvalidReferenceCoordinates); + } let num_motifs = self.motifs.len(); - let bits_per_pos = (num_motifs * 2) as u64; // two strands - let mut end = range.end; - let mut end_w_buffer = std::cmp::min(range.end + buffer_size, ref_end); - let mask = 'fetch_loop: loop { - let seq = - self.reader.get_sequence(contig, range.start, end_w_buffer)?; - let seq = if self.mask { seq } else { seq.to_ascii_uppercase() }; - let mask = self.get_motifs_on_seq( - &seq, + let bits_per_pos = num_motifs * 2; // two strands + if range.is_empty() { + let mask = bitvec![0; 0]; + return Ok(( + FocusPositions2::MotifMask { + mask, + num_motifs: self.num_motifs() as u8, + }, + range.end as u32, + )); + } + + // Discover a little past the requested owner so an ordinary boundary + // pair closes in one scan. If overlapping pairs carry the closure + // farther, scan disjoint lookahead windows that grow geometrically. + // This preserves exact transitive ownership without repeatedly + // fetching and searching the full growing owner prefix. + let mut closure_end = range.end; + let mut scan_start = range.start; + let mut window_size = self.longest_motif_length.max(1); + let mut scan_end = + range.end.saturating_add(window_size).min(owner_limit); + let mut pairs = Vec::new(); + loop { + let mut new_pairs = self.get_combined_motif_pairs( + contig, + scan_start..scan_end, range.start, + owner_limit, tid, stranded_position_filter, - ); - let mut end_idx = ((end - range.start - 1) * bits_per_pos) as usize; - assert!( - (end_idx + num_motifs) < mask.len(), - "off the end {} {}", - end_idx + num_motifs, - mask.len() - ); - - while mask[end_idx..(end_idx + num_motifs)].any() { - debug!("end_idx ({end_idx}) hits a motif.. end={end}",); - end = std::cmp::min( - end.saturating_add(self.longest_motif_length), - ref_end, - ); - end_idx = ((end - range.start - 1) * bits_per_pos) as usize; - debug!( - "moved end_idx to {end_idx}, end={end} ref_end={ref_end}, \ - contig={contig}, range={range:?}" - ); - debug_assert!(end_idx + num_motifs < mask.len()); - if end >= end_w_buffer { - debug!("too close, re-fetching.."); - end_w_buffer = - std::cmp::min(end.saturating_add(buffer_size), ref_end); - continue 'fetch_loop; + )?; + new_pairs.sort_unstable_by_key(|(positive, _, _)| *positive); + for (positive, negative, _) in &new_pairs { + if *positive >= closure_end { + break; + } + if let Some(pair_end) = positive.max(negative).checked_add(1) { + closure_end = closure_end.max(pair_end); } - debug!("re-check.."); } - break mask; - }; + pairs.extend(new_pairs); + + if closure_end <= scan_end { + break; + } + + let previous_scan_end = scan_end; + scan_start = scan_end; + window_size = window_size.saturating_mul(2); + scan_end = scan_end + .saturating_add(window_size) + .max(closure_end) + .min(owner_limit); + debug!( + "growing combined motif lookahead from {previous_scan_end} to \ + {scan_end} for owner closure {closure_end}, contig={contig}, \ + range={range:?}" + ); + debug_assert!(scan_end > previous_scan_end); + } + + pairs.retain(|(positive, _, _)| *positive < closure_end); + let mut mask = + bitvec![0; (closure_end - range.start) as usize * bits_per_pos]; + for (positive, negative, motif_idx) in pairs { + let positive_local = (positive - range.start) as usize; + let negative_local = (negative - range.start) as usize; + mask.set((positive_local * bits_per_pos) + motif_idx, true); + mask.set( + (negative_local * bits_per_pos) + num_motifs + motif_idx, + true, + ); + } Ok(( FocusPositions2::MotifMask { mask, num_motifs: self.num_motifs() as u8, }, - end as u32, + closure_end as u32, )) } @@ -301,15 +465,12 @@ impl MotifLocationsLookup { stranded_position_filter, ) } else { - let seq = - self.reader.get_sequence(contig, range.start, range.end)?; - let seq = if self.mask { seq } else { seq.to_ascii_uppercase() }; - let mask = self.get_motifs_on_seq( - &seq, - range.start, + let mask = self.get_motifs_for_owner( + contig, + range.clone(), tid, stranded_position_filter, - ); + )?; let focus_positions = FocusPositions2::MotifMask { mask, num_motifs: self.num_motifs() as u8, @@ -329,10 +490,477 @@ impl MotifLocationsLookup { #[cfg(test)] mod fasta_mod_tests { - use crate::fasta::HtsFastaHandle; + use std::{ + fs::File, + io::Write, + ops::Range, + path::{Path, PathBuf}, + }; + use rand::prelude::{SeedableRng, StdRng}; + use rust_lapper::Lapper; + use rustc_hash::FxHashMap; use rv::prelude::Rv; + use crate::{ + fasta::{HtsFastaHandle, MotifLocationsLookup}, + interval_chunks::FocusPositions2, + motifs::motif_bed::RegexMotif, + position_filter::{Iv, StrandedPositionFilter}, + util::Strand, + }; + + fn write_fasta(root: &Path, sequence: &str) -> PathBuf { + let fasta_path = root.join("reference.fa"); + File::create(&fasta_path) + .unwrap() + .write_all(format!(">chr1\n{sequence}\n").as_bytes()) + .unwrap(); + File::create(root.join("reference.fa.fai")) + .unwrap() + .write_all( + format!( + "chr1\t{}\t6\t{}\t{}\n", + sequence.len(), + sequence.len(), + sequence.len() + 1 + ) + .as_bytes(), + ) + .unwrap(); + fasta_path + } + + fn decode_motif_mask( + focus_positions: FocusPositions2, + range: Range, + num_motifs: usize, + ) -> Vec<(u64, Strand, usize)> { + let FocusPositions2::MotifMask { mask, num_motifs: observed } = + focus_positions + else { + panic!("expected motif mask") + }; + assert_eq!(observed as usize, num_motifs); + let bits_per_pos = num_motifs * 2; + assert_eq!( + mask.len(), + (range.end - range.start) as usize * bits_per_pos, + "motif mask must cover exactly its owner interval" + ); + let mut hits = Vec::new(); + for local_pos in 0..(range.end - range.start) as usize { + for motif_idx in 0..num_motifs { + if mask[(local_pos * bits_per_pos) + motif_idx] { + hits.push(( + range.start + local_pos as u64, + Strand::Positive, + motif_idx, + )); + } + if mask[(local_pos * bits_per_pos) + num_motifs + motif_idx] { + hits.push(( + range.start + local_pos as u64, + Strand::Negative, + motif_idx, + )); + } + } + } + hits + } + + fn collect_motif_hits( + fasta_path: &PathBuf, + motifs: Vec, + preload: bool, + interval_size: u64, + combine_strands: bool, + position_filter: Option<&StrandedPositionFilter<()>>, + ) -> (Vec<(u64, Strand, usize)>, Vec>) { + let mut lookup = MotifLocationsLookup::from_paths( + fasta_path, false, None, motifs, preload, + ) + .unwrap(); + let contig_end = lookup.reader.contigs["chr1"]; + let num_motifs = lookup.num_motifs(); + let mut start = 0; + let mut hits = Vec::new(); + let mut owners = Vec::new(); + while start < contig_end { + let requested_end = + start.saturating_add(interval_size).min(contig_end); + let (focus_positions, actual_end) = lookup + .get_motif_positions( + "chr1", + 0, + contig_end as u32, + start..requested_end, + position_filter, + combine_strands, + ) + .unwrap(); + let owner = start..actual_end as u64; + assert!(owner.end > owner.start); + hits.extend(decode_motif_mask( + focus_positions, + owner.clone(), + num_motifs, + )); + owners.push(owner.clone()); + start = owner.end; + } + hits.sort_by_key(|(position, strand, motif_idx)| { + (*position, *strand, *motif_idx) + }); + (hits, owners) + } + + #[test] + fn fasta_subsequences_are_half_open_with_or_without_preload() { + let temp_dir = tempfile::tempdir().unwrap(); + let fasta_path = write_fasta(temp_dir.path(), "ACGT"); + + for preload in [false, true] { + let reader = + HtsFastaHandle::from_file(&fasta_path, preload).unwrap(); + for (start, end, expected) in [ + (0, 0, ""), + (0, 1, "A"), + (1, 3, "CG"), + (3, 4, "T"), + (4, 4, ""), + (0, 4, "ACGT"), + ] { + assert_eq!( + reader.get_sequence("chr1", start, end).unwrap(), + expected, + "preload={preload}, range={start}..{end}" + ); + } + assert!(reader.get_sequence("chr1", 3, 2).is_err()); + assert!(reader.get_sequence("chr1", 0, 5).is_err()); + } + } + + #[test] + fn motif_hits_are_interval_and_preload_invariant() { + struct Case<'a> { + sequence: &'a str, + motif: &'a str, + offset: usize, + expected: Vec<(u64, Strand, usize)>, + } + + let cases = [ + Case { + sequence: "ACGTCG", + motif: "CG", + offset: 0, + expected: vec![ + (1, Strand::Positive, 0), + (2, Strand::Negative, 0), + (4, Strand::Positive, 0), + (5, Strand::Negative, 0), + ], + }, + Case { + sequence: "GATCNNGATC", + motif: "GATC", + offset: 1, + expected: vec![ + (1, Strand::Positive, 0), + (2, Strand::Negative, 0), + (7, Strand::Positive, 0), + (8, Strand::Negative, 0), + ], + }, + Case { + sequence: "CGTACG", + motif: "CGT", + offset: 0, + expected: vec![ + (0, Strand::Positive, 0), + (5, Strand::Negative, 0), + ], + }, + ]; + + for case in cases { + let temp_dir = tempfile::tempdir().unwrap(); + let fasta_path = write_fasta(temp_dir.path(), case.sequence); + for preload in [false, true] { + for interval_size in [1, 2, 3, case.sequence.len() as u64, 100] + { + let motif = + RegexMotif::parse_string(case.motif, case.offset) + .unwrap(); + let (observed, _) = collect_motif_hits( + &fasta_path, + vec![motif], + preload, + interval_size, + false, + None, + ); + assert_eq!( + observed, case.expected, + "sequence={}, motif={} {}, preload={preload}, \ + interval_size={interval_size}", + case.sequence, case.motif, case.offset + ); + } + } + } + } + + #[test] + fn motif_position_filter_uses_global_stranded_anchors() { + let temp_dir = tempfile::tempdir().unwrap(); + let fasta_path = write_fasta(temp_dir.path(), "NNGATCNN"); + let mut pos_positions = FxHashMap::default(); + pos_positions + .insert(0, Lapper::new(vec![Iv { start: 3, stop: 4, val: () }])); + let mut neg_positions = FxHashMap::default(); + neg_positions + .insert(0, Lapper::new(vec![Iv { start: 4, stop: 5, val: () }])); + let position_filter = + StrandedPositionFilter { pos_positions, neg_positions }; + let motif = RegexMotif::parse_string("GATC", 1).unwrap(); + + for preload in [false, true] { + let (observed, _) = collect_motif_hits( + &fasta_path, + vec![motif.clone()], + preload, + 1, + false, + Some(&position_filter), + ); + assert_eq!( + observed, + vec![(3, Strand::Positive, 0), (4, Strand::Negative, 0),] + ); + } + } + + #[test] + fn combined_strand_lookup_preserves_boundary_extension() { + let temp_dir = tempfile::tempdir().unwrap(); + let fasta_path = write_fasta(temp_dir.path(), "GATCNNNN"); + let motif = RegexMotif::parse_string("GATC", 1).unwrap(); + + for preload in [false, true] { + let (observed, owners) = collect_motif_hits( + &fasta_path, + vec![motif.clone()], + preload, + 2, + true, + None, + ); + assert_eq!(owners, vec![0..3, 3..5, 5..7, 7..8]); + assert_eq!( + observed, + vec![(1, Strand::Positive, 0), (2, Strand::Negative, 0),] + ); + } + } + + #[test] + fn combined_strand_lookup_extends_to_the_actual_cgcg_pair() { + let temp_dir = tempfile::tempdir().unwrap(); + let fasta_path = write_fasta(temp_dir.path(), "CGCGNN"); + let motif = RegexMotif::parse_string("CGCG", 0).unwrap(); + + for preload in [false, true] { + for (interval_size, expected_owners) in [ + (1, vec![0..4, 4..5, 5..6]), + (2, vec![0..4, 4..6]), + (3, vec![0..4, 4..6]), + (6, vec![0..6]), + ] { + let (observed, owners) = collect_motif_hits( + &fasta_path, + vec![motif.clone()], + preload, + interval_size, + true, + None, + ); + assert_eq!(owners, expected_owners); + assert_eq!( + observed, + vec![(0, Strand::Positive, 0), (3, Strand::Negative, 0),], + "preload={preload}, interval_size={interval_size}" + ); + } + } + } + + #[test] + fn combined_strand_lookup_closes_over_every_new_boundary_pair() { + let temp_dir = tempfile::tempdir().unwrap(); + let fasta_path = write_fasta(temp_dir.path(), "CGCGCGNN"); + let motif = RegexMotif::parse_string("CGCG", 0).unwrap(); + + for preload in [false, true] { + let (observed, owners) = collect_motif_hits( + &fasta_path, + vec![motif.clone()], + preload, + 1, + true, + None, + ); + assert_eq!(owners, vec![0..6, 6..7, 7..8]); + assert_eq!( + observed, + vec![ + (0, Strand::Positive, 0), + (2, Strand::Positive, 0), + (3, Strand::Negative, 0), + (5, Strand::Negative, 0), + ] + ); + } + } + + #[test] + fn combined_strand_dense_overlap_scan_work_is_bounded() { + let temp_dir = tempfile::tempdir().unwrap(); + let repeat_count = 256usize; + let sequence = format!("{}NN", "CG".repeat(repeat_count)); + let fasta_path = write_fasta(temp_dir.path(), &sequence); + let motif = RegexMotif::parse_string("CGCG", 0).unwrap(); + let mut expected = (0..repeat_count - 1) + .flat_map(|i| { + [ + (2 * i as u64, Strand::Positive, 0), + (2 * i as u64 + 3, Strand::Negative, 0), + ] + }) + .collect::>(); + expected.sort_by_key(|(position, strand, motif_idx)| { + (*position, *strand, *motif_idx) + }); + + for preload in [false, true] { + let mut lookup = MotifLocationsLookup::from_paths( + &fasta_path, + false, + None, + vec![motif.clone()], + preload, + ) + .unwrap(); + let (focus_positions, actual_end) = lookup + .get_motif_positions( + "chr1", + 0, + sequence.len() as u32, + 0..1, + None, + true, + ) + .unwrap(); + + assert_eq!(actual_end, (repeat_count * 2) as u32); + assert_eq!( + decode_motif_mask(focus_positions, 0..actual_end as u64, 1,), + expected, + "preload={preload}" + ); + + let (fetches, bases_fetched) = lookup.reader.sequence_fetch_stats(); + assert!( + fetches <= 12, + "dense overlap required {fetches} sequence fetches with \ + preload={preload}" + ); + assert!( + bases_fetched <= 3 * sequence.len() as u64, + "dense overlap fetched {bases_fetched} bases for a {}-base \ + reference with preload={preload}", + sequence.len() + ); + } + } + + #[test] + fn combined_strand_lookup_owns_and_filters_pairs_by_positive_anchor() { + let temp_dir = tempfile::tempdir().unwrap(); + let sequence = format!("{}CGCG{}", "N".repeat(38), "N".repeat(28)); + let fasta_path = write_fasta(temp_dir.path(), &sequence); + let motif = RegexMotif::parse_string("CGCG", 0).unwrap(); + + for preload in [false, true] { + let mut lookup = MotifLocationsLookup::from_paths( + &fasta_path, + false, + None, + vec![motif.clone()], + preload, + ) + .unwrap(); + + // The reverse anchor at 41 is not an independently owned hit when + // its positive pair anchor at 38 is before the selected region. + let (focus_positions, actual_end) = lookup + .get_motif_positions("chr1", 0, 62, 40..62, None, true) + .unwrap(); + assert_eq!(actual_end, 62); + assert!(decode_motif_mask(focus_positions, 40..62, 1).is_empty()); + + let mut pos_positions = FxHashMap::default(); + pos_positions.insert( + 0, + Lapper::new(vec![Iv { start: 38, stop: 39, val: () }]), + ); + let positive_only = StrandedPositionFilter { + pos_positions, + neg_positions: FxHashMap::default(), + }; + let (focus_positions, actual_end) = lookup + .get_motif_positions( + "chr1", + 0, + 62, + 38..62, + Some(&positive_only), + true, + ) + .unwrap(); + assert_eq!(actual_end, 62); + assert_eq!( + decode_motif_mask(focus_positions, 38..62, 1), + vec![(38, Strand::Positive, 0), (41, Strand::Negative, 0),] + ); + + let mut neg_positions = FxHashMap::default(); + neg_positions.insert( + 0, + Lapper::new(vec![Iv { start: 41, stop: 42, val: () }]), + ); + let negative_only = StrandedPositionFilter { + pos_positions: FxHashMap::default(), + neg_positions, + }; + let (focus_positions, actual_end) = lookup + .get_motif_positions( + "chr1", + 0, + 62, + 38..62, + Some(&negative_only), + true, + ) + .unwrap(); + assert_eq!(actual_end, 62); + assert!(decode_motif_mask(focus_positions, 38..62, 1).is_empty()); + } + } + #[test] fn test_hts_fasta_reader() { let compressed_fp = std::path::Path::new( diff --git a/modkit-core/src/pileup/subcommand.rs b/modkit-core/src/pileup/subcommand.rs index 4503f99..c566990 100644 --- a/modkit-core/src/pileup/subcommand.rs +++ b/modkit-core/src/pileup/subcommand.rs @@ -57,6 +57,18 @@ use crate::writers::{ PhasedBedMethylWriter, PileupWriter, }; +const MAX_PILEUP_MOTIFS: usize = u8::BITS as usize; + +fn validate_motif_capacity(motif_count: usize) -> anyhow::Result<()> { + if motif_count > MAX_PILEUP_MOTIFS { + bail!( + "pileup supports at most {MAX_PILEUP_MOTIFS} motifs because motif \ + membership is stored in an 8-bit mask; received {motif_count}" + ) + } + Ok(()) +} + #[derive(Args)] #[command(arg_required_else_help = true)] pub struct ModBamPileup { @@ -304,9 +316,11 @@ pub struct ModBamPileup { /// used, the resulting output BED file will indicate the motif in the /// "name" field as ,,. For example, given /// `--motif CGCG 2 --motif CG 0` there will be output lines with name - /// fields such as "m,CG,0" and "m,CGCG,2". To use this option with - /// `--combine-strands`, all motifs must be reverse-complement - /// palindromic or an error will be raised. + /// fields such as "m,CG,0" and "m,CGCG,2". At most eight motifs are + /// supported, including motifs added by `--cpg` or `--modified-bases`. + /// To use `--combine-strands`, exactly one motif must be supplied. It must + /// be reverse-complement palindromic, and the reverse-strand anchor must + /// not precede the forward-strand anchor, or an error will be raised. #[clap(help_heading = "Modified Base Options")] #[arg(long, action = clap::ArgAction::Append, num_args = 2, requires = "reference_fasta")] motif: Option>, @@ -438,6 +452,48 @@ impl ModBamPileup { motifs == &[RegexMotif::parse_string("CG", 0).unwrap()] } + fn validate_combine_strands_motifs( + combine_strands: bool, + motifs: &[RegexMotif], + ) -> anyhow::Result<()> { + if !combine_strands { + return Ok(()); + } + + let motif = match motifs { + [motif] => motif, + [] => bail!( + "need to provide one reverse-complement palindromic motif to \ + combine strands" + ), + _ => bail!( + "multiple motifs and combine-strands not currently supported" + ), + }; + if !motif.is_palendrome() { + bail!( + "cannot combine strands for motif '{} {}': motif must be \ + reverse-complement palindromic", + motif.raw_motif, + motif.motif_info.forward_offset, + ); + } + + let motif_info = motif.motif_info; + if motif_info.reverse_offset < motif_info.forward_offset { + bail!( + "cannot combine strands for motif '{} {}': reverse-strand \ + anchor offset {} precedes forward-strand anchor offset {}; \ + the reverse anchor must not precede the forward anchor", + motif.raw_motif, + motif_info.forward_offset, + motif_info.reverse_offset, + motif_info.forward_offset, + ); + } + Ok(()) + } + fn is_5mc_5hmc_cpg( modified_bases: &[(DnaBase, ModCodeRepr)], ) -> anyhow::Result { @@ -543,6 +599,11 @@ impl ModBamPileup { ) -> anyhow::Result<(Option, Option>)> { let mut regex_motifs = self.parse_user_motifs().transpose()?.unwrap_or_else(Vec::new); + validate_motif_capacity(regex_motifs.len())?; + Self::validate_combine_strands_motifs( + self.combine_strands, + ®ex_motifs, + )?; if regex_motifs.len() > 1 { if self.combine_strands { bail!( @@ -553,7 +614,7 @@ impl ModBamPileup { info!( "more than one motif requires use of general pileup processor" ); - let motif_primary_bases = regex_motifs + let mut motif_primary_bases = regex_motifs .iter() .map(|mot| mot.motif_info.primary_base) .collect::>(); @@ -561,7 +622,7 @@ impl ModBamPileup { for primary_base in modified_bases.iter().map(|x| x.primary_base) { - if !motif_primary_bases.contains(&primary_base) { + if motif_primary_bases.insert(primary_base) { info!( "adding single-base motif: '{} 0'", primary_base.char() @@ -576,6 +637,7 @@ impl ModBamPileup { } } } + validate_motif_capacity(regex_motifs.len())?; return Ok((None, Some(regex_motifs))); } if let Some(modified_bases) = self.modified_bases.as_ref() { @@ -640,11 +702,7 @@ impl ModBamPileup { motif to combine strands" ) } - let motif_offset = regex_motifs[0].motif_info.offset(); - if motif_offset < 0 { - bail!("invalid palindromic motif"); - } - let motif_offset = motif_offset as u32; + let motif_offset = regex_motifs[0].motif_info.offset() as u32; let motif_bases = [ regex_motifs[0].motif_info.primary_base, DnaBase::C, @@ -1029,6 +1087,9 @@ impl ModBamPileup { Presets::DnaCpGCombineStrands { .. } => { (PileupNumericOptions::Passthrough, true) } + Presets::DynamicAllContext { + motif_offset: Some(_), .. + } => (PileupNumericOptions::Passthrough, true), _ => (PileupNumericOptions::Passthrough, false), }, None => { @@ -1664,6 +1725,20 @@ impl ModBamPileup { } } +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn motif_capacity_accepts_eight_and_rejects_nine() { + assert!(validate_motif_capacity(8).is_ok()); + + let error = validate_motif_capacity(9).unwrap_err().to_string(); + assert!(error.contains("at most 8 motifs")); + assert!(error.contains("received 9")); + } +} + #[derive(Clone, Debug)] enum Presets { /// CpG-special, combine strands, maybe combine mods diff --git a/modkit/tests/test_pileup.rs b/modkit/tests/test_pileup.rs index 68f2cac..de3a139 100644 --- a/modkit/tests/test_pileup.rs +++ b/modkit/tests/test_pileup.rs @@ -1,11 +1,15 @@ 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::{BTreeSet, 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 std::process::{Command, Output}; use common::{check_against_expected_text_file, run_modkit}; use mod_kit::dmr::bedmethyl::BedMethylLine; @@ -13,6 +17,359 @@ use mod_kit::mod_base_code::{ModCodeRepr, METHYL_CYTOSINE}; mod common; +fn write_motif_boundary_fixture(root: &Path) -> (PathBuf, PathBuf) { + let bam_path = root.join("motif-boundaries.bam"); + let fasta_path = root.join("motif-boundaries.fa"); + + let mut header = Header::new(); + let mut sq = HeaderRecord::new(b"SQ"); + sq.push_tag(b"SN", "chr1"); + sq.push_tag(b"LN", 6); + header.push_record(&sq); + + let mut writer = + BamWriter::from_path(&bam_path, &header, Format::Bam).unwrap(); + for (name, reverse) in [("forward", false), ("reverse", true)] { + let cigar = CigarString(vec![Cigar::Match(6)]); + let mut record = Record::new(); + record.set(name.as_bytes(), Some(&cigar), b"ACGTCG", &[30; 6]); + record.set_tid(0); + record.set_pos(0); + record.set_mapq(60); + if reverse { + record.set_flags(16); + } + record.push_aux(b"MM", Aux::String("C+m?,0,0;")).unwrap(); + record.push_aux(b"ML", Aux::ArrayU8((&[255, 255][..]).into())).unwrap(); + record.push_aux(b"MN", Aux::U32(6)).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\nACGTCG\n").unwrap(); + File::create(root.join("motif-boundaries.fa.fai")) + .unwrap() + .write_all(b"chr1\t6\t6\t6\t7\n") + .unwrap(); + + (bam_path, fasta_path) +} + +fn run_motif_boundary_pileup( + root: &Path, + bam_path: &Path, + fasta_path: &Path, + preload: bool, + interval_size: usize, +) -> Vec { + let preload_name = if preload { "preload" } else { "faidx" }; + let output_path = + root.join(format!("motif-{preload_name}-{interval_size}.bed")); + let interval_size = interval_size.to_string(); + let mut args = vec![ + "pileup", + bam_path.to_str().unwrap(), + output_path.to_str().unwrap(), + "--ref", + fasta_path.to_str().unwrap(), + "--motif", + "CG", + "0", + "--modified-bases", + "C:m", + "--no-filtering", + "--interval-size", + &interval_size, + "--threads", + "1", + "--suppress-progress", + ]; + if preload { + args.push("--preload-references"); + } + run_modkit(&args).unwrap(); + std::fs::read(output_path).unwrap() +} + +fn run_cgcg0_combined_pileup( + root: &Path, + optimized: bool, + region: &str, + interval_size: usize, +) -> Vec { + let processor = if optimized { "optimized" } else { "generic" }; + let output_path = root.join(format!( + "cgcg0-{processor}-{}-{interval_size}.bed", + region.replace(':', "-") + )); + let interval_size = interval_size.to_string(); + let mut args = vec![ + "pileup", + "../tests/resources/CG_5mC_20230207_1700_6A_PAG66026_3c0abf27_oligo_741_adapters_modcalls_0th_sort_10_reads.bam", + output_path.to_str().unwrap(), + "--motif", + "CGCG", + "0", + "--combine-strands", + "--no-filtering", + "--ref", + "../tests/resources/CGI_ladder_3.6kb_ref.fa", + "--region", + region, + "--interval-size", + &interval_size, + "--threads", + "1", + "--suppress-progress", + ]; + if optimized { + args.extend(["--modified-bases", "C:m"]); + } + run_modkit(&args).unwrap(); + std::fs::read(output_path).unwrap() +} + +fn run_combine_motif_validation_pileup( + output_path: &Path, + motifs: &[(&str, &str)], + optimized: bool, +) -> std::process::Output { + let mut command = std::process::Command::new(env!("CARGO_BIN_EXE_modkit")); + command.args([ + "pileup", + "../tests/resources/CG_5mC_20230207_1700_6A_PAG66026_3c0abf27_oligo_741_adapters_modcalls_0th_sort_10_reads.bam", + output_path.to_str().unwrap(), + "--combine-strands", + "--no-filtering", + "--ref", + "../tests/resources/CGI_ladder_3.6kb_ref.fa", + "--region", + "oligo_741_adapters:22-62", + "--interval-size", + "40", + "--threads", + "1", + "--suppress-progress", + ]); + for (motif, offset) in motifs { + command.args(["--motif", motif, offset]); + } + if optimized { + command.args(["--modified-bases", "C:m"]); + } + command.output().unwrap() +} + +#[test] +fn test_pileup_motif_boundaries_are_preload_and_interval_invariant() { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + let (bam_path, fasta_path) = write_motif_boundary_fixture(root); + let expected = concat!( + "chr1\t1\t2\tm\t1\t+\t1\t2\t255,0,0\t1\t100.00\t1\t0\t0\t0\t0\t0\t0\n", + "chr1\t2\t3\tm\t1\t-\t2\t3\t255,0,0\t1\t100.00\t1\t0\t0\t0\t0\t0\t0\n", + "chr1\t4\t5\tm\t1\t+\t4\t5\t255,0,0\t1\t100.00\t1\t0\t0\t0\t0\t0\t0\n", + "chr1\t5\t6\tm\t1\t-\t5\t6\t255,0,0\t1\t100.00\t1\t0\t0\t0\t0\t0\t0\n", + ); + + for preload in [false, true] { + for interval_size in [1, 2, 3, 6, 100] { + let observed = run_motif_boundary_pileup( + root, + &bam_path, + &fasta_path, + preload, + interval_size, + ); + assert_eq!( + observed, + expected.as_bytes(), + "preload={preload}, interval_size={interval_size}" + ); + } + } +} + +#[test] +fn test_cgcg0_combined_optimized_and_generic_are_interval_invariant() { + let temp_dir = tempfile::tempdir().unwrap(); + let expected = b"oligo_741_adapters\t38\t39\tm\t11\t.\t38\t39\t255,0,0\t11\t100.00\t11\t0\t0\t1\t0\t0\t0\n"; + + for interval_size in [1, 2, 3, 40] { + let optimized = run_cgcg0_combined_pileup( + temp_dir.path(), + true, + "oligo_741_adapters:22-62", + interval_size, + ); + let generic = run_cgcg0_combined_pileup( + temp_dir.path(), + false, + "oligo_741_adapters:22-62", + interval_size, + ); + assert_eq!(optimized, expected, "optimized interval={interval_size}"); + assert_eq!(generic, expected, "generic interval={interval_size}"); + } + + // The only CGCG0 pair in this span is +38/-41. Its positive owner is + // before the region, so neither processor may fabricate a row at 40. + for optimized in [false, true] { + let observed = run_cgcg0_combined_pileup( + temp_dir.path(), + optimized, + "oligo_741_adapters:40-62", + 2, + ); + assert!(observed.is_empty(), "optimized={optimized}"); + } +} + +#[test] +fn test_negative_combine_anchor_is_rejected_before_opening_output() { + let temp_dir = tempfile::tempdir().unwrap(); + let sentinel = b"existing output must remain unchanged\n"; + + for motif_offset in ["2", "3"] { + for optimized in [false, true] { + let processor = if optimized { "optimized" } else { "generic" }; + let output_path = temp_dir + .path() + .join(format!("cgcg{motif_offset}-{processor}.bed")); + + let output = run_combine_motif_validation_pileup( + &output_path, + &[("CGCG", motif_offset)], + optimized, + ); + let diagnostics = format!( + "{}{}", + String::from_utf8_lossy(&output.stdout), + String::from_utf8_lossy(&output.stderr) + ); + assert!( + !output.status.success(), + "CGCG {motif_offset} unexpectedly succeeded with {processor}" + ); + assert!( + diagnostics.contains(&format!( + "cannot combine strands for motif 'CGCG {motif_offset}'" + )), + "unexpected {processor} diagnostics: {diagnostics}" + ); + assert!( + diagnostics.contains( + "the reverse anchor must not precede the forward anchor" + ), + "unexpected {processor} diagnostics: {diagnostics}" + ); + assert!( + !output_path.exists(), + "invalid {processor} request created its output" + ); + + std::fs::write(&output_path, sentinel).unwrap(); + let output = run_combine_motif_validation_pileup( + &output_path, + &[("CGCG", motif_offset)], + optimized, + ); + assert!( + !output.status.success(), + "CGCG {motif_offset} unexpectedly succeeded with {processor}" + ); + assert_eq!( + std::fs::read(&output_path).unwrap(), + sentinel, + "invalid {processor} request truncated its output" + ); + } + } +} + +fn assert_invalid_combine_motif_configuration_preserves_output( + label: &str, + motifs: &[(&str, &str)], + expected_diagnostic: &str, +) { + let temp_dir = tempfile::tempdir().unwrap(); + let sentinel = b"existing output must remain unchanged\n"; + for optimized in [false, true] { + let processor = if optimized { "optimized" } else { "generic" }; + let output_path = + temp_dir.path().join(format!("{label}-{processor}.bed")); + + let output = run_combine_motif_validation_pileup( + &output_path, + motifs, + optimized, + ); + let diagnostics = format!( + "{}{}", + String::from_utf8_lossy(&output.stdout), + String::from_utf8_lossy(&output.stderr) + ); + assert!( + !output.status.success(), + "{label} unexpectedly succeeded with {processor}" + ); + assert!( + diagnostics.contains(expected_diagnostic), + "unexpected {processor} diagnostics for {label}: {diagnostics}" + ); + assert!( + !output_path.exists(), + "invalid {label} {processor} request created its output" + ); + + std::fs::write(&output_path, sentinel).unwrap(); + let output = run_combine_motif_validation_pileup( + &output_path, + motifs, + optimized, + ); + assert!( + !output.status.success(), + "{label} unexpectedly succeeded with {processor}" + ); + assert_eq!( + std::fs::read(&output_path).unwrap(), + sentinel, + "invalid {label} {processor} request truncated its output" + ); + } +} + +#[test] +fn test_missing_combine_motif_is_rejected_before_output() { + assert_invalid_combine_motif_configuration_preserves_output( + "no-motif", + &[], + "combine strands", + ); +} + +#[test] +fn test_non_palindromic_combine_motif_is_rejected_before_output() { + assert_invalid_combine_motif_configuration_preserves_output( + "non-palindrome", + &[("CAT", "0")], + "combine strands", + ); +} + +#[test] +fn test_multiple_combine_motifs_are_rejected_before_output() { + assert_invalid_combine_motif_configuration_preserves_output( + "multiple-motifs", + &[("CG", "0"), ("GATC", "1")], + "multiple motifs and combine-strands not currently supported", + ); +} + #[test] fn test_pileup_help() { let pileup_help_args = ["pileup", "--help"]; @@ -839,6 +1196,168 @@ fn test_pileup_motifs_cg0_cgcg2() { ); } +fn run_pileup_with_overlapping_motifs( + output_path: &Path, + motif_count: usize, + modified_base: &str, + extra_modified_base: Option<&str>, + add_cpg: bool, + extra_motif: Option<(&str, &str)>, +) -> Output { + let mut command = Command::new(env!("CARGO_BIN_EXE_modkit")); + command.args([ + "pileup", + "../tests/resources/CG_5mC_20230207_1700_6A_PAG66026_3c0abf27_oligo_741_adapters_modcalls_0th_sort_10_reads.bam", + output_path.to_str().unwrap(), + "--no-filtering", + "--ref", + "../tests/resources/CGI_ladder_3.6kb_ref.fa", + "--region", + "oligo_741_adapters:22-62", + "--modified-bases", + modified_base, + ]); + if let Some(extra_modified_base) = extra_modified_base { + command.arg(extra_modified_base); + } + command.args(["--threads", "1", "--io-threads", "1"]); + + for (motif, offset) in [ + ("CGG", "0"), + ("CGGG", "0"), + ("TCG", "1"), + ("CTCG", "2"), + ("GCTCG", "3"), + ("TGCTCG", "4"), + ("TTGCTCG", "5"), + ("ATTGCTCG", "6"), + ] + .into_iter() + .take(motif_count) + { + command.args(["--motif", motif, offset]); + } + if let Some((motif, offset)) = extra_motif { + command.args(["--motif", motif, offset]); + } + if add_cpg { + command.arg("--cpg"); + } + + command.output().unwrap() +} + +#[test] +fn test_pileup_rejects_more_than_eight_motifs() { + let temp_dir = tempfile::tempdir().unwrap(); + + let eight_motifs_path = temp_dir.path().join("eight.bed"); + let output = run_pileup_with_overlapping_motifs( + &eight_motifs_path, + 8, + "5mC", + None, + false, + None, + ); + assert!( + output.status.success(), + "eight motifs should be accepted: {}", + String::from_utf8_lossy(&output.stderr) + ); + let observed_labels = + BufReader::new(File::open(&eight_motifs_path).unwrap()) + .lines() + .map(|line| line.unwrap().split('\t').nth(3).unwrap().to_string()) + .collect::>(); + let expected_labels = [ + "m,CGG,0", + "m,CGGG,0", + "m,TCG,1", + "m,CTCG,2", + "m,GCTCG,3", + "m,TGCTCG,4", + "m,TTGCTCG,5", + "m,ATTGCTCG,6", + ] + .into_iter() + .map(str::to_string) + .collect::>(); + assert_eq!(observed_labels, expected_labels); + + let deduplicated_auto_motif_path = + temp_dir.path().join("deduplicated-auto-motif.bed"); + let output = run_pileup_with_overlapping_motifs( + &deduplicated_auto_motif_path, + 7, + "m6A", + Some("2OmeA"), + false, + None, + ); + assert!( + output.status.success(), + "two codes for one missing primary base should add one motif: {}", + String::from_utf8_lossy(&output.stderr) + ); + let automatically_added_labels = + BufReader::new(File::open(&deduplicated_auto_motif_path).unwrap()) + .lines() + .map(|line| line.unwrap().split('\t').nth(3).unwrap().to_string()) + .collect::>(); + assert!(automatically_added_labels.contains("a,A,0")); + + for (filename, modified_base, add_cpg, extra_motif) in [ + ("ninth-explicit.bed", "5mC", false, Some(("CGCG", "0"))), + ("ninth-cpg.bed", "5mC", true, None), + ("ninth-added-base.bed", "m6A", false, None), + ] { + let output_path = temp_dir.path().join(filename); + let output = run_pileup_with_overlapping_motifs( + &output_path, + 8, + modified_base, + None, + add_cpg, + extra_motif, + ); + let stderr = String::from_utf8_lossy(&output.stderr); + assert!(!output.status.success(), "nine motifs should fail"); + assert!( + stderr.contains("pileup supports at most 8 motifs") + && stderr.contains("received 9"), + "expected a clear motif-capacity diagnostic, got: {stderr}" + ); + assert!( + !output_path.exists(), + "motif validation should happen before creating output" + ); + + let sentinel = b"existing output\n"; + std::fs::write(&output_path, sentinel).unwrap(); + let output = run_pileup_with_overlapping_motifs( + &output_path, + 8, + modified_base, + None, + add_cpg, + extra_motif, + ); + let stderr = String::from_utf8_lossy(&output.stderr); + assert!(!output.status.success(), "nine motifs should fail"); + assert!( + stderr.contains("pileup supports at most 8 motifs") + && stderr.contains("received 9"), + "expected a clear motif-capacity diagnostic, got: {stderr}" + ); + assert_eq!( + std::fs::read(&output_path).unwrap(), + sentinel, + "motif validation should happen before changing output" + ); + } +} + #[test] #[ignore = "multiple motifs and combine strands not supported in v0.6.0"] fn test_pileup_motifs_cg0_cgcg2_combined() {