feat(eval): truth-based variant-representation and methylation-correlation metrics

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 `<prefix>.variants.tsv` and `<prefix>.meth.tsv` alongside the
+  existing `<prefix>.eval.txt`.
+
 ### Fixed
 
 - `simulate` and `methylate` now tolerate VCFs that redeclare a header ID

README.md

@@ -439,21 +439,39 @@ 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 `<prefix>.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:**
 
 | 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 |
+|------|---------------|---------|
+| `<prefix>.eval.txt` | always | placement accuracy per MAPQ bin |
+| `<prefix>.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 |
+| `<prefix>.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)