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/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 119aa3ed..ffd0b806 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; @@ -38,6 +40,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 { @@ -80,18 +86,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 +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 { @@ -136,13 +154,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 +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 }) } @@ -165,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()) @@ -185,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 @@ -192,28 +230,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()); @@ -276,6 +377,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, @@ -297,6 +407,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 +451,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; @@ -371,6 +493,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 @@ -382,16 +505,51 @@ 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) .push(f32::NAN); continue; }; - let mut base = DnaBase::parse(q_seq[*q_pos as usize] as char)?; + 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) + .push(f32::NAN); + continue; + }; 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)) @@ -690,11 +848,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)); @@ -713,6 +873,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, @@ -762,11 +930,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( @@ -898,10 +1063,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,18 +1124,20 @@ 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 { - 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, @@ -1015,6 +1178,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"); @@ -1085,13 +1300,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, @@ -1115,7 +1343,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 @@ -1132,3 +1361,358 @@ impl ValidateFromModBam { Ok(()) } } + +#[cfg(test)] +mod tests { + use std::io::Write; + + use rust_htslib::bam::record::{Aux, Cigar, CigarString}; + use tempfile::NamedTempFile; + + 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(); + 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 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, + call: BaseStatus, + ) -> usize { + status_probs.get(&(truth, call)).map(Vec::len).unwrap_or(0) + } + + fn classify_record(record: &Record) -> anyhow::Result { + 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( + record, + &truth, + &tid_to_chrom, + DnaBase::C, + None, + None, + ) + } + + fn classify(gap: Cigar) -> StatusProbs { + classify_record(&make_record(gap)).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); + } + + #[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); + } + + #[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-core/src/validate/unseeded_observations_tests.rs b/modkit-core/src/validate/unseeded_observations_tests.rs new file mode 100644 index 00000000..8e9d61bc --- /dev/null +++ b/modkit-core/src/validate/unseeded_observations_tests.rs @@ -0,0 +1,481 @@ +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_bam::EdgeFilter; +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 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, +) -> 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(); + for record in records { + let observed = process_bam_record( + &record, + truth, + &tid_to_chrom, + DnaBase::C, + None, + edge_filter, + ) + .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()); +} + +#[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), + ], + ); +} + +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", + "CACNCC", + vec![Cigar::Match(6)], + "G-m?;", + &[], + ); + record.set_reverse(); + let truth = positive_truth_range(6); + + 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] +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)]); +} diff --git a/modkit/tests/test_validate.rs b/modkit/tests/test_validate.rs index aba45a9d..74e26172 100644 --- a/modkit/tests/test_validate.rs +++ b/modkit/tests/test_validate.rs @@ -1,10 +1,50 @@ use crate::common::run_modkit; use anyhow::Context; -use std::fs::File; -use std::io::{BufRead, BufReader}; +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, Write}; +use std::path::Path; +use std::process::Command; +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 +92,489 @@ 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); +} + +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]]", + ); +} + +#[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(); + let stderr = String::from_utf8_lossy(&output.stderr); + assert!(output.status.success(), "validate failed: {}", stderr); + assert!( + !stderr.contains("NaN%"), + "zero-total table cell must be NA:\n{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", + ) + ); +} + +#[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}" + ); + } +}