From e6774349e133f545d14067dd0706163ded1cbe90 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Fri, 26 Jun 2026 08:00:10 -0700 Subject: [PATCH 1/9] refactor(eval): move eval.rs into eval/ module directory --- src/commands/{eval.rs => eval/mod.rs} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/commands/{eval.rs => eval/mod.rs} (100%) diff --git a/src/commands/eval.rs b/src/commands/eval/mod.rs similarity index 100% rename from src/commands/eval.rs rename to src/commands/eval/mod.rs From 2e46bb928d95b12b8feda0f06535b79066fe9491 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Fri, 26 Jun 2026 08:03:22 -0700 Subject: [PATCH 2/9] refactor(eval): extract placement scoring into a submodule The eval command is about to grow optional truth-aware metrics (variant representation, methylation-level correlation, golden-BAM placement), so the existing MAPQ-binned placement scoring moves out of the command struct into a focused `placement` submodule. This leaves `eval/mod.rs` as a thin orchestrator over the shared `Eval` options and gives each future metric its own file. Pure code movement: the placement loop, MAPQ binning, and TSV row formatting are unchanged, and `placement::run` reproduces the previous `execute` behavior. The only edit is renaming the final log line to "Placement results written to" so it reads correctly once sibling metrics write their own tables. --- src/commands/eval/mod.rs | 220 ++---------------------------- src/commands/eval/placement.rs | 235 +++++++++++++++++++++++++++++++++ 2 files changed, 246 insertions(+), 209 deletions(-) create mode 100644 src/commands/eval/placement.rs diff --git a/src/commands/eval/mod.rs b/src/commands/eval/mod.rs index 97dd742..10dd7b8 100644 --- a/src/commands/eval/mod.rs +++ b/src/commands/eval/mod.rs @@ -1,17 +1,20 @@ //! Alignment accuracy evaluation command. +//! +//! `holodeck eval` scores an aligner's BAM against holodeck's own truth. +//! Placement accuracy (true vs mapped position, by MAPQ bin) is always +//! reported from encoded read names. Optional truth inputs unlock further +//! metrics, each written to its own TSV alongside `.eval.txt`: +//! - [`placement`] — placement accuracy (always). + +mod placement; -use std::collections::BTreeMap; -use std::io::Write; use std::path::PathBuf; -use anyhow::{Context, Result}; -use bstr::ByteSlice; +use anyhow::Result; use clap::Parser; -use noodles::bam; -use super::command::{Command, output_path}; +use super::command::Command; use super::common::OutputPrefixOptions; -use crate::read_naming::{parse_encoded_pe_name, parse_encoded_se_name}; /// Evaluate alignment accuracy of simulated reads. /// @@ -43,213 +46,12 @@ pub struct Eval { pub wiggle: u32, } -/// Accuracy counts for a single MAPQ bin. -#[derive(Debug, Default, Clone)] -struct BinCounts { - /// Reads mapped to the correct position (within wiggle). - correct: u64, - /// Reads mapped to a wrong position or wrong contig. - mismapped: u64, - /// Reads that are unmapped. - unmapped: u64, - /// Total reads in this bin. - total: u64, -} - impl Command for Eval { fn execute(&self) -> Result<()> { if self.truth.is_some() { log::warn!("--truth (golden BAM) is not yet implemented; using read names"); } - let mut reader = bam::io::reader::Builder - .build_from_path(&self.mapped) - .with_context(|| format!("Failed to open BAM: {}", self.mapped.display()))?; - - let header = reader.read_header()?; - - // MAPQ bins: 0, 1-9, 10-19, 20-29, 30-39, 40-49, 50-59, 60+ - let mut bins: BTreeMap = BTreeMap::new(); - let mut total_reads: u64 = 0; - let mut parse_failures: u64 = 0; - - for result in reader.records() { - let record = result.with_context(|| "Failed to read BAM record")?; - - // Skip secondary and supplementary alignments before counting. - let flags = record.flags(); - if flags.is_secondary() || flags.is_supplementary() { - continue; - } - total_reads += 1; - - // Get read name. - let name_bytes = record.name().map_or(&b""[..], |n| n.as_bytes()); - let name = name_bytes.to_str().unwrap_or(""); - - // Parse truth from encoded read name. For PE names, pick R1 or R2 - // based on the record's segment flag; mis-selecting here caused R2 - // alignments to be scored against the R1 truth position. - let truth = if let Some((_, r1, r2)) = parse_encoded_pe_name(name) { - if flags.is_last_segment() { Some(r2) } else { Some(r1) } - } else { - parse_encoded_se_name(name).map(|(_, truth)| truth) - }; - - let Some(truth) = truth else { - parse_failures += 1; - continue; - }; - - let mapq = record.mapping_quality().map_or(0, u8::from); - let bin_key = mapq_bin(mapq); - let counts = bins.entry(bin_key).or_default(); - counts.total += 1; - - if flags.is_unmapped() { - counts.unmapped += 1; - continue; - } - - // Get mapped position. - let mapped_contig_idx = record.reference_sequence_id().and_then(Result::ok); - let mapped_pos = record.alignment_start().and_then(Result::ok).map(usize::from); - - let (Some(mapped_contig_idx), Some(mapped_pos_1based)) = - (mapped_contig_idx, mapped_pos) - else { - counts.mismapped += 1; - continue; - }; - - // Resolve mapped contig name. - let Some((contig_name, _)) = header.reference_sequences().get_index(mapped_contig_idx) - else { - counts.mismapped += 1; - continue; - }; - let mapped_contig = String::from_utf8_lossy(contig_name.as_ref()); - - // Compare truth vs mapped. - #[expect(clippy::cast_possible_truncation, reason = "mapped positions fit u32")] - let mapped_pos_u32 = mapped_pos_1based as u32; - let is_correct = mapped_contig == truth.contig - && mapped_pos_u32.abs_diff(truth.position) <= self.wiggle; - - if is_correct { - counts.correct += 1; - } else { - counts.mismapped += 1; - } - } - - if parse_failures > 0 { - log::warn!("{parse_failures} reads had unparseable names; skipped"); - } - - // Write results. - let output_file = output_path(&self.output.output, ".eval.txt"); - let mut out = std::fs::File::create(&output_file) - .with_context(|| format!("Failed to create {}", output_file.display()))?; - - writeln!( - out, - "mapq_bin\ttotal\tcorrect\tmismapped\tunmapped\tpct_correct\tpct_mismapped\tpct_unmapped" - )?; - - let mut grand_total = BinCounts::default(); - for (&bin, counts) in &bins { - write_bin_row(&mut out, &format_bin_label(bin), counts)?; - grand_total.correct += counts.correct; - grand_total.mismapped += counts.mismapped; - grand_total.unmapped += counts.unmapped; - grand_total.total += counts.total; - } - write_bin_row(&mut out, "ALL", &grand_total)?; - - log::info!( - "Evaluated {total_reads} reads: {} correct, {} mismapped, {} unmapped", - grand_total.correct, - grand_total.mismapped, - grand_total.unmapped - ); - log::info!("Results written to: {}", output_file.display()); - - Ok(()) - } -} - -/// Map a MAPQ value to a bin key. -fn mapq_bin(mapq: u8) -> u8 { - match mapq { - 0 => 0, - 1..=9 => 1, - 10..=19 => 10, - 20..=29 => 20, - 30..=39 => 30, - 40..=49 => 40, - 50..=59 => 50, - _ => 60, - } -} - -/// Format a bin key as a label string. -fn format_bin_label(bin: u8) -> String { - match bin { - 0 => "0".to_string(), - 1 => "1-9".to_string(), - 60 => "60+".to_string(), - _ => format!("{}-{}", bin, bin + 9), - } -} - -/// Write one row of the evaluation results table. -fn write_bin_row(out: &mut impl Write, label: &str, counts: &BinCounts) -> Result<()> { - let total = counts.total.max(1) as f64; - writeln!( - out, - "{label}\t{}\t{}\t{}\t{}\t{:.2}\t{:.2}\t{:.2}", - counts.total, - counts.correct, - counts.mismapped, - counts.unmapped, - counts.correct as f64 / total * 100.0, - counts.mismapped as f64 / total * 100.0, - counts.unmapped as f64 / total * 100.0, - )?; - Ok(()) -} - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - fn test_mapq_bin() { - assert_eq!(mapq_bin(0), 0); - assert_eq!(mapq_bin(5), 1); - assert_eq!(mapq_bin(10), 10); - assert_eq!(mapq_bin(15), 10); - assert_eq!(mapq_bin(30), 30); - assert_eq!(mapq_bin(60), 60); - assert_eq!(mapq_bin(255), 60); - } - - #[test] - fn test_format_bin_label() { - assert_eq!(format_bin_label(0), "0"); - assert_eq!(format_bin_label(1), "1-9"); - assert_eq!(format_bin_label(10), "10-19"); - assert_eq!(format_bin_label(60), "60+"); - } - - #[test] - fn test_write_bin_row() { - let counts = BinCounts { correct: 90, mismapped: 8, unmapped: 2, total: 100 }; - let mut buf = Vec::new(); - write_bin_row(&mut buf, "30-39", &counts).unwrap(); - let line = String::from_utf8(buf).unwrap(); - assert!(line.starts_with("30-39\t100\t90\t8\t2\t")); - assert!(line.contains("90.00")); + placement::run(&self.mapped, &self.output.output, self.wiggle) } } diff --git a/src/commands/eval/placement.rs b/src/commands/eval/placement.rs new file mode 100644 index 0000000..829d898 --- /dev/null +++ b/src/commands/eval/placement.rs @@ -0,0 +1,235 @@ +//! Placement accuracy: true vs mapped position, stratified by MAPQ bin. +//! +//! Truth positions are parsed from encoded holodeck read names. For each +//! primary mapped record, the mapped start is compared to the true start on +//! the same contig within a wiggle tolerance; reads are tallied as correct, +//! mismapped, or unmapped within their MAPQ bin. + +use std::collections::BTreeMap; +use std::io::Write; +use std::path::Path; + +use anyhow::{Context, Result}; +use bstr::ByteSlice; +use noodles::bam; + +use crate::commands::command::output_path; +use crate::read_naming::{parse_encoded_pe_name, parse_encoded_se_name}; + +/// Accuracy counts for a single MAPQ bin. +#[derive(Debug, Default, Clone)] +struct BinCounts { + /// Reads mapped to the correct position (within wiggle). + correct: u64, + /// Reads mapped to a wrong position or wrong contig. + mismapped: u64, + /// Reads that are unmapped. + unmapped: u64, + /// Total reads in this bin. + total: u64, +} + +/// Evaluate placement accuracy of `mapped` and write `.eval.txt`. +/// +/// # Errors +/// Returns an error if the BAM cannot be read or the output cannot be written. +pub fn run(mapped: &Path, output_prefix: &Path, wiggle: u32) -> Result<()> { + let mut reader = bam::io::reader::Builder + .build_from_path(mapped) + .with_context(|| format!("Failed to open BAM: {}", mapped.display()))?; + + let header = reader.read_header()?; + + // MAPQ bins: 0, 1-9, 10-19, 20-29, 30-39, 40-49, 50-59, 60+, NA (255) + let mut bins: BTreeMap = BTreeMap::new(); + let mut total_reads: u64 = 0; + let mut parse_failures: u64 = 0; + + for result in reader.records() { + let record = result.with_context(|| "Failed to read BAM record")?; + + // Skip secondary and supplementary alignments before counting. + let flags = record.flags(); + if flags.is_secondary() || flags.is_supplementary() { + continue; + } + total_reads += 1; + + // Get read name. + let name_bytes = record.name().map_or(&b""[..], |n| n.as_bytes()); + let name = name_bytes.to_str().unwrap_or(""); + + // Parse truth from encoded read name. For PE names, pick R1 or R2 + // based on the record's segment flag; mis-selecting here caused R2 + // alignments to be scored against the R1 truth position. + let truth = if let Some((_, r1, r2)) = parse_encoded_pe_name(name) { + if flags.is_last_segment() { Some(r2) } else { Some(r1) } + } else { + parse_encoded_se_name(name).map(|(_, truth)| truth) + }; + + let Some(truth) = truth else { + parse_failures += 1; + continue; + }; + + let mapq = record.mapping_quality().map_or(0, u8::from); + let bin_key = mapq_bin(mapq); + let counts = bins.entry(bin_key).or_default(); + counts.total += 1; + + if flags.is_unmapped() { + counts.unmapped += 1; + continue; + } + + // Get mapped position. + let mapped_contig_idx = record.reference_sequence_id().and_then(Result::ok); + let mapped_pos = record.alignment_start().and_then(Result::ok).map(usize::from); + + let (Some(mapped_contig_idx), Some(mapped_pos_1based)) = (mapped_contig_idx, mapped_pos) + else { + counts.mismapped += 1; + continue; + }; + + // Resolve mapped contig name. + let Some((contig_name, _)) = header.reference_sequences().get_index(mapped_contig_idx) + else { + counts.mismapped += 1; + continue; + }; + let mapped_contig = String::from_utf8_lossy(contig_name.as_ref()); + + // Compare truth vs mapped. + #[expect(clippy::cast_possible_truncation, reason = "mapped positions fit u32")] + let mapped_pos_u32 = mapped_pos_1based as u32; + let is_correct = + mapped_contig == truth.contig && mapped_pos_u32.abs_diff(truth.position) <= wiggle; + + if is_correct { + counts.correct += 1; + } else { + counts.mismapped += 1; + } + } + + if parse_failures > 0 { + log::warn!("{parse_failures} reads had unparseable names; skipped"); + } + + // Write results. + let output_file = output_path(output_prefix, ".eval.txt"); + let mut out = std::fs::File::create(&output_file) + .with_context(|| format!("Failed to create {}", output_file.display()))?; + + writeln!( + out, + "mapq_bin\ttotal\tcorrect\tmismapped\tunmapped\tpct_correct\tpct_mismapped\tpct_unmapped" + )?; + + let mut grand_total = BinCounts::default(); + for (&bin, counts) in &bins { + write_bin_row(&mut out, &format_bin_label(bin), counts)?; + grand_total.correct += counts.correct; + grand_total.mismapped += counts.mismapped; + grand_total.unmapped += counts.unmapped; + grand_total.total += counts.total; + } + write_bin_row(&mut out, "ALL", &grand_total)?; + + log::info!( + "Evaluated {total_reads} reads: {} correct, {} mismapped, {} unmapped", + grand_total.correct, + grand_total.mismapped, + grand_total.unmapped + ); + log::info!("Placement results written to: {}", output_file.display()); + + Ok(()) +} + +/// Map a MAPQ value to a bin key. MAPQ `255` is reserved by the SAM spec for +/// "mapping quality unavailable" and is binned separately so unknown-quality +/// reads are not counted with the highest-confidence `60+` alignments. +fn mapq_bin(mapq: u8) -> u8 { + match mapq { + 0 => 0, + 1..=9 => 1, + 10..=19 => 10, + 20..=29 => 20, + 30..=39 => 30, + 40..=49 => 40, + 50..=59 => 50, + 255 => 255, + _ => 60, + } +} + +/// Format a bin key as a label string. The `255` key is the SAM "unavailable" +/// MAPQ and is labelled `NA`. +fn format_bin_label(bin: u8) -> String { + match bin { + 0 => "0".to_string(), + 1 => "1-9".to_string(), + 60 => "60+".to_string(), + 255 => "NA".to_string(), + _ => format!("{}-{}", bin, bin + 9), + } +} + +/// Write one row of the evaluation results table. +fn write_bin_row(out: &mut impl Write, label: &str, counts: &BinCounts) -> Result<()> { + let total = counts.total.max(1) as f64; + writeln!( + out, + "{label}\t{}\t{}\t{}\t{}\t{:.2}\t{:.2}\t{:.2}", + counts.total, + counts.correct, + counts.mismapped, + counts.unmapped, + counts.correct as f64 / total * 100.0, + counts.mismapped as f64 / total * 100.0, + counts.unmapped as f64 / total * 100.0, + )?; + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_mapq_bin() { + assert_eq!(mapq_bin(0), 0); + assert_eq!(mapq_bin(5), 1); + assert_eq!(mapq_bin(10), 10); + assert_eq!(mapq_bin(15), 10); + assert_eq!(mapq_bin(30), 30); + assert_eq!(mapq_bin(60), 60); + // 61..=254 are valid high-confidence MAPQs and stay in the 60+ bin. + assert_eq!(mapq_bin(254), 60); + // 255 is reserved by the SAM spec for "mapping quality unavailable" and + // must not be counted as the highest-confidence bin. + assert_eq!(mapq_bin(255), 255); + } + + #[test] + fn test_format_bin_label() { + assert_eq!(format_bin_label(0), "0"); + assert_eq!(format_bin_label(1), "1-9"); + assert_eq!(format_bin_label(10), "10-19"); + assert_eq!(format_bin_label(60), "60+"); + assert_eq!(format_bin_label(255), "NA"); + } + + #[test] + fn test_write_bin_row() { + let counts = BinCounts { correct: 90, mismapped: 8, unmapped: 2, total: 100 }; + let mut buf = Vec::new(); + write_bin_row(&mut buf, "30-39", &counts).unwrap(); + let line = String::from_utf8(buf).unwrap(); + assert!(line.starts_with("30-39\t100\t90\t8\t2\t")); + assert!(line.contains("90.00")); + } +} From 12f1f27abd4707a59fde5a0e7a4cd71a893a5981 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Fri, 26 Jun 2026 08:14:13 -0700 Subject: [PATCH 3/9] feat(eval): score variant representation against golden truth `holodeck eval --variants truth.vcf --truth golden.bam` adds an accuracy axis beyond placement: for every simulated substitution a read should carry, does the aligned read actually represent the alternate base, and with what MAPQ and alignment score? Truth is taken entirely from holodeck's own outputs. The golden BAM supplies each read's true span and source haplotype (hp:i); the truth VCF's phased genotypes say which single-base substitutions that haplotype carries within the span. For each expected substitution the pass walks the *mapped* read's CIGAR to the variant's reference position and compares the observed base to the alternate allele. Reads mapped to the wrong locus simply fail to represent their variants, so mismapping is captured without special-casing. Per-read MD:Z / NM:i tags are compared against the golden tags to report tag-level concordance over variant-bearing reads. Results land in `.variants.tsv`. With `--meth`, results break down by substitution class relative to the read's bisulfite conversion direction (XG, falling back to XR). The C->T cell on a CT-strand read (G->A on GA) is intrinsically confounded with conversion and is labelled as such rather than scored as a real signal; the discriminating classes are the mirror (T->C / A->G) and the transversions. This is the axis on which a methylation-aware scoring mode should match the genomic truth without over- or under-penalizing. Adds eval submodules cigar (CIGAR geometry, unit-tested), golden (truth index), and variants (classifier, per-haplotype expected SNVs, the scoring pass). The classifier and expected-SNV queries are covered by unit tests with programmatically built records; the end-to-end pass is exercised by a later integration test. --- src/commands/eval/cigar.rs | 149 ++++++++++ src/commands/eval/golden.rs | 107 +++++++ src/commands/eval/mod.rs | 62 +++- src/commands/eval/variants.rs | 536 ++++++++++++++++++++++++++++++++++ 4 files changed, 842 insertions(+), 12 deletions(-) create mode 100644 src/commands/eval/cigar.rs create mode 100644 src/commands/eval/golden.rs create mode 100644 src/commands/eval/variants.rs diff --git a/src/commands/eval/cigar.rs b/src/commands/eval/cigar.rs new file mode 100644 index 0000000..0bb81d7 --- /dev/null +++ b/src/commands/eval/cigar.rs @@ -0,0 +1,149 @@ +//! CIGAR geometry helpers shared by the truth-aware eval metrics. +//! +//! Both operate on a [`Cigar`] paired with the alignment's 0-based reference +//! start. [`reference_len`] gives the reference span consumed; +//! [`ref_pos_to_read_offset`] maps a reference position to the read offset of +//! the base aligned there, or `None` when that position is deleted, skipped, +//! or otherwise not covered by an aligned (`M`/`=`/`X`) operation. + +use noodles::sam::alignment::record::cigar::op::Kind; +use noodles::sam::alignment::record_buf::Cigar; + +/// Whether a CIGAR operation consumes reference bases. +fn consumes_reference(kind: Kind) -> bool { + matches!( + kind, + Kind::Match | Kind::Deletion | Kind::Skip | Kind::SequenceMatch | Kind::SequenceMismatch + ) +} + +/// Whether a CIGAR operation consumes query (read) bases. +fn consumes_query(kind: Kind) -> bool { + matches!( + kind, + Kind::Match + | Kind::Insertion + | Kind::SoftClip + | Kind::SequenceMatch + | Kind::SequenceMismatch + ) +} + +/// Number of reference bases the alignment spans. +#[must_use] +pub fn reference_len(cigar: &Cigar) -> u32 { + let mut len: u32 = 0; + for op in cigar.as_ref() { + if consumes_reference(op.kind()) { + len += u32::try_from(op.len()).unwrap_or(0); + } + } + len +} + +/// Map a 0-based reference position to the 0-based read offset aligned there. +/// +/// `aln_start0` is the alignment's 0-based reference start. Returns `None` when +/// `target_ref0` falls outside the alignment, inside a deletion/skip, or is +/// otherwise not covered by an aligned base. +#[must_use] +pub fn ref_pos_to_read_offset(cigar: &Cigar, aln_start0: u32, target_ref0: u32) -> Option { + if target_ref0 < aln_start0 { + return None; + } + let mut ref_pos = aln_start0; + let mut read_pos: usize = 0; + for op in cigar.as_ref() { + let kind = op.kind(); + let span = u32::try_from(op.len()).unwrap_or(0); + let consumes_ref = consumes_reference(kind); + let consumes_q = consumes_query(kind); + + if consumes_ref && consumes_q { + // Aligned run: target may land within it. + if target_ref0 < ref_pos + span { + let within = (target_ref0 - ref_pos) as usize; + return Some(read_pos + within); + } + ref_pos += span; + read_pos += op.len(); + } else if consumes_ref { + // Deletion/skip: target inside this run is not represented by a base. + if target_ref0 < ref_pos + span { + return None; + } + ref_pos += span; + } else if consumes_q { + // Insertion/soft-clip: advances the read only. + read_pos += op.len(); + } + } + None +} + +#[cfg(test)] +mod tests { + use super::*; + use noodles::sam::alignment::record::cigar::op::Op; + + fn cigar(ops: &[(Kind, usize)]) -> Cigar { + Cigar::from(ops.iter().map(|&(k, n)| Op::new(k, n)).collect::>()) + } + + #[test] + fn reference_len_sums_ref_consuming_ops() { + // 10M2I5M3D4M -> reference span 10 + 5 + 3 + 4 = 22 (insertion excluded). + let c = cigar(&[ + (Kind::Match, 10), + (Kind::Insertion, 2), + (Kind::Match, 5), + (Kind::Deletion, 3), + (Kind::Match, 4), + ]); + assert_eq!(reference_len(&c), 22); + } + + #[test] + fn ref_offset_simple_match() { + // 100M starting at ref 1000: ref 1005 -> read offset 5. + let c = cigar(&[(Kind::Match, 100)]); + assert_eq!(ref_pos_to_read_offset(&c, 1000, 1005), Some(5)); + assert_eq!(ref_pos_to_read_offset(&c, 1000, 1000), Some(0)); + } + + #[test] + fn ref_offset_after_insertion_shifts_read() { + // 5M2I5M at ref 1000: ref 1006 is in the second M; read offset = 5 (M) + 2 (I) + 1. + let c = cigar(&[(Kind::Match, 5), (Kind::Insertion, 2), (Kind::Match, 5)]); + assert_eq!(ref_pos_to_read_offset(&c, 1000, 1006), Some(8)); + } + + #[test] + fn ref_offset_after_deletion_shifts_ref() { + // 5M3D5M at ref 1000: ref 1008 is the first base after the deletion; + // read offset = 5 (no read bases consumed by D). + let c = cigar(&[(Kind::Match, 5), (Kind::Deletion, 3), (Kind::Match, 5)]); + assert_eq!(ref_pos_to_read_offset(&c, 1000, 1008), Some(5)); + } + + #[test] + fn ref_offset_inside_deletion_is_none() { + let c = cigar(&[(Kind::Match, 5), (Kind::Deletion, 3), (Kind::Match, 5)]); + // ref 1006 is inside the 3bp deletion (1005..1008). + assert_eq!(ref_pos_to_read_offset(&c, 1000, 1006), None); + } + + #[test] + fn ref_offset_softclip_offsets_read() { + // 4S10M at ref 1000: ref 1000 -> read offset 4 (soft-clip consumes read). + let c = cigar(&[(Kind::SoftClip, 4), (Kind::Match, 10)]); + assert_eq!(ref_pos_to_read_offset(&c, 1000, 1000), Some(4)); + } + + #[test] + fn ref_offset_out_of_span_is_none() { + let c = cigar(&[(Kind::Match, 10)]); + assert_eq!(ref_pos_to_read_offset(&c, 1000, 1010), None); // 1000..1010 only + assert_eq!(ref_pos_to_read_offset(&c, 1000, 999), None); + } +} diff --git a/src/commands/eval/golden.rs b/src/commands/eval/golden.rs new file mode 100644 index 0000000..58b764c --- /dev/null +++ b/src/commands/eval/golden.rs @@ -0,0 +1,107 @@ +//! Golden-BAM truth: per-read true alignment loaded for placement and for +//! variant / MD-tag concordance scoring. +//! +//! The golden BAM written by `holodeck simulate --golden-bam` carries the +//! true alignment of every read (MAPQ 60, correct CIGAR), the source +//! haplotype in the `hp:i` tag, and — for methylation runs — Bismark-style +//! `NM:i` / `MD:Z` call tags. This module indexes those records by read end so +//! the eval pass can look up each mapped read's truth in O(1). + +use std::collections::HashMap; +use std::path::Path; + +use anyhow::{Context, Result}; +use bstr::ByteSlice; +use noodles::bam; +use noodles::sam::alignment::record::data::field::Tag; +use noodles::sam::alignment::record_buf::data::field::Value; + +use super::cigar; + +/// Key identifying one read end: read name plus whether it is the last segment +/// (R2). R1 and single-end reads use `false`. +pub type ReadKey = (Vec, bool); + +/// True alignment for one read end, taken from the golden BAM. +#[derive(Debug, Clone)] +pub struct GoldenInfo { + /// True reference contig name. + pub contig: String, + /// 0-based true start position. + pub start0: u32, + /// Reference bases consumed by the true alignment. + pub ref_len: u32, + /// Source haplotype index (`hp:i` tag; `0` if absent). + pub haplotype: usize, + /// `NM:i` edit distance, if present. + pub nm: Option, + /// `MD:Z` string, if present. + pub md: Option, +} + +impl GoldenInfo { + /// 0-based exclusive end of the true alignment. + #[must_use] + pub fn end0(&self) -> u32 { + self.start0 + self.ref_len + } +} + +/// Load every primary, mapped golden record keyed by `(name, is_last_segment)`. +/// +/// # Errors +/// Returns an error if the golden BAM cannot be read. +pub fn load(path: &Path) -> Result> { + let mut reader = bam::io::reader::Builder + .build_from_path(path) + .with_context(|| format!("Failed to open golden BAM: {}", path.display()))?; + let header = reader.read_header()?; + + let mut map = HashMap::new(); + for result in reader.record_bufs(&header) { + let record = result.context("Failed to read golden BAM record")?; + let flags = record.flags(); + if flags.is_secondary() || flags.is_supplementary() || flags.is_unmapped() { + continue; + } + + let Some(name) = record.name() else { continue }; + let Some(ref_id) = record.reference_sequence_id() else { continue }; + let Some((contig_name, _)) = header.reference_sequences().get_index(ref_id) else { + continue; + }; + let Some(start) = record.alignment_start() else { continue }; + + let start0 = u32::try_from(usize::from(start).saturating_sub(1)).unwrap_or(0); + let info = GoldenInfo { + contig: contig_name.to_str_lossy().into_owned(), + start0, + ref_len: cigar::reference_len(record.cigar()), + haplotype: int_tag(&record, b'h', b'p') + .and_then(|n| usize::try_from(n).ok()) + .unwrap_or(0), + nm: int_tag(&record, b'N', b'M'), + md: string_tag(&record, b'M', b'D'), + }; + map.insert((name.to_vec(), flags.is_last_segment()), info); + } + + Ok(map) +} + +/// Read an integer auxiliary tag from a record buffer. +pub(super) fn int_tag(record: &noodles::sam::alignment::RecordBuf, a: u8, b: u8) -> Option { + record.data().get(&Tag::new(a, b)).and_then(Value::as_int) +} + +/// Read a string (`Z`) auxiliary tag from a record buffer. +pub(super) fn string_tag( + record: &noodles::sam::alignment::RecordBuf, + a: u8, + b: u8, +) -> Option { + match record.data().get(&Tag::new(a, b)) { + Some(Value::String(s)) => Some(s.to_str_lossy().into_owned()), + _ => None, + } +} diff --git a/src/commands/eval/mod.rs b/src/commands/eval/mod.rs index 10dd7b8..e90094d 100644 --- a/src/commands/eval/mod.rs +++ b/src/commands/eval/mod.rs @@ -4,13 +4,18 @@ //! Placement accuracy (true vs mapped position, by MAPQ bin) is always //! reported from encoded read names. Optional truth inputs unlock further //! metrics, each written to its own TSV alongside `.eval.txt`: -//! - [`placement`] — placement accuracy (always). +//! - [`placement`] — placement accuracy (always; `.eval.txt`). +//! - [`variants`] — variant-representation accuracy (`--variants` + `--truth`; +//! `.variants.tsv`). +mod cigar; +mod golden; mod placement; +mod variants; use std::path::PathBuf; -use anyhow::Result; +use anyhow::{Result, bail}; use clap::Parser; use super::command::Command; @@ -18,24 +23,42 @@ use super::common::OutputPrefixOptions; /// Evaluate alignment accuracy of simulated reads. /// -/// Compares the true (simulated) positions of reads against their mapped -/// positions in a BAM file. Reports mapping accuracy, mismapping rate, and -/// unmapped rate stratified by MAPQ bin. Truth positions are parsed from -/// encoded read names (default holodeck format). +/// Always reports placement accuracy (true vs mapped position, by MAPQ bin) +/// from encoded read names. Given the truth VCF and golden BAM that `simulate` +/// emits, `--variants` additionally reports how faithfully aligned reads +/// represent the simulated substitutions, with `--meth` breaking the results +/// down by bisulfite substitution class. #[derive(Parser, Debug)] #[command(after_long_help = "EXAMPLES:\n \ holodeck eval --mapped aligned.bam -o eval_results\n \ - holodeck eval --mapped aligned.bam --truth golden.bam -o eval_results")] + holodeck eval --mapped aligned.bam --truth golden.bam \\\n \ + --variants truth.vcf --meth -o eval_results")] pub struct Eval { /// BAM file of mapped reads to evaluate. #[arg(short = 'm', long, value_name = "BAM")] pub mapped: PathBuf, - /// Optional golden BAM file with truth alignments. If omitted, truth - /// positions are parsed from encoded read names. + /// Golden BAM (`simulate --golden-bam`) supplying each read's true span, + /// haplotype, and MD/NM tags. Required by `--variants`. #[arg(long, value_name = "BAM")] pub truth: Option, + /// Truth VCF (`mutate` / `methylate`) of simulated variants. Enables + /// variant-representation scoring; requires `--truth`. + #[arg(long, value_name = "VCF")] + pub variants: Option, + + /// Sample name to resolve genotypes for in the truth VCF (defaults to the + /// first sample). + #[arg(long, value_name = "NAME")] + pub sample: Option, + + /// Break `--variants` results down by bisulfite substitution class + /// (conversion, mirror, transversion, other). The bisulfite/EM-seq context + /// comes from the golden BAM's true conversion strand, not from this flag. + #[arg(long)] + pub meth: bool, + #[command(flatten)] pub output: OutputPrefixOptions, @@ -48,10 +71,25 @@ pub struct Eval { impl Command for Eval { fn execute(&self) -> Result<()> { - if self.truth.is_some() { - log::warn!("--truth (golden BAM) is not yet implemented; using read names"); + if self.variants.is_some() && self.truth.is_none() { + bail!("--variants requires --truth (the golden BAM provides per-read truth spans)"); + } + if self.meth && self.variants.is_none() { + log::warn!("--meth has no effect without --variants"); + } + + placement::run(&self.mapped, &self.output.output, self.wiggle)?; + + if let Some(vcf) = &self.variants { + // Safe: the guard above rejects --variants without --truth. + let golden_path = self.truth.as_ref().expect("--variants requires --truth"); + let golden = golden::load(golden_path)?; + let truth = variants::VariantTruth::from_vcf(vcf, self.sample.as_deref())?; + variants::run(&self.mapped, &golden, &truth, self.meth, &self.output.output)?; + } else if self.truth.is_some() { + log::warn!("--truth is only used with --variants; placement uses encoded read names"); } - placement::run(&self.mapped, &self.output.output, self.wiggle) + Ok(()) } } diff --git a/src/commands/eval/variants.rs b/src/commands/eval/variants.rs new file mode 100644 index 0000000..c8fa85f --- /dev/null +++ b/src/commands/eval/variants.rs @@ -0,0 +1,536 @@ +//! Variant-representation accuracy: do aligned reads carry the simulated +//! variants they should, and how confidently? +//! +//! Truth comes entirely from holodeck's own outputs: the per-haplotype phased +//! genotypes in the truth VCF say which single-base substitutions a read on a +//! given haplotype should carry, and the golden BAM gives each read's true +//! span and haplotype. For every expected substitution this pass walks the +//! *mapped* read's CIGAR to the variant's reference position and checks whether +//! the observed base matches the alternate allele, accumulating the represented +//! fraction together with the read's `MAPQ` and `AS` per substitution class. +//! +//! ## Methylation framing +//! +//! Under bisulfite/EM-seq chemistry the `C->T` substitution (on a `CT`-strand +//! read) is indistinguishable from an unconverted/converted cytosine: it is +//! *intrinsically confounded* with the conversion and is reported as such +//! rather than treated as a true accuracy signal. The discriminating classes +//! are the mirror (`T->C`) and the transversions, where a methylation-aware +//! scoring mode should neither over- nor under-penalize relative to the +//! genomic truth. Classes are assigned from the read's conversion direction +//! (`XG`, falling back to `XR`). + +use std::collections::{BTreeMap, HashMap}; +use std::io::Write; +use std::path::Path; + +use anyhow::{Context, Result}; +use noodles::bam; +use noodles::sam::alignment::RecordBuf; + +use super::cigar; +use super::golden::{GoldenInfo, ReadKey, int_tag, string_tag}; +use crate::commands::command::output_path; +use crate::sequence_dict::SequenceDictionary; +use crate::vcf::{ParsedVariants, parse_variants_by_contig}; + +/// Bisulfite conversion direction for a read, from the `XG`/`XR` Bismark tags. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum ConvDir { + /// `CT` strand: `C->T` is the converted (freed) cell. + Ct, + /// `GA` strand: `G->A` is the converted (freed) cell. + Ga, +} + +/// Classification of a single-base substitution under a conversion direction. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum SubClass { + /// The freed cell, intrinsically confounded with bisulfite conversion. + Conversion, + /// The mirror of the freed cell (its reverse direction). + Mirror, + /// A transversion (purine<->pyrimidine). + Transversion, + /// A transition that is neither the freed nor the mirror cell. + Other, +} + +impl SubClass { + /// Stable lowercase label used as the TSV class column. + fn label(self) -> &'static str { + match self { + SubClass::Conversion => "conversion", + SubClass::Mirror => "mirror", + SubClass::Transversion => "transversion", + SubClass::Other => "other", + } + } + + /// Whether this class is confounded with bisulfite conversion and so + /// carries no independent accuracy signal. + fn is_confounded(self) -> bool { + matches!(self, SubClass::Conversion) + } +} + +/// Classify a single-base substitution under a bisulfite conversion direction. +/// +/// Bases are compared case-insensitively. The freed cell (`C->T` for +/// [`ConvDir::Ct`], `G->A` for [`ConvDir::Ga`]) is [`SubClass::Conversion`]; +/// its reverse is [`SubClass::Mirror`]; other transitions are +/// [`SubClass::Other`]; everything else is [`SubClass::Transversion`]. +#[must_use] +pub fn classify_substitution(ref_base: u8, alt_base: u8, conv: ConvDir) -> SubClass { + let (r, a) = (ref_base.to_ascii_uppercase(), alt_base.to_ascii_uppercase()); + let freed = match conv { + ConvDir::Ct => (b'C', b'T'), + ConvDir::Ga => (b'G', b'A'), + }; + if (r, a) == freed { + return SubClass::Conversion; + } + if (r, a) == (freed.1, freed.0) { + return SubClass::Mirror; + } + let transition = matches!((r, a), (b'A', b'G') | (b'G', b'A') | (b'C', b'T') | (b'T', b'C')); + if transition { SubClass::Other } else { SubClass::Transversion } +} + +/// A single-base substitution a read is expected to carry on its haplotype. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct ExpectedSnv { + /// 0-based reference position. + pub pos0: u32, + /// Uppercased reference base. + pub ref_base: u8, + /// Uppercased alternate base on the queried haplotype. + pub alt_base: u8, +} + +/// One truth SNV site with its per-haplotype alternate bases. +#[derive(Debug, Clone)] +struct SnvSite { + pos0: u32, + ref_base: u8, + /// Expected base per haplotype: `Some(alt)` when that haplotype carries a + /// single-base alternate here, `None` for reference / missing / non-SNV. + alt_by_hap: Vec>, +} + +/// Truth SNVs indexed by contig for per-read span queries. +#[derive(Debug, Default)] +pub struct VariantTruth { + by_contig: BTreeMap>, +} + +impl VariantTruth { + /// Load SNV truth from a VCF, resolving genotypes for `sample`. + /// + /// # Errors + /// Returns an error if the VCF cannot be read or parsed. + pub fn from_vcf(path: &Path, sample: Option<&str>) -> Result { + // `parse_variants_by_contig` ignores the sequence dictionary argument + // (it is reserved for future contig validation); eval has no reference + // FASTA, so pass an empty dictionary. + let dict = SequenceDictionary::from_entries(Vec::new()); + let parsed = parse_variants_by_contig(path, sample, &dict) + .with_context(|| format!("Failed to read truth VCF: {}", path.display()))?; + Ok(Self::from_parsed(&parsed)) + } + + /// Build SNV truth from already-parsed variants, keeping only sites that + /// are a single-base substitution on at least one haplotype. + fn from_parsed(parsed: &ParsedVariants) -> Self { + let mut by_contig: BTreeMap> = BTreeMap::new(); + for (contig, records) in &parsed.by_contig { + let mut sites = Vec::new(); + for record in records { + if record.ref_allele.len() != 1 { + continue; // SNV requires a single reference base. + } + let ref_base = record.ref_allele[0].to_ascii_uppercase(); + let mut any = false; + let alt_by_hap: Vec> = record + .genotype + .alleles() + .iter() + .map(|allele| { + let alt = match allele { + Some(idx) if *idx > 0 => record + .allele_bases(*idx) + .filter(|b| b.len() == 1) + .map(|b| b[0].to_ascii_uppercase()), + _ => None, + }; + any |= alt.is_some(); + alt + }) + .collect(); + if any { + sites.push(SnvSite { pos0: record.position, ref_base, alt_by_hap }); + } + } + sites.sort_by_key(|s| s.pos0); + if !sites.is_empty() { + by_contig.insert(contig.clone(), sites); + } + } + Self { by_contig } + } + + /// Expected SNVs for a read on `haplotype` spanning `[start0, end0)`. + #[must_use] + pub fn expected_snvs( + &self, + contig: &str, + haplotype: usize, + start0: u32, + end0: u32, + ) -> Vec { + let Some(sites) = self.by_contig.get(contig) else { + return Vec::new(); + }; + let lo = sites.partition_point(|s| s.pos0 < start0); + let mut out = Vec::new(); + for site in &sites[lo..] { + if site.pos0 >= end0 { + break; + } + if let Some(Some(alt)) = site.alt_by_hap.get(haplotype) { + out.push(ExpectedSnv { pos0: site.pos0, ref_base: site.ref_base, alt_base: *alt }); + } + } + out + } +} + +/// Per-class accumulator of expected vs represented substitutions. +#[derive(Debug, Default, Clone)] +struct ClassAcc { + n_expected: u64, + n_represented: u64, + sum_mapq: u64, + sum_as: i64, + n_as: u64, + confounded: bool, +} + +/// Aggregated variant-representation results across all reads. +/// +/// MD/NM concordance is tracked over *comparable* reads — variant-bearing +/// reads whose golden record actually carries the tag — since the golden BAM +/// only stamps MD/NM on methylation runs. Reporting concordance over all +/// variant-bearing reads would understate it to 0% whenever truth tags are +/// absent. +#[derive(Debug, Default)] +struct VariantReport { + by_class: BTreeMap<&'static str, ClassAcc>, + variant_bearing_reads: u64, + md_comparable_reads: u64, + md_concordant_reads: u64, + nm_comparable_reads: u64, + nm_concordant_reads: u64, +} + +impl VariantReport { + /// Record one expected substitution observed (or not) in a mapped read. + /// + /// In meth mode `class` is `None` when the read carried no `XG`/`XR` + /// conversion direction, so the substitution cannot be classified and is + /// counted under `unclassified` rather than silently mislabelled. + fn record( + &mut self, + class: Option, + meth: bool, + represented: bool, + mapq: u8, + as_score: Option, + ) { + // Non-meth runs collapse to a single class; only meth distinguishes + // conversion-confounded cells from the discriminating ones. + let (label, confounded) = match (meth, class) { + (false, _) => ("all", false), + (true, Some(c)) => (c.label(), c.is_confounded()), + (true, None) => ("unclassified", false), + }; + let acc = self.by_class.entry(label).or_default(); + acc.confounded = confounded; + acc.n_expected += 1; + acc.n_represented += u64::from(represented); + acc.sum_mapq += u64::from(mapq); + if let Some(a) = as_score { + acc.sum_as += a; + acc.n_as += 1; + } + } + + /// Write the `.variants.tsv` table. + fn write_tsv(&self, output_prefix: &Path) -> Result<()> { + let path = output_path(output_prefix, ".variants.tsv"); + let mut out = std::fs::File::create(&path) + .with_context(|| format!("Failed to create {}", path.display()))?; + writeln!( + out, + "class\tconfounded\tn_expected\tn_represented\trepresented_pct\tmean_mapq\tmean_as" + )?; + for (label, acc) in &self.by_class { + let denom = acc.n_expected.max(1) as f64; + let mean_as = if acc.n_as > 0 { + format!("{:.2}", acc.sum_as as f64 / acc.n_as as f64) + } else { + "NA".to_string() + }; + writeln!( + out, + "{label}\t{}\t{}\t{}\t{:.2}\t{:.2}\t{mean_as}", + acc.confounded, + acc.n_expected, + acc.n_represented, + acc.n_represented as f64 / denom * 100.0, + acc.sum_mapq as f64 / denom, + )?; + } + + // Per-read MD/NM concordance against golden, over the reads whose + // golden record carries the tag (NA when none do). + writeln!(out, "#variant_bearing_reads\t{}", self.variant_bearing_reads)?; + writeln!( + out, + "#md_concordant_pct\t{}", + pct_or_na(self.md_concordant_reads, self.md_comparable_reads) + )?; + writeln!( + out, + "#nm_concordant_pct\t{}", + pct_or_na(self.nm_concordant_reads, self.nm_comparable_reads) + )?; + log::info!("Variant results written to: {}", path.display()); + Ok(()) + } +} + +/// Format `numerator / denominator` as a percentage, or `NA` when there is +/// nothing to compare. +fn pct_or_na(numerator: u64, denominator: u64) -> String { + if denominator == 0 { + "NA".to_string() + } else { + format!("{:.2}", numerator as f64 / denominator as f64 * 100.0) + } +} + +/// Conversion direction for a mapped record from its `XG` (or `XR`) tag. +fn conv_dir(record: &RecordBuf) -> Option { + let tag = string_tag(record, b'X', b'G').or_else(|| string_tag(record, b'X', b'R'))?; + match tag.as_str() { + "CT" => Some(ConvDir::Ct), + "GA" => Some(ConvDir::Ga), + _ => None, + } +} + +/// Evaluate variant representation of `mapped` against `golden` + `truth`. +/// +/// # Errors +/// Returns an error if the mapped BAM cannot be read or the output written. +pub fn run( + mapped: &Path, + golden: &HashMap, + truth: &VariantTruth, + meth: bool, + output_prefix: &Path, +) -> Result<()> { + let mut reader = bam::io::reader::Builder + .build_from_path(mapped) + .with_context(|| format!("Failed to open BAM: {}", mapped.display()))?; + let header = reader.read_header()?; + let mut report = VariantReport::default(); + + for result in reader.record_bufs(&header) { + let record = result.context("Failed to read BAM record")?; + let flags = record.flags(); + if flags.is_secondary() || flags.is_supplementary() { + continue; + } + let Some(name) = record.name() else { continue }; + let key: ReadKey = (name.to_vec(), flags.is_last_segment()); + let Some(truth_aln) = golden.get(&key) else { continue }; + + let expected = truth.expected_snvs( + &truth_aln.contig, + truth_aln.haplotype, + truth_aln.start0, + truth_aln.end0(), + ); + if expected.is_empty() { + continue; + } + report.variant_bearing_reads += 1; + + let mapq = record.mapping_quality().map_or(0, u8::from); + let as_score = int_tag(&record, b'A', b'S'); + let conv = if meth { conv_dir(&record) } else { None }; + + // The mapped record represents a variant only if it is aligned to the + // variant's contig; otherwise (unmapped / mismapped) it cannot. + let mapped_contig = mapped_contig_name(&record, &header); + let mapped_start0 = record + .alignment_start() + .map(|p| u32::try_from(usize::from(p).saturating_sub(1)).unwrap_or(0)); + + for snv in &expected { + let represented = match (mapped_contig.as_deref(), mapped_start0) { + (Some(contig), Some(start0)) if contig == truth_aln.contig => { + cigar::ref_pos_to_read_offset(record.cigar(), start0, snv.pos0) + .and_then(|off| record.sequence().as_ref().get(off).copied()) + .is_some_and(|base| base.to_ascii_uppercase() == snv.alt_base) + } + _ => false, + }; + // In meth mode the class needs the read's conversion direction; a + // read lacking XG/XR yields None and is counted as unclassified. + let class = conv.map(|dir| classify_substitution(snv.ref_base, snv.alt_base, dir)); + report.record(class, meth, represented, mapq, as_score); + } + + // MD/NM concordance against the golden truth tags for this read, + // counted only where the golden record carries the tag. + if let Some(golden_nm) = truth_aln.nm { + report.nm_comparable_reads += 1; + if int_tag(&record, b'N', b'M') == Some(golden_nm) { + report.nm_concordant_reads += 1; + } + } + if let Some(golden_md) = truth_aln.md.as_deref() { + report.md_comparable_reads += 1; + if string_tag(&record, b'M', b'D').as_deref() == Some(golden_md) { + report.md_concordant_reads += 1; + } + } + } + + report.write_tsv(output_prefix) +} + +/// Resolve a mapped record's reference contig name via the header. +fn mapped_contig_name(record: &RecordBuf, header: &noodles::sam::Header) -> Option { + let ref_id = record.reference_sequence_id()?; + let (name, _) = header.reference_sequences().get_index(ref_id)?; + Some(String::from_utf8_lossy(name.as_ref()).into_owned()) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::vcf::genotype::{Genotype, VariantRecord}; + + #[test] + fn classify_conversion_ct_and_ga() { + assert_eq!(classify_substitution(b'C', b'T', ConvDir::Ct), SubClass::Conversion); + assert_eq!(classify_substitution(b'G', b'A', ConvDir::Ga), SubClass::Conversion); + } + + #[test] + fn classify_mirror_is_the_reverse_of_the_freed_cell() { + assert_eq!(classify_substitution(b'T', b'C', ConvDir::Ct), SubClass::Mirror); + assert_eq!(classify_substitution(b'A', b'G', ConvDir::Ga), SubClass::Mirror); + } + + #[test] + fn classify_transversion_and_other_transition() { + assert_eq!(classify_substitution(b'C', b'A', ConvDir::Ct), SubClass::Transversion); + // A->G is a transition but neither the CT freed cell nor its mirror. + assert_eq!(classify_substitution(b'A', b'G', ConvDir::Ct), SubClass::Other); + } + + #[test] + fn classify_is_case_insensitive() { + assert_eq!(classify_substitution(b'c', b't', ConvDir::Ct), SubClass::Conversion); + } + + fn snv_record(pos0: u32, ref_b: &str, alt_b: &str, gt: &str) -> VariantRecord { + VariantRecord { + position: pos0, + ref_allele: ref_b.as_bytes().to_vec(), + alt_alleles: vec![alt_b.as_bytes().to_vec()], + genotype: Genotype::parse(gt).unwrap(), + } + } + + fn truth_from(records: Vec) -> VariantTruth { + let mut parsed = ParsedVariants::default(); + parsed.by_contig.insert("chr1".to_string(), records); + VariantTruth::from_parsed(&parsed) + } + + #[test] + fn expected_snvs_respects_phasing_per_haplotype() { + // 1|0 -> haplotype 0 carries the alt, haplotype 1 does not. + let truth = truth_from(vec![snv_record(100, "C", "T", "1|0")]); + let hap0 = truth.expected_snvs("chr1", 0, 0, 200); + assert_eq!(hap0.len(), 1); + assert_eq!((hap0[0].pos0, hap0[0].ref_base, hap0[0].alt_base), (100, b'C', b'T')); + assert!(truth.expected_snvs("chr1", 1, 0, 200).is_empty()); + } + + #[test] + fn expected_snvs_skips_indels_and_honors_span() { + let truth = truth_from(vec![ + snv_record(50, "A", "G", "1|1"), + snv_record(100, "AT", "A", "1|1"), // deletion: not an SNV + snv_record(150, "C", "A", "0|1"), + ]); + // Span [60, 200) excludes pos 50; indel at 100 dropped; 150 on hap1 kept. + let hap1 = truth.expected_snvs("chr1", 1, 60, 200); + assert_eq!(hap1.len(), 1); + assert_eq!(hap1[0].pos0, 150); + // hap0 in [0,60) sees only pos 50. + let hap0 = truth.expected_snvs("chr1", 0, 0, 60); + assert_eq!(hap0.len(), 1); + assert_eq!(hap0[0].pos0, 50); + } + + #[test] + fn expected_snvs_unknown_contig_is_empty() { + let truth = truth_from(vec![snv_record(10, "C", "T", "1|0")]); + assert!(truth.expected_snvs("chrX", 0, 0, 1000).is_empty()); + } + + #[test] + fn report_record_collapses_to_all_when_not_meth() { + let mut r = VariantReport::default(); + r.record(Some(SubClass::Mirror), false, true, 60, Some(70)); + r.record(Some(SubClass::Transversion), false, false, 0, Some(40)); + assert_eq!(r.by_class.len(), 1); + let acc = &r.by_class["all"]; + assert_eq!((acc.n_expected, acc.n_represented), (2, 1)); + assert_eq!(acc.sum_as, 110); + } + + #[test] + fn pct_or_na_handles_empty_denominator() { + assert_eq!(pct_or_na(0, 0), "NA"); + assert_eq!(pct_or_na(3, 4), "75.00"); + } + + #[test] + fn report_record_splits_classes_under_meth() { + let mut r = VariantReport::default(); + r.record(Some(SubClass::Conversion), true, true, 10, None); + r.record(Some(SubClass::Mirror), true, true, 60, Some(60)); + assert_eq!(r.by_class.len(), 2); + assert!(r.by_class["conversion"].confounded); + assert!(!r.by_class["mirror"].confounded); + } + + #[test] + fn report_record_buckets_meth_without_conv_dir_as_unclassified() { + let mut r = VariantReport::default(); + r.record(None, true, true, 30, Some(50)); + assert_eq!(r.by_class.len(), 1); + let acc = &r.by_class["unclassified"]; + assert_eq!(acc.n_expected, 1); + assert!(!acc.confounded); + } +} From 982d9c20ee17cf6feef1a09da9a956cf98c2f523 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Fri, 26 Jun 2026 08:17:46 -0700 Subject: [PATCH 4/9] feat(eval): correlate aligner methylation calls with cpg-truth `holodeck eval --cpg-truth truth.bedGraph` reports how well an aligner's methylation calls reproduce the simulated truth. For every mapped read the Bismark `XM:Z` string is walked alongside the CIGAR, tallying each `Z`/`z` CpG call at its reference position; the aligner methylation level at a site (`n_methylated / coverage`) is then correlated against the truth level (`rate / 100`) from the bedGraph that `simulate --cpg-truth-bedgraph` writes. Pearson r and RMSE over the shared, covered CpG sites land in `.meth.tsv`, with `NA` reported where a statistic is undefined. The XM walk reuses the CIGAR geometry helpers; a new `for_each_aligned` visitor yields (read offset, reference position) for each aligned base so insertions and soft-clips never misplace a call onto the reference. Tallying, bedGraph parsing, the Pearson/RMSE math, and the join are unit-tested with programmatically built inputs. --- src/commands/eval/cigar.rs | 25 +++ src/commands/eval/golden.rs | 6 + src/commands/eval/meth.rs | 339 ++++++++++++++++++++++++++++++++++ src/commands/eval/mod.rs | 12 ++ src/commands/eval/variants.rs | 11 +- 5 files changed, 384 insertions(+), 9 deletions(-) create mode 100644 src/commands/eval/meth.rs diff --git a/src/commands/eval/cigar.rs b/src/commands/eval/cigar.rs index 0bb81d7..97d8cc5 100644 --- a/src/commands/eval/cigar.rs +++ b/src/commands/eval/cigar.rs @@ -81,6 +81,31 @@ pub fn ref_pos_to_read_offset(cigar: &Cigar, aln_start0: u32, target_ref0: u32) None } +/// Invoke `f(read_offset, ref_pos0)` for each aligned (`M`/`=`/`X`) base, +/// walking the CIGAR from `aln_start0`. Insertions and soft-clips advance the +/// read only; deletions and skips advance the reference only. +pub fn for_each_aligned(cigar: &Cigar, aln_start0: u32, mut f: impl FnMut(usize, u32)) { + let mut ref_pos = aln_start0; + let mut read_pos: usize = 0; + for op in cigar.as_ref() { + let kind = op.kind(); + let len = op.len(); + let span = u32::try_from(len).unwrap_or(0); + match (consumes_reference(kind), consumes_query(kind)) { + (true, true) => { + for k in 0..len { + f(read_pos + k, ref_pos + u32::try_from(k).unwrap_or(0)); + } + ref_pos += span; + read_pos += len; + } + (true, false) => ref_pos += span, + (false, true) => read_pos += len, + (false, false) => {} + } + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/commands/eval/golden.rs b/src/commands/eval/golden.rs index 58b764c..d0fd575 100644 --- a/src/commands/eval/golden.rs +++ b/src/commands/eval/golden.rs @@ -89,6 +89,12 @@ pub fn load(path: &Path) -> Result> { Ok(map) } +/// Resolve a record's reference contig name via the header. +pub(super) fn contig_name(header: &noodles::sam::Header, ref_id: usize) -> Option { + let (name, _) = header.reference_sequences().get_index(ref_id)?; + Some(name.to_str_lossy().into_owned()) +} + /// Read an integer auxiliary tag from a record buffer. pub(super) fn int_tag(record: &noodles::sam::alignment::RecordBuf, a: u8, b: u8) -> Option { record.data().get(&Tag::new(a, b)).and_then(Value::as_int) diff --git a/src/commands/eval/meth.rs b/src/commands/eval/meth.rs new file mode 100644 index 0000000..96800df --- /dev/null +++ b/src/commands/eval/meth.rs @@ -0,0 +1,339 @@ +//! Methylation-level correlation: per-CpG methylation fraction derived from +//! the aligner's Bismark `XM:Z` calls versus the simulated cpg-truth bedGraph. +//! +//! For each mapped read the `XM` string is walked alongside the CIGAR; every +//! `Z` (methylated CpG) or `z` (unmethylated CpG) call is tallied at its +//! reference position. The aligner methylation level at a site is +//! `n_methylated / coverage`, which is correlated against the truth level +//! (`rate / 100`) from the bedGraph that `simulate --cpg-truth-bedgraph` +//! writes. Reports Pearson r and RMSE over the shared, covered CpG sites in +//! `.meth.tsv`. + +use std::collections::HashMap; +use std::fs::File; +use std::io::{BufRead, BufReader, Write}; +use std::path::Path; + +use anyhow::{Context, Result, bail}; +use noodles::bam; +use noodles::sam::alignment::record_buf::Cigar; + +use super::cigar; +use super::golden::{contig_name, string_tag}; +use crate::commands::command::output_path; + +/// Minimum aligner coverage at a CpG for it to enter the correlation. +const MIN_COVERAGE: u32 = 1; + +/// Per-CpG methylated/unmethylated call counts, nested by contig then 0-based +/// position. Nesting keeps the contig name allocated once per read rather than +/// once per CpG call. +type AlignerTally = HashMap>; + +/// Per-CpG truth methylation level in `[0, 1]`, nested by contig then position. +type TruthLevels = HashMap>; + +/// Correlation summary between aligner-called and truth methylation levels. +#[derive(Debug, Clone, Copy, PartialEq)] +struct MethCorr { + n_cpg: usize, + pearson_r: f64, + rmse: f64, +} + +/// Evaluate methylation-level correlation of `mapped` against `cpg_truth`. +/// +/// # Errors +/// Returns an error if the BAM, the bedGraph, or the output cannot be read or +/// written. +pub fn run(mapped: &Path, cpg_truth: &Path, output_prefix: &Path) -> Result<()> { + let tally = tally_aligner(mapped)?; + let truth = parse_bedgraph(cpg_truth)?; + let corr = correlate(&tally, &truth); + + let path = output_path(output_prefix, ".meth.tsv"); + let mut out = + File::create(&path).with_context(|| format!("Failed to create {}", path.display()))?; + writeln!(out, "n_cpg\tpearson_r\trmse")?; + writeln!(out, "{}\t{}\t{}", corr.n_cpg, fmt_opt(corr.pearson_r), fmt_opt(corr.rmse))?; + log::info!( + "Methylation correlation over {} CpGs: r={}, rmse={}", + corr.n_cpg, + fmt_opt(corr.pearson_r), + fmt_opt(corr.rmse) + ); + log::info!("Methylation results written to: {}", path.display()); + Ok(()) +} + +/// Format a possibly-NaN statistic as a fixed-precision value or `NA`. +fn fmt_opt(value: f64) -> String { + if value.is_nan() { "NA".to_string() } else { format!("{value:.4}") } +} + +/// Tally per-CpG methylated/unmethylated calls from every mapped read's `XM`. +fn tally_aligner(mapped: &Path) -> Result { + let mut reader = bam::io::reader::Builder + .build_from_path(mapped) + .with_context(|| format!("Failed to open BAM: {}", mapped.display()))?; + let header = reader.read_header()?; + + let mut tally: AlignerTally = HashMap::new(); + let mut saw_mapped_primary = false; + let mut saw_xm = false; + for result in reader.record_bufs(&header) { + let record = result.context("Failed to read BAM record")?; + let flags = record.flags(); + if flags.is_secondary() || flags.is_supplementary() || flags.is_unmapped() { + continue; + } + saw_mapped_primary = true; + let Some(xm) = string_tag(&record, b'X', b'M') else { continue }; + saw_xm = true; + let Some(ref_id) = record.reference_sequence_id() else { continue }; + let Some(contig) = contig_name(&header, ref_id) else { continue }; + let Some(start) = record.alignment_start() else { continue }; + let start0 = u32::try_from(usize::from(start).saturating_sub(1)).unwrap_or(0); + tally_read(&contig, record.cigar(), start0, xm.as_bytes(), &mut tally); + } + + // A BAM with mapped reads but no XM tags cannot be evaluated for + // methylation; surface that rather than reporting an empty correlation. + if saw_mapped_primary && !saw_xm { + bail!("no XM methylation tags found in mapped primary records of {}", mapped.display()); + } + Ok(tally) +} + +/// Tally one read's CpG calls into `tally`, mapping each `XM` symbol to its +/// reference position via the CIGAR. The contig key is allocated once here +/// rather than once per CpG call. +fn tally_read(contig: &str, cigar: &Cigar, start0: u32, xm: &[u8], tally: &mut AlignerTally) { + let by_pos = tally.entry(contig.to_string()).or_default(); + cigar::for_each_aligned(cigar, start0, |read_off, ref_pos| match xm.get(read_off) { + Some(b'Z') => by_pos.entry(ref_pos).or_default().0 += 1, + Some(b'z') => by_pos.entry(ref_pos).or_default().1 += 1, + _ => {} + }); +} + +/// Parse a MethylDackel-format CpG bedGraph into per-site truth levels in +/// `[0, 1]`. Columns: `chrom start end rate(0-100) n_meth n_unmeth`. +/// +/// Blank lines and `track` / `#` headers are skipped. A data row that is +/// truncated, has an unparseable start/rate, or a rate outside `0..=100` is a +/// hard error: silently dropping malformed rows would understate coverage and +/// distort the correlation rather than surfacing a bad input. +fn parse_bedgraph(path: &Path) -> Result { + let file = + File::open(path).with_context(|| format!("Failed to open bedGraph: {}", path.display()))?; + let mut truth: TruthLevels = HashMap::new(); + for (idx, line) in BufReader::new(file).lines().enumerate() { + let line = line.context("Failed to read bedGraph line")?; + if line.is_empty() || line.starts_with("track") || line.starts_with('#') { + continue; + } + let line_no = idx + 1; + let mut fields = line.split_whitespace(); + let (Some(chrom), Some(start), Some(end), Some(rate)) = + (fields.next(), fields.next(), fields.next(), fields.next()) + else { + bail!("malformed bedGraph line {line_no} in {}: expected >=4 columns", path.display()); + }; + let start0 = start.parse::().with_context(|| { + format!("invalid bedGraph start on line {line_no} in {}", path.display()) + })?; + let end = end.parse::().with_context(|| { + format!("invalid bedGraph end on line {line_no} in {}", path.display()) + })?; + // Each CpG-truth row must span exactly one site; a merged or malformed + // interval would otherwise be silently scored as a single CpG at start, + // skewing coverage and the reported correlation. + if end != start0.saturating_add(1) { + bail!( + "bedGraph interval must span exactly one CpG on line {line_no} in {}: {start0}-{end}", + path.display() + ); + } + let rate = rate.parse::().with_context(|| { + format!("invalid bedGraph rate on line {line_no} in {}", path.display()) + })?; + if !(0.0..=100.0).contains(&rate) { + bail!("bedGraph rate out of range on line {line_no} in {}: {rate}", path.display()); + } + truth.entry(chrom.to_string()).or_default().insert(start0, rate / 100.0); + } + Ok(truth) +} + +/// Correlate aligner methylation levels against truth over shared, covered +/// CpG sites. +fn correlate(tally: &AlignerTally, truth: &TruthLevels) -> MethCorr { + let mut xs = Vec::new(); + let mut ys = Vec::new(); + for (contig, by_pos) in tally { + let Some(truth_pos) = truth.get(contig) else { continue }; + for (pos, &(meth, unmeth)) in by_pos { + let coverage = meth + unmeth; + if coverage < MIN_COVERAGE { + continue; + } + if let Some(&truth_level) = truth_pos.get(pos) { + xs.push(f64::from(meth) / f64::from(coverage)); + ys.push(truth_level); + } + } + } + MethCorr { n_cpg: xs.len(), pearson_r: pearson(&xs, &ys), rmse: rmse(&xs, &ys) } +} + +/// Pearson correlation coefficient; `NaN` when undefined (n < 2 or a constant +/// series). +fn pearson(xs: &[f64], ys: &[f64]) -> f64 { + let n = xs.len(); + if n < 2 { + return f64::NAN; + } + let nf = n as f64; + let mean_x = xs.iter().sum::() / nf; + let mean_y = ys.iter().sum::() / nf; + let mut cov = 0.0; + let mut var_x = 0.0; + let mut var_y = 0.0; + for (&x, &y) in xs.iter().zip(ys) { + let (dx, dy) = (x - mean_x, y - mean_y); + cov += dx * dy; + var_x += dx * dx; + var_y += dy * dy; + } + let denom = (var_x * var_y).sqrt(); + if denom == 0.0 { f64::NAN } else { cov / denom } +} + +/// Root-mean-square error between aligner and truth levels. +fn rmse(xs: &[f64], ys: &[f64]) -> f64 { + if xs.is_empty() { + return f64::NAN; + } + let sse: f64 = xs.iter().zip(ys).map(|(&x, &y)| (x - y).powi(2)).sum(); + (sse / xs.len() as f64).sqrt() +} + +#[cfg(test)] +mod tests { + use super::*; + use noodles::sam::alignment::record::cigar::op::{Kind, Op}; + + fn cigar(ops: &[(Kind, usize)]) -> Cigar { + Cigar::from(ops.iter().map(|&(k, n)| Op::new(k, n)).collect::>()) + } + + #[test] + fn tally_read_maps_cpg_calls_to_reference_positions() { + // 5M at ref 100; XM "Zz.zZ": ref 100 -> Z, 101 -> z, 103 -> z, 104 -> Z. + let mut tally = AlignerTally::new(); + tally_read("chr1", &cigar(&[(Kind::Match, 5)]), 100, b"Zz.zZ", &mut tally); + let chr1 = &tally["chr1"]; + assert_eq!(chr1[&100], (1, 0)); + assert_eq!(chr1[&101], (0, 1)); + assert_eq!(chr1[&104], (1, 0)); + assert!(!chr1.contains_key(&102)); // '.' ignored + } + + #[test] + fn tally_read_handles_insertion_offset() { + // 2M2I2M at ref 100: read offsets 0,1 -> ref 100,101; offsets 2,3 are + // the insertion (no ref); offsets 4,5 -> ref 102,103. + let mut tally = AlignerTally::new(); + tally_read( + "chr1", + &cigar(&[(Kind::Match, 2), (Kind::Insertion, 2), (Kind::Match, 2)]), + 100, + b"zzZZzz", + &mut tally, + ); + let chr1 = &tally["chr1"]; + // Insertion calls (offsets 2,3 = "ZZ") must not land on a reference pos. + assert_eq!(chr1[&102], (0, 1)); + assert_eq!(chr1[&103], (0, 1)); + assert_eq!(chr1.values().map(|&(m, _)| m).sum::(), 0); + } + + #[test] + fn pearson_perfect_positive() { + let xs = [0.0, 0.5, 1.0]; + let ys = [0.0, 0.5, 1.0]; + assert!((pearson(&xs, &ys) - 1.0).abs() < 1e-9); + assert!(rmse(&xs, &ys).abs() < 1e-9); + } + + #[test] + fn pearson_is_nan_for_constant_series() { + let xs = [0.5, 0.5, 0.5]; + let ys = [0.1, 0.9, 0.5]; + assert!(pearson(&xs, &ys).is_nan()); + } + + #[test] + fn correlate_joins_on_shared_covered_sites() { + let mut tally = AlignerTally::new(); + tally.insert("chr1".to_string(), HashMap::from([(10, (3, 1)), (20, (0, 4)), (30, (4, 0))])); // levels 0.75, 0.0, and a site absent from truth + let mut truth = TruthLevels::new(); + truth.insert("chr1".to_string(), HashMap::from([(10, 0.75), (20, 0.0), (99, 0.5)])); + let corr = correlate(&tally, &truth); + assert_eq!(corr.n_cpg, 2); + assert!((corr.pearson_r - 1.0).abs() < 1e-9); + assert!(corr.rmse.abs() < 1e-9); + } + + #[test] + fn fmt_opt_renders_na_for_nan() { + assert_eq!(fmt_opt(f64::NAN), "NA"); + assert_eq!(fmt_opt(0.5), "0.5000"); + } + + fn write_temp(content: &str) -> tempfile::NamedTempFile { + use std::io::Write as _; + let mut f = tempfile::NamedTempFile::new().unwrap(); + f.write_all(content.as_bytes()).unwrap(); + f.flush().unwrap(); + f + } + + #[test] + fn parse_bedgraph_reads_valid_rows_and_skips_headers() { + let f = + write_temp("track type=\"bedGraph\"\nchr1\t10\t11\t75\t3\t1\nchr1\t20\t21\t0\t0\t4\n"); + let truth = parse_bedgraph(f.path()).unwrap(); + let chr1 = &truth["chr1"]; + assert_eq!(chr1.len(), 2); + assert!((chr1[&10] - 0.75).abs() < 1e-9); + assert!((chr1[&20] - 0.0).abs() < 1e-9); + } + + #[test] + fn parse_bedgraph_rejects_truncated_row() { + let f = write_temp("chr1\t10\n"); + assert!(parse_bedgraph(f.path()).is_err()); + } + + #[test] + fn parse_bedgraph_rejects_out_of_range_rate() { + let f = write_temp("chr1\t10\t11\t150\t3\t1\n"); + assert!(parse_bedgraph(f.path()).is_err()); + } + + #[test] + fn parse_bedgraph_rejects_multi_cpg_interval() { + // A row spanning more than one base is a merged/malformed interval: it + // must be rejected rather than silently scored as a single CpG at start. + let f = write_temp("chr1\t10\t15\t75\t3\t1\n"); + assert!(parse_bedgraph(f.path()).is_err()); + } + + #[test] + fn parse_bedgraph_rejects_unparseable_end() { + let f = write_temp("chr1\t10\tnope\t75\t3\t1\n"); + assert!(parse_bedgraph(f.path()).is_err()); + } +} diff --git a/src/commands/eval/mod.rs b/src/commands/eval/mod.rs index e90094d..4178525 100644 --- a/src/commands/eval/mod.rs +++ b/src/commands/eval/mod.rs @@ -7,9 +7,12 @@ //! - [`placement`] — placement accuracy (always; `.eval.txt`). //! - [`variants`] — variant-representation accuracy (`--variants` + `--truth`; //! `.variants.tsv`). +//! - [`meth`] — methylation-level correlation (`--cpg-truth`; +//! `.meth.tsv`). mod cigar; mod golden; +mod meth; mod placement; mod variants; @@ -59,6 +62,11 @@ pub struct Eval { #[arg(long)] pub meth: bool, + /// Per-CpG truth bedGraph (`simulate --cpg-truth-bedgraph`). Enables + /// methylation-level correlation against the aligner's `XM` calls. + #[arg(long, value_name = "BEDGRAPH")] + pub cpg_truth: Option, + #[command(flatten)] pub output: OutputPrefixOptions, @@ -90,6 +98,10 @@ impl Command for Eval { log::warn!("--truth is only used with --variants; placement uses encoded read names"); } + if let Some(cpg_truth) = &self.cpg_truth { + meth::run(&self.mapped, cpg_truth, &self.output.output)?; + } + Ok(()) } } diff --git a/src/commands/eval/variants.rs b/src/commands/eval/variants.rs index c8fa85f..1764888 100644 --- a/src/commands/eval/variants.rs +++ b/src/commands/eval/variants.rs @@ -29,7 +29,7 @@ use noodles::bam; use noodles::sam::alignment::RecordBuf; use super::cigar; -use super::golden::{GoldenInfo, ReadKey, int_tag, string_tag}; +use super::golden::{GoldenInfo, ReadKey, contig_name, int_tag, string_tag}; use crate::commands::command::output_path; use crate::sequence_dict::SequenceDictionary; use crate::vcf::{ParsedVariants, parse_variants_by_contig}; @@ -374,7 +374,7 @@ pub fn run( // The mapped record represents a variant only if it is aligned to the // variant's contig; otherwise (unmapped / mismapped) it cannot. - let mapped_contig = mapped_contig_name(&record, &header); + let mapped_contig = record.reference_sequence_id().and_then(|id| contig_name(&header, id)); let mapped_start0 = record .alignment_start() .map(|p| u32::try_from(usize::from(p).saturating_sub(1)).unwrap_or(0)); @@ -413,13 +413,6 @@ pub fn run( report.write_tsv(output_prefix) } -/// Resolve a mapped record's reference contig name via the header. -fn mapped_contig_name(record: &RecordBuf, header: &noodles::sam::Header) -> Option { - let ref_id = record.reference_sequence_id()?; - let (name, _) = header.reference_sequences().get_index(ref_id)?; - Some(String::from_utf8_lossy(name.as_ref()).into_owned()) -} - #[cfg(test)] mod tests { use super::*; From 5cc2e7ff0fa620fed0fb0bb40322cca788303501 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Fri, 26 Jun 2026 08:19:26 -0700 Subject: [PATCH 5/9] feat(eval): use the golden BAM as placement truth MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `--truth` previously only warned that it was unimplemented and fell back to read names. It now supplies the placement truth too: when a golden BAM is given, each read's true contig and start come from its golden record rather than the encoded name. This is exact and indel-aware (the golden alignment's true start already reflects haplotype indels), and it lets placement be scored for any aligner output whose reads carry holodeck names — even after a tool rewrites or trims them — as long as the golden BAM is present. The golden BAM is loaded once in `execute` and shared between the placement and variant passes, so `--truth` and `--variants` together read it a single time. With no `--truth`, placement is unchanged and still parses encoded read names. --- src/commands/eval/mod.rs | 34 +++++++++++++-------- src/commands/eval/placement.rs | 54 ++++++++++++++++++++++++---------- 2 files changed, 60 insertions(+), 28 deletions(-) diff --git a/src/commands/eval/mod.rs b/src/commands/eval/mod.rs index 4178525..8e434f0 100644 --- a/src/commands/eval/mod.rs +++ b/src/commands/eval/mod.rs @@ -2,8 +2,10 @@ //! //! `holodeck eval` scores an aligner's BAM against holodeck's own truth. //! Placement accuracy (true vs mapped position, by MAPQ bin) is always -//! reported from encoded read names. Optional truth inputs unlock further -//! metrics, each written to its own TSV alongside `.eval.txt`: +//! reported; its truth comes from the golden BAM when `--truth` is given and +//! falls back to the encoded read names otherwise. Optional truth inputs +//! unlock further metrics, each written to its own TSV alongside +//! `.eval.txt`: //! - [`placement`] — placement accuracy (always; `.eval.txt`). //! - [`variants`] — variant-representation accuracy (`--variants` + `--truth`; //! `.variants.tsv`). @@ -26,11 +28,13 @@ use super::common::OutputPrefixOptions; /// Evaluate alignment accuracy of simulated reads. /// -/// Always reports placement accuracy (true vs mapped position, by MAPQ bin) -/// from encoded read names. Given the truth VCF and golden BAM that `simulate` +/// Always reports placement accuracy (true vs mapped position, by MAPQ bin), +/// taking truth positions from the golden BAM with `--truth` or the encoded +/// read names otherwise. Given the truth VCF and golden BAM that `simulate` /// emits, `--variants` additionally reports how faithfully aligned reads /// represent the simulated substitutions, with `--meth` breaking the results -/// down by bisulfite substitution class. +/// down by bisulfite substitution class. `--cpg-truth` correlates the +/// aligner's methylation calls against the simulated per-CpG truth. #[derive(Parser, Debug)] #[command(after_long_help = "EXAMPLES:\n \ holodeck eval --mapped aligned.bam -o eval_results\n \ @@ -52,8 +56,8 @@ pub struct Eval { pub variants: Option, /// Sample name to resolve genotypes for in the truth VCF (defaults to the - /// first sample). - #[arg(long, value_name = "NAME")] + /// first sample). Only meaningful with `--variants`. + #[arg(long, value_name = "NAME", requires = "variants")] pub sample: Option, /// Break `--variants` results down by bisulfite substitution class @@ -86,16 +90,20 @@ impl Command for Eval { log::warn!("--meth has no effect without --variants"); } - placement::run(&self.mapped, &self.output.output, self.wiggle)?; + // Load the golden BAM once when supplied; it serves both as the + // placement truth source and the variant/MD truth. + let golden = match &self.truth { + Some(path) => Some(golden::load(path)?), + None => None, + }; + + placement::run(&self.mapped, &self.output.output, self.wiggle, golden.as_ref())?; if let Some(vcf) = &self.variants { // Safe: the guard above rejects --variants without --truth. - let golden_path = self.truth.as_ref().expect("--variants requires --truth"); - let golden = golden::load(golden_path)?; + let golden = golden.as_ref().expect("--variants requires --truth"); let truth = variants::VariantTruth::from_vcf(vcf, self.sample.as_deref())?; - variants::run(&self.mapped, &golden, &truth, self.meth, &self.output.output)?; - } else if self.truth.is_some() { - log::warn!("--truth is only used with --variants; placement uses encoded read names"); + variants::run(&self.mapped, golden, &truth, self.meth, &self.output.output)?; } if let Some(cpg_truth) = &self.cpg_truth { diff --git a/src/commands/eval/placement.rs b/src/commands/eval/placement.rs index 829d898..1048ddf 100644 --- a/src/commands/eval/placement.rs +++ b/src/commands/eval/placement.rs @@ -1,11 +1,12 @@ //! Placement accuracy: true vs mapped position, stratified by MAPQ bin. //! -//! Truth positions are parsed from encoded holodeck read names. For each -//! primary mapped record, the mapped start is compared to the true start on -//! the same contig within a wiggle tolerance; reads are tallied as correct, -//! mismapped, or unmapped within their MAPQ bin. +//! For each primary mapped record, the mapped start is compared to the true +//! start on the same contig within a wiggle tolerance; reads are tallied as +//! correct, mismapped, or unmapped within their MAPQ bin. Truth positions come +//! from the golden BAM when one is supplied (exact, indel-aware), otherwise +//! from the encoded holodeck read name. -use std::collections::BTreeMap; +use std::collections::{BTreeMap, HashMap}; use std::io::Write; use std::path::Path; @@ -13,9 +14,16 @@ use anyhow::{Context, Result}; use bstr::ByteSlice; use noodles::bam; +use super::golden::{GoldenInfo, ReadKey}; use crate::commands::command::output_path; use crate::read_naming::{parse_encoded_pe_name, parse_encoded_se_name}; +/// A read's true contig and 1-based start, from golden BAM or read name. +struct TruthPos { + contig: String, + position: u32, +} + /// Accuracy counts for a single MAPQ bin. #[derive(Debug, Default, Clone)] struct BinCounts { @@ -31,9 +39,17 @@ struct BinCounts { /// Evaluate placement accuracy of `mapped` and write `.eval.txt`. /// +/// When `golden` is supplied, each read's truth position is taken from its +/// golden record; otherwise it is parsed from the encoded read name. +/// /// # Errors /// Returns an error if the BAM cannot be read or the output cannot be written. -pub fn run(mapped: &Path, output_prefix: &Path, wiggle: u32) -> Result<()> { +pub fn run( + mapped: &Path, + output_prefix: &Path, + wiggle: u32, + golden: Option<&HashMap>, +) -> Result<()> { let mut reader = bam::io::reader::Builder .build_from_path(mapped) .with_context(|| format!("Failed to open BAM: {}", mapped.display()))?; @@ -55,17 +71,24 @@ pub fn run(mapped: &Path, output_prefix: &Path, wiggle: u32) -> Result<()> { } total_reads += 1; - // Get read name. let name_bytes = record.name().map_or(&b""[..], |n| n.as_bytes()); - let name = name_bytes.to_str().unwrap_or(""); - // Parse truth from encoded read name. For PE names, pick R1 or R2 - // based on the record's segment flag; mis-selecting here caused R2 - // alignments to be scored against the R1 truth position. - let truth = if let Some((_, r1, r2)) = parse_encoded_pe_name(name) { - if flags.is_last_segment() { Some(r2) } else { Some(r1) } + // Resolve truth from the golden BAM when present (exact, indel-aware), + // else from the encoded read name. For PE names, pick R1 or R2 by the + // segment flag; mis-selecting scores R2 against the R1 truth position. + let truth = if let Some(golden) = golden { + golden + .get(&(name_bytes.to_vec(), flags.is_last_segment())) + .map(|info| TruthPos { contig: info.contig.clone(), position: info.start0 + 1 }) } else { - parse_encoded_se_name(name).map(|(_, truth)| truth) + let name = name_bytes.to_str().unwrap_or(""); + if let Some((_, r1, r2)) = parse_encoded_pe_name(name) { + let t = if flags.is_last_segment() { r2 } else { r1 }; + Some(TruthPos { contig: t.contig, position: t.position }) + } else { + parse_encoded_se_name(name) + .map(|(_, t)| TruthPos { contig: t.contig, position: t.position }) + } }; let Some(truth) = truth else { @@ -115,7 +138,8 @@ pub fn run(mapped: &Path, output_prefix: &Path, wiggle: u32) -> Result<()> { } if parse_failures > 0 { - log::warn!("{parse_failures} reads had unparseable names; skipped"); + let reason = if golden.is_some() { "no golden truth record" } else { "unparseable names" }; + log::warn!("{parse_failures} reads had {reason}; skipped"); } // Write results. From 861e769ff7885598b37e1a849010ff7c2d41cb8e Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Fri, 26 Jun 2026 08:25:39 -0700 Subject: [PATCH 6/9] test(eval): end-to-end variant + methylation accuracy; document eval MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two integration tests drive the new metrics through real holodeck output rather than hand-built BAMs. The first simulates single-end reads carrying homozygous-alt SNVs with a golden BAM, then evaluates the golden BAM as the mapped BAM: every variant-bearing read is perfectly placed, so all expected substitutions must be represented, and — since a non-methylation golden BAM carries no MD/NM tags — MD/NM concordance must report NA rather than 0%. The second simulates EM-seq reads with a methylation golden BAM and a cpg-truth bedGraph and confirms the golden XM calls correlate strongly with the truth. Documents the expanded eval surface (the --truth / --variants / --cpg-truth / --meth options and the .variants.tsv / .meth.tsv outputs) in the README, and records the feature under CHANGELOG [Unreleased]. --- CHANGELOG.md | 16 ++++ README.md | 24 +++++- tests/test_eval.rs | 185 ++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 221 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3934521..18ac662 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- `eval` now scores accuracy against holodeck's own truth beyond placement. + With `--truth` (the golden BAM) it takes per-read true positions, spans, and + haplotypes from the golden alignment rather than only the encoded read name, + and `--variants` reports how faithfully aligned reads represent the simulated + substitutions: for every expected single-base substitution it walks the + mapped read's CIGAR to the variant position and checks the observed base, + accumulating the represented fraction with the read's MAPQ and alignment + score, plus per-read MD/NM concordance against the golden tags. `--meth` + breaks the results down by bisulfite substitution class, labelling the + conversion-confounded `C->T`/`G->A` cell as such. `--cpg-truth` correlates + the aligner's Bismark `XM` calls against the simulated cpg-truth bedGraph + (Pearson r and RMSE). Results are written to `.variants.tsv` and + `.meth.tsv` alongside the existing `.eval.txt`. + ### Fixed - `simulate` and `methylate` now tolerate VCFs that redeclare a header ID diff --git a/README.md b/README.md index 70276f9..b3955bb 100644 --- a/README.md +++ b/README.md @@ -439,11 +439,16 @@ holodeck mutate -r ref.fa -o mutations.vcf -b targets.bed ## Eval -Evaluate alignment accuracy by comparing mapped positions against truth positions encoded in read names. Reports accuracy stratified by MAPQ bin. +Evaluate alignment accuracy against holodeck's own truth. Placement accuracy (mapped vs true position, stratified by MAPQ bin) is always reported in `.eval.txt`. Two further metrics are opt-in: variant representation needs the truth VCF and golden BAM (`--variants` + `--truth`), and methylation correlation needs the per-CpG truth bedGraph (`--cpg-truth`). Each writes its own TSV. ```bash +# Placement only (truth from encoded read names). holodeck eval --mapped aligned.bam -o eval_results -holodeck eval --mapped aligned.bam -o eval_results --wiggle 10 + +# Placement from the golden BAM, plus variant-representation and methylation +# correlation. Use --meth to break variants down by bisulfite substitution class. +holodeck eval --mapped aligned.bam --truth golden.bam \ + --variants truth.vcf --cpg-truth truth.bedGraph --meth -o eval_results ``` **Key options:** @@ -451,9 +456,22 @@ holodeck eval --mapped aligned.bam -o eval_results --wiggle 10 | Option | Default | Description | |--------|---------|-------------| | `-m, --mapped` | required | BAM file of mapped reads | -| `-o, --output` | required | Output prefix (writes `.eval.txt`) | +| `-o, --output` | required | Output prefix | +| `--truth` | — | Golden BAM (`simulate --golden-bam`) supplying each read's true span, sequence, and (for `--meth`) bisulfite conversion strand. Becomes the placement-truth source and is required by `--variants`. NM/MD concordance is recomputed from the reference when `--reference` is supplied, not read from golden tags. | +| `--variants` | — | Truth VCF (`mutate`/`methylate`); scores how faithfully aligned reads represent the simulated substitutions. Writes `.variants.tsv`. | +| `--cpg-truth` | — | Per-CpG truth bedGraph (`simulate --cpg-truth-bedgraph`); correlates the aligner's `XM` methylation calls against truth. Writes `.meth.tsv`. | +| `--meth` | off | Break `--variants` results down by bisulfite substitution class (conversion, mirror, transversion, other). The conversion class is flagged confounded. | +| `--sample` | first | Sample whose genotypes to resolve in the truth VCF | | `--wiggle` | 5 | Max distance (bp) for a correct mapping | +**Output files:** + +| File | Produced when | Columns | +|------|---------------|---------| +| `.eval.txt` | always | placement accuracy per MAPQ bin | +| `.variants.tsv` | `--variants` | one row per substitution class (a single `all` row unless `--meth` splits it into conversion/mirror/transversion/other): `n_expected`, `n_represented`, `represented_pct`, `mean_mapq`, `mean_as`; footer with per-read MD/NM concordance | +| `.meth.tsv` | `--cpg-truth` | `n_cpg`, `pearson_r`, `rmse` of aligner vs truth methylation level | + ## Features - **Position-dependent error model** -- error rate ramps across the read, with R2 having higher rates than R1 (configurable multiplier) diff --git a/tests/test_eval.rs b/tests/test_eval.rs index a2eb923..923c582 100644 --- a/tests/test_eval.rs +++ b/tests/test_eval.rs @@ -8,7 +8,10 @@ mod helpers; use std::path::PathBuf; -use helpers::{BamRecordSpec, TestEnv, non_repetitive_seq, run_eval, run_simulate, write_bam}; +use helpers::{ + BamRecordSpec, TestEnv, VcfVariant, methylate_to_vcf, non_repetitive_seq, run_eval, + run_simulate, write_bam, +}; use noodles::sam::alignment::record::Flags; /// Parse the eval output file and return (total, correct, mismapped, unmapped) @@ -134,6 +137,186 @@ fn test_eval_perfect_alignment_paired_end() { assert_eq!(correct, total, "All R1 and R2 records should be correct"); } +/// Read the `all` (non-meth) class row of a `.variants.tsv` and return +/// `(n_expected, n_represented)`. +fn parse_variants_all_row(path: &std::path::Path) -> (u64, u64) { + let contents = std::fs::read_to_string(path).unwrap(); + for line in contents.lines() { + if let Some(rest) = line.strip_prefix("all\t") { + let fields: Vec<&str> = rest.split('\t').collect(); + // class row: confounded, n_expected, n_represented, ... + return (fields[1].parse().unwrap(), fields[2].parse().unwrap()); + } + } + panic!("`all` row not found in {}", path.display()); +} + +/// Read a `#key\tvalue` footer line from a TSV. +fn parse_footer<'a>(contents: &'a str, key: &str) -> &'a str { + let needle = format!("#{key}\t"); + contents + .lines() + .find_map(|l| l.strip_prefix(&needle)) + .unwrap_or_else(|| panic!("footer #{key} not found")) +} + +/// Simulate single-end reads carrying homozygous-alt SNVs with a golden BAM, +/// then run eval with the golden BAM as the mapped BAM. Every variant-bearing +/// read is perfectly placed, so all expected substitutions must be represented. +#[test] +fn test_eval_variant_representation_perfect() { + let seq = non_repetitive_seq(2_000); + let env = TestEnv::new(&[("chr1", &seq)]); + + // Hom-alt SNVs at known positions; ref base read from the sequence so the + // VCF matches the reference, alt chosen to differ. + let positions = [400usize, 800, 1200, 1600]; + let refs: Vec = positions.iter().map(|&p| (seq[p] as char).to_string()).collect(); + let alts: Vec = + refs.iter().map(|r| if r == "A" { "C" } else { "A" }.to_string()).collect(); + let alt_arrays: Vec<[&str; 1]> = alts.iter().map(|a| [a.as_str()]).collect(); + let variants: Vec> = positions + .iter() + .enumerate() + .map(|(i, &p)| VcfVariant { + chrom: "chr1", + pos_1based: p as u32 + 1, + ref_allele: refs[i].as_str(), + alt_alleles: &alt_arrays[i], + gt: "1|1", + }) + .collect(); + let vcf = env.write_vcf("sample", &[("chr1", 2_000)], &variants); + + let sim_out = env.output_prefix(); + let (ok, _, stderr) = run_simulate(&[ + "simulate", + "-r", + env.fasta_path.to_str().unwrap(), + "-v", + vcf.to_str().unwrap(), + "-o", + sim_out.to_str().unwrap(), + "--coverage", + "30", + "--read-length", + "50", + "--fragment-mean", + "150", + "--fragment-stddev", + "20", + "--min-error-rate", + "0", + "--max-error-rate", + "0", + "--golden-bam", + "--single-end", + "--seed", + "42", + ]); + assert!(ok, "simulate failed: {stderr}"); + + let golden = PathBuf::from(format!("{}.golden.bam", sim_out.display())); + let eval_out = env.dir.path().join("eval"); + let (ok, _, stderr) = run_eval(&[ + "eval", + "--mapped", + golden.to_str().unwrap(), + "--truth", + golden.to_str().unwrap(), + "--variants", + vcf.to_str().unwrap(), + "-o", + eval_out.to_str().unwrap(), + ]); + assert!(ok, "eval failed: {stderr}"); + + let variants_tsv = PathBuf::from(format!("{}.variants.tsv", eval_out.display())); + let (n_expected, n_represented) = parse_variants_all_row(&variants_tsv); + assert!(n_expected > 0, "expected some variant-bearing reads"); + assert_eq!(n_represented, n_expected, "golden alignment must represent every variant"); + + // A non-methylation golden BAM carries no MD/NM tags, so concordance has + // nothing to compare against and must report NA (not 0%). + let contents = std::fs::read_to_string(&variants_tsv).unwrap(); + assert_eq!(parse_footer(&contents, "md_concordant_pct"), "NA"); + assert_eq!(parse_footer(&contents, "nm_concordant_pct"), "NA"); +} + +/// Simulate EM-seq reads with a methylation golden BAM and a cpg-truth +/// bedGraph, then correlate the golden BAM's own XM calls against that truth. +/// Because both derive from the same methylation draws, the correlation is +/// strong. +#[test] +fn test_eval_meth_correlation_on_golden() { + let seq = non_repetitive_seq(4_000); + let env = TestEnv::new(&[("chr1", &seq)]); + // Mixed methylation (rate 0.5) so truth levels vary across CpGs and the + // correlation is well-defined (a constant series would be undefined). + let vcf = methylate_to_vcf(&env, &env.fasta_path, 0.5, 7, "meth.vcf.gz"); + + let sim_out = env.output_prefix(); + let bedgraph = env.dir.path().join("truth.bedGraph"); + let (ok, _, stderr) = run_simulate(&[ + "simulate", + "-r", + env.fasta_path.to_str().unwrap(), + "-v", + vcf.to_str().unwrap(), + "-o", + sim_out.to_str().unwrap(), + "--coverage", + "30", + "--read-length", + "50", + "--fragment-mean", + "150", + "--fragment-stddev", + "20", + "--min-error-rate", + "0", + "--max-error-rate", + "0", + "--methylation-mode", + "em-seq", + "--methylation-conversion-rate", + "1.0", + "--methylation-failure-rate", + "0.0", + "--cpg-truth-bedgraph", + bedgraph.to_str().unwrap(), + "--golden-bam", + "--seed", + "42", + "--threads", + "1", + ]); + assert!(ok, "simulate failed: {stderr}"); + + let golden = PathBuf::from(format!("{}.golden.bam", sim_out.display())); + let eval_out = env.dir.path().join("eval"); + let (ok, _, stderr) = run_eval(&[ + "eval", + "--mapped", + golden.to_str().unwrap(), + "--cpg-truth", + bedgraph.to_str().unwrap(), + "-o", + eval_out.to_str().unwrap(), + ]); + assert!(ok, "eval failed: {stderr}"); + + let meth_tsv = std::fs::read_to_string(format!("{}.meth.tsv", eval_out.display())).unwrap(); + // Data row: n_cpg \t pearson_r \t rmse + let row = meth_tsv.lines().nth(1).expect("meth.tsv data row"); + let fields: Vec<&str> = row.split('\t').collect(); + let n_cpg: u64 = fields[0].parse().unwrap(); + assert!(n_cpg > 0, "expected covered CpGs"); + assert_ne!(fields[1], "NA", "pearson_r should be defined"); + let r: f64 = fields[1].parse().unwrap(); + assert!(r > 0.8, "golden XM should track truth strongly, got r={r}"); +} + /// Create a BAM where all reads are unmapped. Eval should report 100% /// unmapped. #[test] From c5d177756abb48eced7d1ce37284e6caa16daadc Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Fri, 26 Jun 2026 13:12:17 -0700 Subject: [PATCH 7/9] fix(eval): read each read's alt allele from the golden sequence, not VCF phasing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Variant representation inferred which haplotype carried the alt from the truth VCF genotype (`alt_by_hap[read's haplotype]`), assigning an unphased het alt to haplotype 1 deterministically. But `simulate` assigns unphased genotypes to haplotypes by a random permutation, so eval's guess disagreed with the simulator on roughly half of heterozygous sites — it checked the reference-copy reads for an alt that was actually on the other copy, scoring them as misses. The effect was large: a byte-identical golden-as-mapped run scored ~74% representation instead of 100%, and every aligner inherited the same artifact. Resolve the allele a read truly carries from the golden read's own sequence instead. The golden BAM is the per-read oracle: at each truth substitution site within a read's true span, the golden read's base is exactly what the simulator placed on the copy that read was sequenced from. A read showing the reference base is (correctly) not expected to carry the alt. This makes representation correct whether or not the truth VCF is phased — golden-as-mapped now scores 100%, and a real aligner's shortfall reflects only its actual mismapping. `GoldenInfo` carries the read sequence and CIGAR (new `base_at` resolver); the now-unused `hp:i` haplotype field is dropped from the eval-side record. --- CHANGELOG.md | 14 +-- src/commands/eval/golden.rs | 39 +++++--- src/commands/eval/variants.rs | 165 ++++++++++++++++++++-------------- 3 files changed, 135 insertions(+), 83 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 18ac662..b163b3c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,12 +11,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - `eval` now scores accuracy against holodeck's own truth beyond placement. With `--truth` (the golden BAM) it takes per-read true positions, spans, and - haplotypes from the golden alignment rather than only the encoded read name, + sequences from the golden alignment rather than only the encoded read name, and `--variants` reports how faithfully aligned reads represent the simulated - substitutions: for every expected single-base substitution it walks the - mapped read's CIGAR to the variant position and checks the observed base, - accumulating the represented fraction with the read's MAPQ and alignment - score, plus per-read MD/NM concordance against the golden tags. `--meth` + substitutions. The allele a read truly carries at each truth site is read + from the golden read's own sequence — the per-read oracle — so scoring is + correct whether or not the truth VCF is phased (a read sequenced from the + reference copy shows the reference base and is not expected to carry the + alt). For every such expected substitution it walks the mapped read's CIGAR + to the variant position and checks the observed base, accumulating the + represented fraction with the read's MAPQ and alignment score, plus per-read + MD/NM concordance against the golden tags. `--meth` breaks the results down by bisulfite substitution class, labelling the conversion-confounded `C->T`/`G->A` cell as such. `--cpg-truth` correlates the aligner's Bismark `XM` calls against the simulated cpg-truth bedGraph diff --git a/src/commands/eval/golden.rs b/src/commands/eval/golden.rs index d0fd575..123a314 100644 --- a/src/commands/eval/golden.rs +++ b/src/commands/eval/golden.rs @@ -2,10 +2,11 @@ //! variant / MD-tag concordance scoring. //! //! The golden BAM written by `holodeck simulate --golden-bam` carries the -//! true alignment of every read (MAPQ 60, correct CIGAR), the source -//! haplotype in the `hp:i` tag, and — for methylation runs — Bismark-style -//! `NM:i` / `MD:Z` call tags. This module indexes those records by read end so -//! the eval pass can look up each mapped read's truth in O(1). +//! true alignment of every read (MAPQ 60, correct CIGAR, the sequence the read +//! was given) and — for methylation runs — Bismark-style `NM:i` / `MD:Z` call +//! tags. This module indexes those records by read end so the eval pass can +//! look up each mapped read's truth in O(1), including the actual allele the +//! read carries at any reference position (see [`GoldenInfo::base_at`]). use std::collections::HashMap; use std::path::Path; @@ -14,6 +15,7 @@ use anyhow::{Context, Result}; use bstr::ByteSlice; use noodles::bam; use noodles::sam::alignment::record::data::field::Tag; +use noodles::sam::alignment::record_buf::Cigar; use noodles::sam::alignment::record_buf::data::field::Value; use super::cigar; @@ -31,12 +33,18 @@ pub struct GoldenInfo { pub start0: u32, /// Reference bases consumed by the true alignment. pub ref_len: u32, - /// Source haplotype index (`hp:i` tag; `0` if absent). - pub haplotype: usize, /// `NM:i` edit distance, if present. pub nm: Option, /// `MD:Z` string, if present. pub md: Option, + /// Uppercased read sequence of the true alignment. Paired with `cigar` and + /// `start0`, this is the per-read oracle for which allele a read actually + /// carries at a variant site — making variant-representation scoring + /// independent of whether the truth VCF is phased. + pub sequence: Vec, + /// CIGAR of the true alignment, for mapping a reference position to a read + /// offset within `sequence`. + pub cigar: Cigar, } impl GoldenInfo { @@ -45,6 +53,17 @@ impl GoldenInfo { pub fn end0(&self) -> u32 { self.start0 + self.ref_len } + + /// The uppercased base this read carries at 0-based reference position + /// `ref_pos0`, or `None` when that position is deleted, clipped, or outside + /// the alignment. This is read straight from the golden sequence, so it + /// reflects exactly what the simulator placed on this read (the alt allele + /// for a read sequenced from the alt copy, the reference base otherwise). + #[must_use] + pub fn base_at(&self, ref_pos0: u32) -> Option { + let offset = cigar::ref_pos_to_read_offset(&self.cigar, self.start0, ref_pos0)?; + self.sequence.get(offset).map(u8::to_ascii_uppercase) + } } /// Load every primary, mapped golden record keyed by `(name, is_last_segment)`. @@ -73,15 +92,15 @@ pub fn load(path: &Path) -> Result> { let Some(start) = record.alignment_start() else { continue }; let start0 = u32::try_from(usize::from(start).saturating_sub(1)).unwrap_or(0); + let cigar = record.cigar().clone(); let info = GoldenInfo { contig: contig_name.to_str_lossy().into_owned(), start0, - ref_len: cigar::reference_len(record.cigar()), - haplotype: int_tag(&record, b'h', b'p') - .and_then(|n| usize::try_from(n).ok()) - .unwrap_or(0), + ref_len: cigar::reference_len(&cigar), nm: int_tag(&record, b'N', b'M'), md: string_tag(&record, b'M', b'D'), + sequence: record.sequence().as_ref().to_vec(), + cigar, }; map.insert((name.to_vec(), flags.is_last_segment()), info); } diff --git a/src/commands/eval/variants.rs b/src/commands/eval/variants.rs index 1764888..b2ebcaf 100644 --- a/src/commands/eval/variants.rs +++ b/src/commands/eval/variants.rs @@ -97,25 +97,25 @@ pub fn classify_substitution(ref_base: u8, alt_base: u8, conv: ConvDir) -> SubCl if transition { SubClass::Other } else { SubClass::Transversion } } -/// A single-base substitution a read is expected to carry on its haplotype. +/// A single-base substitution a read actually carries, with the alternate base +/// read from the golden truth sequence (not inferred from VCF phasing). #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub struct ExpectedSnv { /// 0-based reference position. pub pos0: u32, /// Uppercased reference base. pub ref_base: u8, - /// Uppercased alternate base on the queried haplotype. + /// Uppercased alternate base the golden read carries at this site. pub alt_base: u8, } -/// One truth SNV site with its per-haplotype alternate bases. +/// One truth substitution site. Only the position and reference base are kept: +/// which copy carries the alt — and which base — is read per-read from the +/// golden BAM, so the VCF need not be phased. #[derive(Debug, Clone)] struct SnvSite { pos0: u32, ref_base: u8, - /// Expected base per haplotype: `Some(alt)` when that haplotype carries a - /// single-base alternate here, `None` for reference / missing / non-SNV. - alt_by_hap: Vec>, } /// Truth SNVs indexed by contig for per-read span queries. @@ -139,8 +139,10 @@ impl VariantTruth { Ok(Self::from_parsed(&parsed)) } - /// Build SNV truth from already-parsed variants, keeping only sites that - /// are a single-base substitution on at least one haplotype. + /// Build SNV truth from already-parsed variants, keeping each site that is + /// a single-base substitution the sample carries (single-base reference and + /// at least one single-base ALT in the genotype). The realized allele per + /// read is resolved later from the golden sequence, so no phasing is stored. fn from_parsed(parsed: &ParsedVariants) -> Self { let mut by_contig: BTreeMap> = BTreeMap::new(); for (contig, records) in &parsed.by_contig { @@ -149,26 +151,17 @@ impl VariantTruth { if record.ref_allele.len() != 1 { continue; // SNV requires a single reference base. } - let ref_base = record.ref_allele[0].to_ascii_uppercase(); - let mut any = false; - let alt_by_hap: Vec> = record - .genotype - .alleles() - .iter() - .map(|allele| { - let alt = match allele { - Some(idx) if *idx > 0 => record - .allele_bases(*idx) - .filter(|b| b.len() == 1) - .map(|b| b[0].to_ascii_uppercase()), - _ => None, - }; - any |= alt.is_some(); - alt - }) - .collect(); - if any { - sites.push(SnvSite { pos0: record.position, ref_base, alt_by_hap }); + // Keep only sites the sample actually carries as a substitution: + // some genotype allele indexes a single-base ALT. + let carries_snv_alt = record.genotype.alleles().iter().any(|allele| { + matches!(allele, Some(idx) if *idx > 0 + && record.allele_bases(*idx).is_some_and(|b| b.len() == 1)) + }); + if carries_snv_alt { + sites.push(SnvSite { + pos0: record.position, + ref_base: record.ref_allele[0].to_ascii_uppercase(), + }); } } sites.sort_by_key(|s| s.pos0); @@ -179,26 +172,30 @@ impl VariantTruth { Self { by_contig } } - /// Expected SNVs for a read on `haplotype` spanning `[start0, end0)`. + /// Substitutions a read actually carries, read from the golden truth. + /// + /// For every truth SNV site within the read's true span, the golden read's + /// own base at that site is the per-read oracle: when it differs from the + /// reference, this read carries that alternate (the value the simulator + /// placed on the copy this read was sequenced from). A read sequenced from + /// the reference copy shows the reference base and yields nothing — so the + /// result is correct whether or not the truth VCF is phased. #[must_use] - pub fn expected_snvs( - &self, - contig: &str, - haplotype: usize, - start0: u32, - end0: u32, - ) -> Vec { - let Some(sites) = self.by_contig.get(contig) else { + pub fn expected_for_read(&self, golden: &GoldenInfo) -> Vec { + let Some(sites) = self.by_contig.get(&golden.contig) else { return Vec::new(); }; - let lo = sites.partition_point(|s| s.pos0 < start0); + let end0 = golden.end0(); + let lo = sites.partition_point(|s| s.pos0 < golden.start0); let mut out = Vec::new(); for site in &sites[lo..] { if site.pos0 >= end0 { break; } - if let Some(Some(alt)) = site.alt_by_hap.get(haplotype) { - out.push(ExpectedSnv { pos0: site.pos0, ref_base: site.ref_base, alt_base: *alt }); + if let Some(base) = golden.base_at(site.pos0) + && base != site.ref_base + { + out.push(ExpectedSnv { pos0: site.pos0, ref_base: site.ref_base, alt_base: base }); } } out @@ -357,12 +354,7 @@ pub fn run( let key: ReadKey = (name.to_vec(), flags.is_last_segment()); let Some(truth_aln) = golden.get(&key) else { continue }; - let expected = truth.expected_snvs( - &truth_aln.contig, - truth_aln.haplotype, - truth_aln.start0, - truth_aln.end0(), - ); + let expected = truth.expected_for_read(truth_aln); if expected.is_empty() { continue; } @@ -415,6 +407,9 @@ pub fn run( #[cfg(test)] mod tests { + use noodles::sam::alignment::record::cigar::op::{Kind, Op}; + use noodles::sam::alignment::record_buf::Cigar; + use super::*; use crate::vcf::genotype::{Genotype, VariantRecord}; @@ -457,37 +452,71 @@ mod tests { VariantTruth::from_parsed(&parsed) } + fn golden_read(start0: u32, seq: &[u8], ops: &[(Kind, usize)]) -> GoldenInfo { + let cigar = Cigar::from(ops.iter().map(|&(k, n)| Op::new(k, n)).collect::>()); + GoldenInfo { + contig: "chr1".to_string(), + start0, + ref_len: cigar::reference_len(&cigar), + nm: None, + md: None, + sequence: seq.to_vec(), + cigar, + } + } + #[test] - fn expected_snvs_respects_phasing_per_haplotype() { - // 1|0 -> haplotype 0 carries the alt, haplotype 1 does not. - let truth = truth_from(vec![snv_record(100, "C", "T", "1|0")]); - let hap0 = truth.expected_snvs("chr1", 0, 0, 200); - assert_eq!(hap0.len(), 1); - assert_eq!((hap0[0].pos0, hap0[0].ref_base, hap0[0].alt_base), (100, b'C', b'T')); - assert!(truth.expected_snvs("chr1", 1, 0, 200).is_empty()); + fn expected_reads_alt_from_golden_even_when_vcf_is_unphased() { + // SNV at ref 100 (ref C). The VCF is UNPHASED (0/1) — eval must not guess + // a haplotype; it reads the read's actual base from the golden sequence. + let truth = truth_from(vec![snv_record(100, "C", "T", "0/1")]); + // Golden read covering [90, 200) as 110M; base at ref 100 = offset 10 = T. + let mut seq = vec![b'A'; 110]; + seq[10] = b'T'; + let exp = truth.expected_for_read(&golden_read(90, &seq, &[(Kind::Match, 110)])); + assert_eq!(exp.len(), 1); + assert_eq!((exp[0].pos0, exp[0].ref_base, exp[0].alt_base), (100, b'C', b'T')); } #[test] - fn expected_snvs_skips_indels_and_honors_span() { + fn expected_skips_reads_carrying_the_reference_allele() { + // Same unphased site, but this golden read shows the REFERENCE base C at + // 100 — sequenced from the reference copy, so nothing is expected. + let truth = truth_from(vec![snv_record(100, "C", "T", "0/1")]); + let mut seq = vec![b'A'; 110]; + seq[10] = b'C'; + assert!(truth.expected_for_read(&golden_read(90, &seq, &[(Kind::Match, 110)])).is_empty()); + } + + #[test] + fn expected_honors_span_and_skips_indels() { let truth = truth_from(vec![ - snv_record(50, "A", "G", "1|1"), - snv_record(100, "AT", "A", "1|1"), // deletion: not an SNV - snv_record(150, "C", "A", "0|1"), + snv_record(50, "A", "G", "1/1"), + snv_record(100, "AT", "A", "1/1"), // deletion: not an SNV site + snv_record(150, "C", "A", "0/1"), ]); - // Span [60, 200) excludes pos 50; indel at 100 dropped; 150 on hap1 kept. - let hap1 = truth.expected_snvs("chr1", 1, 60, 200); - assert_eq!(hap1.len(), 1); - assert_eq!(hap1[0].pos0, 150); - // hap0 in [0,60) sees only pos 50. - let hap0 = truth.expected_snvs("chr1", 0, 0, 60); - assert_eq!(hap0.len(), 1); - assert_eq!(hap0[0].pos0, 50); + // Golden read [60, 200) as 140M: pos 50 is before the span, the indel at + // 100 is dropped, only site 150 (offset 90 = alt A) remains. + let mut seq = vec![b'C'; 140]; + seq[90] = b'A'; + let exp = truth.expected_for_read(&golden_read(60, &seq, &[(Kind::Match, 140)])); + assert_eq!(exp.len(), 1); + assert_eq!((exp[0].pos0, exp[0].alt_base), (150, b'A')); } #[test] - fn expected_snvs_unknown_contig_is_empty() { - let truth = truth_from(vec![snv_record(10, "C", "T", "1|0")]); - assert!(truth.expected_snvs("chrX", 0, 0, 1000).is_empty()); + fn expected_for_read_on_unknown_contig_is_empty() { + let truth = truth_from(vec![snv_record(10, "C", "T", "0/1")]); + let g = GoldenInfo { + contig: "chrX".to_string(), + start0: 0, + ref_len: 100, + nm: None, + md: None, + sequence: vec![b'T'; 100], + cigar: Cigar::from(vec![Op::new(Kind::Match, 100)]), + }; + assert!(truth.expected_for_read(&g).is_empty()); } #[test] From 4e730066826f956de8c0eac0d92f204d35cd5174 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Fri, 26 Jun 2026 13:17:17 -0700 Subject: [PATCH 8/9] fix(eval): tolerate aligners without XM tags instead of failing --meth MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `--meth` methylation correlation aborted the whole eval when a mapped BAM had no XM tags. But a plain bisulfite aligner (e.g. bwameth) emits no per-base methylation calls at all — calling is a separate extractor (MethylDackel) step — so its BAM legitimately has no XM, and bailing made it impossible to score such an aligner's placement and variant representation alongside callers that do emit XM. Warn and report NA methylation correlation (n_cpg 0) instead of erroring, so the placement and variant-representation axes still complete for these BAMs. --- src/commands/eval/meth.rs | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/commands/eval/meth.rs b/src/commands/eval/meth.rs index 96800df..6d28dc6 100644 --- a/src/commands/eval/meth.rs +++ b/src/commands/eval/meth.rs @@ -97,10 +97,17 @@ fn tally_aligner(mapped: &Path) -> Result { tally_read(&contig, record.cigar(), start0, xm.as_bytes(), &mut tally); } - // A BAM with mapped reads but no XM tags cannot be evaluated for - // methylation; surface that rather than reporting an empty correlation. + // A BAM with mapped reads but no XM tags carries no methylation calls — for + // example a plain bisulfite aligner like bwameth, where calling is a + // separate extractor (MethylDackel) step. Warn and return the empty tally so + // the correlation is reported as NA rather than failing the whole eval; + // placement and variant representation are still meaningful for such a BAM. if saw_mapped_primary && !saw_xm { - bail!("no XM methylation tags found in mapped primary records of {}", mapped.display()); + log::warn!( + "no XM methylation tags in mapped primaries of {}; \ + reporting NA methylation correlation", + mapped.display() + ); } Ok(tally) } From 6d2fa7fe9ae3948aa8e4cda2d0930f0cab510c1e Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Fri, 26 Jun 2026 13:50:34 -0700 Subject: [PATCH 9/9] feat(eval): bisulfite-aware genomic NM/MD concordance via --reference MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit NM/MD concordance compared the aligner's raw NM:i/MD:Z tags against the golden tags. Those tags are convention-dependent for bisulfite data: an aligner may score edits against the original 4-letter reference (every unmethylated C->T is a "mismatch"), against a C->T-converted reference (conversions match), or in a bisulfite-aware convention. Two correctly-placed reads therefore disagree on their tags purely by convention — D3 (scoring vs the original reference) matched the golden ~1% of the time while bwameth (vs the converted reference) matched ~75%, measuring which convention each picked rather than alignment quality. Add `--reference ` and recompute concordance as a convention-independent genomic edit distance: walk each read against the reference, exclude bisulfite conversions using the read's TRUE strand (from the golden truth, so it works for aligners like bwameth that emit no XG), and compare the resulting non-conversion mismatch + indel profiles of the aligned and golden reads. D3 and bwameth now both land at ~96-100% on the same scale. The new `edits` module holds the genomic-edit walk and an on-demand reference-contig cache; the golden record now carries its true conversion strand instead of the raw NM/MD tags. Without `--reference`, NM/MD concordance is reported as NA. --- CHANGELOG.md | 22 ++-- src/commands/eval/edits.rs | 213 ++++++++++++++++++++++++++++++++++ src/commands/eval/golden.rs | 74 +++++++++--- src/commands/eval/mod.rs | 35 +++++- src/commands/eval/variants.rs | 154 +++++++++++++----------- tests/test_eval.rs | 80 +++++++++++++ 6 files changed, 487 insertions(+), 91 deletions(-) create mode 100644 src/commands/eval/edits.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index b163b3c..e4ae6b3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,13 +19,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 reference copy shows the reference base and is not expected to carry the alt). For every such expected substitution it walks the mapped read's CIGAR to the variant position and checks the observed base, accumulating the - represented fraction with the read's MAPQ and alignment score, plus per-read - MD/NM concordance against the golden tags. `--meth` - breaks the results down by bisulfite substitution class, labelling the - conversion-confounded `C->T`/`G->A` cell as such. `--cpg-truth` correlates - the aligner's Bismark `XM` calls against the simulated cpg-truth bedGraph - (Pearson r and RMSE). Results are written to `.variants.tsv` and - `.meth.tsv` alongside the existing `.eval.txt`. + represented fraction with the read's MAPQ and alignment score. With + `--reference` it also reports per-read NM/MD concordance as a **bisulfite-aware + genomic edit distance**: rather than comparing raw `NM:i`/`MD:Z` tags (which + are convention-dependent — a bisulfite aligner may score against the original + or the converted reference, so the tags differ even when both placed the read + correctly), it recomputes each read's edits against the reference and excludes + conversions using the read's TRUE strand (taken from the golden truth, so it + works even for aligners such as bwameth that emit no `XG`). The result is + comparable across aligners; without `--reference` NM/MD concordance is `NA`. + `--meth` breaks the variant results down by bisulfite substitution class, + labelling the conversion-confounded `C->T`/`G->A` cell as such. `--cpg-truth` + correlates the aligner's Bismark `XM` calls against the simulated cpg-truth + bedGraph (Pearson r and RMSE; `NA` for aligners that emit no `XM`). Results + are written to `.variants.tsv` and `.meth.tsv` alongside the + existing `.eval.txt`. ### Fixed diff --git a/src/commands/eval/edits.rs b/src/commands/eval/edits.rs new file mode 100644 index 0000000..1782c59 --- /dev/null +++ b/src/commands/eval/edits.rs @@ -0,0 +1,213 @@ +//! Bisulfite-aware genomic edit distance for truth-vs-aligner concordance. +//! +//! `NM:i` / `MD:Z` tags are reference- and convention-dependent: a bisulfite +//! aligner may report edits against the original 4-letter reference (every +//! unmethylated `C->T` is a "mismatch"), against a `C->T`-converted reference +//! (conversions match), or in a bisulfite-aware convention. Comparing two +//! aligners' raw tags therefore measures which convention each picked, not +//! whether they aligned correctly. +//! +//! Instead, [`genomic_edits`] recomputes a read's edits directly against the +//! reference and excludes bisulfite conversions using the read's TRUE strand +//! (taken from the golden truth, so it works even for aligners that emit no +//! `XG`). The result — non-conversion mismatches plus indels — is a +//! convention-independent genomic edit distance that is comparable across +//! aligners. [`RefCache`] loads reference contigs on demand to support it. + +use std::collections::BTreeSet; +use std::collections::HashMap; + +use noodles::sam::alignment::record::cigar::op::Kind; +use noodles::sam::alignment::record_buf::Cigar; +use rand::SeedableRng; +use rand::rngs::SmallRng; + +use super::golden::ConvDir; +use crate::fasta::Fasta; +use crate::seed::compute_seed; + +/// A read's genomic edits against the reference: mismatches and indels that are +/// NOT explained by bisulfite conversion. Independent of the aligner's NM/MD +/// tagging convention. +#[derive(Debug, Default, Clone, PartialEq, Eq)] +pub struct GenomicEdits { + /// Edit distance: non-conversion mismatches + inserted + deleted bases + /// (the bisulfite-aware analogue of `NM:i`). + pub nm: u32, + /// Reference positions of non-conversion mismatches and deletions (the + /// bisulfite-aware analogue of the `MD:Z`-encoded edit set). Insertions + /// carry no reference position and so contribute to `nm` only. + pub positions: BTreeSet, +} + +/// Whether `read_base` at a reference `ref_base` is a valid bisulfite +/// conversion under the read's true strand `conv` (and therefore not a genomic +/// edit): `C->T` on the CT strand, `G->A` on the GA strand. +fn is_conversion(ref_base: u8, read_base: u8, conv: Option) -> bool { + match conv { + Some(ConvDir::Ct) => ref_base == b'C' && read_base == b'T', + Some(ConvDir::Ga) => ref_base == b'G' && read_base == b'A', + None => false, + } +} + +/// Compute a read's genomic edits against `ref_seq` (the full contig the read +/// aligns to), walking its `cigar` from 0-based `start0`. `conv` is the read's +/// TRUE bisulfite strand; conversions consistent with it are excluded so the +/// result reflects only genuine genomic differences. +#[must_use] +pub fn genomic_edits( + seq: &[u8], + cigar: &Cigar, + start0: u32, + ref_seq: &[u8], + conv: Option, +) -> GenomicEdits { + let mut edits = GenomicEdits::default(); + let mut ref_pos = start0 as usize; + let mut read_pos = 0usize; + for op in cigar.as_ref() { + let len = op.len(); + match op.kind() { + Kind::Match | Kind::SequenceMatch | Kind::SequenceMismatch => { + for k in 0..len { + let r = ref_seq.get(ref_pos + k).map(u8::to_ascii_uppercase); + let q = seq.get(read_pos + k).map(u8::to_ascii_uppercase); + if let (Some(r), Some(q)) = (r, q) + && r != q + && !is_conversion(r, q, conv) + { + edits.nm += 1; + edits.positions.insert(u32::try_from(ref_pos + k).unwrap_or(0)); + } + } + ref_pos += len; + read_pos += len; + } + Kind::Insertion => { + edits.nm += u32::try_from(len).unwrap_or(0); + read_pos += len; + } + Kind::Deletion => { + edits.nm += u32::try_from(len).unwrap_or(0); + for k in 0..len { + edits.positions.insert(u32::try_from(ref_pos + k).unwrap_or(0)); + } + ref_pos += len; + } + Kind::Skip => ref_pos += len, + Kind::SoftClip => read_pos += len, + Kind::HardClip | Kind::Pad => {} + } + } + edits +} + +/// On-demand cache of uppercased reference contig sequences over an indexed +/// FASTA, so the eval pass loads only the contigs its reads actually touch. +pub struct RefCache { + fasta: Fasta, + cache: HashMap>>, +} + +impl RefCache { + /// Wrap an opened reference FASTA. + #[must_use] + pub fn new(fasta: Fasta) -> Self { + Self { fasta, cache: HashMap::new() } + } + + /// The uppercased sequence of `contig`, loaded and cached on first use. + /// Returns `None` (cached) when the contig is absent from the reference. + pub fn contig(&mut self, contig: &str) -> Option<&[u8]> { + if !self.cache.contains_key(contig) { + // IUPAC ambiguity codes (rare, and not at variant sites) are + // resolved randomly; seed a fresh RNG per contig from a deterministic + // FNV-1a hash of the contig name so the resolution does not depend on + // the order contigs are first requested. + let mut rng = SmallRng::seed_from_u64(compute_seed(contig)); + // Cache the attempt either way (so a failure is not retried). Warn + // on a genuine load error rather than silently dropping the contig's + // reads from NM/MD concordance — it is indistinguishable from an + // absent contig at the call site otherwise. + let loaded = match self.fasta.load_contig(contig, &mut rng) { + Ok(seq) => Some(seq), + Err(e) => { + log::warn!( + "reference contig {contig:?} unavailable ({e:#}); \ + its reads are excluded from NM/MD concordance" + ); + None + } + }; + self.cache.insert(contig.to_string(), loaded); + } + self.cache.get(contig).and_then(Option::as_deref) + } +} + +#[cfg(test)] +mod tests { + use noodles::sam::alignment::record::cigar::op::Op; + + use super::*; + + fn cigar(ops: &[(Kind, usize)]) -> Cigar { + Cigar::from(ops.iter().map(|&(k, n)| Op::new(k, n)).collect::>()) + } + + #[test] + fn counts_a_plain_mismatch_as_a_genomic_edit() { + // read T vs ref A at ref pos 12 (offset 2 in a 10M from start 10). + let edits = genomic_edits( + b"AATAAAAAAA", + &cigar(&[(Kind::Match, 10)]), + 10, + b"AAAAAAAAAAAAAAA", + None, + ); + assert_eq!(edits.nm, 1); + assert_eq!(edits.positions.iter().copied().collect::>(), vec![12]); + } + + #[test] + fn excludes_bisulfite_conversion_on_its_strand() { + // ref C, read T at every position. On the CT strand these are all + // conversions (no genomic edits); with no strand they are all edits. + let seq = b"TTTT"; + let cig = cigar(&[(Kind::Match, 4)]); + let reference = b"CCCC"; + assert_eq!(genomic_edits(seq, &cig, 0, reference, Some(ConvDir::Ct)).nm, 0); + assert_eq!(genomic_edits(seq, &cig, 0, reference, None).nm, 4); + // On the GA strand a C->T is NOT the freed cell, so it stays an edit. + assert_eq!(genomic_edits(seq, &cig, 0, reference, Some(ConvDir::Ga)).nm, 4); + } + + #[test] + fn a_real_variant_survives_conversion_masking() { + // ref C, read A (a transversion variant) on the CT strand: A is not the + // conversion product T, so it is a genuine genomic edit. + let edits = + genomic_edits(b"A", &cigar(&[(Kind::Match, 1)]), 5, b"CCCCCCC", Some(ConvDir::Ct)); + assert_eq!(edits.nm, 1); + assert!(edits.positions.contains(&5)); + } + + #[test] + fn counts_indels_in_nm_and_deletions_in_positions() { + // 2M1I2M1D2M over ref AAAAAAA: insertion adds to nm only; deletion adds + // to nm and contributes its reference position. + let cig = cigar(&[ + (Kind::Match, 2), + (Kind::Insertion, 1), + (Kind::Match, 2), + (Kind::Deletion, 1), + (Kind::Match, 2), + ]); + // read consumes 2+1+2+0+2 = 7 bases, all matching ref where aligned. + let edits = genomic_edits(b"AAAAAAA", &cig, 0, b"AAAAAAAAAA", None); + assert_eq!(edits.nm, 2); // 1 inserted + 1 deleted + // deletion reference position: 2M(0,1) 2M(2,3) D at ref 4. + assert!(edits.positions.contains(&4)); + } +} diff --git a/src/commands/eval/golden.rs b/src/commands/eval/golden.rs index 123a314..48b309b 100644 --- a/src/commands/eval/golden.rs +++ b/src/commands/eval/golden.rs @@ -3,15 +3,16 @@ //! //! The golden BAM written by `holodeck simulate --golden-bam` carries the //! true alignment of every read (MAPQ 60, correct CIGAR, the sequence the read -//! was given) and — for methylation runs — Bismark-style `NM:i` / `MD:Z` call -//! tags. This module indexes those records by read end so the eval pass can +//! was given) and — for methylation runs — the Bismark `XG` genome-conversion +//! strand. This module indexes those records by read end so the eval pass can //! look up each mapped read's truth in O(1), including the actual allele the -//! read carries at any reference position (see [`GoldenInfo::base_at`]). +//! read carries at any reference position (see [`GoldenInfo::base_at`]) and the +//! read's true bisulfite strand (see [`GoldenInfo::conv_dir`]). use std::collections::HashMap; use std::path::Path; -use anyhow::{Context, Result}; +use anyhow::{Context, Result, bail}; use bstr::ByteSlice; use noodles::bam; use noodles::sam::alignment::record::data::field::Tag; @@ -24,6 +25,27 @@ use super::cigar; /// (R2). R1 and single-end reads use `false`. pub type ReadKey = (Vec, bool); +/// Bisulfite conversion direction of a read, from the Bismark `XG`/`XR` tags. +/// On the golden record this is ground truth (the strand the simulator drew); +/// on an aligner record it is the aligner's own call. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum ConvDir { + /// `CT` strand: `C->T` is the converted (freed) cell. + Ct, + /// `GA` strand: `G->A` is the converted (freed) cell. + Ga, +} + +/// Parse a record's bisulfite strand from its `XG` tag, falling back to `XR`. +pub(super) fn conv_dir_from_tags(record: &noodles::sam::alignment::RecordBuf) -> Option { + let tag = string_tag(record, b'X', b'G').or_else(|| string_tag(record, b'X', b'R'))?; + match tag.as_str() { + "CT" => Some(ConvDir::Ct), + "GA" => Some(ConvDir::Ga), + _ => None, + } +} + /// True alignment for one read end, taken from the golden BAM. #[derive(Debug, Clone)] pub struct GoldenInfo { @@ -33,10 +55,6 @@ pub struct GoldenInfo { pub start0: u32, /// Reference bases consumed by the true alignment. pub ref_len: u32, - /// `NM:i` edit distance, if present. - pub nm: Option, - /// `MD:Z` string, if present. - pub md: Option, /// Uppercased read sequence of the true alignment. Paired with `cigar` and /// `start0`, this is the per-read oracle for which allele a read actually /// carries at a variant site — making variant-representation scoring @@ -45,6 +63,9 @@ pub struct GoldenInfo { /// CIGAR of the true alignment, for mapping a reference position to a read /// offset within `sequence`. pub cigar: Cigar, + /// True bisulfite strand (golden `XG`), or `None` for non-meth reads. Used + /// to exclude bisulfite conversions when computing genomic edit distance. + pub conv_dir: Option, } impl GoldenInfo { @@ -84,12 +105,25 @@ pub fn load(path: &Path) -> Result> { continue; } - let Some(name) = record.name() else { continue }; - let Some(ref_id) = record.reference_sequence_id() else { continue }; + // The golden BAM is holodeck's own truth; a primary mapped record + // missing a required field, or a duplicated read-end, means corrupt or + // non-holodeck input. Fail fast rather than silently skipping or + // overwriting, which would quietly mis-score the eval. + let Some(name) = record.name() else { + bail!("Golden BAM has a primary mapped record with no read name"); + }; + let Some(ref_id) = record.reference_sequence_id() else { + bail!("Golden BAM record {} has no reference sequence id", name.to_str_lossy()); + }; let Some((contig_name, _)) = header.reference_sequences().get_index(ref_id) else { - continue; + bail!( + "Golden BAM record {} references unknown sequence id {ref_id}", + name.to_str_lossy() + ); + }; + let Some(start) = record.alignment_start() else { + bail!("Golden BAM record {} has no alignment start", name.to_str_lossy()); }; - let Some(start) = record.alignment_start() else { continue }; let start0 = u32::try_from(usize::from(start).saturating_sub(1)).unwrap_or(0); let cigar = record.cigar().clone(); @@ -97,12 +131,20 @@ pub fn load(path: &Path) -> Result> { contig: contig_name.to_str_lossy().into_owned(), start0, ref_len: cigar::reference_len(&cigar), - nm: int_tag(&record, b'N', b'M'), - md: string_tag(&record, b'M', b'D'), - sequence: record.sequence().as_ref().to_vec(), + // Uppercase at ingestion so `sequence` honors its documented + // contract for every consumer (noodles usually decodes uppercase, + // but normalize defensively rather than rely on it). + sequence: record.sequence().as_ref().iter().map(u8::to_ascii_uppercase).collect(), cigar, + conv_dir: conv_dir_from_tags(&record), }; - map.insert((name.to_vec(), flags.is_last_segment()), info); + let segment = if flags.is_last_segment() { "R2" } else { "R1" }; + if map.insert((name.to_vec(), flags.is_last_segment()), info).is_some() { + bail!( + "Golden BAM has a duplicate primary {segment} record for read {}", + name.to_str_lossy() + ); + } } Ok(map) diff --git a/src/commands/eval/mod.rs b/src/commands/eval/mod.rs index 8e434f0..7deb16b 100644 --- a/src/commands/eval/mod.rs +++ b/src/commands/eval/mod.rs @@ -13,6 +13,7 @@ //! `.meth.tsv`). mod cigar; +mod edits; mod golden; mod meth; mod placement; @@ -20,11 +21,13 @@ mod variants; use std::path::PathBuf; -use anyhow::{Result, bail}; +use anyhow::{Context, Result, bail}; use clap::Parser; use super::command::Command; use super::common::OutputPrefixOptions; +use crate::fasta::Fasta; +use edits::RefCache; /// Evaluate alignment accuracy of simulated reads. /// @@ -46,7 +49,7 @@ pub struct Eval { pub mapped: PathBuf, /// Golden BAM (`simulate --golden-bam`) supplying each read's true span, - /// haplotype, and MD/NM tags. Required by `--variants`. + /// sequence, and (for meth) bisulfite strand. Required by `--variants`. #[arg(long, value_name = "BAM")] pub truth: Option, @@ -71,6 +74,14 @@ pub struct Eval { #[arg(long, value_name = "BEDGRAPH")] pub cpg_truth: Option, + /// Reference FASTA (indexed). Enables bisulfite-aware genomic NM/MD + /// concordance under `--variants`: each read's edits are recomputed against + /// the reference and conversions excluded via its true strand, so the + /// metric is comparable across aligners regardless of their NM/MD + /// convention. Without it NM/MD concordance is reported as NA. + #[arg(short = 'r', long, value_name = "FASTA")] + pub reference: Option, + #[command(flatten)] pub output: OutputPrefixOptions, @@ -103,7 +114,25 @@ impl Command for Eval { // Safe: the guard above rejects --variants without --truth. let golden = golden.as_ref().expect("--variants requires --truth"); let truth = variants::VariantTruth::from_vcf(vcf, self.sample.as_deref())?; - variants::run(&self.mapped, golden, &truth, self.meth, &self.output.output)?; + let mut reference = if let Some(path) = &self.reference { + Some(RefCache::new( + Fasta::from_path(path) + .with_context(|| format!("Failed to open reference {}", path.display()))?, + )) + } else { + log::warn!( + "--reference not given; NM/MD genomic-edit concordance will be reported as NA" + ); + None + }; + variants::run( + &self.mapped, + golden, + &truth, + self.meth, + reference.as_mut(), + &self.output.output, + )?; } if let Some(cpg_truth) = &self.cpg_truth { diff --git a/src/commands/eval/variants.rs b/src/commands/eval/variants.rs index b2ebcaf..3326c93 100644 --- a/src/commands/eval/variants.rs +++ b/src/commands/eval/variants.rs @@ -1,13 +1,15 @@ //! Variant-representation accuracy: do aligned reads carry the simulated //! variants they should, and how confidently? //! -//! Truth comes entirely from holodeck's own outputs: the per-haplotype phased -//! genotypes in the truth VCF say which single-base substitutions a read on a -//! given haplotype should carry, and the golden BAM gives each read's true -//! span and haplotype. For every expected substitution this pass walks the -//! *mapped* read's CIGAR to the variant's reference position and checks whether -//! the observed base matches the alternate allele, accumulating the represented -//! fraction together with the read's `MAPQ` and `AS` per substitution class. +//! Truth comes entirely from holodeck's own outputs: the truth VCF enumerates +//! the simulated single-base substitutions and the golden BAM gives each read's +//! true span and the actual base it carries at each site (so whether a read is +//! expected to show the alt is read from the golden sequence, independent of +//! whether the VCF is phased). For every expected substitution this pass walks +//! the *mapped* read's CIGAR to the variant's reference position and checks +//! whether the observed base matches the alternate allele, accumulating the +//! represented fraction together with the read's `MAPQ` and `AS` per +//! substitution class. //! //! ## Methylation framing //! @@ -17,8 +19,10 @@ //! rather than treated as a true accuracy signal. The discriminating classes //! are the mirror (`T->C`) and the transversions, where a methylation-aware //! scoring mode should neither over- nor under-penalize relative to the -//! genomic truth. Classes are assigned from the read's conversion direction -//! (`XG`, falling back to `XR`). +//! genomic truth. Classes are assigned from the read's TRUE conversion +//! direction, taken from the golden BAM (its `XG`) rather than the mapped +//! record's tags, so an aligner that omits or rewrites `XG`/`XR` cannot shift a +//! variant into the wrong class. use std::collections::{BTreeMap, HashMap}; use std::io::Write; @@ -26,23 +30,14 @@ use std::path::Path; use anyhow::{Context, Result}; use noodles::bam; -use noodles::sam::alignment::RecordBuf; use super::cigar; -use super::golden::{GoldenInfo, ReadKey, contig_name, int_tag, string_tag}; +use super::edits::{self, RefCache}; +use super::golden::{ConvDir, GoldenInfo, ReadKey, contig_name, int_tag}; use crate::commands::command::output_path; use crate::sequence_dict::SequenceDictionary; use crate::vcf::{ParsedVariants, parse_variants_by_contig}; -/// Bisulfite conversion direction for a read, from the `XG`/`XR` Bismark tags. -#[derive(Debug, Clone, Copy, PartialEq, Eq)] -pub enum ConvDir { - /// `CT` strand: `C->T` is the converted (freed) cell. - Ct, - /// `GA` strand: `G->A` is the converted (freed) cell. - Ga, -} - /// Classification of a single-base substitution under a conversion direction. #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub enum SubClass { @@ -109,13 +104,17 @@ pub struct ExpectedSnv { pub alt_base: u8, } -/// One truth substitution site. Only the position and reference base are kept: -/// which copy carries the alt — and which base — is read per-read from the -/// golden BAM, so the VCF need not be phased. +/// One truth substitution site. The position and reference base, plus the +/// single-base ALT base(s) the sample's genotype declares here (uppercased). +/// Which copy carries the alt — and thus which read shows it — is read per-read +/// from the golden BAM, so the VCF need not be phased; the ALT set is retained +/// only so a bisulfite conversion or sequencing error that yields some *other* +/// non-reference base is not mistaken for the variant. #[derive(Debug, Clone)] struct SnvSite { pos0: u32, ref_base: u8, + alts: Vec, } /// Truth SNVs indexed by contig for per-read span queries. @@ -151,16 +150,28 @@ impl VariantTruth { if record.ref_allele.len() != 1 { continue; // SNV requires a single reference base. } - // Keep only sites the sample actually carries as a substitution: - // some genotype allele indexes a single-base ALT. - let carries_snv_alt = record.genotype.alleles().iter().any(|allele| { - matches!(allele, Some(idx) if *idx > 0 - && record.allele_bases(*idx).is_some_and(|b| b.len() == 1)) - }); - if carries_snv_alt { + // The single-base ALT base(s) the genotype carries here. Keep + // the site only if it has at least one (i.e. the sample carries + // a substitution); skip ref-only and non-SNV-ALT records. + let mut alts: Vec = record + .genotype + .alleles() + .iter() + .filter_map(|allele| match allele { + Some(idx) if *idx > 0 => record + .allele_bases(*idx) + .filter(|b| b.len() == 1) + .map(|b| b[0].to_ascii_uppercase()), + _ => None, + }) + .collect(); + alts.sort_unstable(); + alts.dedup(); + if !alts.is_empty() { sites.push(SnvSite { pos0: record.position, ref_base: record.ref_allele[0].to_ascii_uppercase(), + alts, }); } } @@ -175,11 +186,13 @@ impl VariantTruth { /// Substitutions a read actually carries, read from the golden truth. /// /// For every truth SNV site within the read's true span, the golden read's - /// own base at that site is the per-read oracle: when it differs from the - /// reference, this read carries that alternate (the value the simulator - /// placed on the copy this read was sequenced from). A read sequenced from - /// the reference copy shows the reference base and yields nothing — so the - /// result is correct whether or not the truth VCF is phased. + /// own base at that site is the per-read oracle: when it matches one of the + /// site's declared ALT bases, this read carries that variant (the value the + /// simulator placed on the copy it was sequenced from). A read sequenced + /// from the reference copy shows the reference base and yields nothing — so + /// the result is correct whether or not the truth VCF is phased. Requiring + /// an ALT match (rather than merely "non-reference") keeps a bisulfite + /// conversion or sequencing error at the site from posing as the variant. #[must_use] pub fn expected_for_read(&self, golden: &GoldenInfo) -> Vec { let Some(sites) = self.by_contig.get(&golden.contig) else { @@ -193,7 +206,7 @@ impl VariantTruth { break; } if let Some(base) = golden.base_at(site.pos0) - && base != site.ref_base + && site.alts.contains(&base) { out.push(ExpectedSnv { pos0: site.pos0, ref_base: site.ref_base, alt_base: base }); } @@ -317,16 +330,6 @@ fn pct_or_na(numerator: u64, denominator: u64) -> String { } } -/// Conversion direction for a mapped record from its `XG` (or `XR`) tag. -fn conv_dir(record: &RecordBuf) -> Option { - let tag = string_tag(record, b'X', b'G').or_else(|| string_tag(record, b'X', b'R'))?; - match tag.as_str() { - "CT" => Some(ConvDir::Ct), - "GA" => Some(ConvDir::Ga), - _ => None, - } -} - /// Evaluate variant representation of `mapped` against `golden` + `truth`. /// /// # Errors @@ -336,6 +339,7 @@ pub fn run( golden: &HashMap, truth: &VariantTruth, meth: bool, + mut reference: Option<&mut RefCache>, output_prefix: &Path, ) -> Result<()> { let mut reader = bam::io::reader::Builder @@ -362,7 +366,7 @@ pub fn run( let mapq = record.mapping_quality().map_or(0, u8::from); let as_score = int_tag(&record, b'A', b'S'); - let conv = if meth { conv_dir(&record) } else { None }; + let conv = if meth { truth_aln.conv_dir } else { None }; // The mapped record represents a variant only if it is aligned to the // variant's contig; otherwise (unmapped / mismapped) it cannot. @@ -380,24 +384,46 @@ pub fn run( } _ => false, }; - // In meth mode the class needs the read's conversion direction; a - // read lacking XG/XR yields None and is counted as unclassified. + // In meth mode the class needs the read's TRUE conversion direction + // (from the golden truth, so an aligner that omits or rewrites XG/XR + // cannot move a variant into the wrong class); a non-meth read yields + // None and is counted as unclassified. let class = conv.map(|dir| classify_substitution(snv.ref_base, snv.alt_base, dir)); report.record(class, meth, represented, mapq, as_score); } - // MD/NM concordance against the golden truth tags for this read, - // counted only where the golden record carries the tag. - if let Some(golden_nm) = truth_aln.nm { - report.nm_comparable_reads += 1; - if int_tag(&record, b'N', b'M') == Some(golden_nm) { - report.nm_concordant_reads += 1; - } - } - if let Some(golden_md) = truth_aln.md.as_deref() { - report.md_comparable_reads += 1; - if string_tag(&record, b'M', b'D').as_deref() == Some(golden_md) { - report.md_concordant_reads += 1; + // Bisulfite-aware genomic NM/MD concordance against the golden read, + // computed from the reference (not the aligner's tags) using the read's + // TRUE strand, so it is comparable across aligners regardless of their + // NM/MD convention. Requires --reference; NA without it. + if let (Some(ref_cache), Some(m_contig), Some(m_start0)) = + (reference.as_deref_mut(), mapped_contig.as_deref(), mapped_start0) + { + // Borrow each contig in turn: both calls return owned GenomicEdits, + // so the &[u8] borrows do not overlap. + let golden_edits = ref_cache.contig(&truth_aln.contig).map(|r| { + edits::genomic_edits( + &truth_aln.sequence, + &truth_aln.cigar, + truth_aln.start0, + r, + truth_aln.conv_dir, + ) + }); + let aligned_edits = ref_cache.contig(m_contig).map(|r| { + edits::genomic_edits( + record.sequence().as_ref(), + record.cigar(), + m_start0, + r, + truth_aln.conv_dir, + ) + }); + if let (Some(g), Some(a)) = (golden_edits, aligned_edits) { + report.nm_comparable_reads += 1; + report.nm_concordant_reads += u64::from(a.nm == g.nm); + report.md_comparable_reads += 1; + report.md_concordant_reads += u64::from(a.positions == g.positions); } } } @@ -458,10 +484,9 @@ mod tests { contig: "chr1".to_string(), start0, ref_len: cigar::reference_len(&cigar), - nm: None, - md: None, sequence: seq.to_vec(), cigar, + conv_dir: None, } } @@ -511,10 +536,9 @@ mod tests { contig: "chrX".to_string(), start0: 0, ref_len: 100, - nm: None, - md: None, sequence: vec![b'T'; 100], cigar: Cigar::from(vec![Op::new(Kind::Match, 100)]), + conv_dir: None, }; assert!(truth.expected_for_read(&g).is_empty()); } diff --git a/tests/test_eval.rs b/tests/test_eval.rs index 923c582..23d67f1 100644 --- a/tests/test_eval.rs +++ b/tests/test_eval.rs @@ -243,6 +243,86 @@ fn test_eval_variant_representation_perfect() { assert_eq!(parse_footer(&contents, "nm_concordant_pct"), "NA"); } +/// With `--reference`, NM/MD concordance is the bisulfite-aware genomic edit +/// distance recomputed against the reference (not the raw tags). Grading the +/// golden BAM against itself must be perfect — every read's genomic edits match +/// its own — so both report 100%, not NA. +#[test] +fn test_eval_genomic_nm_md_concordance_with_reference() { + let seq = non_repetitive_seq(2_000); + let env = TestEnv::new(&[("chr1", &seq)]); + + let positions = [400usize, 800, 1200, 1600]; + let refs: Vec = positions.iter().map(|&p| (seq[p] as char).to_string()).collect(); + let alts: Vec = + refs.iter().map(|r| if r == "A" { "C" } else { "A" }.to_string()).collect(); + let alt_arrays: Vec<[&str; 1]> = alts.iter().map(|a| [a.as_str()]).collect(); + let variants: Vec> = positions + .iter() + .enumerate() + .map(|(i, &p)| VcfVariant { + chrom: "chr1", + pos_1based: p as u32 + 1, + ref_allele: refs[i].as_str(), + alt_alleles: &alt_arrays[i], + gt: "1|1", + }) + .collect(); + let vcf = env.write_vcf("sample", &[("chr1", 2_000)], &variants); + + let sim_out = env.output_prefix(); + let (ok, _, stderr) = run_simulate(&[ + "simulate", + "-r", + env.fasta_path.to_str().unwrap(), + "-v", + vcf.to_str().unwrap(), + "-o", + sim_out.to_str().unwrap(), + "--coverage", + "30", + "--read-length", + "50", + "--fragment-mean", + "150", + "--fragment-stddev", + "20", + "--min-error-rate", + "0", + "--max-error-rate", + "0", + "--golden-bam", + "--single-end", + "--seed", + "42", + ]); + assert!(ok, "simulate failed: {stderr}"); + + let golden = PathBuf::from(format!("{}.golden.bam", sim_out.display())); + let eval_out = env.dir.path().join("eval_ref"); + let (ok, _, stderr) = run_eval(&[ + "eval", + "--mapped", + golden.to_str().unwrap(), + "--truth", + golden.to_str().unwrap(), + "--variants", + vcf.to_str().unwrap(), + "--reference", + env.fasta_path.to_str().unwrap(), + "-o", + eval_out.to_str().unwrap(), + ]); + assert!(ok, "eval failed: {stderr}"); + + let variants_tsv = PathBuf::from(format!("{}.variants.tsv", eval_out.display())); + let contents = std::fs::read_to_string(&variants_tsv).unwrap(); + // Golden vs itself: genomic edits are identical → 100% concordance, and + // crucially NOT "NA" (which is what the raw-tag path returned here). + assert_eq!(parse_footer(&contents, "md_concordant_pct"), "100.00"); + assert_eq!(parse_footer(&contents, "nm_concordant_pct"), "100.00"); +} + /// Simulate EM-seq reads with a methylation golden BAM and a cpg-truth /// bedGraph, then correlate the golden BAM's own XM calls against that truth. /// Because both derive from the same methylation draws, the correlation is