perf(methylate): methylate contigs in parallel by nh13 · Pull Request #21 · fg-labs/holodeck

nh13 · 2026-06-26T23:00:26Z

What

Parallelize holodeck methylate across contigs. Stacked on #20 (base nh/eval-variants-meth) so the accuracy-bench work can pin one SHA that carries both the truth-aware eval changes and this speedup.

Why

Generating whole-genome simulated methylation datasets calls methylate over the full hs38DH reference, and it was single-threaded — the slowest step in dataset generation. Methylation is drawn per contig by a spatial Markov chain whose correlation never crosses contig boundaries (chromosomes are physically independent), so the per-contig loop is embarrassingly parallel.

How

The per-contig loop is now a work-stealing rayon parallel map. Each contig's RNGs are already seeded purely from (seed, contig) with no state carried between contigs, so the parallel output is byte-identical to the previous single-threaded version regardless of thread count — no seeding change was needed. Each worker serializes its contig's VCF (and optional bedGraph) into its own buffer; the buffers are written in dict order so the VCF stays coordinate-sorted.
with_max_len(1) forces one contig per job. Per-contig cost is wildly uneven (chr1 vs a 50 kb alt) and the large chromosomes sort first, so rayon's default count-balanced contiguous chunking otherwise hands one worker the whole block of big contigs while the rest drain thousands of trivial alt/decoy contigs and starve. This was the difference between ~1× and ~5× on the real reference.
map_init gives each worker one reusable FASTA handle instead of reopening per contig. Fasta::from_path rebuilds the full sequence dictionary, so per-contig opens are O(contigs²) and dominate on references with thousands of contigs (hs38DH has ~3.4k). The handle is only seeked + read by load_contig, never shared across threads, so reuse is safe and deterministic.

Correctness

Verified byte-identical to the single-threaded baseline on the whole hs38DH genome — decompressed VCF content, plain VCF, and the population-fraction bedGraph all match — and stable across RAYON_NUM_THREADS=1/4/8/12. Added an integration test (methylate_output_is_identical_across_thread_counts) that asserts byte-identical output across thread counts on a multi-contig reference.

Performance

Whole human genome on a 12-core host, measured against the byte-identical sequential baseline:

output	threads=1	threads=12	speedup
compute only (plain VCF)	72 s	11 s	6.5×
real-world (BGZF VCF + bedGraph)	~110 s	~25 s	~4.6×

(For reference, methylating chr22 alone — a single contig, so no parallel benefit — is 2.4 s; the win is entirely in spreading the many large chromosomes across cores.)

Notes / tradeoffs

Peak memory now holds every contig's serialized output at once — O(whole reference output), a few GB resident for hg38 — rather than O(one contig) as the old streaming loop did. Fine on a multi-GB host; a bounded ordered-streaming writer would remove the ceiling if a far larger reference ever needs a constrained host. Documented at the call site.
On a multi-contig failure, rayon surfaces an arbitrary failing contig's error rather than deterministically the first in dict order; the job aborts either way, so this only affects which contig the message cites. Documented at the collect.

Reviewed locally with the CodeRabbit CLI and a CodeRabbit-style pass before opening.

Additional fix: allele-specific `methylate --vcf` speedup

Orthogonal follow-up found while benchmarking this PR. methylate --vcf (allele-specific methylation) was ~100× slower than the reference-only path — ≈43 min vs ≈25 s on a whole human genome with a few-million-variant VCF.

The per-haplotype CpG classifier (classify_cpgs_with_haplotypes and per_variant_per_hap_cpg_offsets_with_haplotypes) resolved each CpG's variant-vs-reference source with a linear scan over every alt span on the contig, making it O(CpGs × variants) per haplotype. The reference-only path skips this entirely via the no-variant fast path, which is why plain methylate stayed fast (and why this is a separate axis from the parallelization above, which speeds both paths equally). The Markov walk and haplotype materialization are identical between the two paths, so they were not the cause.

Because the alt spans are sorted by hap_start and disjoint, and the CpG scan walks positions in ascending order, a single monotonic cursor advances past spent ranges and resolves each lookup in O(1) amortized. Classification is now O(haplotype length + variants) per haplotype.

Output is byte-identical to before — the cursor returns exactly what the linear scan did, pinned by an equivalence test against a naive scan under the real (h, h+1) ascending access pattern. Verified on phased-variant inputs (VCF data lines + population-fraction bedGraph) and benchmarked at fixed CpG count: previously linear in variant count (6.4 s → 39.6 s for 2.4k → 20k variants), now flat at ~0.4 s.

Summary by CodeRabbit

Performance
- methylate now runs per-contig generation in parallel for faster throughput.
- methylate --vcf no longer slows down quadratically with higher variant density.
Reliability
- Output remains byte-identical to the prior single-threaded behavior, regardless of configured thread count.
Tests
- Added end-to-end determinism coverage across different Rayon thread counts, validating both VCF and bedGraph outputs.

coderabbitai · 2026-06-26T23:00:34Z

