From d4409906f2c357c9dc707a1399c0c53cb64fa960 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Tue, 4 Aug 2026 23:27:50 -0700 Subject: [PATCH 01/11] Reopen validation BAM for each truth BED --- modkit-core/src/validate/subcommand.rs | 6 +- modkit/tests/test_validate.rs | 110 ++++++++++++++++++++++++- 2 files changed, 112 insertions(+), 4 deletions(-) diff --git a/modkit-core/src/validate/subcommand.rs b/modkit-core/src/validate/subcommand.rs index 119aa3ed..69e23caa 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -966,15 +966,15 @@ impl ValidateFromModBam { let mut all_probs = HashMap::new(); for (bam_path, bed_indices) in bam_path_to_bed_indices { - let mut reader = Reader::from_path(bam_path.as_path())?; - reader.set_threads(self.threads)?; - let tid_to_chrom = get_tid_to_chrom(&reader)?; info!( "Parsing mapping at {}", bam_path.to_str().unwrap_or("invalid-UTF-8") ); for bed_idx in bed_indices { + let mut reader = Reader::from_path(bam_path.as_path())?; + reader.set_threads(self.threads)?; + let tid_to_chrom = get_tid_to_chrom(&reader)?; let status_probs = process_bam_file( &mut reader, &read_filter, diff --git a/modkit/tests/test_validate.rs b/modkit/tests/test_validate.rs index aba45a9d..57ce90a2 100644 --- a/modkit/tests/test_validate.rs +++ b/modkit/tests/test_validate.rs @@ -1,10 +1,46 @@ use crate::common::run_modkit; use anyhow::Context; -use std::fs::File; +use std::collections::{HashMap, HashSet}; +use std::fs::{self, File}; use std::io::{BufRead, BufReader}; +use std::path::Path; +use tempfile::tempdir; mod common; +fn run_validate( + bam_bed_pairs: &[(&Path, &Path)], + output_file: &Path, +) -> (String, HashMap) { + let mut args = vec!["validate".to_string()]; + for (bam, bed) in bam_bed_pairs { + args.push("--bam-and-bed".to_string()); + args.push(bam.to_str().unwrap().to_string()); + args.push(bed.to_str().unwrap().to_string()); + } + args.extend([ + "--out-filepath".to_string(), + output_file.to_str().unwrap().to_string(), + "--filter-threshold".to_string(), + "0".to_string(), + "--suppress-progress".to_string(), + "--threads".to_string(), + "1".to_string(), + ]); + let args = args.iter().map(String::as_str).collect::>(); + run_modkit(&args).unwrap(); + + let output = fs::read_to_string(output_file).unwrap(); + let fields = output + .lines() + .map(|line| { + let (key, value) = line.split_once(": ").unwrap(); + (key.to_string(), value.to_string()) + }) + .collect(); + (output, fields) +} + #[test] fn test_validate_help() { run_modkit(&["validate", "--help"]) @@ -52,3 +88,75 @@ fn test_validate_expected() { } } } + +#[test] +fn test_validate_reopens_bam_for_each_truth_bed() { + let temp_dir = tempdir().unwrap(); + let bam = Path::new("../tests/resources/input_5mC.bam"); + let truth = Path::new("../tests/resources/CGI_ladder_3.6kb_ref_CG_5mC.bed"); + let truth_lines = fs::read_to_string(truth) + .unwrap() + .lines() + .map(str::to_string) + .collect::>(); + assert_eq!(truth_lines.len(), 48); + + let mut split_a_lines = Vec::new(); + let mut split_b_lines = Vec::new(); + for (idx, line) in truth_lines.iter().enumerate() { + if idx % 2 == 0 { + split_a_lines.push(line.clone()); + } else { + split_b_lines.push(line.clone()); + } + } + let truth_set = truth_lines.iter().collect::>(); + let split_a_set = split_a_lines.iter().collect::>(); + let split_b_set = split_b_lines.iter().collect::>(); + assert_eq!(truth_set.len(), truth_lines.len()); + assert!(split_a_set.is_disjoint(&split_b_set)); + assert_eq!( + split_a_set.union(&split_b_set).copied().collect::>(), + truth_set + ); + + let split_a = temp_dir.path().join("truth-a.bed"); + let split_b = temp_dir.path().join("truth-b.bed"); + fs::write(&split_a, format!("{}\n", split_a_lines.join("\n"))).unwrap(); + fs::write(&split_b, format!("{}\n", split_b_lines.join("\n"))).unwrap(); + + let (union_output, union) = run_validate( + &[(bam, truth)], + &temp_dir.path().join("union-output.tsv"), + ); + assert_eq!( + union["full_contingency_table"], + "[[\"ground_truth_label\",\"C\",\"h\",\"m\",\"No Call\",\"A\",\"G\",\"T\",\"Deletion\"],[\"m\",16625,20432,120926,9069,4826,8337,859,4608]]" + ); + assert_eq!( + union["raw_contingency_table"], + "[[\"ground_truth_label\",\"C\",\"h\",\"m\"],[\"m\",16625,20432,120926]]" + ); + + let (split_ab_output, split_ab) = run_validate( + &[(bam, split_a.as_path()), (bam, split_b.as_path())], + &temp_dir.path().join("split-ab-output.tsv"), + ); + assert_eq!( + split_ab["full_contingency_table"], + union["full_contingency_table"] + ); + assert_eq!(split_ab, union); + assert_eq!(split_ab_output, union_output); + + let (split_ba_output, split_ba) = run_validate( + &[(bam, split_b.as_path()), (bam, split_a.as_path())], + &temp_dir.path().join("split-ba-output.tsv"), + ); + assert_eq!( + split_ba["full_contingency_table"], + union["full_contingency_table"] + ); + assert_eq!(split_ba, union); + assert_eq!(split_ba_output, union_output); +} From 56ce4f38e1641b5fcf3bf74f6d9092724f31ff88 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 05:52:38 -0700 Subject: [PATCH 02/11] Count unseeded validate observations --- modkit-core/src/mod_bam.rs | 31 ++ modkit-core/src/validate/subcommand.rs | 41 ++- .../validate/unseeded_observations_tests.rs | 313 ++++++++++++++++++ modkit/tests/test_validate.rs | 248 +++++++++++++- 4 files changed, 622 insertions(+), 11 deletions(-) create mode 100644 modkit-core/src/validate/unseeded_observations_tests.rs diff --git a/modkit-core/src/mod_bam.rs b/modkit-core/src/mod_bam.rs index ef6a80e7..fa1e0599 100644 --- a/modkit-core/src/mod_bam.rs +++ b/modkit-core/src/mod_bam.rs @@ -1506,6 +1506,7 @@ fn parse_raw_mod_tags(record: &bam::Record) -> MkResult { pub struct ModBaseInfo { pub pos_seq_base_mod_probs: HashMap, pub neg_seq_base_mod_probs: HashMap, + modified_primary_base_strands: HashMap>, converters: HashMap, pub mm_style: &'static str, pub ml_style: &'static str, @@ -1531,6 +1532,24 @@ impl ModBaseInfo { HashMap::::new(); let mut neg_seq_base_mod_probs = HashMap::::new(); + let modified_primary_base_strands = tag_infos.iter().fold( + HashMap::>::new(), + |mut strands_by_base, tag_info| { + for &raw_base in tag_info.fundamental_base.expand_bases() { + let modified_primary_base = + if tag_info.strand == Strand::Negative { + raw_base.complement() + } else { + raw_base + }; + strands_by_base + .entry(modified_primary_base) + .or_default() + .insert(tag_info.strand); + } + strands_by_base + }, + ); let mut converters = HashMap::new(); let mut pointer = 0usize; @@ -1604,12 +1623,24 @@ impl ModBaseInfo { Ok(Self { pos_seq_base_mod_probs, neg_seq_base_mod_probs, + modified_primary_base_strands, converters, mm_style: raw_mod_tags.mm_style, ml_style: raw_mod_tags.ml_style, }) } + pub(crate) fn mod_strands_for_modified_primary_base( + &self, + canonical_base: DnaBase, + ) -> impl Iterator + '_ { + self.modified_primary_base_strands + .get(&canonical_base) + .into_iter() + .flatten() + .copied() + } + pub fn into_iter_base_mod_probs( self, ) -> ( diff --git a/modkit-core/src/validate/subcommand.rs b/modkit-core/src/validate/subcommand.rs index 69e23caa..bcceb730 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -38,6 +38,10 @@ use crate::util::{ record_is_not_primary, Strand, }; +#[cfg(test)] +#[path = "unseeded_observations_tests.rs"] +mod unseeded_observations_tests; + /// todo investigate using this type in BaseModCall #[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)] enum BaseStatus { @@ -297,6 +301,17 @@ fn process_bam_record( let cgt_mod_pos = mod_positions .get(chrom) .ok_or_else(|| anyhow!("No ground truth on this contig.",))?; + let alignment_strand = + if record.is_reverse() { Strand::Negative } else { Strand::Positive }; + let mut called_ref_pos = mbi + .mod_strands_for_modified_primary_base(can_base) + .map(|mod_strand| { + ( + get_reference_mod_strand(mod_strand, alignment_strand), + HashSet::new(), + ) + }) + .collect::>(); let mbp = ReadBaseModProfile::process_record( &record, &record_name, @@ -330,20 +345,21 @@ fn process_bam_record( .map(|gt_code| (mod_call, ref_pos, ref_strand, gt_code)) }); - let mut called_ref_pos = HashMap::new(); let mut result = HashMap::new(); for (mod_call, ref_pos, ref_mod_strand, gt_code) in mod_call_iter { - called_ref_pos - .entry(ref_mod_strand) - .or_insert_with(HashSet::new) - .insert(ref_pos); + let Some(positions) = called_ref_pos.get_mut(&ref_mod_strand) else { + continue; + }; + positions.insert(ref_pos); - if mod_call.canonical_base != can_base { + let modified_primary_base = if mod_call.mod_strand == Strand::Negative { + mod_call.canonical_base.complement() + } else { + mod_call.canonical_base + }; + if modified_primary_base != can_base { result - .entry(( - *gt_code, - BaseStatus::Mismatch(mod_call.canonical_base), - )) + .entry((*gt_code, BaseStatus::Mismatch(modified_primary_base))) .or_insert_with(Vec::new) .push(f32::NAN); continue; @@ -392,6 +408,11 @@ fn process_bam_record( if record.is_reverse() { base = base.complement(); } + let mod_strand = + get_reference_mod_strand(*strand, alignment_strand); + if mod_strand == Strand::Negative { + base = base.complement(); + } if base == can_base { result .entry((*gt_code, BaseStatus::NoCall)) diff --git a/modkit-core/src/validate/unseeded_observations_tests.rs b/modkit-core/src/validate/unseeded_observations_tests.rs new file mode 100644 index 00000000..93ac9ddc --- /dev/null +++ b/modkit-core/src/validate/unseeded_observations_tests.rs @@ -0,0 +1,313 @@ +use std::collections::{BTreeMap, HashMap}; + +use rust_htslib::bam::record::{Aux, Cigar, CigarString}; +use rust_htslib::bam::Record; + +use super::{ + machine_parseable_table, process_bam_record, BaseStatus, + ChromStrandPositionNames, StatusProbs, TidToChrom, +}; +use crate::mod_base_code::{DnaBase, ModCodeRepr}; +use crate::util::Strand; + +const M: BaseStatus = BaseStatus::Modified(ModCodeRepr::Code('m')); + +fn make_record( + name: &str, + sequence: &str, + cigar: Vec, + mm: &str, + ml: &[u8], +) -> Record { + let mut record = Record::new(); + let cigar = CigarString(cigar); + record.set( + name.as_bytes(), + Some(&cigar), + sequence.as_bytes(), + &vec![30; sequence.len()], + ); + record.set_tid(0); + record.set_pos(0); + record.set_mapq(60); + record.push_aux(b"MM", Aux::String(mm)).unwrap(); + record.push_aux(b"ML", Aux::ArrayU8(ml.into())).unwrap(); + record.push_aux(b"MN", Aux::U32(sequence.len() as u32)).unwrap(); + record +} + +fn truth( + strands: impl IntoIterator, + status: BaseStatus, +) -> ChromStrandPositionNames { + let strand_positions = strands + .into_iter() + .map(|strand| (strand, BTreeMap::from([(1, status)]))) + .collect(); + HashMap::from([("chr1".to_string(), strand_positions)]) +} + +fn classify( + records: impl IntoIterator, + truth: &ChromStrandPositionNames, +) -> StatusProbs { + let tid_to_chrom: TidToChrom = HashMap::from([(0, "chr1".to_string())]); + let mut combined = StatusProbs::new(); + for record in records { + let observed = process_bam_record( + &record, + truth, + &tid_to_chrom, + DnaBase::C, + None, + None, + ) + .unwrap(); + for (status, probs) in observed { + combined.entry(status).or_default().extend(probs); + } + } + combined +} + +fn called_record(name: &str, mod_strand: char) -> Record { + let (sequence, fundamental_base) = match mod_strand { + '+' => ("CC", 'C'), + '-' => ("GG", 'G'), + _ => panic!("invalid modification strand"), + }; + make_record( + name, + sequence, + vec![Cigar::Match(2)], + &format!("{fundamental_base}{mod_strand}m?,1;"), + &[255], + ) +} + +fn uncalled_record(name: &str, sequence: &str, mod_strand: char) -> Record { + let fundamental_base = match mod_strand { + '+' => 'C', + '-' => 'G', + _ => panic!("invalid modification strand"), + }; + make_record( + name, + sequence, + vec![Cigar::Match(2)], + &format!("{fundamental_base}{mod_strand}m?,0;"), + &[255], + ) +} + +fn deletion_record(name: &str) -> Record { + make_record( + name, + "C", + vec![Cigar::Match(1), Cigar::Del(1)], + "C+m?,0;", + &[255], + ) +} + +fn empty_descriptor_record( + name: &str, + sequence: &str, + mod_strand: char, +) -> Record { + let fundamental_base = match mod_strand { + '+' => 'C', + '-' => 'G', + _ => panic!("invalid modification strand"), + }; + make_record( + name, + sequence, + vec![Cigar::Match(sequence.len() as u32)], + &format!("{fundamental_base}{mod_strand}m?;"), + &[], + ) +} + +fn assert_exact_counts( + observed: &StatusProbs, + expected: &[((BaseStatus, BaseStatus), usize)], +) { + let observed = observed + .iter() + .map(|(statuses, probs)| (*statuses, probs.len())) + .collect::>(); + let expected = expected.iter().copied().collect::>(); + assert_eq!(observed, expected); +} + +#[test] +fn explicit_mode_counts_no_calls_without_truth_overlapping_seed_call() { + let records = + (0..3).map(|i| called_record(&format!("called-{i}"), '+')).chain( + (0..6) + .map(|i| uncalled_record(&format!("uncalled-{i}"), "CC", '+')), + ); + let observed = classify(records, &truth([Strand::Positive], M)); + + assert_exact_counts( + &observed, + &[((M, M), 3), ((M, BaseStatus::NoCall), 6)], + ); + assert_eq!( + machine_parseable_table(DnaBase::C, &observed), + "[[\"ground_truth_label\",\"m\",\"No Call\"],[\"m\",3,6]]" + ); +} + +#[test] +fn explicit_mode_counts_unseeded_no_calls_mismatches_and_deletions() { + let records = (0..2) + .map(|i| called_record(&format!("called-{i}"), '+')) + .chain( + (0..6) + .map(|i| uncalled_record(&format!("uncalled-{i}"), "CC", '+')), + ) + .chain( + (0..2) + .map(|i| uncalled_record(&format!("mismatch-{i}"), "CA", '+')), + ) + .chain((0..2).map(|i| deletion_record(&format!("deletion-{i}")))); + let observed = classify(records, &truth([Strand::Positive], M)); + + assert_exact_counts( + &observed, + &[ + ((M, M), 2), + ((M, BaseStatus::NoCall), 6), + ((M, BaseStatus::Mismatch(DnaBase::A)), 2), + ((M, BaseStatus::Deletion), 2), + ], + ); + assert_eq!( + machine_parseable_table(DnaBase::C, &observed), + "[[\"ground_truth_label\",\"m\",\"No Call\",\"A\",\"Deletion\"],[\"m\",2,6,2,2]]" + ); +} + +#[test] +fn explicit_mode_counts_only_the_descriptor_reference_strand() { + let records = std::iter::once(called_record("positive-called", '+')) + .chain((0..5).map(|i| { + uncalled_record(&format!("positive-uncalled-{i}"), "CC", '+') + })) + .chain(std::iter::once(called_record("negative-called", '-'))) + .chain((0..5).map(|i| { + uncalled_record(&format!("negative-uncalled-{i}"), "GG", '-') + })); + let observed = + classify(records, &truth([Strand::Positive, Strand::Negative], M)); + + assert_exact_counts( + &observed, + &[((M, M), 2), ((M, BaseStatus::NoCall), 10)], + ); + assert_eq!( + machine_parseable_table(DnaBase::C, &observed), + "[[\"ground_truth_label\",\"m\",\"No Call\"],[\"m\",2,10]]" + ); +} + +#[test] +fn implicit_mode_canonical_calls_are_not_reclassified_as_no_calls() { + let records = (0..3).map(|i| { + make_record( + &format!("implicit-{i}"), + "CC", + vec![Cigar::Match(2)], + "C+m.,0;", + &[255], + ) + }); + let observed = + classify(records, &truth([Strand::Positive], BaseStatus::Canonical)); + + assert_exact_counts( + &observed, + &[((BaseStatus::Canonical, BaseStatus::Canonical), 3)], + ); + assert_eq!( + machine_parseable_table(DnaBase::C, &observed), + "[[\"ground_truth_label\",\"C\"],[\"C\",3]]" + ); +} + +#[test] +fn explicit_empty_positive_descriptor_counts_no_call() { + let observed = classify( + [empty_descriptor_record("empty-positive", "CC", '+')], + &truth([Strand::Positive], M), + ); + + assert_exact_counts(&observed, &[((M, BaseStatus::NoCall), 1)]); +} + +#[test] +fn explicit_empty_negative_descriptor_counts_no_call() { + let observed = classify( + [empty_descriptor_record("empty-negative", "GG", '-')], + &truth([Strand::Negative], M), + ); + + assert_exact_counts(&observed, &[((M, BaseStatus::NoCall), 1)]); +} + +#[test] +fn explicit_empty_n_descriptor_expands_to_target_base() { + let record = + make_record("empty-n", "CC", vec![Cigar::Match(2)], "N+m?;", &[]); + let observed = classify([record], &truth([Strand::Positive], M)); + + assert_exact_counts(&observed, &[((M, BaseStatus::NoCall), 1)]); +} + +#[test] +fn unrelated_canonical_base_descriptor_does_not_create_observations() { + let unrelated = make_record( + "adenine-only", + "AA", + vec![Cigar::Match(2)], + "A+a?,1;", + &[255], + ); + let observed = + classify([unrelated], &truth([Strand::Positive, Strand::Negative], M)); + + assert!(observed.is_empty()); +} + +#[test] +fn minus_strand_call_on_reverse_alignment_uses_modified_primary_base() { + let mut record = make_record( + "reverse-g-minus", + "CC", + vec![Cigar::Match(2)], + "G-m?,0;", + &[255], + ); + record.set_reverse(); + + let observed = classify([record], &truth([Strand::Positive], M)); + + assert_exact_counts(&observed, &[((M, M), 1)]); +} + +#[test] +fn unrelated_empty_descriptor_does_not_create_observations() { + let unrelated = make_record( + "empty-adenine-only", + "AA", + vec![Cigar::Match(2)], + "A+a?;", + &[], + ); + let observed = + classify([unrelated], &truth([Strand::Positive, Strand::Negative], M)); + + assert!(observed.is_empty()); +} diff --git a/modkit/tests/test_validate.rs b/modkit/tests/test_validate.rs index 57ce90a2..0763e575 100644 --- a/modkit/tests/test_validate.rs +++ b/modkit/tests/test_validate.rs @@ -1,8 +1,11 @@ use crate::common::run_modkit; use anyhow::Context; +use rust_htslib::bam::header::HeaderRecord; +use rust_htslib::bam::record::{Aux, Cigar, CigarString}; +use rust_htslib::bam::{Format, Header, Record, Writer}; use std::collections::{HashMap, HashSet}; use std::fs::{self, File}; -use std::io::{BufRead, BufReader}; +use std::io::{BufRead, BufReader, Write}; use std::path::Path; use tempfile::tempdir; @@ -160,3 +163,246 @@ fn test_validate_reopens_bam_for_each_truth_bed() { assert_eq!(split_ba, union); assert_eq!(split_ba_output, union_output); } + +fn synthetic_record( + name: &str, + sequence: &str, + cigar: Vec, + mm: &str, + ml: &[u8], + nm: u32, +) -> Record { + let mut record = Record::new(); + let cigar = CigarString(cigar); + record.set( + name.as_bytes(), + Some(&cigar), + sequence.as_bytes(), + &vec![30; sequence.len()], + ); + record.set_tid(0); + record.set_pos(0); + record.set_mapq(60); + record.push_aux(b"MM", Aux::String(mm)).unwrap(); + record.push_aux(b"ML", Aux::ArrayU8(ml.into())).unwrap(); + record.push_aux(b"MN", Aux::U32(sequence.len() as u32)).unwrap(); + record.push_aux(b"NM", Aux::U32(nm)).unwrap(); + record +} + +fn synthetic_called(name: &str, mod_strand: char) -> Record { + let (sequence, fundamental_base) = match mod_strand { + '+' => ("CC", 'C'), + '-' => ("GG", 'G'), + _ => panic!("invalid modification strand"), + }; + synthetic_record( + name, + sequence, + vec![Cigar::Match(2)], + &format!("{fundamental_base}{mod_strand}m?,1;"), + &[255], + 0, + ) +} + +fn synthetic_uncalled(name: &str, sequence: &str, mod_strand: char) -> Record { + let (expected_sequence, fundamental_base) = match mod_strand { + '+' => ("CC", 'C'), + '-' => ("GG", 'G'), + _ => panic!("invalid modification strand"), + }; + synthetic_record( + name, + sequence, + vec![Cigar::Match(2)], + &format!("{fundamental_base}{mod_strand}m?,0;"), + &[255], + u32::from(sequence != expected_sequence), + ) +} + +fn write_synthetic_bam(path: &Path, records: Vec) { + let mut header = Header::new(); + let mut sq = HeaderRecord::new(b"SQ"); + sq.push_tag(b"SN", "chr1"); + sq.push_tag(b"LN", 2); + header.push_record(&sq); + let mut writer = Writer::from_path(path, &header, Format::Bam).unwrap(); + for record in records { + writer.write(&record).unwrap(); + } +} + +fn assert_validate_fixture( + root: &Path, + name: &str, + records: Vec, + bed: &str, + expected_full_table: &str, +) { + let bam = root.join(format!("{name}.bam")); + let bed_path = root.join(format!("{name}.bed")); + let output = root.join(format!("{name}.tsv")); + let check_tags_dir = root.join(format!("{name}-check-tags")); + write_synthetic_bam(&bam, records); + let mut bed_file = File::create(&bed_path).unwrap(); + bed_file.write_all(bed.as_bytes()).unwrap(); + + run_modkit(&[ + "modbam", + "check-tags", + bam.to_str().unwrap(), + "--ignore-index", + "--suppress-progress", + "--out-dir", + check_tags_dir.to_str().unwrap(), + ]) + .with_context(|| { + format!("{name}: synthetic records should pass check-tags") + }) + .unwrap(); + + run_modkit(&[ + "validate", + "--bam-and-bed", + bam.to_str().unwrap(), + bed_path.to_str().unwrap(), + "--canonical-base", + "C", + "--filter-threshold", + "0", + "--threads", + "1", + "--suppress-progress", + "--out-filepath", + output.to_str().unwrap(), + ]) + .with_context(|| format!("{name}: validate should succeed")) + .unwrap(); + + let first_line = BufReader::new(File::open(output).unwrap()) + .lines() + .next() + .unwrap() + .unwrap(); + assert_eq!( + first_line, + format!("full_contingency_table: {expected_full_table}"), + "{name}" + ); +} + +#[test] +fn test_validate_counts_observations_without_truth_overlapping_seed_calls() { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + + let gt_m_records = + (0..3) + .map(|i| synthetic_called(&format!("called-{i}"), '+')) + .chain((0..6).map(|i| { + synthetic_uncalled(&format!("uncalled-{i}"), "CC", '+') + })) + .collect(); + assert_validate_fixture( + root, + "gt-m", + gt_m_records, + "chr1\t1\t2\tm\t0\t+\n", + "[[\"ground_truth_label\",\"m\",\"No Call\"],[\"m\",3,6]]", + ); + + let mirror_records = + (0..2) + .map(|i| synthetic_called(&format!("called-{i}"), '+')) + .chain((0..6).map(|i| { + synthetic_uncalled(&format!("uncalled-{i}"), "CC", '+') + })) + .chain((0..2).map(|i| { + synthetic_uncalled(&format!("mismatch-{i}"), "CA", '+') + })) + .chain((0..2).map(|i| { + synthetic_record( + &format!("deletion-{i}"), + "C", + vec![Cigar::Match(1), Cigar::Del(1)], + "C+m?,0;", + &[255], + 1, + ) + })) + .collect(); + assert_validate_fixture( + root, + "mirror", + mirror_records, + "chr1\t1\t2\tm\t0\t+\n", + "[[\"ground_truth_label\",\"m\",\"No Call\",\"A\",\"Deletion\"],[\"m\",2,6,2,2]]", + ); + + let both_strands_records = + std::iter::once(synthetic_called("positive-called", '+')) + .chain((0..5).map(|i| { + synthetic_uncalled(&format!("positive-uncalled-{i}"), "CC", '+') + })) + .chain(std::iter::once(synthetic_called("negative-called", '-'))) + .chain((0..5).map(|i| { + synthetic_uncalled(&format!("negative-uncalled-{i}"), "GG", '-') + })) + .collect(); + assert_validate_fixture( + root, + "both-strands", + both_strands_records, + "chr1\t1\t2\tm\t0\t+\nchr1\t1\t2\tm\t0\t-\n", + "[[\"ground_truth_label\",\"m\",\"No Call\"],[\"m\",2,10]]", + ); + + let empty_descriptor_records = vec![ + synthetic_called("called-positive", '+'), + synthetic_record( + "empty-positive", + "CC", + vec![Cigar::Match(2)], + "C+m?;", + &[], + 0, + ), + synthetic_record( + "empty-negative", + "GG", + vec![Cigar::Match(2)], + "G-m?;", + &[], + 0, + ), + ]; + assert_validate_fixture( + root, + "empty-descriptors", + empty_descriptor_records, + "chr1\t1\t2\tm\t0\t+\nchr1\t1\t2\tm\t0\t-\n", + "[[\"ground_truth_label\",\"m\",\"No Call\"],[\"m\",1,2]]", + ); + + let implicit_records = (0..3) + .map(|i| { + synthetic_record( + &format!("implicit-{i}"), + "CC", + vec![Cigar::Match(2)], + "C+m.,0;", + &[255], + 0, + ) + }) + .collect(); + assert_validate_fixture( + root, + "implicit", + implicit_records, + "chr1\t1\t2\t-\t0\t+\n", + "[[\"ground_truth_label\",\"C\"],[\"C\",3]]", + ); +} From 97fda63e08beb8205d93c8e0e42d7739131c5eeb Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 05:17:13 -0700 Subject: [PATCH 03/11] Fix validate reference-skip classification --- modkit-core/src/validate/subcommand.rs | 90 ++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/modkit-core/src/validate/subcommand.rs b/modkit-core/src/validate/subcommand.rs index bcceb730..0f8dd0c3 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -280,6 +280,15 @@ fn derive_canonical_base( // ground truth and observed base status pointing to vector of mod probabilities type StatusProbs = HashMap<(BaseStatus, BaseStatus), Vec>; +fn position_is_reference_skip( + position: i64, + reference_skips: &[[i64; 2]], +) -> bool { + let insertion_idx = + reference_skips.partition_point(|[start, _]| *start <= position); + insertion_idx > 0 && position < reference_skips[insertion_idx - 1][1] +} + fn process_bam_record( record: &Record, mod_positions: &ChromStrandPositionNames, @@ -387,6 +396,7 @@ fn process_bam_record( let q_seq = record.seq(); let ref_to_query: FxHashMap = record.aligned_pairs().map(|pos| (pos[1], pos[0])).collect(); + let reference_skips = record.introns().collect::>(); for (strand, positions) in called_ref_pos.iter() { let Some(cs_mod_pos) = cgt_mod_pos.get(&strand) else { // should be unnecessary @@ -398,6 +408,9 @@ fn process_bam_record( continue; }; let Some(q_pos) = ref_to_query.get(pos) else { + if position_is_reference_skip(*pos, &reference_skips) { + continue; + } result .entry((*gt_code, BaseStatus::Deletion)) .or_insert_with(Vec::new) @@ -1153,3 +1166,80 @@ impl ValidateFromModBam { Ok(()) } } + +#[cfg(test)] +mod tests { + use rust_htslib::bam::record::{Aux, Cigar, CigarString}; + + use super::*; + + fn make_record(gap: Cigar) -> Record { + let cigar = CigarString(vec![Cigar::Match(1), gap, Cigar::Match(1)]); + let mut record = Record::new(); + record.set(b"read", Some(&cigar), b"CC", &[255, 255]); + record.set_tid(0); + record.set_pos(0); + record.push_aux(b"MM", Aux::String("C+m?,0;")).unwrap(); + record.push_aux(b"ML", Aux::ArrayU8((&[255][..]).into())).unwrap(); + record + } + + fn category_count( + status_probs: &StatusProbs, + truth: BaseStatus, + call: BaseStatus, + ) -> usize { + status_probs.get(&(truth, call)).map(Vec::len).unwrap_or(0) + } + + fn classify(gap: Cigar) -> StatusProbs { + let truth_status = BaseStatus::Modified(ModCodeRepr::Code('m')); + let truth = HashMap::from([( + "chr1".to_string(), + HashMap::from([( + Strand::Positive, + (0..5).map(|pos| (pos, truth_status)).collect(), + )]), + )]); + let tid_to_chrom = HashMap::from([(0, "chr1".to_string())]); + + process_bam_record( + &make_record(gap), + &truth, + &tid_to_chrom, + DnaBase::C, + None, + None, + ) + .unwrap() + } + + #[test] + fn reference_skips_are_not_classified_as_deletions() { + let truth_status = BaseStatus::Modified(ModCodeRepr::Code('m')); + + let skipped = classify(Cigar::RefSkip(3)); + assert_eq!(category_count(&skipped, truth_status, truth_status), 1); + assert_eq!( + category_count(&skipped, truth_status, BaseStatus::NoCall), + 1 + ); + assert_eq!( + category_count(&skipped, truth_status, BaseStatus::Deletion), + 0 + ); + assert_eq!(skipped.values().map(Vec::len).sum::(), 2); + + let deleted = classify(Cigar::Del(3)); + assert_eq!(category_count(&deleted, truth_status, truth_status), 1); + assert_eq!( + category_count(&deleted, truth_status, BaseStatus::NoCall), + 1 + ); + assert_eq!( + category_count(&deleted, truth_status, BaseStatus::Deletion), + 3 + ); + assert_eq!(deleted.values().map(Vec::len).sum::(), 5); + } +} From d636dd3412034c7a48e42a7e7e758ebd9868ef61 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 05:24:07 -0700 Subject: [PATCH 04/11] Handle ambiguous validate bases site-locally --- modkit-core/src/validate/subcommand.rs | 57 ++++++++++++++++++++++++-- 1 file changed, 53 insertions(+), 4 deletions(-) diff --git a/modkit-core/src/validate/subcommand.rs b/modkit-core/src/validate/subcommand.rs index 0f8dd0c3..f43070e2 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -417,7 +417,14 @@ fn process_bam_record( .push(f32::NAN); continue; }; - let mut base = DnaBase::parse(q_seq[*q_pos as usize] as char)?; + let Ok(mut base) = DnaBase::parse(q_seq[*q_pos as usize] as char) + else { + result + .entry((*gt_code, BaseStatus::NoCall)) + .or_insert_with(Vec::new) + .push(f32::NAN); + continue; + }; if record.is_reverse() { base = base.complement(); } @@ -1184,6 +1191,17 @@ mod tests { record } + fn make_record_with_ambiguous_bases() -> Record { + let cigar = CigarString(vec![Cigar::Match(5)]); + let mut record = Record::new(); + record.set(b"ambiguous", Some(&cigar), b"CNRAC", &[255; 5]); + record.set_tid(0); + record.set_pos(0); + record.push_aux(b"MM", Aux::String("C+m?,0,0;")).unwrap(); + record.push_aux(b"ML", Aux::ArrayU8((&[255, 0][..]).into())).unwrap(); + record + } + fn category_count( status_probs: &StatusProbs, truth: BaseStatus, @@ -1192,7 +1210,7 @@ mod tests { status_probs.get(&(truth, call)).map(Vec::len).unwrap_or(0) } - fn classify(gap: Cigar) -> StatusProbs { + fn classify_record(record: &Record) -> anyhow::Result { let truth_status = BaseStatus::Modified(ModCodeRepr::Code('m')); let truth = HashMap::from([( "chr1".to_string(), @@ -1204,14 +1222,17 @@ mod tests { let tid_to_chrom = HashMap::from([(0, "chr1".to_string())]); process_bam_record( - &make_record(gap), + record, &truth, &tid_to_chrom, DnaBase::C, None, None, ) - .unwrap() + } + + fn classify(gap: Cigar) -> StatusProbs { + classify_record(&make_record(gap)).unwrap() } #[test] @@ -1242,4 +1263,32 @@ mod tests { ); assert_eq!(deleted.values().map(Vec::len).sum::(), 5); } + + #[test] + fn ambiguous_aligned_bases_are_site_local_no_calls() { + let truth_status = BaseStatus::Modified(ModCodeRepr::Code('m')); + let classified = classify_record(&make_record_with_ambiguous_bases()) + .expect( + "N and IUPAC query bases must not discard the whole record", + ); + + assert_eq!(category_count(&classified, truth_status, truth_status), 1); + assert_eq!( + category_count(&classified, truth_status, BaseStatus::Canonical), + 1 + ); + assert_eq!( + category_count(&classified, truth_status, BaseStatus::NoCall), + 2 + ); + assert_eq!( + category_count( + &classified, + truth_status, + BaseStatus::Mismatch(DnaBase::A), + ), + 1 + ); + assert_eq!(classified.values().map(Vec::len).sum::(), 5); + } } From dbe1eb2303608c2b37df9a50decafcc39318aa7b Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 06:35:29 -0700 Subject: [PATCH 05/11] Cover composed validate classifications --- .../validate/unseeded_observations_tests.rs | 59 +++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/modkit-core/src/validate/unseeded_observations_tests.rs b/modkit-core/src/validate/unseeded_observations_tests.rs index 93ac9ddc..4efe83d1 100644 --- a/modkit-core/src/validate/unseeded_observations_tests.rs +++ b/modkit-core/src/validate/unseeded_observations_tests.rs @@ -47,6 +47,20 @@ fn truth( HashMap::from([("chr1".to_string(), strand_positions)]) } +fn truth_positions( + strand: Strand, + positions: impl IntoIterator, + status: BaseStatus, +) -> ChromStrandPositionNames { + HashMap::from([( + "chr1".to_string(), + HashMap::from([( + strand, + positions.into_iter().map(|pos| (pos, status)).collect(), + )]), + )]) +} + fn classify( records: impl IntoIterator, truth: &ChromStrandPositionNames, @@ -311,3 +325,48 @@ fn unrelated_empty_descriptor_does_not_create_observations() { assert!(observed.is_empty()); } + +#[test] +fn explicit_unseeded_sites_compose_reference_skips_and_iupac_no_calls() { + let record = make_record( + "unseeded-refskip-iupac", + "CNRAC", + vec![Cigar::Match(1), Cigar::RefSkip(3), Cigar::Match(4)], + "C+m?,0;", + &[255], + ); + let observed = + classify([record], &truth_positions(Strand::Positive, 1..8, M)); + + assert_exact_counts( + &observed, + &[ + ((M, BaseStatus::NoCall), 3), + ((M, BaseStatus::Mismatch(DnaBase::A)), 1), + ], + ); +} + +#[test] +fn reverse_minus_fallback_preserves_no_call_mismatch_and_deletion_orientation() +{ + let mut record = make_record( + "reverse-minus-fallback", + "CAC", + vec![Cigar::Match(2), Cigar::Del(1), Cigar::Match(1)], + "G-m?;", + &[], + ); + record.set_reverse(); + let observed = + classify([record], &truth_positions(Strand::Positive, 0..3, M)); + + assert_exact_counts( + &observed, + &[ + ((M, BaseStatus::NoCall), 1), + ((M, BaseStatus::Mismatch(DnaBase::A)), 1), + ((M, BaseStatus::Deletion), 1), + ], + ); +} From 6ea6b0cd043c07dfac9938cb1f57cc2686bab192 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 06:42:07 -0700 Subject: [PATCH 06/11] Respect edge filters in validate fallback calls --- modkit-core/src/validate/subcommand.rs | 24 +++- .../validate/unseeded_observations_tests.rs | 106 +++++++++++++++++- 2 files changed, 127 insertions(+), 3 deletions(-) diff --git a/modkit-core/src/validate/subcommand.rs b/modkit-core/src/validate/subcommand.rs index f43070e2..540cb7ed 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -417,8 +417,28 @@ fn process_bam_record( .push(f32::NAN); continue; }; - let Ok(mut base) = DnaBase::parse(q_seq[*q_pos as usize] as char) - else { + let q_pos = *q_pos as usize; + if let Some(edge_filter) = edge_filter { + // EdgeFilter uses forward/as-sequenced query coordinates, + // while aligned_pairs follows the alignment orientation. + let forward_q_pos = if record.is_reverse() { + q_pos + .checked_add(1) + .and_then(|p| record.seq_len().checked_sub(p)) + } else { + Some(q_pos) + }; + let Some(forward_q_pos) = forward_q_pos else { + continue; + }; + if !edge_filter + .keep_position(forward_q_pos, record.seq_len()) + .unwrap_or(false) + { + continue; + } + } + let Ok(mut base) = DnaBase::parse(q_seq[q_pos] as char) else { result .entry((*gt_code, BaseStatus::NoCall)) .or_insert_with(Vec::new) diff --git a/modkit-core/src/validate/unseeded_observations_tests.rs b/modkit-core/src/validate/unseeded_observations_tests.rs index 4efe83d1..a348854f 100644 --- a/modkit-core/src/validate/unseeded_observations_tests.rs +++ b/modkit-core/src/validate/unseeded_observations_tests.rs @@ -7,6 +7,7 @@ use super::{ machine_parseable_table, process_bam_record, BaseStatus, ChromStrandPositionNames, StatusProbs, TidToChrom, }; +use crate::mod_bam::EdgeFilter; use crate::mod_base_code::{DnaBase, ModCodeRepr}; use crate::util::Strand; @@ -64,6 +65,14 @@ fn truth_positions( fn classify( records: impl IntoIterator, truth: &ChromStrandPositionNames, +) -> StatusProbs { + classify_with_edge_filter(records, truth, None) +} + +fn classify_with_edge_filter( + records: impl IntoIterator, + truth: &ChromStrandPositionNames, + edge_filter: Option<&EdgeFilter>, ) -> StatusProbs { let tid_to_chrom: TidToChrom = HashMap::from([(0, "chr1".to_string())]); let mut combined = StatusProbs::new(); @@ -74,7 +83,7 @@ fn classify( &tid_to_chrom, DnaBase::C, None, - None, + edge_filter, ) .unwrap(); for (status, probs) in observed { @@ -370,3 +379,98 @@ fn reverse_minus_fallback_preserves_no_call_mismatch_and_deletion_orientation() ], ); } + +fn positive_truth_range(end: i64) -> ChromStrandPositionNames { + HashMap::from([( + "chr1".to_string(), + HashMap::from([( + Strand::Positive, + (0..end).map(|position| (position, M)).collect(), + )]), + )]) +} + +fn forward_edge_filter_record() -> Record { + make_record( + "forward-edge-filter", + "CCCACCC", + vec![Cigar::Match(7)], + "C+m?,0,0,1,1;", + &[255; 4], + ) +} + +#[test] +fn ordinary_edge_filter_excludes_filtered_fallback_sites() { + let edge_filter = EdgeFilter::new(1, 1, false); + let observed = classify_with_edge_filter( + [forward_edge_filter_record()], + &positive_truth_range(7), + Some(&edge_filter), + ); + + assert_exact_counts( + &observed, + &[ + ((M, M), 2), + ((M, BaseStatus::NoCall), 2), + ((M, BaseStatus::Mismatch(DnaBase::A)), 1), + ], + ); +} + +#[test] +fn inverted_edge_filter_excludes_filtered_fallback_sites() { + let edge_filter = EdgeFilter::new(1, 1, true); + let observed = classify_with_edge_filter( + [forward_edge_filter_record()], + &positive_truth_range(7), + Some(&edge_filter), + ); + + assert_exact_counts(&observed, &[((M, M), 2)]); +} + +#[test] +fn reverse_alignment_edge_filter_uses_forward_query_coordinates() { + let mut record = make_record( + "reverse-edge-filter", + "CCCCCC", + vec![Cigar::Match(6)], + "G-m?,0,0,0,0,0,0;", + &[255; 6], + ); + record.set_reverse(); + let truth = positive_truth_range(6); + + for inverted in [false, true] { + let edge_filter = EdgeFilter::new(1, 2, inverted); + let observed = classify_with_edge_filter( + [record.clone()], + &truth, + Some(&edge_filter), + ); + assert_exact_counts(&observed, &[((M, M), 3)]); + } +} + +#[test] +fn edge_filter_does_not_reclassify_deletion_without_query_position() { + let record = make_record( + "edge-filter-deletion", + "CCC", + vec![Cigar::Match(1), Cigar::Del(1), Cigar::Match(2)], + "C+m?,0;", + &[255], + ); + let truth = HashMap::from([( + "chr1".to_string(), + HashMap::from([(Strand::Positive, BTreeMap::from([(1, M)]))]), + )]); + let edge_filter = EdgeFilter::new(1, 1, false); + + let observed = + classify_with_edge_filter([record], &truth, Some(&edge_filter)); + + assert_exact_counts(&observed, &[((M, BaseStatus::Deletion), 1)]); +} From 74191180161a72e94be3965ef96c5bf202eee607 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 06:41:14 -0700 Subject: [PATCH 07/11] Handle validate runs without confident calls --- modkit-core/src/validate/subcommand.rs | 94 ++++++++++++++++++--- modkit/tests/test_validate.rs | 112 +++++++++++++++++++++++++ 2 files changed, 193 insertions(+), 13 deletions(-) diff --git a/modkit-core/src/validate/subcommand.rs b/modkit-core/src/validate/subcommand.rs index 540cb7ed..3dae2e32 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -751,11 +751,13 @@ fn machine_parseable_table( gt_codes.iter().chain(call_codes.iter()).unique().collect(); all_codes.sort(); - let mut out_str = "[[\"ground_truth_label\",\"".to_string(); - out_str.push_str( - &all_codes.iter().map(|x| x.human_display(validate_base)).join("\",\""), - ); - out_str.push_str("\"]"); + let mut out_str = "[[\"ground_truth_label\"".to_string(); + for &call_code in &all_codes { + out_str.push_str(",\""); + out_str.push_str(&call_code.human_display(validate_base)); + out_str.push_str("\""); + } + out_str.push_str("]"); for gt_code in >_codes { out_str.push_str(",[\""); out_str.push_str(>_code.human_display(validate_base)); @@ -1076,6 +1078,58 @@ impl ValidateFromModBam { _ => false, }); + // Accuracy and filtering rates have no denominator when every + // observation is a non-call status. + if all_probs.values().all(Vec::is_empty) { + info!("Balancing ground truth call totals"); + print_table(can_base, &all_probs, false, "Balanced counts summary"); + info!("Raw accuracy: NA"); + print_table( + can_base, + &all_probs, + true, + "Raw modified base calls contingency table", + ); + + let threshold = self + .filter_threshold + .map(|threshold| threshold.to_string()) + .unwrap_or_else(|| "NA".to_string()); + info!("Call probability threshold: {threshold}"); + info!("Percent of modified base calls removed: NA"); + info!("{}", Style::new().bold().paint("Filtered accuracy: NA")); + print_table( + can_base, + &all_probs, + true, + "Filtered modified base calls contingency table", + ); + + if let Some(valid_out_handle) = &mut out_handle { + let empty_table = machine_parseable_table(can_base, &all_probs); + valid_out_handle + .write_all( + format!( + concat!( + "raw_accuracy: NA\n", + "raw_contingency_table: {}\n", + "filter_threshold: {}\n", + "percent_of_mod_called_removed: NA\n", + "filtered_accuracy: NA\n", + "filtered_contingency_table: {}\n", + ), + empty_table, threshold, empty_table, + ) + .as_bytes(), + ) + .map_err(|e| { + anyhow::anyhow!("Error writing to file: {}", e) + })?; + } + + return Ok(()); + } + info!("Balancing ground truth call totals"); balance_ground_truth(&mut all_probs)?; print_table(can_base, &all_probs, false, "Balanced counts summary"); @@ -1146,13 +1200,26 @@ impl ValidateFromModBam { .filter(|&((gt_code, call_code), _)| gt_code == call_code) .map(|(_, values)| values.len()) .sum::(); - let filt_acc = 100.0 * correct_filt_calls as f32 / filt_calls as f32; - info!( - "{}", - Style::new() - .bold() - .paint(format!("Filtered accuracy: {:.2}%", filt_acc)), - ); + let filt_acc = if filt_calls == 0 { + None + } else { + Some(100.0 * correct_filt_calls as f32 / filt_calls as f32) + }; + let filtered_accuracy = match filt_acc { + Some(accuracy) => { + info!( + "{}", + Style::new() + .bold() + .paint(format!("Filtered accuracy: {:.2}%", accuracy)), + ); + accuracy.to_string() + } + None => { + info!("{}", Style::new().bold().paint("Filtered accuracy: NA")); + "NA".to_string() + } + }; print_table( can_base, &all_probs, @@ -1176,7 +1243,8 @@ impl ValidateFromModBam { .map_err(|e| anyhow::anyhow!("Error writing to file: {}", e))?; valid_out_handle .write_all( - &format!("filtered_accuracy: {}\n", filt_acc).into_bytes(), + &format!("filtered_accuracy: {}\n", filtered_accuracy) + .into_bytes(), ) .map_err(|e| anyhow::anyhow!("Error writing to file: {}", e))?; valid_out_handle diff --git a/modkit/tests/test_validate.rs b/modkit/tests/test_validate.rs index 0763e575..b1a975a5 100644 --- a/modkit/tests/test_validate.rs +++ b/modkit/tests/test_validate.rs @@ -406,3 +406,115 @@ fn test_validate_counts_observations_without_truth_overlapping_seed_calls() { "[[\"ground_truth_label\",\"C\"],[\"C\",3]]", ); } + +#[test] +fn test_validate_all_no_calls_reports_counts_and_unavailable_accuracies() { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + let bam = root.join("all-no-calls.bam"); + let bed = root.join("all-no-calls.bed"); + let report = root.join("all-no-calls.tsv"); + + let records = (0..3) + .map(|i| { + synthetic_record( + &format!("no-call-{i}"), + "CC", + vec![Cigar::Match(2)], + "C+m?;", + &[], + 0, + ) + }) + .collect(); + write_synthetic_bam(&bam, records); + File::create(&bed).unwrap().write_all(b"chr1\t1\t2\tm\t0\t+\n").unwrap(); + + let output = std::process::Command::new(env!("CARGO_BIN_EXE_modkit")) + .args([ + "validate", + "--bam-and-bed", + bam.to_str().unwrap(), + bed.to_str().unwrap(), + "--canonical-base", + "C", + "--filter-threshold", + "0", + "--threads", + "1", + "--suppress-progress", + "--out-filepath", + report.to_str().unwrap(), + ]) + .output() + .unwrap(); + assert!( + output.status.success(), + "validate failed: {}", + String::from_utf8_lossy(&output.stderr) + ); + + let observed = std::fs::read_to_string(report).unwrap(); + assert_eq!( + observed, + concat!( + "full_contingency_table: [[\"ground_truth_label\",\"m\",\"No Call\"],[\"m\",0,3]]\n", + "raw_accuracy: NA\n", + "raw_contingency_table: [[\"ground_truth_label\"]]\n", + "filter_threshold: 0\n", + "percent_of_mod_called_removed: NA\n", + "filtered_accuracy: NA\n", + "filtered_contingency_table: [[\"ground_truth_label\"]]\n", + ) + ); +} + +#[test] +fn test_validate_all_calls_filtered_reports_unavailable_filtered_accuracy() { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + let bam = root.join("all-filtered.bam"); + let bed = root.join("all-filtered.bed"); + let report = root.join("all-filtered.tsv"); + + write_synthetic_bam(&bam, vec![synthetic_called("called", '+')]); + File::create(&bed).unwrap().write_all(b"chr1\t1\t2\tm\t0\t+\n").unwrap(); + + let output = std::process::Command::new(env!("CARGO_BIN_EXE_modkit")) + .args([ + "validate", + "--bam-and-bed", + bam.to_str().unwrap(), + bed.to_str().unwrap(), + "--canonical-base", + "C", + "--filter-threshold", + "1", + "--threads", + "1", + "--suppress-progress", + "--out-filepath", + report.to_str().unwrap(), + ]) + .output() + .unwrap(); + assert!( + output.status.success(), + "validate failed: {}", + String::from_utf8_lossy(&output.stderr) + ); + + let observed = std::fs::read_to_string(report).unwrap(); + assert_eq!( + observed, + concat!( + "full_contingency_table: [[\"ground_truth_label\",\"m\"],[\"m\",1]]\n", + "raw_accuracy: 100\n", + "raw_contingency_table: [[\"ground_truth_label\",\"m\"],[\"m\",1]]\n", + "filter_threshold: 1\n", + "percent_of_mod_called_removed: 100\n", + "filtered_accuracy: NA\n", + "filtered_contingency_table: [[\"ground_truth_label\",\"m\"],[\"m\",0]]\n", + ) + ); +} From 03bd539ed62244855d216dbbe6f3bb77161e54a7 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 06:56:50 -0700 Subject: [PATCH 08/11] Render undefined validate percentages as unavailable --- modkit-core/src/validate/subcommand.rs | 21 ++++++++++++++++----- modkit/tests/test_validate.rs | 7 ++++--- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/modkit-core/src/validate/subcommand.rs b/modkit-core/src/validate/subcommand.rs index 3dae2e32..57d4bc71 100644 --- a/modkit-core/src/validate/subcommand.rs +++ b/modkit-core/src/validate/subcommand.rs @@ -776,6 +776,14 @@ fn machine_parseable_table( out_str } +fn format_percentage_cell(count: usize, total: usize) -> String { + if total == 0 { + "NA".to_string() + } else { + format!("{:.2}%", 100.0 * count as f32 / total as f32) + } +} + fn print_table( validate_base: DnaBase, status_probs: &StatusProbs, @@ -825,11 +833,8 @@ fn print_table( if show_percentages { let gt_total = gt_totals.get(gt_code).unwrap(); row.add_cell( - cell!(&format!( - "{:.2}%", - 100.0 * vector_length as f32 / *gt_total as f32 - )) - .style_spec("r"), + cell!(&format_percentage_cell(vector_length, *gt_total)) + .style_spec("r"), ); } else { row.add_cell( @@ -1268,6 +1273,12 @@ mod tests { use super::*; + #[test] + fn percentage_cells_are_unavailable_without_a_denominator() { + assert_eq!(format_percentage_cell(0, 0), "NA"); + assert_eq!(format_percentage_cell(1, 4), "25.00%"); + } + fn make_record(gap: Cigar) -> Record { let cigar = CigarString(vec![Cigar::Match(1), gap, Cigar::Match(1)]); let mut record = Record::new(); diff --git a/modkit/tests/test_validate.rs b/modkit/tests/test_validate.rs index b1a975a5..b6a0d899 100644 --- a/modkit/tests/test_validate.rs +++ b/modkit/tests/test_validate.rs @@ -498,10 +498,11 @@ fn test_validate_all_calls_filtered_reports_unavailable_filtered_accuracy() { ]) .output() .unwrap(); + let stderr = String::from_utf8_lossy(&output.stderr); + assert!(output.status.success(), "validate failed: {}", stderr); assert!( - output.status.success(), - "validate failed: {}", - String::from_utf8_lossy(&output.stderr) + !stderr.contains("NaN%"), + "zero-total table cell must be NA:\n{stderr}" ); let observed = std::fs::read_to_string(report).unwrap(); From 384ce82738818837852cd971f2eeae0780621217 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 06:57:49 -0700 Subject: [PATCH 09/11] Strengthen reverse edge-filter validation oracle --- .../validate/unseeded_observations_tests.rs | 29 +++++++++++-------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/modkit-core/src/validate/unseeded_observations_tests.rs b/modkit-core/src/validate/unseeded_observations_tests.rs index a348854f..8e9d61bc 100644 --- a/modkit-core/src/validate/unseeded_observations_tests.rs +++ b/modkit-core/src/validate/unseeded_observations_tests.rs @@ -435,23 +435,28 @@ fn inverted_edge_filter_excludes_filtered_fallback_sites() { fn reverse_alignment_edge_filter_uses_forward_query_coordinates() { let mut record = make_record( "reverse-edge-filter", - "CCCCCC", + "CACNCC", vec![Cigar::Match(6)], - "G-m?,0,0,0,0,0,0;", - &[255; 6], + "G-m?;", + &[], ); record.set_reverse(); let truth = positive_truth_range(6); - for inverted in [false, true] { - let edge_filter = EdgeFilter::new(1, 2, inverted); - let observed = classify_with_edge_filter( - [record.clone()], - &truth, - Some(&edge_filter), - ); - assert_exact_counts(&observed, &[((M, M), 3)]); - } + let ordinary = EdgeFilter::new(1, 2, false); + let observed = + classify_with_edge_filter([record.clone()], &truth, Some(&ordinary)); + assert_exact_counts(&observed, &[((M, BaseStatus::NoCall), 3)]); + + let inverted = EdgeFilter::new(1, 2, true); + let observed = classify_with_edge_filter([record], &truth, Some(&inverted)); + assert_exact_counts( + &observed, + &[ + ((M, BaseStatus::NoCall), 2), + ((M, BaseStatus::Mismatch(DnaBase::A)), 1), + ], + ); } #[test] From 93195bf0ac9dc317a2753db13c08f145e6b14e63 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 10/11] Return validate truth BED errors --- book/src/intro_validate.md | 6 +- modkit-core/src/validate/subcommand.rs | 253 +++++++++++++++++++++---- modkit/tests/test_validate.rs | 54 ++++++ 3 files changed, 280 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 57d4bc71..7116e73c 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; @@ -84,18 +84,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" + )), } } } @@ -140,13 +152,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 { @@ -196,16 +216,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) @@ -217,7 +250,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()); @@ -966,10 +1002,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)?; @@ -1031,6 +1063,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 { @@ -1269,7 +1303,10 @@ impl ValidateFromModBam { #[cfg(test)] mod tests { + use std::io::Write; + use rust_htslib::bam::record::{Aux, Cigar, CigarString}; + use tempfile::NamedTempFile; use super::*; @@ -1390,4 +1427,156 @@ mod tests { ); assert_eq!(classified.values().map(Vec::len).sum::(), 5); } + + #[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 b6a0d899..c0ec0156 100644 --- a/modkit/tests/test_validate.rs +++ b/modkit/tests/test_validate.rs @@ -7,6 +7,7 @@ use std::collections::{HashMap, HashSet}; use std::fs::{self, File}; use std::io::{BufRead, BufReader, Write}; use std::path::Path; +use std::process::Command; use tempfile::tempdir; mod common; @@ -519,3 +520,56 @@ fn test_validate_all_calls_filtered_reports_unavailable_filtered_accuracy() { ) ); } + +#[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 e90e6ff8a2cdc5196d35a5364c4c992c6b37da91 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 11/11] 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 7116e73c..ffd0b806 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; @@ -138,7 +140,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 { @@ -179,7 +181,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 }) } @@ -189,6 +191,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()) @@ -209,6 +222,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 @@ -239,14 +253,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() { @@ -1562,6 +1623,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 c0ec0156..74e26172 100644 --- a/modkit/tests/test_validate.rs +++ b/modkit/tests/test_validate.rs @@ -533,6 +533,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 {