diff --git a/book/src/intro_validate.md b/book/src/intro_validate.md index b8eea706..9d93e59b 100644 --- a/book/src/intro_validate.md +++ b/book/src/intro_validate.md @@ -76,10 +76,17 @@ 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`. +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. 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..19335e0b 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -1,13 +1,15 @@ 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; 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 +82,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" + )), } } } @@ -122,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 { @@ -136,13 +150,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 { @@ -155,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 }) } @@ -165,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()) @@ -185,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 @@ -192,28 +226,91 @@ 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 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() { - 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 +995,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 +1056,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 +1227,239 @@ 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_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(); + 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..7496c7af 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,61 @@ 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"), + ( + "conflicting-label", + "chr1\t4\t6\th\t.\t+", + "conflicting ground truth labels", + ), + ]; + + 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}" + ); + } +}