📝 Walkthrough

Walkthrough

methylate now processes contigs in parallel with deterministic per-contig seeding, uses monotonic cursors for allele-specific CpG lookup, and verifies identical output across thread counts. Changelog entries document these changes.

Changes

Methylation processing

Layer / File(s)	Summary
Cursor-based CpG resolution `src/vcf/methylation.rs`	CpG classification and variant ownership use monotonic cursors over sorted variant ranges, with equivalence tests against linear lookup.
Parallel contig methylation `Cargo.toml`, `src/commands/methylate.rs`	Rayon workers process contigs with deterministic seeds and worker-local FASTA handles, then write serialized results in dictionary order.
Determinism and behavior validation `tests/test_methylate.rs`, `CHANGELOG.md`	Tests compare single-threaded and multithreaded VCF and bedGraph output; changelog entries document parallel processing and cursor-based lookup.

Estimated code review effort: 4 (Complex) | ~45 minutes

Sequence Diagram(s)

sequenceDiagram
  participant User
  participant MethylateExecute
  participant RayonWorkers
  participant FASTA
  participant ContigOutput
  User->>MethylateExecute: run methylate
  MethylateExecute->>RayonWorkers: process contigs with par_iter
  RayonWorkers->>FASTA: load contig with worker-local handle
  RayonWorkers->>ContigOutput: serialize VCF and optional bedGraph buffers
  MethylateExecute->>ContigOutput: write buffers in dictionary order

Suggested reviewers: tfenne

Poem

I hop through contigs, parallel and bright,
Seeded paws keep each output right.
Cursors glide where variants grew,
VCFs and bedGraphs emerge byte-true.
This bunny approves with a joyful chew.

🚥 Pre-merge checks | ✅ 5

✅ Passed checks (5 passed)

Check name	Status	Explanation
Docstring Coverage	✅ Passed	No functions found in the changed files to evaluate docstring coverage. Skipping docstring coverage check.
Linked Issues check	✅ Passed	Check skipped because no linked issues were found for this pull request.
Out of Scope Changes check	✅ Passed	Check skipped because no linked issues were found for this pull request.
Description Check	✅ Passed	Check skipped - CodeRabbit’s high-level summary is enabled.
Title check	✅ Passed	The title clearly summarizes the main change: parallelizing methylate across contigs for performance.

✨ Finishing Touches

📝 Generate docstrings

Create stacked PR
Commit on current branch

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.

❤️ Share

_{Comment @coderabbitai help to get the list of available commands.}

nh13 · 2026-06-26T23:13:01Z

@coderabbitai review

coderabbitai · 2026-06-26T23:13:06Z

✅ Action performed

Review finished.

Note: CodeRabbit is an incremental review system and does not re-review already reviewed commits. This command is applicable only when automatic reviews are paused.

…r22 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).

* 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/.

nh13 · 2026-07-15T01:09:03Z

@coderabbitai resume

nh13 · 2026-07-15T01:09:04Z

@coderabbitai review

coderabbitai · 2026-07-15T01:09:10Z

✅ Action performed

Review finished.

Note: CodeRabbit is an incremental review system and does not re-review already reviewed commits. This command is applicable only when automatic reviews are paused.

coderabbitai · 2026-07-15T01:09:11Z

✅ Action performed

Reviews resumed.

nh13 · 2026-07-15T02:09:05Z

@coderabbitai resume

nh13 · 2026-07-15T02:09:06Z

@coderabbitai review

coderabbitai · 2026-07-15T02:09:12Z

✅ Action performed

Review finished.

Note: CodeRabbit is an incremental review system and does not re-review already reviewed commits. This command is applicable only when automatic reviews are paused.

coderabbitai · 2026-07-15T02:09:12Z

✅ Action performed

Reviews resumed.

coderabbitai

Actionable comments posted: 1

🤖 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 `@src/commands/methylate.rs`:
- Around line 204-273: Replace the map_init-based FASTA handle initialization in
the contig processing pipeline with a true per-worker-thread cache, so each
Rayon worker creates and reuses one Fasta handle across all assigned contigs
despite with_max_len(1). Preserve the existing methylate_contig inputs, error
propagation, deterministic output ordering, and per-thread handle isolation.

🪄 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: b7610b72-1ef4-4143-b34e-a9f01fac8ac2

📥 Commits

Reviewing files that changed from the base of the PR and between 6d2fa7f and efae5ce.

⛔ Files ignored due to path filters (1)

