From b6eb50751e1786627dd9c6cfd8bf852fea253a37 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 03:38:26 -0700 Subject: [PATCH 1/9] Validate localize region files --- modkit-core/src/localise/subcommand.rs | 705 ++++++++++++++++++++++--- modkit/tests/test_localize.rs | 106 ++++ 2 files changed, 740 insertions(+), 71 deletions(-) diff --git a/modkit-core/src/localise/subcommand.rs b/modkit-core/src/localise/subcommand.rs index 0a36c143..6f2b55ef 100644 --- a/modkit-core/src/localise/subcommand.rs +++ b/modkit-core/src/localise/subcommand.rs @@ -1,4 +1,3 @@ -use std::collections::HashMap; use std::fs::File; use std::io::{stdout, BufRead, BufReader, BufWriter, Write}; use std::path::PathBuf; @@ -7,7 +6,6 @@ use std::sync::Arc; use anyhow::{anyhow, bail}; use clap::Args; use indicatif::{MultiProgress, ParallelProgressIterator, ProgressIterator}; -use itertools::Itertools; use log::{debug, info}; use rayon::prelude::*; use rustc_hash::FxHashMap; @@ -102,6 +100,94 @@ pub struct EntryLocalize { batch_size_bp: u64, } +#[derive(Debug, Copy, Clone, Eq, PartialEq)] +enum BedDelimiter { + Tabs, + Whitespace, +} + +impl BedDelimiter { + fn from_first_data_line(line: &str) -> Self { + if line.contains('\t') { + Self::Tabs + } else { + Self::Whitespace + } + } + + fn matches(self, line: &str) -> bool { + match self { + Self::Tabs => line.contains('\t'), + Self::Whitespace => !line.contains('\t'), + } + } + + fn split<'a>(self, line: &'a str) -> Vec<&'a str> { + match self { + Self::Tabs => line.split('\t').collect(), + Self::Whitespace => line.split_whitespace().collect(), + } + } +} + +fn parse_region_core_fields(fields: &[&str]) -> anyhow::Result { + if fields.len() < 3 { + bail!("expected at least 3 BED core columns, found {}", fields.len()) + } + let chrom = fields[0]; + if chrom.is_empty() || chrom.chars().any(char::is_whitespace) { + bail!("invalid chromosome name {chrom:?}") + } + let start = fields[1].parse::().map_err(|_| { + anyhow!( + "invalid start coordinate {:?}, expected a complete u64", + fields[1] + ) + })?; + let end = fields[2].parse::().map_err(|_| { + anyhow!( + "invalid end coordinate {:?}, expected a complete u64", + fields[2] + ) + })?; + let name = fields + .get(3) + .map(|name| { + if name.is_empty() { + bail!("BED name must not be empty") + } + Ok((*name).to_string()) + }) + .transpose()?; + if let Some(score) = fields.get(4) { + if *score != "." { + let score_value = score.parse::().map_err(|_| { + anyhow!( + "invalid BED score {score:?}, expected '.' or a complete \ + finite number" + ) + })?; + if !score_value.is_finite() { + bail!("expected a finite BED score, found {score:?}") + } + } + } + let strand = match fields.get(5) { + None => StrandRule::Both, + Some(&"+") => StrandRule::Positive, + Some(&"-") => StrandRule::Negative, + Some(&".") => StrandRule::Both, + Some(strand) => bail!( + "invalid strand {strand:?}, expected exactly '+', '-', or '.'" + ), + }; + if start > end { + bail!("start {start} is greater than end {end}") + } + + Ok(GenomeRegion { chrom: chrom.to_string(), start, end, strand, name }) +} + impl EntryLocalize { fn load_focus_regions( &self, @@ -111,68 +197,107 @@ impl EntryLocalize { ) -> anyhow::Result> { let pb = multi_progress.add(get_ticker()); pb.set_message("regions parsed"); - let mut reader = BufReader::new(File::open(&self.regions)?) - .lines() - .skip_while(|r| { - r.as_ref().map(|l| l.starts_with('#')).unwrap_or(true) - }) - .peekable(); - let parser = match reader.peek() { - Some(Ok(l)) => { - let num_fields = l.split_whitespace().count(); - if num_fields <= 5 { - |l: &str| GenomeRegion::parse_unstranded_bed_line(l) - } else { - |l: &str| GenomeRegion::parse_stranded_bed_line(l) - } - } - Some(Err(e)) => bail!("failed to inspect regions BED, {e}"), - None => bail!("failed to inspect regions BED, no valid lines"), - }; - - let (regions, errs) = BufReader::new(File::open(&self.regions)?) + let mut data_lines = BufReader::new(File::open(&self.regions)?) .lines() + .enumerate() .progress_with(pb) - .skip_while(|r| { - r.as_ref().map(|l| l.starts_with('#')).unwrap_or(false) - }) - .map(|r| { - r.map_err(|e| { - anyhow!("failed to read from regions BED file, {e}") - }) - .and_then(|l| parser(&l)) - }) - .fold((Vec::new(), HashMap::new()), |(mut acc, mut errs), next| { - match next { - Ok(gr) => { - acc.push(gr); - (acc, errs) - } - Err(e) => { - *errs.entry(e.to_string()).or_insert(0u32) += 1u32; - (acc, errs) + .filter_map(|(line_idx, result)| { + let line_number = line_idx + 1; + match result { + Ok(mut line) => { + if line.trim().is_empty() + || line.trim_start().starts_with('#') + { + None + } else { + let trimmed_len = line + .trim_end_matches(|c: char| { + c.is_ascii_whitespace() + }) + .len(); + line.truncate(trimmed_len); + Some(Ok((line_number, line))) + } } + Err(e) => Some(Err(anyhow!( + "failed to read regions BED file {} at line \ + {line_number}: {e}", + self.regions.display() + ))), } }); - if regions.is_empty() { - let mut err_msg = "failure reasons: ".to_string(); - for (err, count) in - errs.into_iter().sorted_by(|(_, x), (_, y)| x.cmp(y)) - { - err_msg.push_str(&format!("\t{err}: {count}\n")); - } - bail!("failed to load any regions, {err_msg}") - } + let Some(first_line) = data_lines.next().transpose()? else { + bail!( + "failed to inspect regions BED {}: no data lines", + self.regions.display() + ) + }; + let first_line_number = first_line.0; + let delimiter = BedDelimiter::from_first_data_line(&first_line.1); + let expected_fields = delimiter.split(&first_line.1).len(); + let regions = std::iter::once(Ok(first_line)) + .chain(data_lines) + .map(|result| { + let (line_number, line) = result?; + if !delimiter.matches(&line) { + bail!( + "failed to parse regions BED file {} at line \ + {line_number}: delimiter mode changed from \ + {delimiter:?}, selected at line {first_line_number}", + self.regions.display() + ) + } + let fields = delimiter.split(&line); + let actual_fields = fields.len(); + if actual_fields != expected_fields { + bail!( + "failed to parse regions BED file {} at line \ + {line_number}: expected {expected_fields} BED fields \ + to match line {first_line_number}, found \ + {actual_fields}", + self.regions.display() + ) + } + parse_region_core_fields(&fields) + .map(|region| (line_number, region)) + .map_err(|e| { + anyhow!( + "failed to parse regions BED file {} at line \ + {line_number}: {e}", + self.regions.display() + ) + }) + }) + .collect::>>()?; + let (regions, missing_from_sizes, missing_from_index) = - regions.into_iter().fold( + regions.into_iter().try_fold( (Vec::new(), 0usize, 0usize), - |(mut acc, no_size, no_index), mut next| { - if !sequence_lengths.contains_key(&next.chrom) { - (acc, no_size + 1, no_index) + |(mut acc, no_size, no_index), (line_number, mut next)| { + let contig_length = + sequence_lengths.get(&next.chrom).copied(); + if let Some(contig_length) = contig_length { + if next.start > contig_length + || next.end > contig_length + { + return Err(anyhow!( + "invalid region in {} at line {line_number}: \ + {}:{}-{} exceeds contig length \ + {contig_length}", + self.regions.display(), + next.chrom, + next.start, + next.end + )); + } + } + if contig_length.is_none() { + Ok((acc, no_size + 1, no_index)) } else if !index.has_contig(&next.chrom) { - (acc, no_size, no_index + 1) + Ok((acc, no_size, no_index + 1)) } else { + let contig_length = contig_length.unwrap(); let mp = next.midpoint(); let start = mp .checked_sub(self.expand_window + 1) @@ -180,15 +305,15 @@ impl EntryLocalize { // safe because of conditional above let end = std::cmp::min( mp.saturating_add(self.expand_window), - *sequence_lengths.get(&next.chrom).unwrap(), + contig_length, ); next.start = start; next.end = end; acc.push(next); - (acc, no_size, no_index) + Ok((acc, no_size, no_index)) } }, - ); + )?; if missing_from_sizes > 0usize { debug!( @@ -232,20 +357,6 @@ impl EntryLocalize { let pool = rayon::ThreadPoolBuilder::new() .num_threads(self.threads) .build()?; - let writer: Box = - if let Some(out_fp) = self.out_file.as_ref() { - if self.force { - Box::new(BufWriter::new(File::create(out_fp)?)) - } else { - Box::new(BufWriter::new(File::create_new(out_fp)?)) - } - } else { - Box::new(BufWriter::new(stdout())) - }; - - let writer = - csv::WriterBuilder::new().delimiter('\t' as u8).from_writer(writer); - info!("loading sequence lengths from {:?}", &self.genome_sizes); let sequence_lengths = read_sequence_lengths_file(&self.genome_sizes)? @@ -258,6 +369,18 @@ impl EntryLocalize { &tabix_index, &multi_progress, )?; + let writer: Box = + if let Some(out_fp) = self.out_file.as_ref() { + if self.force { + Box::new(BufWriter::new(File::create(out_fp)?)) + } else { + Box::new(BufWriter::new(File::create_new(out_fp)?)) + } + } else { + Box::new(BufWriter::new(stdout())) + }; + let writer = + csv::WriterBuilder::new().delimiter('\t' as u8).from_writer(writer); info!("loaded {} regions", genome_regions.len()); let successes = @@ -308,3 +431,443 @@ impl EntryLocalize { Ok(()) } } + +#[cfg(test)] +mod tests { + use std::fs::write; + use std::path::PathBuf; + + use indicatif::MultiProgress; + use rustc_hash::FxHashMap; + + use super::EntryLocalize; + use crate::dmr::bedmethyl::BedMethylLine; + use crate::tabix::HtsTabixHandler; + use crate::util::{GenomeRegion, StrandRule}; + + fn test_bedmethyl() -> PathBuf { + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join( + "../tests/resources/\ + lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz", + ) + } + + fn load_regions_with_sequence_lengths( + regions: PathBuf, + sequence_lengths: FxHashMap, + ) -> anyhow::Result> { + let entry = EntryLocalize { + in_bedmethyl: test_bedmethyl(), + regions, + chart_filepath: None, + chart_name: None, + expand_window: 10, + stranded: None, + stranded_features: None, + min_coverage: 1, + genome_sizes: PathBuf::new(), + out_file: None, + log_filepath: None, + threads: 1, + io_threads: 1, + force: false, + batch_size_bp: 1, + }; + let index = + HtsTabixHandler::::from_path(&entry.in_bedmethyl)?; + + entry.load_focus_regions( + &sequence_lengths, + &index, + &MultiProgress::new(), + ) + } + + fn load_regions(regions: PathBuf) -> anyhow::Result> { + load_regions_with_sequence_lengths( + regions, + FxHashMap::from_iter([("chr20".to_string(), 100_000_000)]), + ) + } + + #[test] + fn load_focus_regions_rejects_malformed_row_after_valid_row() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("mixed-regions.bed"); + write( + ®ions, + "# test regions\n\ + chr20\t9681998\t9681999\n\ + not-a-valid-bed-row\n\ + chr20\t9682013\t9682014\n", + ) + .unwrap(); + + let error = load_regions(regions.clone()).unwrap_err(); + let message = error.to_string(); + + assert!( + message.contains(regions.to_string_lossy().as_ref()), + "missing region filepath in error: {message}" + ); + assert!(message.contains("line 3"), "missing line number: {message}"); + assert!( + message.contains("delimiter mode changed"), + "missing parse reason: {message}" + ); + } + + #[test] + fn load_focus_regions_accepts_valid_rows() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("valid-regions.bed"); + write( + ®ions, + "# test regions\n\ + chr20\t9681998\t9681999\n\ + chr20\t9682013\t9682014\n", + ) + .unwrap(); + + assert_eq!(load_regions(regions).unwrap().len(), 2); + } + + #[test] + fn load_focus_regions_preserves_bed4_names_with_spaces() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("bed4-regions.bed"); + write(®ions, "chr20\t9681998\t9681999\tregion name with spaces\n") + .unwrap(); + + let regions = load_regions(regions).unwrap(); + + assert_eq!(regions.len(), 1); + assert_eq!(regions[0].name.as_deref(), Some("region name with spaces")); + assert_eq!(regions[0].strand, StrandRule::Both); + } + + #[test] + fn load_focus_regions_accepts_valid_bed5_score_and_opaque_extensions() { + for (name, row, expected_strand) in [ + ( + "bed5", + "chr20\t9681998\t9681999\tregion name\t500.5", + StrandRule::Both, + ), + ( + "bed8", + "chr20\t9681998\t9681999\tregion name\t.\t-\t\topaque", + StrandRule::Negative, + ), + ] { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join(format!("{name}-regions.bed")); + write(®ions, format!("{row}\n")).unwrap(); + + let parsed = load_regions(regions).unwrap(); + + assert_eq!(parsed.len(), 1); + assert_eq!(parsed[0].name.as_deref(), Some("region name")); + assert_eq!(parsed[0].strand, expected_strand); + } + } + + #[test] + fn load_focus_regions_ignores_spacing_and_preserves_exact_strands() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("stranded-regions.bed"); + write( + ®ions, + " \n\ + # leading comment\n\ + chr20\t9681998\t9681999\tplus\t0\t+\n\ + \t\n\ + # interleaved comment\n\ + chr20\t9682013\t9682014\tminus\t0\t-\n\ + chr20\t9682030\t9682030\tzero-length\t0\t.\n", + ) + .unwrap(); + + let regions = load_regions(regions).unwrap(); + + assert_eq!(regions.len(), 3); + assert_eq!(regions[0].strand, StrandRule::Positive); + assert_eq!(regions[1].strand, StrandRule::Negative); + assert_eq!(regions[2].strand, StrandRule::Both); + } + + #[test] + fn load_focus_regions_rejects_mixed_bed_field_counts() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("mixed-schema-regions.bed"); + write( + ®ions, + "chr20\t9681998\t9681999\n\ + chr20\t9682013\t9682014\tminus\t0\t-\n", + ) + .unwrap(); + + let error = load_regions(regions.clone()).unwrap_err(); + let message = error.to_string(); + + assert!( + message.contains(regions.to_string_lossy().as_ref()), + "missing region filepath in error: {message}" + ); + assert!(message.contains("line 2"), "missing line number: {message}"); + assert!( + message.contains("expected 3 BED fields") + && message.contains("found 6"), + "missing field-count reason: {message}" + ); + } + + #[test] + fn load_focus_regions_accepts_legacy_space_delimited_bed6() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("space-delimited-regions.bed"); + write(®ions, "chr20 9681998 9681999 name 0 +\n").unwrap(); + + let parsed = load_regions(regions).unwrap(); + + assert_eq!(parsed.len(), 1); + assert_eq!(parsed[0].name.as_deref(), Some("name")); + assert_eq!(parsed[0].strand, StrandRule::Positive); + } + + #[test] + fn load_focus_regions_rejects_delimiter_mode_changes() { + for (name, rows) in [ + ( + "space-then-tabs", + "chr20 9681998 9681999 name 0 +\n\ + chr20\t9682013\t9682014\tname\t0\t-\n", + ), + ( + "tabs-then-space", + "chr20\t9681998\t9681999\tname\t0\t+\n\ + chr20 9682013 9682014 name 0 -\n", + ), + ] { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join(format!("{name}.bed")); + write(®ions, rows).unwrap(); + + let error = load_regions(regions.clone()).unwrap_err(); + let message = error.to_string(); + + assert!( + message.contains(regions.to_string_lossy().as_ref()), + "missing region filepath in error: {message}" + ); + assert!(message.contains("line 2"), "{message}"); + assert!(message.contains("delimiter mode"), "{message}"); + } + } + + #[test] + fn load_focus_regions_rejects_invalid_consumed_bed_fields() { + for (name, row, reason) in [ + ( + "partial-end", + "chr20\t9681998\t9681999x", + "invalid end coordinate", + ), + ( + "invalid-score", + "chr20\t9681998\t9681999\tname\t500x", + "invalid BED score", + ), + ( + "non-finite-score", + "chr20\t9681998\t9681999\tname\tNaN", + "finite BED score", + ), + ( + "partial-strand", + "chr20\t9681998\t9681999\tname\t0\t+garbage", + "invalid strand", + ), + ] { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join(format!("{name}-regions.bed")); + write(®ions, format!("{row}\n")).unwrap(); + + let error = load_regions(regions.clone()).unwrap_err(); + let message = error.to_string(); + + assert!( + message.contains(regions.to_string_lossy().as_ref()), + "missing region filepath in error: {message}" + ); + assert!(message.contains("line 1"), "{message}"); + assert!(message.contains(reason), "{message}"); + } + } + + #[test] + fn load_focus_regions_rejects_interior_empty_bed_name() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("empty-name-regions.bed"); + write(®ions, "chr20\t9681998\t9681999\t\t500\n").unwrap(); + + let error = load_regions(regions.clone()).unwrap_err(); + let message = error.to_string(); + + assert!( + message.contains(regions.to_string_lossy().as_ref()), + "missing region filepath in error: {message}" + ); + assert!(message.contains("line 1"), "{message}"); + assert!(message.contains("BED name must not be empty"), "{message}"); + } + + #[test] + fn load_focus_regions_preserves_trailing_delimiter_compatibility() { + for (name, row) in [ + ("tab", "chr20\t9681998\t9681999\t\t "), + ("space", "chr20 9681998 9681999 "), + ] { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join(format!("{name}-trailing.bed")); + write(®ions, format!("{row}\n")).unwrap(); + + let parsed = load_regions(regions).unwrap(); + + assert_eq!(parsed.len(), 1); + assert_eq!(parsed[0].name, None); + assert_eq!(parsed[0].strand, StrandRule::Both); + } + } + + #[test] + fn load_focus_regions_rejects_missing_strand_after_bed6_schema() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("missing-strand-regions.bed"); + write( + ®ions, + "chr20\t9681998\t9681999\tplus\t0\t+\n\ + chr20\t9682013\t9682014\tmissing\t0\t\n", + ) + .unwrap(); + + let error = load_regions(regions.clone()).unwrap_err(); + let message = error.to_string(); + + assert!( + message.contains(regions.to_string_lossy().as_ref()), + "missing region filepath in error: {message}" + ); + assert!(message.contains("line 2"), "{message}"); + assert!( + message.contains("expected 6 BED fields") + && message.contains("found 5"), + "{message}" + ); + } + + #[test] + fn load_focus_regions_rejects_headers_and_browser_directives() { + for (name, first_line) in [ + ("header", "chrom\tstart\tend"), + ("track", "track name=test"), + ("browser", "browser position chr20:9681998-9682014"), + ] { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join(format!("{name}-regions.bed")); + write( + ®ions, + format!( + "{first_line}\n\ + chr20\t9681998\t9681999\n" + ), + ) + .unwrap(); + + let error = load_regions(regions.clone()).unwrap_err(); + let message = error.to_string(); + + assert!( + message.contains(regions.to_string_lossy().as_ref()), + "missing region filepath in error: {message}" + ); + assert!(message.contains("line 1"), "{message}"); + assert!(message.contains("failed to parse"), "{message}"); + } + } + + #[test] + fn load_focus_regions_rejects_invalid_coordinate_ranges() { + for (name, row, reason) in [ + ( + "reversed", + "chr20\t9682014\t9682013", + "start 9682014 is greater than end 9682013", + ), + ( + "past-contig-end", + "chr20\t100000001\t100000001", + "exceeds contig length 100000000", + ), + ] { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join(format!("{name}-regions.bed")); + write(®ions, format!("{row}\n")).unwrap(); + + let error = load_regions(regions.clone()).unwrap_err(); + let message = error.to_string(); + + assert!( + message.contains(regions.to_string_lossy().as_ref()), + "missing region filepath in error: {message}" + ); + assert!(message.contains("line 1"), "{message}"); + assert!(message.contains(reason), "{message}"); + } + } + + #[test] + fn load_focus_regions_accepts_zero_length_at_contig_end() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("boundary-region.bed"); + write(®ions, "chr20\t100000000\t100000000\n").unwrap(); + + assert_eq!(load_regions(regions).unwrap().len(), 1); + } + + #[test] + fn load_focus_regions_validates_coordinates_before_missing_contigs() { + for (name, row, sequence_lengths, reason) in [ + ( + "reversed-unknown-contig", + "unknown\t10\t5", + FxHashMap::default(), + "start 10 is greater than end 5", + ), + ( + "past-known-unindexed-contig", + "known-unindexed\t101\t101", + FxHashMap::from_iter([("known-unindexed".to_string(), 100)]), + "exceeds contig length 100", + ), + ] { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join(format!("{name}.bed")); + write(®ions, format!("{row}\n")).unwrap(); + + let error = load_regions_with_sequence_lengths( + regions.clone(), + sequence_lengths, + ) + .unwrap_err(); + let message = error.to_string(); + + assert!( + message.contains(regions.to_string_lossy().as_ref()), + "missing region filepath in error: {message}" + ); + assert!(message.contains("line 1"), "{message}"); + assert!(message.contains(reason), "{message}"); + } + } +} diff --git a/modkit/tests/test_localize.rs b/modkit/tests/test_localize.rs index 22257737..14016d77 100644 --- a/modkit/tests/test_localize.rs +++ b/modkit/tests/test_localize.rs @@ -1,7 +1,42 @@ use crate::common::run_modkit; +use std::fs::{read_to_string, write}; +use std::path::Path; +use std::process::{Command, Output}; mod common; +fn run_localize( + regions: &Path, + genome_sizes: &Path, + output: &Path, + force: bool, +) -> Output { + let mut command = Command::new(env!("CARGO_BIN_EXE_modkit")); + command.args([ + "localize", + "../tests/resources/\ + lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz", + "--regions", + regions.to_str().unwrap(), + "--genome-sizes", + genome_sizes.to_str().unwrap(), + "--window", + "10", + "--min-coverage", + "1", + "--threads", + "1", + "--io-threads", + "1", + "--out-file", + output.to_str().unwrap(), + ]); + if force { + command.arg("--force"); + } + command.output().unwrap() +} + #[test] fn test_localise_helps() { let _ = run_modkit(&["localize", "--help"]) @@ -9,3 +44,74 @@ fn test_localise_helps() { let _ = run_modkit(&["localise", "--help"]) .expect("failed to run modkit localise help"); } + +#[test] +fn test_localize_rejects_mixed_regions_without_mutating_output() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("mixed-regions.bed"); + let genome_sizes = temp_dir.path().join("genome-sizes.tsv"); + let new_output = temp_dir.path().join("new-localize.tsv"); + let existing_output = temp_dir.path().join("existing-localize.tsv"); + write( + ®ions, + "# test regions\n\ + chr20\t9681998\t9681999\n\ + not-a-valid-bed-row\n", + ) + .unwrap(); + write(&genome_sizes, "chr20\t100000000\n").unwrap(); + + let create_result = + run_localize(®ions, &genome_sizes, &new_output, false); + let sentinel = "keep existing output content\n"; + write(&existing_output, sentinel).unwrap(); + let force_result = + run_localize(®ions, &genome_sizes, &existing_output, true); + + let create_stderr = String::from_utf8_lossy(&create_result.stderr); + assert!(!create_result.status.success(), "{create_stderr}"); + assert!( + create_stderr.contains(regions.to_string_lossy().as_ref()), + "failure did not reach region parsing: {create_stderr}" + ); + assert!(create_stderr.contains("line 3"), "{create_stderr}"); + assert!( + create_stderr.contains("delimiter mode changed"), + "{create_stderr}" + ); + assert!(!force_result.status.success()); + assert!( + !new_output.exists() + && read_to_string(&existing_output).unwrap() == sentinel, + "failed region parsing mutated output paths: new_exists={}, \ + existing_content={:?}", + new_output.exists(), + read_to_string(existing_output).unwrap() + ); +} + +#[test] +fn test_localize_accepts_valid_regions() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("valid-regions.bed"); + let genome_sizes = temp_dir.path().join("genome-sizes.tsv"); + let output = temp_dir.path().join("localize.tsv"); + write( + ®ions, + "# test regions\n\ + chr20\t9681998\t9681999\n\ + chr20\t9682013\t9682014\n", + ) + .unwrap(); + write(&genome_sizes, "chr20\t100000000\n").unwrap(); + + let result = run_localize(®ions, &genome_sizes, &output, false); + assert!( + result.status.success(), + "{}", + String::from_utf8_lossy(&result.stderr) + ); + assert!(read_to_string(output) + .unwrap() + .starts_with("mod_code\toffset\tn_valid\tn_mod\tpercent_modified\n")); +} From 076f0fd9a69d1661134bd6087d52043ddb67f7d6 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 17:33:24 -0700 Subject: [PATCH 2/9] test: expose swallowed regional query errors --- modkit/tests/common/mod.rs | 2 + modkit/tests/common/regional_query.rs | 50 +++++++++++ modkit/tests/test_localize.rs | 121 ++++++++++++++++++++++++++ modkit/tests/test_stats.rs | 102 ++++++++++++++++++++++ 4 files changed, 275 insertions(+) create mode 100644 modkit/tests/common/regional_query.rs create mode 100644 modkit/tests/test_stats.rs diff --git a/modkit/tests/common/mod.rs b/modkit/tests/common/mod.rs index b05c06e4..d3869a1b 100644 --- a/modkit/tests/common/mod.rs +++ b/modkit/tests/common/mod.rs @@ -12,6 +12,8 @@ use std::io::Read; use std::path::{Path, PathBuf}; use std::process::Output; +pub mod regional_query; + pub fn run_modkit(args: &[&str]) -> AnyhowResult { let exe = Path::new(env!("CARGO_BIN_EXE_modkit")); assert!(exe.exists()); diff --git a/modkit/tests/common/regional_query.rs b/modkit/tests/common/regional_query.rs new file mode 100644 index 00000000..96f78e8b --- /dev/null +++ b/modkit/tests/common/regional_query.rs @@ -0,0 +1,50 @@ +use std::ffi::CString; +use std::io::Write; +use std::path::{Path, PathBuf}; + +const VALID_ROWS: [&str; 4] = [ + "chr1\t100\t101\tm\t4\t+\t100\t101\t255,0,0\t4 25.00 1 3 0 0 0 0 0\n", + "chr1\t200\t201\tm\t6\t+\t200\t201\t255,0,0\t6 50.00 3 3 0 0 0 0 0\n", + "chr1\t300\t301\tm\t8\t-\t300\t301\t255,0,0\t8 50.00 4 4 0 0 0 0 0\n", + "chr1\t400\t401\tm\t10\t-\t400\t401\t255,0,0\t10 20.00 2 8 0 0 0 0 0\n", +]; + +const INVALID_ROW: &str = + "chr1\t200\t201\tm\tnot-a-number\t+\t200\t201\t255,0,0\t6 50.00 3 3 0 0 0 0 0\n"; + +pub const REGIONS: &str = "chr1\t100\t101\n\ + chr1\t200\t201\n\ + chr1\t300\t301\n\ + chr1\t400\t401\n\ + absent\t0\t1\n"; + +pub const GENOME_SIZES: &str = "chr1\t1000\nabsent\t10\n"; + +pub const OUTPUT_SENTINEL: &[u8] = b"existing output must survive\n"; + +pub fn make_indexed_bedmethyl(directory: &Path, malformed: bool) -> PathBuf { + let path = directory.join(if malformed { + "malformed.bed.gz" + } else { + "valid.bed.gz" + }); + let mut writer = rust_htslib::bgzf::Writer::from_path(&path).unwrap(); + for (index, row) in VALID_ROWS.iter().enumerate() { + let row = if malformed && index == 1 { INVALID_ROW } else { row }; + writer.write_all(row.as_bytes()).unwrap(); + } + writer.flush().unwrap(); + drop(writer); + + let path_cstr = CString::new(path.to_str().unwrap()).unwrap(); + let result = unsafe { + rust_htslib::htslib::tbx_index_build( + path_cstr.as_ptr(), + 0, + &rust_htslib::htslib::tbx_conf_bed, + ) + }; + assert_eq!(result, 0, "failed to index {}", path.display()); + assert!(PathBuf::from(format!("{}.tbi", path.display())).is_file()); + path +} diff --git a/modkit/tests/test_localize.rs b/modkit/tests/test_localize.rs index 14016d77..bb7c4a0c 100644 --- a/modkit/tests/test_localize.rs +++ b/modkit/tests/test_localize.rs @@ -1,3 +1,6 @@ +use crate::common::regional_query::{ + make_indexed_bedmethyl, GENOME_SIZES, OUTPUT_SENTINEL, REGIONS, +}; use crate::common::run_modkit; use std::fs::{read_to_string, write}; use std::path::Path; @@ -37,6 +40,53 @@ fn run_localize( command.output().unwrap() } +fn run_localize_query( + bedmethyl: &Path, + regions: &Path, + genome_sizes: &Path, + output: &Path, + threads: usize, + force: bool, +) -> Output { + let mut command = Command::new(env!("CARGO_BIN_EXE_modkit")); + command + .arg("localize") + .arg(bedmethyl) + .arg("--regions") + .arg(regions) + .arg("--genome-sizes") + .arg(genome_sizes) + .arg("--out-file") + .arg(output) + .args([ + "--window", + "1", + "--min-coverage", + "1", + "--threads", + &threads.to_string(), + "--io-threads", + "1", + ]); + if force { + command.arg("--force"); + } + command.output().unwrap() +} + +fn assert_query_failure(output: &Output) { + let stderr = String::from_utf8_lossy(&output.stderr); + assert!( + !output.status.success(), + "localize unexpectedly succeeded: {stderr}" + ); + assert!( + stderr.contains("failed to localize region chr1:") + && stderr.contains("invalid-bedmethyl-data"), + "localize error lacks region/decode context: {stderr}" + ); +} + #[test] fn test_localise_helps() { let _ = run_modkit(&["localize", "--help"]) @@ -115,3 +165,74 @@ fn test_localize_accepts_valid_regions() { .unwrap() .starts_with("mod_code\toffset\tn_valid\tn_mod\tpercent_modified\n")); } + +#[test] +fn localize_propagates_regional_query_errors_without_mutating_output() { + for threads in [1usize, 4] { + let temp_dir = tempfile::tempdir().unwrap(); + let bedmethyl = make_indexed_bedmethyl(temp_dir.path(), true); + let regions = temp_dir.path().join("regions.bed"); + let genome_sizes = temp_dir.path().join("genome-sizes.tsv"); + write(®ions, REGIONS).unwrap(); + write(&genome_sizes, GENOME_SIZES).unwrap(); + + let absent_output = temp_dir.path().join("absent.tsv"); + assert_query_failure(&run_localize_query( + &bedmethyl, + ®ions, + &genome_sizes, + &absent_output, + threads, + false, + )); + assert!( + !absent_output.exists(), + "failed localize run created output with {threads} threads" + ); + + let existing_output = temp_dir.path().join("existing.tsv"); + write(&existing_output, OUTPUT_SENTINEL).unwrap(); + assert_query_failure(&run_localize_query( + &bedmethyl, + ®ions, + &genome_sizes, + &existing_output, + threads, + true, + )); + assert_eq!(std::fs::read(existing_output).unwrap(), OUTPUT_SENTINEL); + } +} + +#[test] +fn localize_valid_output_is_byte_identical_across_thread_counts() { + let temp_dir = tempfile::tempdir().unwrap(); + let bedmethyl = make_indexed_bedmethyl(temp_dir.path(), false); + let regions = temp_dir.path().join("regions.bed"); + let genome_sizes = temp_dir.path().join("genome-sizes.tsv"); + write(®ions, REGIONS).unwrap(); + write(&genome_sizes, GENOME_SIZES).unwrap(); + + let one_thread = temp_dir.path().join("one.tsv"); + let many_threads = temp_dir.path().join("many.tsv"); + for (threads, output) in [(1usize, &one_thread), (4, &many_threads)] { + let result = run_localize_query( + &bedmethyl, + ®ions, + &genome_sizes, + output, + threads, + false, + ); + assert!( + result.status.success(), + "valid localize run failed: {}", + String::from_utf8_lossy(&result.stderr) + ); + } + + let expected = b"mod_code\toffset\tn_valid\tn_mod\tpercent_modified\n\ +m\t-1\t28\t10\t35.714287\n"; + assert_eq!(std::fs::read(one_thread).unwrap(), expected); + assert_eq!(std::fs::read(many_threads).unwrap(), expected); +} diff --git a/modkit/tests/test_stats.rs b/modkit/tests/test_stats.rs new file mode 100644 index 00000000..4a6eb47e --- /dev/null +++ b/modkit/tests/test_stats.rs @@ -0,0 +1,102 @@ +use std::fs; +use std::path::Path; +use std::process::{Command, Output}; + +use common::regional_query::{ + make_indexed_bedmethyl, OUTPUT_SENTINEL, REGIONS, +}; + +mod common; + +fn run_stats( + bedmethyl: &Path, + regions: &Path, + output: &Path, + threads: usize, + force: bool, +) -> Output { + let mut command = Command::new(env!("CARGO_BIN_EXE_modkit")); + command + .arg("stats") + .arg(bedmethyl) + .arg("--regions") + .arg(regions) + .arg("--out-table") + .arg(output) + .args(["--threads", &threads.to_string(), "--io-threads", "1"]); + if force { + command.arg("--force"); + } + command.output().unwrap() +} + +fn assert_query_failure(output: &Output) { + let stderr = String::from_utf8_lossy(&output.stderr); + assert!(!output.status.success(), "stats unexpectedly succeeded: {stderr}"); + assert!( + stderr.contains("failed to calculate stats for region chr1:") + && stderr.contains("invalid-bedmethyl-data"), + "stats error lacks region/decode context: {stderr}" + ); +} + +#[test] +fn stats_propagates_regional_query_errors_without_mutating_output() { + for threads in [1usize, 4] { + let temp_dir = tempfile::tempdir().unwrap(); + let bedmethyl = make_indexed_bedmethyl(temp_dir.path(), true); + let regions = temp_dir.path().join("regions.bed"); + fs::write(®ions, REGIONS).unwrap(); + + let absent_output = temp_dir.path().join("absent.tsv"); + assert_query_failure(&run_stats( + &bedmethyl, + ®ions, + &absent_output, + threads, + false, + )); + assert!( + !absent_output.exists(), + "failed stats run created output with {threads} threads" + ); + + let existing_output = temp_dir.path().join("existing.tsv"); + fs::write(&existing_output, OUTPUT_SENTINEL).unwrap(); + assert_query_failure(&run_stats( + &bedmethyl, + ®ions, + &existing_output, + threads, + true, + )); + assert_eq!(fs::read(existing_output).unwrap(), OUTPUT_SENTINEL); + } +} + +#[test] +fn stats_valid_output_is_byte_identical_across_thread_counts() { + let temp_dir = tempfile::tempdir().unwrap(); + let bedmethyl = make_indexed_bedmethyl(temp_dir.path(), false); + let regions = temp_dir.path().join("regions.bed"); + fs::write(®ions, REGIONS).unwrap(); + + let one_thread = temp_dir.path().join("one.tsv"); + let many_threads = temp_dir.path().join("many.tsv"); + for (threads, output) in [(1usize, &one_thread), (4, &many_threads)] { + let result = run_stats(&bedmethyl, ®ions, output, threads, false); + assert!( + result.status.success(), + "valid stats run failed: {}", + String::from_utf8_lossy(&result.stderr) + ); + } + + let expected = b"#chrom\tstart\tend\tname\tstrand\tcount_m\tcount_valid_m\tpercent_m\n\ +chr1\t100\t101\t.\t.\t1\t4\t25\n\ +chr1\t200\t201\t.\t.\t3\t6\t50\n\ +chr1\t300\t301\t.\t.\t4\t8\t50\n\ +chr1\t400\t401\t.\t.\t2\t10\t20\n"; + assert_eq!(fs::read(one_thread).unwrap(), expected); + assert_eq!(fs::read(many_threads).unwrap(), expected); +} From 77cfd4981a99771addde90bea91db4f79b2f780c Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 17:39:15 -0700 Subject: [PATCH 3/9] fix: propagate regional query failures --- modkit-core/src/localise/subcommand.rs | 45 ++++++++++----------- modkit-core/src/stats/subcommand.rs | 54 +++++++++++++------------- 2 files changed, 50 insertions(+), 49 deletions(-) diff --git a/modkit-core/src/localise/subcommand.rs b/modkit-core/src/localise/subcommand.rs index 6f2b55ef..61b8a66b 100644 --- a/modkit-core/src/localise/subcommand.rs +++ b/modkit-core/src/localise/subcommand.rs @@ -3,7 +3,7 @@ use std::io::{stdout, BufRead, BufReader, BufWriter, Write}; use std::path::PathBuf; use std::sync::Arc; -use anyhow::{anyhow, bail}; +use anyhow::{anyhow, bail, Context}; use clap::Args; use indicatif::{MultiProgress, ParallelProgressIterator, ProgressIterator}; use log::{debug, info}; @@ -369,18 +369,6 @@ impl EntryLocalize { &tabix_index, &multi_progress, )?; - let writer: Box = - if let Some(out_fp) = self.out_file.as_ref() { - if self.force { - Box::new(BufWriter::new(File::create(out_fp)?)) - } else { - Box::new(BufWriter::new(File::create_new(out_fp)?)) - } - } else { - Box::new(BufWriter::new(stdout())) - }; - let writer = - csv::WriterBuilder::new().delimiter('\t' as u8).from_writer(writer); info!("loaded {} regions", genome_regions.len()); let successes = @@ -392,27 +380,40 @@ impl EntryLocalize { .into_par_iter() .progress_with(successes) .map(|gr| { + let region = + format!("{}:{}-{}", gr.chrom, gr.start, gr.end); gr.into_localized_mod_counts( &tabix_index, self.stranded_features, stranded_features, self.io_threads, ) + .with_context(|| { + format!("failed to localize region {region}") + }) }) - .fold( + .try_fold( || LocalizedModCounts::zero(), - |counts, next| match next { - Ok(lc) => counts.op(lc), - Err(e) => { - debug!("region failed, {e}"); - counts - } + |counts, next| -> anyhow::Result<_> { + Ok(counts.op(next?)) }, ) - .reduce(|| LocalizedModCounts::zero(), |a, b| a.op(b)) - }); + .try_reduce(|| LocalizedModCounts::zero(), |a, b| Ok(a.op(b))) + })?; let table = counts.get_table(&multi_progress); + let writer: Box = + if let Some(out_fp) = self.out_file.as_ref() { + if self.force { + Box::new(BufWriter::new(File::create(out_fp)?)) + } else { + Box::new(BufWriter::new(File::create_new(out_fp)?)) + } + } else { + Box::new(BufWriter::new(stdout())) + }; + let writer = + csv::WriterBuilder::new().delimiter('\t' as u8).from_writer(writer); table.to_csv_writer(writer)?; if let Some(p) = self.chart_filepath.as_ref() { diff --git a/modkit-core/src/stats/subcommand.rs b/modkit-core/src/stats/subcommand.rs index 5a35c9d8..99c2f522 100644 --- a/modkit-core/src/stats/subcommand.rs +++ b/modkit-core/src/stats/subcommand.rs @@ -2,7 +2,7 @@ use std::fs::File; use std::io::{stdout, BufRead, BufReader, BufWriter, Write}; use std::path::PathBuf; -use anyhow::{anyhow, bail}; +use anyhow::{anyhow, bail, Context}; use clap::Args; use indicatif::{ParallelProgressIterator, ProgressIterator}; use itertools::Itertools; @@ -69,17 +69,6 @@ impl EntryStats { let _ = init_logging(self.log_filepath.as_ref()); let index: HtsTabixHandler = HtsTabixHandler::from_path(&self.in_bedmethyl)?; - let handle: Box = match self.out_table.as_str() { - "-" | "stdout" => Box::new(BufWriter::new(stdout())), - p @ _ => { - let fp = std::path::Path::new(p); - if self.force { - Box::new(BufWriter::new(File::create(fp)?)) - } else { - Box::new(BufWriter::new(File::create_new(fp)?)) - } - } - }; let pool = rayon::ThreadPoolBuilder::new() .num_threads(self.threads) .build()?; @@ -159,37 +148,37 @@ impl EntryStats { genome_regions .into_par_iter() .progress_with(stats_pb) - .filter_map(|gr| { - match gr.into_stats( + .map(|gr| { + let region = + format!("{}:{}-{}", gr.chrom, gr.start, gr.end); + gr.into_stats( &index, self.min_coverage, mod_codes.as_ref(), self.io_threads, - ) { - Ok(stats) => Some(stats), - Err(e) => { - debug!("failed to get stats, {e}"); - None - } - } + ) + .with_context(|| { + format!("failed to calculate stats for region {region}") + }) }) - .fold( + .try_fold( || (Vec::new(), FxHashSet::default()), - |(mut agg, mut codes), next| { + |(mut agg, mut codes), next| -> anyhow::Result<_> { + let next = next?; if mod_codes.is_none() { codes.extend( next.per_mod_methylation.keys().copied(), ); } agg.push(next); - (agg, codes) + Ok((agg, codes)) }, ) - .reduce( + .try_reduce( || (Vec::new(), FxHashSet::default()), - |(a, b), (c, d)| (a.op(c), b.op(d)), + |(a, b), (c, d)| Ok((a.op(c), b.op(d))), ) - }); + })?; let mod_codes = if let Some(codes) = mod_codes { codes } else { obs_codes } @@ -209,6 +198,17 @@ impl EntryStats { table.add_row(x.into_row(&mod_codes)); }); + let handle: Box = match self.out_table.as_str() { + "-" | "stdout" => Box::new(BufWriter::new(stdout())), + p @ _ => { + let fp = std::path::Path::new(p); + if self.force { + Box::new(BufWriter::new(File::create(fp)?)) + } else { + Box::new(BufWriter::new(File::create_new(fp)?)) + } + } + }; let csv_writer = csv::WriterBuilder::new().delimiter('\t' as u8).from_writer(handle); table.to_csv_writer(csv_writer)?; From 3e7f0f80383ed9dd128def7fdf60a35cdaa71976 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 02:33:46 -0700 Subject: [PATCH 4/9] Fix localize window geometry --- modkit-core/src/localise/subcommand.rs | 34 +++++------ modkit-core/src/localise/util.rs | 79 ++++++++++++++++++++++++-- modkit/tests/test_localize.rs | 32 +++++++++++ 3 files changed, 122 insertions(+), 23 deletions(-) diff --git a/modkit-core/src/localise/subcommand.rs b/modkit-core/src/localise/subcommand.rs index 61b8a66b..1721db69 100644 --- a/modkit-core/src/localise/subcommand.rs +++ b/modkit-core/src/localise/subcommand.rs @@ -13,7 +13,9 @@ use rustc_hash::FxHashMap; use modkit_logging::init_logging; use crate::dmr::bedmethyl::BedMethylLine; -use crate::localise::util::{LocalizedModCounts, StrandedFeatures}; +use crate::localise::util::{ + FocusRegion, LocalizedModCounts, StrandedFeatures, +}; use crate::monoid::Moniod; use crate::tabix::HtsTabixHandler; use crate::util::{ @@ -194,7 +196,7 @@ impl EntryLocalize { sequence_lengths: &FxHashMap, index: &HtsTabixHandler, multi_progress: &MultiProgress, - ) -> anyhow::Result> { + ) -> anyhow::Result> { let pb = multi_progress.add(get_ticker()); pb.set_message("regions parsed"); let mut data_lines = BufReader::new(File::open(&self.regions)?) @@ -274,7 +276,7 @@ impl EntryLocalize { let (regions, missing_from_sizes, missing_from_index) = regions.into_iter().try_fold( (Vec::new(), 0usize, 0usize), - |(mut acc, no_size, no_index), (line_number, mut next)| { + |(mut acc, no_size, no_index), (line_number, next)| { let contig_length = sequence_lengths.get(&next.chrom).copied(); if let Some(contig_length) = contig_length { @@ -297,19 +299,12 @@ impl EntryLocalize { } else if !index.has_contig(&next.chrom) { Ok((acc, no_size, no_index + 1)) } else { - let contig_length = contig_length.unwrap(); - let mp = next.midpoint(); - let start = mp - .checked_sub(self.expand_window + 1) - .unwrap_or(0u64); // safe because of conditional above - let end = std::cmp::min( - mp.saturating_add(self.expand_window), - contig_length, - ); - next.start = start; - next.end = end; - acc.push(next); + acc.push(FocusRegion::from_feature( + next, + self.expand_window, + contig_length.unwrap(), + )); Ok((acc, no_size, no_index)) } }, @@ -380,8 +375,13 @@ impl EntryLocalize { .into_par_iter() .progress_with(successes) .map(|gr| { - let region = - format!("{}:{}-{}", gr.chrom, gr.start, gr.end); + let query_region = gr.query_region(); + let region = format!( + "{}:{}-{}", + query_region.chrom, + query_region.start, + query_region.end + ); gr.into_localized_mod_counts( &tabix_index, self.stranded_features, diff --git a/modkit-core/src/localise/util.rs b/modkit-core/src/localise/util.rs index 54fcffbe..84b54998 100644 --- a/modkit-core/src/localise/util.rs +++ b/modkit-core/src/localise/util.rs @@ -27,6 +27,32 @@ pub(super) struct LocalizedModCounts { offsets: FxHashMap>>, } +#[derive(Debug, Clone, PartialEq, Eq)] +pub(super) struct FocusRegion { + region: GenomeRegion, + anchor_point: u64, +} + +impl FocusRegion { + pub(super) fn from_feature( + mut region: GenomeRegion, + window: u64, + chrom_length: u64, + ) -> Self { + let anchor_point = region.midpoint(); + region.start = anchor_point.saturating_sub(window).min(chrom_length); + region.end = anchor_point + .saturating_add(window) + .saturating_add(1) + .min(chrom_length); + Self { region, anchor_point } + } + + pub(super) fn query_region(&self) -> &GenomeRegion { + &self.region + } +} + impl LocalizedModCounts { fn add_bedmethyl_record( &mut self, @@ -186,7 +212,7 @@ impl Moniod for LocalizedModCounts { } } -impl GenomeRegion { +impl FocusRegion { pub(super) fn into_localized_mod_counts( self, index: &HtsTabixHandler, @@ -194,19 +220,19 @@ impl GenomeRegion { stranded_features: Option, io_threads: usize, ) -> anyhow::Result { + let Self { region, anchor_point } = self; let bedmethyl_records = index.fetch_region( - &self.chrom, - &(self.start..self.end), - strand_rule.unwrap_or(self.strand), + ®ion.chrom, + &(region.start..region.end), + strand_rule.unwrap_or(region.strand), io_threads, )?; - let anchor_point = self.midpoint(); let loc_counts = bedmethyl_records .into_par_iter() .filter(|bm| { stranded_features .map(|f| { - let overlaps = self.strand.overlaps(&bm.strand); + let overlaps = region.strand.overlaps(&bm.strand); match f { StrandedFeatures::Same => overlaps, StrandedFeatures::Opposite => !overlaps, @@ -233,3 +259,44 @@ pub(super) enum StrandedFeatures { #[clap(name = "opposite")] Opposite, } + +#[cfg(test)] +mod tests { + use super::FocusRegion; + use crate::util::{GenomeRegion, StrandRule}; + + #[test] + fn focus_region_uses_inclusive_window_and_preserves_anchor() { + let cases = [ + (10, 11, 0, 100, 10, 11, 10), + (10, 14, 2, 100, 10, 15, 12), + (0, 1, 2, 100, 0, 3, 0), + (99, 100, 2, 100, 97, 100, 99), + ]; + + for ( + feature_start, + feature_end, + window, + chrom_length, + expected_start, + expected_end, + expected_anchor, + ) in cases + { + let feature = GenomeRegion::new( + "chr1".to_string(), + feature_start, + feature_end, + StrandRule::Both, + None, + ); + let focus = + FocusRegion::from_feature(feature, window, chrom_length); + + assert_eq!(focus.region.start, expected_start); + assert_eq!(focus.region.end, expected_end); + assert_eq!(focus.anchor_point, expected_anchor); + } + } +} diff --git a/modkit/tests/test_localize.rs b/modkit/tests/test_localize.rs index bb7c4a0c..998de384 100644 --- a/modkit/tests/test_localize.rs +++ b/modkit/tests/test_localize.rs @@ -236,3 +236,35 @@ m\t-1\t28\t10\t35.714287\n"; assert_eq!(std::fs::read(one_thread).unwrap(), expected); assert_eq!(std::fs::read(many_threads).unwrap(), expected); } + +#[test] +fn test_localize_zero_window_includes_anchor() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("regions.bed"); + let genome_sizes = temp_dir.path().join("genome-sizes.tsv"); + let output = temp_dir.path().join("localize.tsv"); + + write(®ions, "chr20\t9681998\t9681999\n").unwrap(); + write(&genome_sizes, "chr20\t100000000\n").unwrap(); + + run_modkit(&[ + "localize", + "../tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz", + "--regions", + regions.to_str().unwrap(), + "--genome-sizes", + genome_sizes.to_str().unwrap(), + "--window", + "0", + "--min-coverage", + "1", + "--out-file", + output.to_str().unwrap(), + ]) + .expect("failed to run modkit localize"); + + assert_eq!( + read_to_string(output).unwrap(), + "mod_code\toffset\tn_valid\tn_mod\tpercent_modified\nC\t0\t1\t1\t100\n" + ); +} From d9c03c9f928c5579c10dd00e49dda4d9c8668e54 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 17:25:18 -0700 Subject: [PATCH 5/9] test: define localize reference-axis offsets --- modkit/tests/test_localize.rs | 44 +++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/modkit/tests/test_localize.rs b/modkit/tests/test_localize.rs index 998de384..2cf2af15 100644 --- a/modkit/tests/test_localize.rs +++ b/modkit/tests/test_localize.rs @@ -268,3 +268,47 @@ fn test_localize_zero_window_includes_anchor() { "mod_code\toffset\tn_valid\tn_mod\tpercent_modified\nC\t0\t1\t1\t100\n" ); } + +#[test] +fn test_localize_offsets_use_reference_axis_for_all_feature_strands() { + let temp_dir = tempfile::tempdir().unwrap(); + let genome_sizes = temp_dir.path().join("genome-sizes.tsv"); + write(&genome_sizes, "chr20\t100000000\n").unwrap(); + + let expected = concat!( + "mod_code\toffset\tn_valid\tn_mod\tpercent_modified\n", + "C\t-7\t1\t1\t100\n", + "C\t8\t1\t1\t100\n", + ); + let feature_strands = [ + ("positive", "chr20\t9682005\t9682006\tfeature\t0\t+\n"), + ("negative", "chr20\t9682005\t9682006\tfeature\t0\t-\n"), + ("both", "chr20\t9682005\t9682006\n"), + ]; + + for (label, region_line) in feature_strands { + let regions = temp_dir.path().join(format!("regions-{label}.bed")); + let output = temp_dir.path().join(format!("localize-{label}.tsv")); + write(®ions, region_line).unwrap(); + + run_modkit(&[ + "localize", + "../tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz", + "--regions", + regions.to_str().unwrap(), + "--genome-sizes", + genome_sizes.to_str().unwrap(), + "--window", + "8", + "--stranded-features", + "both", + "--min-coverage", + "1", + "--out-file", + output.to_str().unwrap(), + ]) + .expect("failed to run modkit localize"); + + assert_eq!(read_to_string(output).unwrap(), expected, "{label}"); + } +} From fa09c007638c8d24f483a0b9ad3657214afaca66 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 17:27:19 -0700 Subject: [PATCH 6/9] fix: use reference-axis localize offsets --- book/src/advanced_usage.md | 10 +++++---- book/src/intro_localize.md | 18 ++++++++++------- modkit-core/src/localise/subcommand.rs | 28 +++++++++++++++++--------- modkit-core/src/localise/util.rs | 2 +- modkit/tests/test_localize.rs | 2 +- 5 files changed, 37 insertions(+), 23 deletions(-) diff --git a/book/src/advanced_usage.md b/book/src/advanced_usage.md index 6745d3ae..e648307f 100644 --- a/book/src/advanced_usage.md +++ b/book/src/advanced_usage.md @@ -1493,10 +1493,12 @@ Options: centered -w, --window - Number of base pairs to search around, for example if your BED region - records are single positions, a window of 500 will look 500 base pairs - upstream and downstream of that position. If your region BED records - are larger regions, this will expand from the midpoint of that region + Number of base pairs to search around. The original BED feature + midpoint remains offset zero. Output offsets are the bedMethyl + position minus that midpoint: negative is toward lower reference + coordinates and positive toward higher reference coordinates, + independent of feature strand. Earlier releases used the opposite + sign. For example, a window of 500 searches 500 bases on each side [default: 2000] diff --git a/book/src/intro_localize.md b/book/src/intro_localize.md index ecbe42cc..df46a1ac 100644 --- a/book/src/intro_localize.md +++ b/book/src/intro_localize.md @@ -20,12 +20,16 @@ modkit localise ${bedmethyl} --regions ${ctcf} --genome-sizes ${sizes} The output table has the following schema: -| column | Name | Description | type | -|--------|------------------|---------------------------------------------------------------------------------------------------------------------|-------| -| 1 | mod code | modification code as present in the bedmethyl | str | -| 2 | offset | distance in base pairs from the center of the genome features, negative values reflect towards the 5' of the genome | int | -| 3 | n_valid | number of valid calls at this offset for this modification code | int | -| 4 | n_mod | number of calls for this modification code at this offset | int | -| 5 | percent_modified | `n_mod` / `n_valid` * 100 | float | +| column | Name | Description | type | +|--------|------------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|-------| +| 1 | mod code | modification code as present in the bedmethyl | str | +| 2 | offset | bedMethyl position minus the original feature midpoint; negative is toward lower reference coordinates and positive toward higher coordinates, independent of the BED feature strand | int | +| 3 | n_valid | number of valid calls at this offset for this modification code | int | +| 4 | n_mod | number of calls for this modification code at this offset | int | +| 5 | percent_modified | `n_mod` / `n_valid` * 100 | float | + +Earlier modkit releases reported the opposite offset sign. Negate offsets from +those releases before comparing them with new output. The reference-coordinate +axis is not reversed for negative-strand features. Optionally the `--chart` argument can be used to create HTML charts of the modification patterns. diff --git a/modkit-core/src/localise/subcommand.rs b/modkit-core/src/localise/subcommand.rs index 1721db69..d536efda 100644 --- a/modkit-core/src/localise/subcommand.rs +++ b/modkit-core/src/localise/subcommand.rs @@ -44,11 +44,12 @@ pub struct EntryLocalize { #[clap(help_heading = "Output Options")] #[arg(long = "name", requires = "chart_filepath")] chart_name: Option, - /// Number of base pairs to search around, for example if your BED region - /// records are single positions, a window of 500 will look 500 base - /// pairs upstream and downstream of that position. If your region BED - /// records are larger regions, this will expand from the midpoint of - /// that region. + /// Number of base pairs to search around. The original BED feature + /// midpoint remains offset zero. Output offsets are the bedMethyl + /// position minus that midpoint: negative is toward lower reference + /// coordinates and positive toward higher reference coordinates, + /// independent of feature strand. Earlier releases used the opposite + /// sign. For example, a window of 500 searches 500 bases on each side. #[arg(short = 'w', long = "window", default_value_t = 2000)] expand_window: u64, // todo @@ -477,11 +478,18 @@ mod tests { let index = HtsTabixHandler::::from_path(&entry.in_bedmethyl)?; - entry.load_focus_regions( - &sequence_lengths, - &index, - &MultiProgress::new(), - ) + entry + .load_focus_regions( + &sequence_lengths, + &index, + &MultiProgress::new(), + ) + .map(|regions| { + regions + .into_iter() + .map(|region| region.query_region().clone()) + .collect() + }) } fn load_regions(regions: PathBuf) -> anyhow::Result> { diff --git a/modkit-core/src/localise/util.rs b/modkit-core/src/localise/util.rs index 84b54998..9e0973ed 100644 --- a/modkit-core/src/localise/util.rs +++ b/modkit-core/src/localise/util.rs @@ -60,7 +60,7 @@ impl LocalizedModCounts { anchor_point: u64, ) { let pos = bed_methyl_line.start() as i64; - let offset = (anchor_point as i64).saturating_sub(pos); + let offset = pos.saturating_sub(anchor_point as i64); let mod_pos_info = self .offsets .entry(bed_methyl_line.raw_mod_code) diff --git a/modkit/tests/test_localize.rs b/modkit/tests/test_localize.rs index 2cf2af15..9c859792 100644 --- a/modkit/tests/test_localize.rs +++ b/modkit/tests/test_localize.rs @@ -232,7 +232,7 @@ fn localize_valid_output_is_byte_identical_across_thread_counts() { } let expected = b"mod_code\toffset\tn_valid\tn_mod\tpercent_modified\n\ -m\t-1\t28\t10\t35.714287\n"; +m\t0\t28\t10\t35.714287\n"; assert_eq!(std::fs::read(one_thread).unwrap(), expected); assert_eq!(std::fs::read(many_threads).unwrap(), expected); } From d83f09a8467678a4bcfa67de3d6f6f8f7fdf3609 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Thu, 6 Aug 2026 03:56:51 -0700 Subject: [PATCH 7/9] docs: clarify localize offset migration --- book/src/advanced_usage.md | 6 ++++-- book/src/intro_localize.md | 6 +++--- modkit-core/src/localise/subcommand.rs | 6 ++++-- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/book/src/advanced_usage.md b/book/src/advanced_usage.md index e648307f..d2e476a0 100644 --- a/book/src/advanced_usage.md +++ b/book/src/advanced_usage.md @@ -1497,8 +1497,10 @@ Options: midpoint remains offset zero. Output offsets are the bedMethyl position minus that midpoint: negative is toward lower reference coordinates and positive toward higher reference coordinates, - independent of feature strand. Earlier releases used the opposite - sign. For example, a window of 500 searches 500 bases on each side + independent of feature strand. Profiles produced by affected earlier + releases must be regenerated because both queried membership and + offset bins can differ. For example, a window of 500 searches 500 bases + on each side [default: 2000] diff --git a/book/src/intro_localize.md b/book/src/intro_localize.md index df46a1ac..1d7b951e 100644 --- a/book/src/intro_localize.md +++ b/book/src/intro_localize.md @@ -28,8 +28,8 @@ The output table has the following schema: | 4 | n_mod | number of calls for this modification code at this offset | int | | 5 | percent_modified | `n_mod` / `n_valid` * 100 | float | -Earlier modkit releases reported the opposite offset sign. Negate offsets from -those releases before comparing them with new output. The reference-coordinate -axis is not reversed for negative-strand features. +Profiles produced by affected earlier modkit releases must be regenerated before +comparison because both queried membership and offset bins can differ. The +reference-coordinate axis is not reversed for negative-strand features. Optionally the `--chart` argument can be used to create HTML charts of the modification patterns. diff --git a/modkit-core/src/localise/subcommand.rs b/modkit-core/src/localise/subcommand.rs index d536efda..6ed81fdb 100644 --- a/modkit-core/src/localise/subcommand.rs +++ b/modkit-core/src/localise/subcommand.rs @@ -48,8 +48,10 @@ pub struct EntryLocalize { /// midpoint remains offset zero. Output offsets are the bedMethyl /// position minus that midpoint: negative is toward lower reference /// coordinates and positive toward higher reference coordinates, - /// independent of feature strand. Earlier releases used the opposite - /// sign. For example, a window of 500 searches 500 bases on each side. + /// independent of feature strand. Profiles produced by affected earlier + /// releases must be regenerated because both queried membership and offset + /// bins can differ. For example, a window of 500 searches 500 bases on each + /// side. #[arg(short = 'w', long = "window", default_value_t = 2000)] expand_window: u64, // todo From 2d68c7e1881705456b312ffb77a5c32bdd27fa8b Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 03:14:33 -0700 Subject: [PATCH 8/9] Fix localize chart rendering --- modkit-core/src/localise/subcommand.rs | 2 +- modkit-core/src/localise/util.rs | 175 ++++++++++++++++++++++--- modkit/tests/test_localize.rs | 55 ++++++++ 3 files changed, 213 insertions(+), 19 deletions(-) diff --git a/modkit-core/src/localise/subcommand.rs b/modkit-core/src/localise/subcommand.rs index 6ed81fdb..c55aabdc 100644 --- a/modkit-core/src/localise/subcommand.rs +++ b/modkit-core/src/localise/subcommand.rs @@ -420,13 +420,13 @@ impl EntryLocalize { table.to_csv_writer(writer)?; if let Some(p) = self.chart_filepath.as_ref() { + let blob = counts.get_plot(self.chart_name.as_ref())?; let fh = if self.force { File::create(p)? } else { File::create_new(p)? }; let mut writer = BufWriter::new(fh); - let blob = counts.get_plot(self.chart_name.as_ref())?; writer.write(blob.as_bytes())?; } diff --git a/modkit-core/src/localise/util.rs b/modkit-core/src/localise/util.rs index 9e0973ed..3823d7f8 100644 --- a/modkit-core/src/localise/util.rs +++ b/modkit-core/src/localise/util.rs @@ -118,10 +118,27 @@ impl LocalizedModCounts { .flat_map(|counts| counts.keys().copied()) .sorted() .collect::>(); - let (left, right) = match xs.iter().minmax() { - MinMaxResult::MinMax(x, y) => (*x, *y), - _ => bail!("should be at least one offset"), + let (left, right, has_single_offset) = match xs.iter().minmax() { + MinMaxResult::NoElements => { + bail!("cannot create localize chart: no offsets available") + } + MinMaxResult::OneElement(x) => { + let x = *x; + (x.saturating_sub(1), x.saturating_add(1), true) + } + MinMaxResult::MinMax(x, y) if x == y => { + let x = *x; + (x.saturating_sub(1), x.saturating_add(1), true) + } + MinMaxResult::MinMax(x, y) => (*x, *y, false), }; + let x_axis = Axis::new() + .type_(AxisType::Value) + .min(left) + .max(right) + .name("offset"); + let x_axis = + if has_single_offset { x_axis.min_interval(1) } else { x_axis }; let mut chart = Chart::new() .data_zoom( DataZoom::new() @@ -143,13 +160,7 @@ impl LocalizedModCounts { .save_as_image(SaveAsImage::new()), ), ) - .x_axis( - Axis::new() - .type_(AxisType::Value) - .min(left) - .max(right) - .name("offset"), - ) + .x_axis(x_axis) .y_axis( Axis::new().type_(AxisType::Value).name("percent modified"), ); @@ -161,18 +172,22 @@ impl LocalizedModCounts { DataPoint::Value(CompositeValue::Array(vec![ CompositeValue::Number(NumericValue::Integer(*offset)), CompositeValue::Number(NumericValue::Float( - info.frac_modified() as f64, + info.percent_modified() as f64, )), ])) }) .collect::>(); - chart = chart.series( - Line::new() - .name(format!("{mod_code}")) - .data(dat) - .symbol(Symbol::None) - .line_style(LineStyle::new().width(1.5)), - ); + let is_singleton = dat.len() == 1; + let line = Line::new() + .name(format!("{mod_code}")) + .data(dat) + .line_style(LineStyle::new().width(1.5)); + let line = if is_singleton { + line.symbol(Symbol::Circle).show_symbol(true) + } else { + line.symbol(Symbol::None) + }; + chart = chart.series(line); } HtmlRenderer::new(chart_name.unwrap_or(&default_name), 800, 800) @@ -299,4 +314,128 @@ mod tests { assert_eq!(focus.anchor_point, expected_anchor); } } + use serde_json::Value; + + use super::LocalizedModCounts; + use crate::mod_base_code::{ + ModCodeRepr, HYDROXY_METHYL_CYTOSINE, METHYL_CYTOSINE, + }; + use crate::util::ModPositionInfo; + + fn counts( + series: impl IntoIterator)>, + ) -> LocalizedModCounts { + let offsets = series + .into_iter() + .map(|(code, values)| { + let values = values + .into_iter() + .map(|(offset, n_valid, n_mod)| { + (offset, ModPositionInfo::new(n_valid, n_mod)) + }) + .collect(); + (code, values) + }) + .collect(); + LocalizedModCounts { offsets } + } + + fn chart_json(counts: &LocalizedModCounts) -> Value { + let html = counts.get_plot(None).unwrap(); + let raw_chart = html + .split_once("var option = ") + .unwrap() + .1 + .split_once(";\n chart.setOption") + .unwrap() + .0; + serde_json::from_str(raw_chart).unwrap() + } + + fn series_named<'a>(chart: &'a Value, name: &str) -> &'a Value { + chart["series"] + .as_array() + .unwrap() + .iter() + .find(|series| series["name"] == name) + .unwrap() + } + + #[test] + fn plot_uses_percent_modified() { + let counts = counts([(METHYL_CYTOSINE, vec![(-1, 4, 1), (1, 4, 3)])]); + let chart = chart_json(&counts); + let data = chart["series"][0]["data"].as_array().unwrap(); + + assert_eq!(chart["yAxis"][0]["name"], "percent modified"); + assert_eq!(data[0][1].as_f64(), Some(25.0)); + assert_eq!(data[1][1].as_f64(), Some(75.0)); + } + + #[test] + fn plot_one_point_has_padded_axis_and_visible_symbol() { + let counts = counts([(METHYL_CYTOSINE, vec![(0, 4, 1)])]); + let chart = chart_json(&counts); + let x_axis = &chart["xAxis"][0]; + let series = series_named(&chart, "m"); + + assert_eq!(x_axis["min"].as_i64(), Some(-1)); + assert_eq!(x_axis["max"].as_i64(), Some(1)); + assert_eq!(x_axis["minInterval"].as_f64(), Some(1.0)); + assert_eq!(series["symbol"], "circle"); + assert_eq!(series["showSymbol"], true); + assert_eq!(series["data"].as_array().unwrap().len(), 1); + } + + #[test] + fn plot_multiple_singleton_codes_share_padded_axis() { + let counts = counts([ + (METHYL_CYTOSINE, vec![(7, 4, 1)]), + (HYDROXY_METHYL_CYTOSINE, vec![(7, 5, 2)]), + ]); + let chart = chart_json(&counts); + let x_axis = &chart["xAxis"][0]; + + assert_eq!(x_axis["min"].as_i64(), Some(6)); + assert_eq!(x_axis["max"].as_i64(), Some(8)); + assert_eq!(x_axis["minInterval"].as_f64(), Some(1.0)); + for name in ["m", "h"] { + let series = series_named(&chart, name); + assert_eq!(series["symbol"], "circle"); + assert_eq!(series["showSymbol"], true); + } + } + + #[test] + fn plot_singleton_symbol_does_not_change_multi_point_series() { + let counts = counts([ + (METHYL_CYTOSINE, vec![(-2, 4, 1), (2, 4, 3)]), + (HYDROXY_METHYL_CYTOSINE, vec![(0, 5, 2)]), + ]); + let chart = chart_json(&counts); + let multi = series_named(&chart, "m"); + let singleton = series_named(&chart, "h"); + + assert_eq!(multi["symbol"], "none"); + assert!(multi.get("showSymbol").is_none()); + assert_eq!(singleton["symbol"], "circle"); + assert_eq!(singleton["showSymbol"], true); + } + + #[test] + fn plot_preserves_multi_offset_bounds() { + let counts = counts([(METHYL_CYTOSINE, vec![(-5, 4, 1), (9, 4, 3)])]); + let chart = chart_json(&counts); + let x_axis = &chart["xAxis"][0]; + + assert_eq!(x_axis["min"].as_i64(), Some(-5)); + assert_eq!(x_axis["max"].as_i64(), Some(9)); + assert!(x_axis.get("minInterval").is_none()); + } + + #[test] + fn plot_empty_counts_returns_clear_error() { + let error = LocalizedModCounts::default().get_plot(None).unwrap_err(); + assert!(error.to_string().contains("no offsets")); + } } diff --git a/modkit/tests/test_localize.rs b/modkit/tests/test_localize.rs index 9c859792..1c355937 100644 --- a/modkit/tests/test_localize.rs +++ b/modkit/tests/test_localize.rs @@ -312,3 +312,58 @@ fn test_localize_offsets_use_reference_axis_for_all_feature_strands() { assert_eq!(read_to_string(output).unwrap(), expected, "{label}"); } } + +#[test] +fn test_localize_empty_chart_does_not_create_or_truncate_file() { + const EMPTY_TABLE: &str = + "mod_code\toffset\tn_valid\tn_mod\tpercent_modified\n"; + + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("regions.bed"); + let genome_sizes = temp_dir.path().join("genome-sizes.tsv"); + let table = temp_dir.path().join("localize.tsv"); + let new_chart = temp_dir.path().join("new-localize.html"); + let existing_chart = temp_dir.path().join("existing-localize.html"); + + write(®ions, "chr20\t1\t2\n").unwrap(); + write(&genome_sizes, "chr20\t100000000\n").unwrap(); + + let run = |chart: &std::path::Path, force: bool| { + let mut args = vec![ + "localize", + "../tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz", + "--regions", + regions.to_str().unwrap(), + "--genome-sizes", + genome_sizes.to_str().unwrap(), + "--window", + "1", + "--out-file", + table.to_str().unwrap(), + "--chart", + chart.to_str().unwrap(), + ]; + if force { + args.push("--force"); + } + run_modkit(&args) + }; + + let create_result = run(&new_chart, false); + assert!(create_result.is_err()); + assert_eq!(read_to_string(&table).unwrap(), EMPTY_TABLE); + + let sentinel = "keep existing chart content\n"; + write(&existing_chart, sentinel).unwrap(); + let truncate_result = run(&existing_chart, true); + + assert!(truncate_result.is_err()); + assert_eq!(read_to_string(&table).unwrap(), EMPTY_TABLE); + let existing_content = read_to_string(existing_chart).unwrap(); + assert!( + !new_chart.exists() && existing_content == sentinel, + "failed chart rendering mutated chart paths: new_exists={}, \ + existing_content={existing_content:?}", + new_chart.exists() + ); +} From a290d201e189c6c4ea375fddea2036bd7601b88c Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Tue, 4 Aug 2026 23:08:12 -0700 Subject: [PATCH 9/9] Apply localize minimum coverage before aggregation --- modkit-core/src/localise/subcommand.rs | 1 + modkit-core/src/localise/util.rs | 2 ++ modkit/tests/test_localize.rs | 46 ++++++++++++++++++++++++++ 3 files changed, 49 insertions(+) diff --git a/modkit-core/src/localise/subcommand.rs b/modkit-core/src/localise/subcommand.rs index c55aabdc..82039282 100644 --- a/modkit-core/src/localise/subcommand.rs +++ b/modkit-core/src/localise/subcommand.rs @@ -387,6 +387,7 @@ impl EntryLocalize { ); gr.into_localized_mod_counts( &tabix_index, + min_cov, self.stranded_features, stranded_features, self.io_threads, diff --git a/modkit-core/src/localise/util.rs b/modkit-core/src/localise/util.rs index 3823d7f8..40981284 100644 --- a/modkit-core/src/localise/util.rs +++ b/modkit-core/src/localise/util.rs @@ -231,6 +231,7 @@ impl FocusRegion { pub(super) fn into_localized_mod_counts( self, index: &HtsTabixHandler, + min_coverage: u64, strand_rule: Option, stranded_features: Option, io_threads: usize, @@ -244,6 +245,7 @@ impl FocusRegion { )?; let loc_counts = bedmethyl_records .into_par_iter() + .filter(|bm| bm.valid_coverage >= min_coverage) .filter(|bm| { stranded_features .map(|f| { diff --git a/modkit/tests/test_localize.rs b/modkit/tests/test_localize.rs index 1c355937..ed138670 100644 --- a/modkit/tests/test_localize.rs +++ b/modkit/tests/test_localize.rs @@ -367,3 +367,49 @@ fn test_localize_empty_chart_does_not_create_or_truncate_file() { new_chart.exists() ); } + +#[test] +fn test_localize_min_coverage_filters_before_aggregation() { + let temp_dir = tempfile::tempdir().unwrap(); + let regions = temp_dir.path().join("regions.bed"); + let genome_sizes = temp_dir.path().join("genome-sizes.tsv"); + write(®ions, "chr20\t9681998\t9681999\nchr20\t9838537\t9838538\n") + .unwrap(); + write(&genome_sizes, "chr20\t64444167\n").unwrap(); + + let run = |min_coverage: u64| { + let output = temp_dir.path().join(format!("min-{min_coverage}.tsv")); + let min_coverage = min_coverage.to_string(); + run_modkit(&[ + "localize", + "../tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz", + "--regions", + regions.to_str().unwrap(), + "--genome-sizes", + genome_sizes.to_str().unwrap(), + "--window", + "0", + "--min-coverage", + &min_coverage, + "--threads", + "1", + "--io-threads", + "1", + "--out-file", + output.to_str().unwrap(), + ]) + .unwrap(); + read_to_string(output).unwrap().replace("\r\n", "\n") + }; + + assert_eq!( + run(1), + "mod_code\toffset\tn_valid\tn_mod\tpercent_modified\n\ + C\t0\t24\t3\t12.5\n" + ); + assert_eq!( + run(3), + "mod_code\toffset\tn_valid\tn_mod\tpercent_modified\n\ + C\t0\t23\t2\t8.695652\n" + ); +}