feat(eval): truth-based variant-representation and methylation-correlation metrics - #20
Open
nh13 wants to merge 9 commits into
Open
feat(eval): truth-based variant-representation and methylation-correlation metrics#20nh13 wants to merge 9 commits into
nh13 wants to merge 9 commits into
Conversation
|
Note Reviews pausedUse the following commands to manage reviews:
Use the checkboxes below for quick actions:
No actionable comments were generated in the recent review. 🎉 ℹ️ Recent review info⚙️ Run configurationConfiguration used: Organization UI Review profile: CHILL Plan: Pro Run ID: 📒 Files selected for processing (10)
✅ Files skipped from review due to trivial changes (1)
🚧 Files skipped from review as they are similar to previous changes (8)
📝 WalkthroughWalkthrough
ChangesEval scoring stack
Sequence Diagram(s)sequenceDiagram
participant Eval
participant GoldenLoad as "golden::load"
participant PlacementRun as "placement::run"
participant VariantsRun as "variants::run"
participant RefCache
participant MethRun as "meth::run"
opt `--truth`
Eval->>GoldenLoad: load golden BAM
end
Eval->>PlacementRun: run placement scoring
opt `--variants`
Eval->>VariantsRun: run variant scoring
opt `--reference`
VariantsRun->>RefCache: contig lookups
end
end
opt `--cpg-truth`
Eval->>MethRun: run methylation correlation
end
Estimated code review effort🎯 4 (Complex) | ⏱️ ~60 minutes Possibly related PRs
Suggested reviewers
Poem
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✨ Finishing Touches📝 Generate docstrings
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
nh13
force-pushed
the
nh/eval-variants-meth
branch
2 times, most recently
from
June 26, 2026 21:59
e632077 to
94c063e
Compare
nh13
marked this pull request as ready for review
June 26, 2026 22:01
Contributor
Author
|
@coderabbitai review |
✅ Action performedReview finished.
|
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 6
🧹 Nitpick comments (1)
src/commands/eval/meth.rs (1)
218-218: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick winAdd
//!docs to the test module.
mod testsis a module, so it should carry the same inner-module documentation as the rest of the Rust modules in this PR. As per coding guidelines,**/*.rs: Add module-level//!documentation on all modules.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/commands/eval/meth.rs` at line 218, The test module declaration is missing inner-module documentation. Add a `//!` doc comment to the `mod tests` module so it follows the same module-level documentation pattern used elsewhere in this Rust file; place the documentation directly above the `mod tests` block and keep it focused on describing the purpose of the tests module.Source: Coding guidelines
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@README.md`:
- Around line 454-466: Update the `--truth` row in the options table so it
matches what `eval` actually consumes from the golden BAM via `golden::load`:
describe it as providing per-read true span, sequence, and bisulfite strand, and
remove the claims about haplotype and MD/NM tags. Keep the wording consistent
with the `eval`/`mod.rs` behavior and mention that MD/NM concordance is derived
from the reference when `--reference` is supplied.
In `@src/commands/eval/edits.rs`:
- Around line 107-130: RefCache currently reuses a single SmallRng seeded with
0, so ambiguity resolution in load_contig can vary based on contig access order.
Update RefCache::contig to create a fresh deterministic SmallRng for each first
load, seeding it from a stable FNV-1a hash of the contig name before calling
Fasta::load_contig. Keep the existing caching behavior, but remove the
order-dependent shared RNG state in RefCache::new and RefCache::contig.
In `@src/commands/eval/meth.rs`:
- Around line 145-159: The bedGraph parsing in meth.rs is ignoring the end
column in the row tuple and treating every record as a single CpG at start,
which allows merged or malformed intervals to slip through. Update the parsing
logic in the bedGraph reader to validate the end value alongside start, and
reject rows unless the interval length is exactly one CpG as expected before
inserting into truth. Use the existing parsing flow around the start0/rate
conversion and the truth.entry insertion to keep the check localized and surface
bad input early.
In `@src/commands/eval/mod.rs`:
- Around line 66-69: Update the doc comment for the meth flag in the eval
command so it no longer says it interprets reads as bisulfite/EM-seq; instead,
describe that --meth only breaks --variants results down by bisulfite
substitution class (conversion, mirror, transversion, other), while the actual
bisulfite/EM-seq context comes from the aligner output or golden BAM tags. Keep
the clarification attached to the meth field in src/commands/eval/mod.rs so the
CLI help matches behavior.
In `@src/commands/eval/placement.rs`:
- Around line 176-187: Update the mapq_bin function so MAPQ 255 is treated as
unknown rather than falling into the 60+ bucket. Add an explicit 255 case in the
match, and either return a distinct NA/unknown bin key or handle it separately
where the MAPQ buckets are aggregated in the placement evaluation flow so it
does not get counted with high-confidence alignments.
In `@src/commands/eval/variants.rs`:
- Line 32: The variant classification logic in eval variants is pulling class
meth labels from mapped BAM tags instead of the truth alignment, which can
misclassify variants as unclassified. Update the classification path in the eval
variants code to use truth_aln.conv_dir, matching the reference-based NM/MD
handling, and remove the dependence on XG/XR tags read from the mapped record in
the affected classification block.
---
Nitpick comments:
In `@src/commands/eval/meth.rs`:
- Line 218: The test module declaration is missing inner-module documentation.
Add a `//!` doc comment to the `mod tests` module so it follows the same
module-level documentation pattern used elsewhere in this Rust file; place the
documentation directly above the `mod tests` block and keep it focused on
describing the purpose of the tests module.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Organization UI
Review profile: CHILL
Plan: Pro
Run ID: 4526ef3f-cecb-4df5-8689-8cf56affb342
📒 Files selected for processing (11)
CHANGELOG.mdREADME.mdsrc/commands/eval.rssrc/commands/eval/cigar.rssrc/commands/eval/edits.rssrc/commands/eval/golden.rssrc/commands/eval/meth.rssrc/commands/eval/mod.rssrc/commands/eval/placement.rssrc/commands/eval/variants.rstests/test_eval.rs
💤 Files with no reviewable changes (1)
- src/commands/eval.rs
The eval command is about to grow optional truth-aware metrics (variant representation, methylation-level correlation, golden-BAM placement), so the existing MAPQ-binned placement scoring moves out of the command struct into a focused `placement` submodule. This leaves `eval/mod.rs` as a thin orchestrator over the shared `Eval` options and gives each future metric its own file. Pure code movement: the placement loop, MAPQ binning, and TSV row formatting are unchanged, and `placement::run` reproduces the previous `execute` behavior. The only edit is renaming the final log line to "Placement results written to" so it reads correctly once sibling metrics write their own tables.
`holodeck eval --variants truth.vcf --truth golden.bam` adds an accuracy axis beyond placement: for every simulated substitution a read should carry, does the aligned read actually represent the alternate base, and with what MAPQ and alignment score? Truth is taken entirely from holodeck's own outputs. The golden BAM supplies each read's true span and source haplotype (hp:i); the truth VCF's phased genotypes say which single-base substitutions that haplotype carries within the span. For each expected substitution the pass walks the *mapped* read's CIGAR to the variant's reference position and compares the observed base to the alternate allele. Reads mapped to the wrong locus simply fail to represent their variants, so mismapping is captured without special-casing. Per-read MD:Z / NM:i tags are compared against the golden tags to report tag-level concordance over variant-bearing reads. Results land in `<prefix>.variants.tsv`. With `--meth`, results break down by substitution class relative to the read's bisulfite conversion direction (XG, falling back to XR). The C->T cell on a CT-strand read (G->A on GA) is intrinsically confounded with conversion and is labelled as such rather than scored as a real signal; the discriminating classes are the mirror (T->C / A->G) and the transversions. This is the axis on which a methylation-aware scoring mode should match the genomic truth without over- or under-penalizing. Adds eval submodules cigar (CIGAR geometry, unit-tested), golden (truth index), and variants (classifier, per-haplotype expected SNVs, the scoring pass). The classifier and expected-SNV queries are covered by unit tests with programmatically built records; the end-to-end pass is exercised by a later integration test.
`holodeck eval --cpg-truth truth.bedGraph` reports how well an aligner's methylation calls reproduce the simulated truth. For every mapped read the Bismark `XM:Z` string is walked alongside the CIGAR, tallying each `Z`/`z` CpG call at its reference position; the aligner methylation level at a site (`n_methylated / coverage`) is then correlated against the truth level (`rate / 100`) from the bedGraph that `simulate --cpg-truth-bedgraph` writes. Pearson r and RMSE over the shared, covered CpG sites land in `<prefix>.meth.tsv`, with `NA` reported where a statistic is undefined. The XM walk reuses the CIGAR geometry helpers; a new `for_each_aligned` visitor yields (read offset, reference position) for each aligned base so insertions and soft-clips never misplace a call onto the reference. Tallying, bedGraph parsing, the Pearson/RMSE math, and the join are unit-tested with programmatically built inputs.
`--truth` previously only warned that it was unimplemented and fell back to read names. It now supplies the placement truth too: when a golden BAM is given, each read's true contig and start come from its golden record rather than the encoded name. This is exact and indel-aware (the golden alignment's true start already reflects haplotype indels), and it lets placement be scored for any aligner output whose reads carry holodeck names — even after a tool rewrites or trims them — as long as the golden BAM is present. The golden BAM is loaded once in `execute` and shared between the placement and variant passes, so `--truth` and `--variants` together read it a single time. With no `--truth`, placement is unchanged and still parses encoded read names.
Two integration tests drive the new metrics through real holodeck output rather than hand-built BAMs. The first simulates single-end reads carrying homozygous-alt SNVs with a golden BAM, then evaluates the golden BAM as the mapped BAM: every variant-bearing read is perfectly placed, so all expected substitutions must be represented, and — since a non-methylation golden BAM carries no MD/NM tags — MD/NM concordance must report NA rather than 0%. The second simulates EM-seq reads with a methylation golden BAM and a cpg-truth bedGraph and confirms the golden XM calls correlate strongly with the truth. Documents the expanded eval surface (the --truth / --variants / --cpg-truth / --meth options and the .variants.tsv / .meth.tsv outputs) in the README, and records the feature under CHANGELOG [Unreleased].
…VCF phasing Variant representation inferred which haplotype carried the alt from the truth VCF genotype (`alt_by_hap[read's haplotype]`), assigning an unphased het alt to haplotype 1 deterministically. But `simulate` assigns unphased genotypes to haplotypes by a random permutation, so eval's guess disagreed with the simulator on roughly half of heterozygous sites — it checked the reference-copy reads for an alt that was actually on the other copy, scoring them as misses. The effect was large: a byte-identical golden-as-mapped run scored ~74% representation instead of 100%, and every aligner inherited the same artifact. Resolve the allele a read truly carries from the golden read's own sequence instead. The golden BAM is the per-read oracle: at each truth substitution site within a read's true span, the golden read's base is exactly what the simulator placed on the copy that read was sequenced from. A read showing the reference base is (correctly) not expected to carry the alt. This makes representation correct whether or not the truth VCF is phased — golden-as-mapped now scores 100%, and a real aligner's shortfall reflects only its actual mismapping. `GoldenInfo` carries the read sequence and CIGAR (new `base_at` resolver); the now-unused `hp:i` haplotype field is dropped from the eval-side record.
`--meth` methylation correlation aborted the whole eval when a mapped BAM had no XM tags. But a plain bisulfite aligner (e.g. bwameth) emits no per-base methylation calls at all — calling is a separate extractor (MethylDackel) step — so its BAM legitimately has no XM, and bailing made it impossible to score such an aligner's placement and variant representation alongside callers that do emit XM. Warn and report NA methylation correlation (n_cpg 0) instead of erroring, so the placement and variant-representation axes still complete for these BAMs.
NM/MD concordance compared the aligner's raw NM:i/MD:Z tags against the golden tags. Those tags are convention-dependent for bisulfite data: an aligner may score edits against the original 4-letter reference (every unmethylated C->T is a "mismatch"), against a C->T-converted reference (conversions match), or in a bisulfite-aware convention. Two correctly-placed reads therefore disagree on their tags purely by convention — D3 (scoring vs the original reference) matched the golden ~1% of the time while bwameth (vs the converted reference) matched ~75%, measuring which convention each picked rather than alignment quality. Add `--reference <fasta>` and recompute concordance as a convention-independent genomic edit distance: walk each read against the reference, exclude bisulfite conversions using the read's TRUE strand (from the golden truth, so it works for aligners like bwameth that emit no XG), and compare the resulting non-conversion mismatch + indel profiles of the aligned and golden reads. D3 and bwameth now both land at ~96-100% on the same scale. The new `edits` module holds the genomic-edit walk and an on-demand reference-contig cache; the golden record now carries its true conversion strand instead of the raw NM/MD tags. Without `--reference`, NM/MD concordance is reported as NA.
nh13
force-pushed
the
nh/eval-variants-meth
branch
2 times, most recently
from
June 26, 2026 23:17
94c063e to
6d2fa7f
Compare
Contributor
Author
|
@coderabbitai review |
✅ Action performedReview finished.
|
nh13
added a commit
to fg-labs/bwa-mem3-bench
that referenced
this pull request
Jun 27, 2026
The accuracy benchmark needs holodeck to generate simulated truth datasets and to score aligner output via `holodeck eval`. holodeck is a public fg-labs Rust crate, so the builder stage cargo-installs it from a pinned git ref rather than vendoring it (the minibwa submodule pattern exists only because lh3/minibwa is private). It is built with the +stable toolchain already installed for tricord — the +toolchain override takes precedence over holodeck's rust-toolchain.toml pin, so the image isn't coupled to holodeck's exact channel. holodeck's `built` build-dependency links libgit2, so libgit2-dev joins the builder apt set (build-only; the resulting binary does not link it). HOLODECK_REPO / HOLODECK_REF join docker/build-arg-defaults.env, exposed via bwa_mem3_bench.holodeck_ref() (factored alongside minibwa_sha() onto a shared _build_arg_default helper) and passed through `cli build` with an optional --holodeck-ref override. The ref is pinned to the fg-labs/holodeck#20 draft-PR SHA while that PR is validated end-to-end here; it moves to a tagged release once holodeck merges.
nh13
added a commit
to fg-labs/bwa-mem3-bench
that referenced
this pull request
Jun 27, 2026
The accuracy benchmark needs holodeck to generate simulated truth datasets and to score aligner output via `holodeck eval`. holodeck is a public fg-labs Rust crate, so the builder stage cargo-installs it from a pinned git ref rather than vendoring it (the minibwa submodule pattern exists only because lh3/minibwa is private). It is built with the +stable toolchain already installed for tricord — the +toolchain override takes precedence over holodeck's rust-toolchain.toml pin, so the image isn't coupled to holodeck's exact channel. holodeck's `built` build-dependency links libgit2, so libgit2-dev joins the builder apt set (build-only; the resulting binary does not link it). HOLODECK_REPO / HOLODECK_REF join docker/build-arg-defaults.env, exposed via bwa_mem3_bench.holodeck_ref() (factored alongside minibwa_sha() onto a shared _build_arg_default helper) and passed through `cli build` with an optional --holodeck-ref override. The ref is pinned to the fg-labs/holodeck#20 draft-PR SHA while that PR is validated end-to-end here; it moves to a tagged release once holodeck merges.
nh13
added a commit
to fg-labs/bwa-mem3-bench
that referenced
this pull request
Jun 27, 2026
The accuracy benchmark needs holodeck to generate simulated truth datasets and to score aligner output via `holodeck eval`. holodeck is a public fg-labs Rust crate, so the builder stage cargo-installs it from a pinned git ref rather than vendoring it (the minibwa submodule pattern exists only because lh3/minibwa is private). It is built with the +stable toolchain already installed for tricord — the +toolchain override takes precedence over holodeck's rust-toolchain.toml pin, so the image isn't coupled to holodeck's exact channel. holodeck's `built` build-dependency links libgit2, so libgit2-dev joins the builder apt set (build-only; the resulting binary does not link it). HOLODECK_REPO / HOLODECK_REF join docker/build-arg-defaults.env, exposed via bwa_mem3_bench.holodeck_ref() (factored alongside minibwa_sha() onto a shared _build_arg_default helper) and passed through `cli build` with an optional --holodeck-ref override. The ref is pinned to the fg-labs/holodeck#20 draft-PR SHA while that PR is validated end-to-end here; it moves to a tagged release once holodeck merges.
nh13
added a commit
to fg-labs/bwa-mem3-bench
that referenced
this pull request
Jun 28, 2026
* feat(docker): install holodeck for the truth-based accuracy benchmark The accuracy benchmark needs holodeck to generate simulated truth datasets and to score aligner output via `holodeck eval`. holodeck is a public fg-labs Rust crate, so the builder stage cargo-installs it from a pinned git ref rather than vendoring it (the minibwa submodule pattern exists only because lh3/minibwa is private). It is built with the +stable toolchain already installed for tricord — the +toolchain override takes precedence over holodeck's rust-toolchain.toml pin, so the image isn't coupled to holodeck's exact channel. holodeck's `built` build-dependency links libgit2, so libgit2-dev joins the builder apt set (build-only; the resulting binary does not link it). HOLODECK_REPO / HOLODECK_REF join docker/build-arg-defaults.env, exposed via bwa_mem3_bench.holodeck_ref() (factored alongside minibwa_sha() onto a shared _build_arg_default helper) and passed through `cli build` with an optional --holodeck-ref override. The ref is pinned to the fg-labs/holodeck#20 draft-PR SHA while that PR is validated end-to-end here; it moves to a tagged release once holodeck merges. * feat(sim): deterministic holodeck truth-dataset generator scripts/gen_holodeck_dataset.sh produces one accuracy-benchmark truth dataset from a reference FASTA: mutate -> [methylate] -> simulate, fully seeded. It emits the truth VCF (the SNVs eval scores), the simulated FASTQs, the golden BAM (eval's --truth), and — for the meth kinds — the coverage-weighted cpg-truth bedGraph (eval's --cpg-truth). Four kinds cover the matrix: wgs-place / meth-place sweep the genome at low coverage for placement and MAPQ calibration; wgs-vars / meth-vars put depth over a target BED to exercise variant representation. The reference is always full hg38 — depth is bounded by the BED, never by a reduced reference, so off-target mismapping stays observable. Verified with shellcheck. The generated datasets are large and live in S3 (staged under each sim sample's source prefix), never committed; this is the reproducer that builds them. * feat(workflow): sim sample class + eval.smk accuracy rule + targets Wire the truth-based accuracy axis into the workflow. A new `truth: true` flag on a sample marks its S3 source prefix as also holding the holodeck truth artifacts (golden.bam, truth.vcf, and cpg-truth.bedGraph for meth); the new eval_accuracy rule grades an aligner's BAM directly against that truth via `holodeck eval`, emitting placement (.eval.txt), per-read variant representation (.variants.tsv), and methylation-level correlation (.meth.tsv) in one invocation — no variant caller or methylation extractor in the path. The eval rule's `tool` wildcard routes to each arm's BAM cache (fg-labs under runs/, the bwa-mem2/bwameth baseline under baseline/, minibwa under minibwa/). The `--meth-scoring genomic` D3 arm is expressed as a separate sample sharing its FASTQs and truth with its collapsed sibling, since the align rules key flags off the sample name; those `*-genomic` samples are fg-labs-arm-only. Accuracy is a property of the aligner build and is ~arch-invariant, so the accuracy targets pin each chemistry to a single arch (non-meth -> c6a, meth -> m7i) instead of sweeping SIMD tiers like the timing matrix. The new `accuracy` (full hg38) and `accuracy_smoke` (chr22-slice) targets enumerate the arm matrix per dataset. Truth samples are excluded from `rule all` and `baseline_all` — they have no upstream-concordance question. * feat(storage): accuracy table + ingest_accuracy Add a v4 schema migration introducing the `accuracy` table — truth-based holodeck eval results, one row per (run, sim-sample, arch, rep, aligner arm). Distinct from `comparisons` (tool-vs-tool agreement), this is graded against simulation truth: placement + MAPQ calibration (the per-bin table as JSON plus the ALL-row rates as headline columns), per-read variant representation (the per-class accumulators as JSON plus the MD/NM concordance footers), and methylation-level correlation (NULL for non-meth). All arms of a run share the run's fg_labs_sha — the eval outputs live under runs/<sha>/ — so the `tool` column disambiguates fg-labs / baseline / minibwa. Being a new table, the migration needs no ALTER: connect()'s unconditional executescript creates it. ingest_accuracy walks runs/<sha>/<sample>/<arch>/rep-N/eval/<tool>.{eval.txt, variants.tsv,meth.tsv}, parsing holodeck's exact TSV formats (the literal "NA" maps to NULL; the empty .meth.tsv placeholder the eval rule writes for non-meth samples parses to no correlation). cli collect calls it after ingest_run, since accuracy outputs share the runs/ tree. * feat(report): bench accuracy (placement + variant-class honesty + meth-r) Add the `bench accuracy` report over the `accuracy` table. It renders three markdown sections per sim dataset, one row per aligner arm: placement + MAPQ calibration (the cross-tool bwa-mem3 / minibwa / baseline comparison), variant representation + methylation correlation (variant-bearing reads, MD/NM concordance, per-CpG Pearson r / RMSE), and per-class AS/MAPQ honesty — the genomic-vs-collapsed headline, with the conversion-direction class flagged `confounded` since no mode resolves it from a single read. Headline metrics are averaged across reps; the per-class table takes the lowest rep per cell (the JSON blobs aren't meaningfully averaged). NULL cells (non-meth MD/NM, non-meth methylation, NA mean_as) render as an em-dash. Wired as `cli bench accuracy --fg-labs-sha <sha> [--out file.md]`, mirroring `bench speedup`. * chore(accuracy): pin holodeck #21, align full -place to ~0.5x, add chr22 vars target - HOLODECK_REF -> 35e4b47 (fg-labs/holodeck#21, per-contig parallel methylate, stacked on #20 so it carries the truth-aware eval work too). - gen_holodeck_dataset.sh: -place coverage 1 -> 0.5 to hit the spec's ~5M read-pair genome-wide placement scale (no -place smoke uses this script). - scripts/sim-targets/chr22.bed: whole-chr22 target for the full -vars datasets (single region, ~34Mb non-N => ~3.9M pairs at 30x; larger than the smoke's 3x100kb). * docs(accuracy): committed wrapper + runbook for simulated truth datasets scripts/gen_all_sim_datasets.sh is the single source of truth for the dataset matrix the accuracy benchmark consumes (name -> kind -> coverage -> target -> chemistry), driving the per-dataset gen_holodeck_dataset.sh and optionally staging the canonical files to S3. docs/data-setup.md gains a 'Simulated truth datasets (holodeck)' section documenting the matrix, canonical filenames, the S3 layout, and seed/holodeck-SHA provenance for reproducibility. * chore(accuracy): bump holodeck pin to efae5ce (PR #21 meth perf fix) PR #21 (per-contig parallel methylate) also carries efae5ce, an O(1) per-CpG variant-lookup fix that makes methylate --vcf and em-seq simulate fast on a whole genome (methylate --vcf 43 min -> 31 s; em-seq simulate chr1-stall -> ~9 s/contig; output byte-identical). Pin HOLODECK_REF to that tip so the bench image and eval use it. * feat(accuracy): add minibwa as the 4th meth arm (bisulfite-aware, --meth) The meth accuracy matrix had three arms (bwameth, D3-collapsed, D3-genomic); add minibwa so all four aligners are graded on the same meth truth. minibwa has a directional BS-seq mode: `map --meth` maps the EM-seq reads bisulfite-aware (read1 C->T, read2 G->A) against the `.meth.mbw` index built by `index --meth` (which lives alongside the plain DNA index, not under hg38-meth, so we strip the -meth suffix to reach it and stage .l2b + .meth.mbw). It is scored on placement + variant representation; methylation-level correlation is NA (minibwa emits no XM tags even in --meth mode — same as the bwameth arm; only bwa-mem3 calls methylation natively). _accuracy_targets now requests the minibwa arm for every primary (non-genomic) sim sample, meth included; the genomic D3 sibling stays fg-labs-only (its baseline/minibwa arms would be byte-identical to the collapsed sample's). minibwa runs on the meth arch (m7i) so the 4-way wall-time comparison is arch-matched. The .meth.mbw sidecar is built once (index --meth, in-memory path) and staged to references/hg38/.
Contributor
Author
|
@coderabbitai pause |
✅ Action performedReviews paused. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Extends
holodeck evalto score an aligner's BAM against holodeck's own truth along three axes, all reading only holodeck-produced outputs (golden BAM, truth VCF, cpg-truth bedGraph) — no external variant or methylation callers.<prefix>.eval.txt): unchanged accuracy-by-MAPQ-bin, but--truthnow makes the golden BAM the truth source (exact and indel-aware) instead of only the encoded read name. This implements the previously-stubbed--truthflag.--variants truth.vcf+--truth golden.bam,<prefix>.variants.tsv): for every single-base substitution a read should carry — derived from the truth VCF's phased per-haplotype genotypes within the read's golden span — walk the mapped read's CIGAR to the variant position and check whether the observed base matches the alternate allele. Reports the represented fraction with the read's MAPQ and alignment score per substitution class, plus per-read MD/NM concordance against the golden tags. Reads mapped to the wrong locus simply fail to represent their variants, so mismapping is captured without special-casing.--cpg-truth truth.bedGraph,<prefix>.meth.tsv): walk each mapped read's BismarkXM:Zstring alongside the CIGAR, tally per-CpG methylated/unmethylated calls, and correlate the aligner methylation level against the truth bedGraph (Pearson r + RMSE).Methylation framing
With
--meth, variant results break down by substitution class relative to the read's bisulfite conversion direction (XG, falling back toXR). TheC->Tcell on aCT-strand read (G->AonGA) is intrinsically confounded with bisulfite conversion and is labelled as such rather than scored as a real signal; the discriminating classes are the mirror (T->C/A->G) and the transversions. A read lacking a conversion-direction tag is counted asunclassifiedinstead of being silently mislabelled.Structure
The
evalcommand moves into aneval/module:placement(existing, extracted),cigar(CIGAR geometry helpers),golden(golden-BAM truth index),variants(substitution classifier + per-haplotype expected SNVs + scoring pass), andmeth(XM tally + correlation). The two refactor commits (move, extract) are behavior-preserving and worth reading first.Testing
CIGAR geometry, the substitution classifier, per-haplotype expected-SNV queries, the bedGraph parser, and the Pearson/RMSE math are unit-tested with programmatically built inputs. Two end-to-end integration tests drive the new metrics through real
simulateoutput: one confirms a golden BAM scored as its own mapped BAM represents every variant (and reports MD/NM asNA, since non-methylation golden BAMs carry no such tags), the other confirms goldenXMcalls correlate strongly with the cpg-truth bedGraph.cargo ci-test(456 tests),cargo ci-lint(clippy pedantic), andcargo fmt --checkall pass.Status
Draft. The metrics are being exercised end-to-end in the bwa-mem3-bench accuracy benchmark before this is marked ready; expect possible small follow-ups from that integration.
Summary by CodeRabbit
holodeck evalto provide placement accuracy plus opt-in variant representation, methylation correlation, and reference-based NM/MD concordance.*.variants.tsv(per substitution class) and*.meth.tsv(per-CpG correlation), alongside the always-produced*.eval.txt.--referenceis set.