diff --git a/CHANGELOG.md b/CHANGELOG.md index 3934521..e4ae6b3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,34 @@ 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 + sequences from the golden alignment rather than only the encoded read name, + and `--variants` reports how faithfully aligned reads represent the simulated + 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. 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 - `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/src/commands/eval.rs b/src/commands/eval.rs deleted file mode 100644 index 97dd742..0000000 --- a/src/commands/eval.rs +++ /dev/null @@ -1,255 +0,0 @@ -//! Alignment accuracy evaluation command. - -use std::collections::BTreeMap; -use std::io::Write; -use std::path::PathBuf; - -use anyhow::{Context, Result}; -use bstr::ByteSlice; -use clap::Parser; -use noodles::bam; - -use super::command::{Command, output_path}; -use super::common::OutputPrefixOptions; -use crate::read_naming::{parse_encoded_pe_name, parse_encoded_se_name}; - -/// 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). -#[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")] -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. - #[arg(long, value_name = "BAM")] - pub truth: Option, - - #[command(flatten)] - pub output: OutputPrefixOptions, - - /// Maximum distance (in bases) between the true and mapped start positions - /// of a read for it to be considered correctly mapped. Uses - /// `|mapped_start - true_start| <= wiggle` on the same contig. - #[arg(long, default_value_t = 5, value_name = "INT")] - 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")); - } -} diff --git a/src/commands/eval/cigar.rs b/src/commands/eval/cigar.rs new file mode 100644 index 0000000..97d8cc5 --- /dev/null +++ b/src/commands/eval/cigar.rs @@ -0,0 +1,174 @@ +//! 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 +} + +/// 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::*; + 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/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 new file mode 100644 index 0000000..48b309b --- /dev/null +++ b/src/commands/eval/golden.rs @@ -0,0 +1,174 @@ +//! 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 sequence the read +//! 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`]) and the +//! read's true bisulfite strand (see [`GoldenInfo::conv_dir`]). + +use std::collections::HashMap; +use std::path::Path; + +use anyhow::{Context, Result, bail}; +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; + +/// 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); + +/// 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 { + /// 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, + /// 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, + /// 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 { + /// 0-based exclusive end of the true alignment. + #[must_use] + 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)`. +/// +/// # 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; + } + + // 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 { + 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 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(&cigar), + // 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), + }; + 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) +} + +/// 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) +} + +/// 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/meth.rs b/src/commands/eval/meth.rs new file mode 100644 index 0000000..6d28dc6 --- /dev/null +++ b/src/commands/eval/meth.rs @@ -0,0 +1,346 @@ +//! 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 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 { + log::warn!( + "no XM methylation tags in mapped primaries of {}; \ + reporting NA methylation correlation", + 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 new file mode 100644 index 0000000..7deb16b --- /dev/null +++ b/src/commands/eval/mod.rs @@ -0,0 +1,144 @@ +//! 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; 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`). +//! - [`meth`] — methylation-level correlation (`--cpg-truth`; +//! `.meth.tsv`). + +mod cigar; +mod edits; +mod golden; +mod meth; +mod placement; +mod variants; + +use std::path::PathBuf; + +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. +/// +/// 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. `--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 \ + 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, + + /// Golden BAM (`simulate --golden-bam`) supplying each read's true span, + /// sequence, and (for meth) bisulfite strand. 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). Only meaningful with `--variants`. + #[arg(long, value_name = "NAME", requires = "variants")] + 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, + + /// 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, + + /// 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, + + /// Maximum distance (in bases) between the true and mapped start positions + /// of a read for it to be considered correctly mapped. Uses + /// `|mapped_start - true_start| <= wiggle` on the same contig. + #[arg(long, default_value_t = 5, value_name = "INT")] + pub wiggle: u32, +} + +impl Command for Eval { + fn execute(&self) -> Result<()> { + 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"); + } + + // 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 = golden.as_ref().expect("--variants requires --truth"); + let truth = variants::VariantTruth::from_vcf(vcf, self.sample.as_deref())?; + 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 { + meth::run(&self.mapped, cpg_truth, &self.output.output)?; + } + + Ok(()) + } +} diff --git a/src/commands/eval/placement.rs b/src/commands/eval/placement.rs new file mode 100644 index 0000000..1048ddf --- /dev/null +++ b/src/commands/eval/placement.rs @@ -0,0 +1,259 @@ +//! Placement accuracy: true vs mapped position, stratified by 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, HashMap}; +use std::io::Write; +use std::path::Path; + +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 { + /// 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`. +/// +/// 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, + golden: Option<&HashMap>, +) -> 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; + + let name_bytes = record.name().map_or(&b""[..], |n| n.as_bytes()); + + // 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 { + 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 { + 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 { + let reason = if golden.is_some() { "no golden truth record" } else { "unparseable names" }; + log::warn!("{parse_failures} reads had {reason}; 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")); + } +} diff --git a/src/commands/eval/variants.rs b/src/commands/eval/variants.rs new file mode 100644 index 0000000..3326c93 --- /dev/null +++ b/src/commands/eval/variants.rs @@ -0,0 +1,582 @@ +//! Variant-representation accuracy: do aligned reads carry the simulated +//! variants they should, and how confidently? +//! +//! 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 +//! +//! 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 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; +use std::path::Path; + +use anyhow::{Context, Result}; +use noodles::bam; + +use super::cigar; +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}; + +/// 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 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 the golden read carries at this site. + pub alt_base: u8, +} + +/// 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. +#[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 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 { + let mut sites = Vec::new(); + for record in records { + if record.ref_allele.len() != 1 { + continue; // SNV requires a single reference base. + } + // 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, + }); + } + } + sites.sort_by_key(|s| s.pos0); + if !sites.is_empty() { + by_contig.insert(contig.clone(), sites); + } + } + Self { by_contig } + } + + /// 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 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 { + return Vec::new(); + }; + 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(base) = golden.base_at(site.pos0) + && site.alts.contains(&base) + { + out.push(ExpectedSnv { pos0: site.pos0, ref_base: site.ref_base, alt_base: base }); + } + } + 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) + } +} + +/// 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, + mut reference: Option<&mut RefCache>, + 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_for_read(truth_aln); + 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 { 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. + 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)); + + 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 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); + } + + // 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); + } + } + } + + report.write_tsv(output_prefix) +} + +#[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}; + + #[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) + } + + 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), + sequence: seq.to_vec(), + cigar, + conv_dir: None, + } + } + + #[test] + 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_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 site + snv_record(150, "C", "A", "0/1"), + ]); + // 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_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, + 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()); + } + + #[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); + } +} diff --git a/tests/test_eval.rs b/tests/test_eval.rs index a2eb923..23d67f1 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,266 @@ 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"); +} + +/// 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 +/// 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]