Cargo.lock is excluded by !**/*.lock

📒 Files selected for processing (5)

CHANGELOG.md
Cargo.toml
src/commands/methylate.rs
src/vcf/methylation.rs
tests/test_methylate.rs

Each contig is methylated independently — its RNGs are seeded purely from `(seed, contig)` with no state carried between contigs (the methylation Markov chain runs within a single contig) — so the per-contig loop now runs as a work-stealing rayon parallel map. Each worker serializes its contig's VCF (and optional bedGraph) into its own in-memory buffer; the buffers are written out in dict order so the VCF stays coordinate-sorted and the output is byte-identical to the previous single-threaded version regardless of thread count. Two details make it actually scale on a real human reference: - `with_max_len(1)` forces one contig per job. Per-contig cost is wildly uneven (chr1 vs a 50 kb alt) and the large chromosomes sort first, so rayon's default count-balanced contiguous chunking otherwise hands one worker the whole block of big contigs while the rest drain thousands of trivial alts and starve. - `map_init` gives each worker ONE reusable FASTA handle instead of reopening per contig; `Fasta::from_path` rebuilds the sequence dictionary, so per-contig opens are O(contigs^2) and dominate on references with thousands of contigs (hs38DH has ~3.4k). The handle is only seeked+read, never shared across threads, so reuse is safe and deterministic. On a whole human genome this is ~5x faster on a 12-core host (~110s -> ~25s including BGZF output). Adds an integration test asserting byte-identical output across RAYON_NUM_THREADS=1 and 8.

methylate --vcf classifies every CpG on each haplotype as variant- or reference-sourced. base_source did this with a linear scan over every alt span on the contig, so classify_cpgs_with_haplotypes and per_variant_per_hap_cpg_offsets_with_haplotypes were O(CpGs × variants) per haplotype. On a whole human genome with a few-million-variant VCF this was ~100× the work of the reference-only path (~43 min vs ~25 s); the no-variant fast path skips the scan entirely, which is why plain methylate stayed fast. The alt spans are sorted by hap_start and disjoint, and the CpG scan walks positions in ascending order, so a single monotonic cursor advances past spent ranges and resolves each lookup in O(1) amortized. Classification is now O(haplotype length + variants) per haplotype. Output is byte-identical: the cursor returns exactly what the linear scan did, pinned by an equivalence test against a naive scan under the real (h, h+1) ascending access pattern. Orthogonal to the per-contig parallelization in this PR, which speeds both paths equally.

nh13 · 2026-07-15T14:25:15Z

@coderabbitai resume

nh13 · 2026-07-15T14:25:17Z

@coderabbitai review

coderabbitai · 2026-07-15T14:25:24Z

✅ Action performed

Review finished.

Note: CodeRabbit is an incremental review system and does not re-review already reviewed commits. This command is applicable only when automatic reviews are paused.

coderabbitai · 2026-07-15T14:27:57Z

✅ Action performed

Reviews resumed.

nh13 requested a review from tfenne June 26, 2026 23:00

nh13 assigned tfenne Jun 26, 2026

nh13 force-pushed the nh/eval-variants-meth branch from 94c063e to 6d2fa7f Compare June 26, 2026 23:17

nh13 force-pushed the nh/methylate-parallel branch from 35e4b47 to 2722c36 Compare June 26, 2026 23:21

nh13 mentioned this pull request Jun 27, 2026

feat(accuracy): holodeck truth-based alignment-accuracy benchmark fg-labs/bwa-mem3-bench#20

Merged

coderabbitai Bot reviewed Jul 15, 2026

View reviewed changes

Comment thread src/commands/methylate.rs

nh13 added 2 commits July 15, 2026 09:38

nh13 force-pushed the nh/methylate-parallel branch from efae5ce to e1e2604 Compare July 15, 2026 13:39

Uh oh!

Conversation

nh13 commented Jun 26, 2026 • edited by coderabbitai Bot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

What

Why

How

Correctness

Performance

Notes / tradeoffs

Additional fix: allele-specific methylate --vcf speedup

Summary by CodeRabbit

Uh oh!

coderabbitai Bot commented Jun 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Walkthrough

Changes

Sequence Diagram(s)

Uh oh!

nh13 commented Jun 26, 2026

Uh oh!

coderabbitai Bot commented Jun 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

nh13 commented Jul 15, 2026

Uh oh!

nh13 commented Jul 15, 2026

Uh oh!

coderabbitai Bot commented Jul 15, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

coderabbitai Bot commented Jul 15, 2026

Uh oh!

nh13 commented Jul 15, 2026

Uh oh!

nh13 commented Jul 15, 2026

Uh oh!

coderabbitai Bot commented Jul 15, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

coderabbitai Bot commented Jul 15, 2026

Uh oh!

coderabbitai Bot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

nh13 commented Jul 15, 2026

Uh oh!

nh13 commented Jul 15, 2026

Uh oh!

coderabbitai Bot commented Jul 15, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

coderabbitai Bot commented Jul 15, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

nh13 commented Jun 26, 2026 •

edited by coderabbitai Bot

Loading

Additional fix: allele-specific `methylate --vcf` speedup

coderabbitai Bot commented Jun 26, 2026 •

edited

Loading

coderabbitai Bot commented Jun 26, 2026 •

edited

Loading

coderabbitai Bot commented Jul 15, 2026 •

edited

Loading

coderabbitai Bot commented Jul 15, 2026 •

edited

Loading

coderabbitai Bot commented Jul 15, 2026 •

edited

Loading