From ac05e20ed68fd263ea1041046b75e29d630e3b7e Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 05:04:03 -0700 Subject: [PATCH 1/2] Return validate truth BED errors --- book/src/intro_validate.md | 6 +- modkit-core/src/validate/subcommand.rs | 259 ++++++++++++++++++++++--- modkit/tests/test_validate.rs | 55 ++++++ 3 files changed, 287 insertions(+), 33 deletions(-) diff --git a/book/src/intro_validate.md b/book/src/intro_validate.md index b8eea706..a37af135 100644 --- a/book/src/intro_validate.md +++ b/book/src/intro_validate.md @@ -76,10 +76,14 @@ These fields are as follows: | 2 | start position | 0-based start position | int | | 3 | end position | 0-based exclusive end position | int | | 4 | mod code | modified base code | str | -| 6 | strand | strand (e.g. +,-,.) | str | +| 6 | strand | strand (`+` or `-`) | str | The 5th column is ignored in the validate command. +The start position must be non-negative and the end position must be greater +than the start position. Zero-length BED features do not annotate a reference +base and are rejected by `validate`. + The 4th column represents the modified base code annotating the status at this reference position (or range of reference positions). This value can be `-` representing a canonical base (note that this differs from the `remora validate` annotation), a single letter code as defined in the modBAM tag specification, or any ChEBI code. The validate command will assume that any base from the associated modBAM file overlapping these positions should match this annotation. diff --git a/modkit-core/src/validate/subcommand.rs b/modkit-core/src/validate/subcommand.rs index 119aa3ed..13145f95 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -7,7 +7,7 @@ use std::path::PathBuf; use std::string::FromUtf8Error; use ansi_term::Style; -use anyhow::{anyhow, bail}; +use anyhow::{anyhow, bail, Context}; use clap::Args; use derive_new::new; use itertools::Itertools; @@ -80,18 +80,30 @@ impl BaseStatus { impl BaseStatus { pub fn parse(raw: &str) -> anyhow::Result { - if let Ok(code) = raw.parse::() { - if code == '-' { - Ok(Self::Canonical) - } else { + if raw == "-" { + return Ok(Self::Canonical); + } + + if !raw.is_empty() && raw.bytes().all(|b| b.is_ascii_digit()) { + let chebi = raw.parse::().map_err(|e| { + anyhow!("failed to parse ChEBI code {raw}: {e}") + })?; + return Ok(Self::Modified(ModCodeRepr::ChEbi(chebi))); + } + + let mut chars = raw.chars(); + match (chars.next(), chars.next()) { + (Some(code), None) + if code.is_ascii_lowercase() + || matches!(code, 'A' | 'C' | 'G' | 'T' | 'U' | 'N') => + { Ok(Self::Modified(ModCodeRepr::Code(code))) } - } else { - if let Ok(chebi) = raw.parse::() { - Ok(Self::Modified(ModCodeRepr::ChEbi(chebi))) - } else { - Err(anyhow!("failed to parse mod code {raw}")) - } + _ => Err(anyhow!( + "failed to parse mod code {raw}: expected `-`, a numeric \ + ChEBI code, one lowercase ASCII letter, or one of \ + A/C/G/T/U/N" + )), } } } @@ -136,13 +148,21 @@ fn parse_ground_truth_bed_line(line: &str) -> anyhow::Result { fields[1].parse().map_err(|e| anyhow!("Error parsing start: {}", e))?; let end: i64 = fields[2].parse().map_err(|e| anyhow!("Error parsing end: {}", e))?; + if start < 0 { + bail!("BED start must be non-negative, found {start}"); + } + if end <= start { + bail!( + "BED end must be greater than start, found start {start} and end \ + {end}" + ); + } let raw_mod_code = fields[3]; - let strand_char = fields[5] - .chars() - .next() - .ok_or_else(|| anyhow!("Error parsing strand {}", fields[5]))?; - - let strand = Strand::parse_char(strand_char)?; + let strand = match fields[5] { + "+" => Strand::Positive, + "-" => Strand::Negative, + raw => bail!("Error parsing strand {raw}: expected `+` or `-`"), + }; let base_status = BaseStatus::parse(&raw_mod_code) .map_err(|e| anyhow!("Error parsing base status code: {}", e))?; if let BaseStatus::Modified(mod_code) = base_status { @@ -192,16 +212,29 @@ fn parse_ground_truth_bed_file( } lines_processed.set_message("rows processed"); - let reader = BufReader::new(File::open(file_path)?); - for ground_truth_site in reader - .lines() - .skip_while(|r| r.as_ref().map(|l| l.starts_with('#')).unwrap_or(false)) - .filter_map(|r| { - r.map_err(|e| anyhow!("failed to read, {}", e.to_string())) - .and_then(|line| parse_ground_truth_bed_line(&line)) - .ok() - }) - { + let reader = BufReader::new(File::open(file_path).with_context(|| { + format!("failed to open ground truth BED {}", file_path.display()) + })?); + for (line_idx, raw_line) in reader.lines().enumerate() { + let line_number = line_idx + 1; + let line = raw_line.with_context(|| { + format!( + "failed to read ground truth BED {} at line {line_number}", + file_path.display() + ) + })?; + let line = line.trim(); + if line.is_empty() || line.starts_with('#') { + continue; + } + let ground_truth_site = parse_ground_truth_bed_line(line) + .with_context(|| { + format!( + "failed to parse ground truth BED {} at line \ + {line_number}", + file_path.display() + ) + })?; let cs_res = result .entry(ground_truth_site.chrom) .or_insert_with(HashMap::new) @@ -213,7 +246,10 @@ fn parse_ground_truth_bed_file( lines_processed.inc(1); } if result.is_empty() { - bail!("zero valid positions parsed from BED file".to_string()); + bail!( + "zero valid positions parsed from BED file {}", + file_path.display() + ); } lines_processed.finish_and_clear(); info!("Processed {} BED lines", lines_processed.position()); @@ -898,10 +934,6 @@ pub struct ValidateFromModBam { impl ValidateFromModBam { pub fn run(&self) -> anyhow::Result<()> { let _handle = init_logging(self.log_filepath.as_ref()); - let mut out_handle: Option = None; - if let Some(file_path) = self.out_filepath.clone() { - out_handle = Some(File::create(&file_path)?); - } let collapse_method = match &self.ignore { Some(raw_mod_code) => { let mod_code = ModCodeRepr::parse(raw_mod_code)?; @@ -963,6 +995,8 @@ impl ValidateFromModBam { let can_base = derive_canonical_base(>_positions, self.canonical_base)?; info!("Canonical base: {}", can_base); + let mut out_handle = + self.out_filepath.as_ref().map(File::create).transpose()?; let mut all_probs = HashMap::new(); for (bam_path, bed_indices) in bam_path_to_bed_indices { @@ -1132,3 +1166,164 @@ impl ValidateFromModBam { Ok(()) } } + +#[cfg(test)] +mod tests { + use std::io::Write; + + use tempfile::NamedTempFile; + + use super::{parse_ground_truth_bed_file, BaseStatus, ModCodeRepr, Strand}; + + #[test] + fn ground_truth_bed_reports_mixed_row_error_with_physical_line() { + let mut bed = NamedTempFile::new().unwrap(); + writeln!(bed, "chr1\t4\t5\tm\t.\t+").unwrap(); + writeln!(bed, "not-a-bed-row").unwrap(); + + let error = + parse_ground_truth_bed_file(&bed.path().to_path_buf(), true) + .unwrap_err(); + let message = format!("{error:#}"); + assert!( + message.contains(&bed.path().display().to_string()), + "{message}" + ); + assert!(message.contains("line 2"), "{message}"); + assert!(message.contains("Invalid number of fields"), "{message}"); + } + + #[test] + fn ground_truth_bed_allows_blank_and_comment_rows_anywhere() { + let mut bed = NamedTempFile::new().unwrap(); + writeln!(bed, "# header").unwrap(); + writeln!(bed, "chr1\t4\t5\tm\t.\t+").unwrap(); + writeln!(bed).unwrap(); + writeln!(bed, " # another comment").unwrap(); + writeln!(bed, "chr1\t8\t9\t-\t.\t-").unwrap(); + + let parsed = + parse_ground_truth_bed_file(&bed.path().to_path_buf(), true) + .unwrap(); + assert_eq!( + parsed["chr1"][&Strand::Positive][&4], + BaseStatus::Modified('m'.into()) + ); + assert_eq!( + parsed["chr1"][&Strand::Negative][&8], + BaseStatus::Canonical + ); + } + + #[test] + fn ground_truth_bed_rejects_invalid_truth_fields_with_physical_line() { + let cases = [ + ( + "negative start", + "chr1\t-1\t1\tm\t.\t+", + "BED start must be non-negative", + ), + ( + "reversed interval", + "chr1\t5\t4\tm\t.\t+", + "BED end must be greater than start", + ), + ( + "zero-width interval", + "chr1\t5\t5\tm\t.\t+", + "BED end must be greater than start", + ), + ("strand suffix", "chr1\t5\t6\tm\t.\t+junk", "expected `+` or `-`"), + ( + "punctuation mod code", + "chr1\t5\t6\t.\t.\t+", + "failed to parse mod code", + ), + ( + "multi-letter mod code", + "chr1\t5\t6\tmm\t.\t+", + "failed to parse mod code", + ), + ( + "invalid uppercase mod code", + "chr1\t5\t6\tM\t.\t+", + "failed to parse mod code", + ), + ( + "non-ASCII mod code", + "chr1\t5\t6\té\t.\t+", + "failed to parse mod code", + ), + ]; + + for (case_name, invalid_line, expected_error) in cases { + let mut bed = NamedTempFile::new().unwrap(); + writeln!(bed, "# header").unwrap(); + writeln!(bed, "chr1\t4\t5\tm\t.\t+").unwrap(); + writeln!(bed).unwrap(); + writeln!(bed, "{invalid_line}").unwrap(); + + let error = + parse_ground_truth_bed_file(&bed.path().to_path_buf(), true) + .unwrap_err(); + let message = format!("{error:#}"); + assert!( + message.contains(&bed.path().display().to_string()), + "{case_name}: {message}" + ); + assert!(message.contains("line 4"), "{case_name}: {message}"); + assert!(message.contains(expected_error), "{case_name}: {message}"); + } + } + + #[test] + fn ground_truth_bed_accepts_documented_mod_tokens_and_extra_columns() { + let cases = [ + ("m", BaseStatus::Modified(ModCodeRepr::Code('m'))), + ("C", BaseStatus::Modified(ModCodeRepr::Code('C'))), + ("1", BaseStatus::Modified(ModCodeRepr::ChEbi(1))), + ("21839", BaseStatus::Modified(ModCodeRepr::ChEbi(21839))), + ("-", BaseStatus::Canonical), + ]; + let mut bed = NamedTempFile::new().unwrap(); + for (offset, (raw_code, _)) in cases.iter().enumerate() { + let start = 10 + offset as i64; + writeln!( + bed, + "chr1\t{start}\t{}\t{raw_code}\t.\t+\textra\tcolumns", + start + 1 + ) + .unwrap(); + } + + let parsed = + parse_ground_truth_bed_file(&bed.path().to_path_buf(), true) + .unwrap(); + for (offset, (raw_code, expected_status)) in cases.iter().enumerate() { + let position = 10 + offset as i64; + assert_eq!( + parsed["chr1"][&Strand::Positive][&position], + *expected_status, + "token {raw_code}" + ); + } + } + + #[test] + fn ground_truth_bed_reports_read_error_with_physical_line() { + let mut bed = NamedTempFile::new().unwrap(); + bed.write_all(b"chr1\t4\t5\tm\t.\t+\n").unwrap(); + bed.write_all(&[0xff, b'\n']).unwrap(); + + let error = + parse_ground_truth_bed_file(&bed.path().to_path_buf(), true) + .unwrap_err(); + let message = format!("{error:#}"); + assert!( + message.contains(&bed.path().display().to_string()), + "{message}" + ); + assert!(message.contains("line 2"), "{message}"); + assert!(message.contains("failed to read"), "{message}"); + } +} diff --git a/modkit/tests/test_validate.rs b/modkit/tests/test_validate.rs index aba45a9d..e404a059 100644 --- a/modkit/tests/test_validate.rs +++ b/modkit/tests/test_validate.rs @@ -2,6 +2,8 @@ use crate::common::run_modkit; use anyhow::Context; use std::fs::File; use std::io::{BufRead, BufReader}; +use std::path::Path; +use std::process::Command; mod common; @@ -52,3 +54,56 @@ fn test_validate_expected() { } } } + +#[test] +fn test_validate_bed_errors_preserve_output_and_report_line() { + let temp_dir = tempfile::tempdir().unwrap(); + let cases = [ + ("field-count", "not-a-bed-row", "Invalid number of fields"), + ( + "coordinate", + "chr1\t5\t5\tm\t.\t+", + "BED end must be greater than start", + ), + ("strand", "chr1\t5\t6\tm\t.\t+junk", "expected `+` or `-`"), + ("mod-code", "chr1\t5\t6\t.\t.\t+", "failed to parse mod code"), + ]; + + for (case_name, invalid_line, expected_error) in cases { + let bed_path = temp_dir.path().join(format!("{case_name}.bed")); + let output_path = temp_dir.path().join(format!("{case_name}.tsv")); + std::fs::write( + &bed_path, + format!("chr1\t4\t5\tm\t.\t+\n{invalid_line}\n"), + ) + .unwrap(); + std::fs::write(&output_path, "sentinel\n").unwrap(); + + let output = Command::new(Path::new(env!("CARGO_BIN_EXE_modkit"))) + .arg("validate") + .arg("--bam-and-bed") + .arg("../tests/resources/input_5mC.bam") + .arg(&bed_path) + .arg("--canonical-base") + .arg("C") + .arg("--out-filepath") + .arg(&output_path) + .arg("--suppress-progress") + .output() + .unwrap(); + let stderr = String::from_utf8_lossy(&output.stderr); + + assert!(!output.status.success(), "{case_name}: {stderr}"); + assert!( + stderr.contains(&bed_path.display().to_string()), + "{case_name}: {stderr}" + ); + assert!(stderr.contains("line 2"), "{case_name}: {stderr}"); + assert!(stderr.contains(expected_error), "{case_name}: {stderr}"); + assert_eq!( + std::fs::read(&output_path).unwrap(), + b"sentinel\n", + "{case_name}" + ); + } +} From 585a9791b711e2bd2fb18cd9ed08bdad131b3f14 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 05:27:58 -0700 Subject: [PATCH 2/2] Reject conflicting validate truth labels --- book/src/intro_validate.md | 3 + modkit-core/src/validate/subcommand.rs | 144 ++++++++++++++++++++++++- modkit/tests/test_validate.rs | 5 + 3 files changed, 148 insertions(+), 4 deletions(-) diff --git a/book/src/intro_validate.md b/book/src/intro_validate.md index a37af135..9d93e59b 100644 --- a/book/src/intro_validate.md +++ b/book/src/intro_validate.md @@ -83,6 +83,9 @@ The 5th column is ignored in the validate command. The start position must be non-negative and the end position must be greater than the start position. Zero-length BED features do not annotate a reference base and are rejected by `validate`. +Repeated or overlapping annotations may assign the same status more than once, +but conflicting statuses at the same reference position and strand are +rejected. The 4th column represents the modified base code annotating the status at this reference position (or range of reference positions). This value can be `-` representing a canonical base (note that this differs from the `remora validate` annotation), a single letter code as defined in the modBAM tag specification, or any ChEBI code. diff --git a/modkit-core/src/validate/subcommand.rs b/modkit-core/src/validate/subcommand.rs index 13145f95..19335e0b 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -1,8 +1,10 @@ use std::cmp::Ordering; +use std::collections::btree_map::Entry as BTreeMapEntry; use std::collections::{BTreeMap, HashMap, HashSet}; use std::fmt::{Debug, Display, Formatter}; use std::fs::File; use std::io::{BufRead, BufReader, Write}; +use std::ops::Range; use std::path::PathBuf; use std::string::FromUtf8Error; @@ -134,7 +136,7 @@ struct GroundTruthSite { pub chrom: String, pub strand: Strand, pub base_status: BaseStatus, - pub positions: Vec, + pub positions: Range, } fn parse_ground_truth_bed_line(line: &str) -> anyhow::Result { @@ -175,7 +177,7 @@ fn parse_ground_truth_bed_line(line: &str) -> anyhow::Result { ) } } - let positions = (start..end).collect(); + let positions = start..end; Ok(GroundTruthSite { chrom, strand, base_status, positions }) } @@ -185,6 +187,17 @@ type TidToChrom = HashMap; type ChromStrandPositionNames = HashMap>>; +#[derive(Debug, Copy, Clone)] +struct GroundTruthIntervalProvenance { + start: i64, + end: i64, + base_status: BaseStatus, + line_number: usize, +} + +type ChromStrandIntervalProvenance = + HashMap>>; + fn get_tid_to_chrom(reader: &Reader) -> anyhow::Result { let header = reader.header().to_owned(); let tid_to_chrom_result = (0..header.target_count()) @@ -205,6 +218,7 @@ fn parse_ground_truth_bed_file( ) -> anyhow::Result { info!("Parsing BED at {}", file_path.to_str().unwrap_or("invalid-UTF-8")); let mut result = HashMap::new(); + let mut provenance: ChromStrandIntervalProvenance = HashMap::new(); let lines_processed = get_ticker(); if suppress_pb { lines_processed @@ -235,14 +249,61 @@ fn parse_ground_truth_bed_file( file_path.display() ) })?; + let interval_start = ground_truth_site.positions.start; + let interval_end = ground_truth_site.positions.end; let cs_res = result - .entry(ground_truth_site.chrom) + .entry(ground_truth_site.chrom.clone()) .or_insert_with(HashMap::new) .entry(ground_truth_site.strand) .or_insert_with(BTreeMap::new); + let prior_intervals = provenance + .entry(ground_truth_site.chrom.clone()) + .or_insert_with(HashMap::new) + .entry(ground_truth_site.strand) + .or_insert_with(Vec::new); for pos in ground_truth_site.positions { - cs_res.insert(pos, ground_truth_site.base_status); + match cs_res.entry(pos) { + BTreeMapEntry::Vacant(entry) => { + entry.insert(ground_truth_site.base_status); + } + BTreeMapEntry::Occupied(entry) => { + let existing_status = *entry.get(); + if existing_status == ground_truth_site.base_status { + continue; + } + let first_line = prior_intervals + .iter() + .find(|assignment| { + assignment.start <= pos + && pos < assignment.end + && assignment.base_status == existing_status + }) + .map(|assignment| assignment.line_number) + .ok_or_else(|| { + anyhow!( + "missing provenance for existing ground truth \ + assignment" + ) + })?; + bail!( + "conflicting ground truth labels in BED {} at line \ + {line_number} for {}:{pos} strand {}: existing \ + `{existing_status}` from line {first_line}, new \ + `{}`", + file_path.display(), + ground_truth_site.chrom, + ground_truth_site.strand, + ground_truth_site.base_status + ); + } + } } + prior_intervals.push(GroundTruthIntervalProvenance { + start: interval_start, + end: interval_end, + base_status: ground_truth_site.base_status, + line_number, + }); lines_processed.inc(1); } if result.is_empty() { @@ -1309,6 +1370,81 @@ mod tests { } } + #[test] + fn ground_truth_bed_rejects_conflicting_overlap_in_both_orders() { + let cases = [("m", "h"), ("h", "m")]; + + for (existing_label, new_label) in cases { + let mut bed = NamedTempFile::new().unwrap(); + writeln!(bed, "chr1\t4\t7\t{existing_label}\t.\t+").unwrap(); + writeln!(bed, "chr1\t6\t8\t{new_label}\t.\t+").unwrap(); + + let error = + parse_ground_truth_bed_file(&bed.path().to_path_buf(), true) + .unwrap_err(); + let message = format!("{error:#}"); + assert!( + message.contains(&bed.path().display().to_string()), + "{message}" + ); + assert!(message.contains("line 2"), "{message}"); + assert!(message.contains("chr1:6 strand +"), "{message}"); + assert!( + message.contains(&format!( + "existing `{existing_label}` from line 1" + )), + "{message}" + ); + assert!( + message.contains(&format!("new `{new_label}`")), + "{message}" + ); + } + } + + #[test] + fn ground_truth_bed_retains_first_line_for_later_conflict() { + let mut bed = NamedTempFile::new().unwrap(); + writeln!(bed, "chr1\t4\t8\tm\t.\t+").unwrap(); + writeln!(bed, "chr1\t6\t9\tm\t.\t+").unwrap(); + writeln!(bed, "chr1\t7\t8\th\t.\t+").unwrap(); + + let error = + parse_ground_truth_bed_file(&bed.path().to_path_buf(), true) + .unwrap_err(); + let message = format!("{error:#}"); + assert!(message.contains("line 3"), "{message}"); + assert!(message.contains("chr1:7 strand +"), "{message}"); + assert!(message.contains("existing `m` from line 1"), "{message}"); + assert!(message.contains("new `h`"), "{message}"); + } + + #[test] + fn ground_truth_bed_accepts_identical_duplicate_and_overlap() { + let mut bed = NamedTempFile::new().unwrap(); + writeln!(bed, "chr1\t4\t8\tm\t.\t+").unwrap(); + writeln!(bed, "chr1\t4\t8\tm\t.\t+").unwrap(); + writeln!(bed, "chr1\t6\t9\tm\t.\t+").unwrap(); + writeln!(bed, "chr1\t6\t7\th\t.\t-").unwrap(); + + let parsed = + parse_ground_truth_bed_file(&bed.path().to_path_buf(), true) + .unwrap(); + let positive = &parsed["chr1"][&Strand::Positive]; + assert_eq!(positive.len(), 5); + assert_eq!( + positive.keys().copied().collect::>(), + vec![4, 5, 6, 7, 8] + ); + assert!(positive + .values() + .all(|status| *status == BaseStatus::Modified('m'.into()))); + assert_eq!( + parsed["chr1"][&Strand::Negative][&6], + BaseStatus::Modified('h'.into()) + ); + } + #[test] fn ground_truth_bed_reports_read_error_with_physical_line() { let mut bed = NamedTempFile::new().unwrap(); diff --git a/modkit/tests/test_validate.rs b/modkit/tests/test_validate.rs index e404a059..7496c7af 100644 --- a/modkit/tests/test_validate.rs +++ b/modkit/tests/test_validate.rs @@ -67,6 +67,11 @@ fn test_validate_bed_errors_preserve_output_and_report_line() { ), ("strand", "chr1\t5\t6\tm\t.\t+junk", "expected `+` or `-`"), ("mod-code", "chr1\t5\t6\t.\t.\t+", "failed to parse mod code"), + ( + "conflicting-label", + "chr1\t4\t6\th\t.\t+", + "conflicting ground truth labels", + ), ]; for (case_name, invalid_line, expected_error) in cases {