diff --git a/CHANGELOG.md b/CHANGELOG.md index e4ae6b3..6c7aa74 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -35,8 +35,31 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 are written to `.variants.tsv` and `.meth.tsv` alongside the existing `.eval.txt`. +### Changed + +- `methylate` now methylates contigs in parallel. Each contig is independent — + 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 runs as a work-stealing parallel map and the output is + byte-identical to the previous single-threaded version regardless of thread + count. Single-item jobs spread the wildly-uneven per-contig cost (chr1 ≫ a + 50 kb alt) evenly across the pool, and one reused FASTA handle per worker + avoids re-parsing the sequence dictionary per contig. On a whole human genome + this is roughly 5× faster on a 12-core host (~110 s → ~25 s including BGZF + output); the thread count honors `RAYON_NUM_THREADS`. + ### Fixed +- `methylate --vcf` (allele-specific methylation) no longer slows down + quadratically with variant density. The per-haplotype CpG classifier looked + up each CpG's variant/reference source with a linear scan over every variant + on the contig, making it 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). Because the alt spans are sorted and + disjoint and the CpG scan is ascending, a single monotonic cursor now resolves + each lookup in O(1) amortized, making classification O(haplotype length + + variants); output is byte-identical. The reference-only path was already fast + and is unchanged. - `simulate` and `methylate` now tolerate VCFs that redeclare a header ID (e.g. `duplicate INFO ID: BREAKSIMLENGTH`). Such duplicates are common in files from upstream tools and are accepted by bcftools; holodeck drops the diff --git a/Cargo.lock b/Cargo.lock index bc1d95d..76d8cca 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -293,6 +293,25 @@ dependencies = [ "crossbeam-utils", ] +[[package]] +name = "crossbeam-deque" +version = "0.8.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51" +dependencies = [ + "crossbeam-epoch", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.9.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" +dependencies = [ + "crossbeam-utils", +] + [[package]] name = "crossbeam-utils" version = "0.8.21" @@ -310,6 +329,12 @@ dependencies = [ "syn", ] +[[package]] +name = "either" +version = "1.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91622ff5e7162018101f2fea40d6ebf4a78bbe5a49736a2020649edf9693679e" + [[package]] name = "enum_dispatch" version = "0.3.13" @@ -543,6 +568,7 @@ dependencies = [ "pooled-writer", "rand", "rand_distr", + "rayon", "tempfile", "thiserror 2.0.18", ] @@ -1287,6 +1313,26 @@ dependencies = [ "rand", ] +[[package]] +name = "rayon" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fb39b166781f92d482534ef4b4b1b2568f42613b53e5b6c160e24cfbfa30926d" +dependencies = [ + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91" +dependencies = [ + "crossbeam-deque", + "crossbeam-utils", +] + [[package]] name = "redox_syscall" version = "0.5.18" diff --git a/Cargo.toml b/Cargo.toml index b128e18..883cde0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -38,6 +38,11 @@ chrono = "0.4" # Vec would cost 8x the memory at whole-chromosome scale. See # MethylationTable in src/meth.rs. bitvec = "1" +# Data-parallel `methylate`: every contig is methylated independently (its RNGs +# are seeded from `(seed, contig)` with no state carried between contigs), so +# the per-contig loop runs as a work-stealing `par_iter`. See the per-contig +# map in `Methylate::execute` (src/commands/methylate.rs). +rayon = "1" [dev-dependencies] tempfile = "3" diff --git a/src/commands/methylate.rs b/src/commands/methylate.rs index dae7104..e2b7d3f 100644 --- a/src/commands/methylate.rs +++ b/src/commands/methylate.rs @@ -5,6 +5,7 @@ use std::fs::File; use std::io::{BufWriter, Write as _}; use std::path::PathBuf; +use std::sync::Mutex; use anyhow::Result; use clap::Parser; @@ -13,6 +14,12 @@ use rand::SeedableRng as _; use super::command::Command; use super::common::{ReferenceOptions, SeedOptions, VcfOptions}; +/// Serialized methylation output for one contig: its VCF body, and an optional +/// bedGraph body (present only when a bedGraph was requested). Produced +/// independently per contig so the per-contig work can run in parallel and the +/// buffers be concatenated in dict order. See [`methylate_contig`]. +type ContigOutput = (Vec, Option>); + /// Generate a methylation-annotated VCF. /// /// Scans every CpG dinucleotide in the reference (optionally after applying @@ -112,13 +119,11 @@ impl Methylate { impl Command for Methylate { fn execute(&self) -> Result<()> { use crate::fasta::Fasta; - use crate::haplotype::build_haplotypes; - use crate::meth::ContigMethylation; - use crate::seed::{derive_seed, resolve_seed}; - use crate::vcf::methylation::write_contig; + use crate::seed::resolve_seed; use crate::vcf::methylation::write_vcf_header; use crate::vcf::writer::VcfWriter; use crate::version::VERSION; + use rayon::prelude::*; // 1. Build and validate the methylation model from the per-context flags. let model = self.methylation_model(); @@ -141,8 +146,11 @@ impl Command for Methylate { let seed = resolve_seed(self.seed.seed, &seed_desc); log::info!("Using random seed: {seed}"); - // 3. Open the reference and collect contig names. - let mut fasta = Fasta::from_path(&self.reference.reference)?; + // 3. Open the reference and collect contig names. Only the dictionary + // is needed here; the per-contig sequence is loaded per worker below + // (each opens its own handle, since `load_contig` seeks a shared + // reader). + let fasta = Fasta::from_path(&self.reference.reference)?; let dict = fasta.dict().clone(); let contig_names: Vec = dict.names().into_iter().map(String::from).collect(); @@ -194,59 +202,96 @@ impl Command for Methylate { }; let sample_ploidy = parsed_variants.as_ref().map_or(2, |p| p.sample_ploidy); - // 7. For each contig, build haplotypes, compute methylation, emit rows. - for contig_name in &contig_names { - // Per-contig deterministic seed for reference normalization. - // Use the bare contig name (no prefix) to match `simulate`'s - // `derive_seed(main_seed, contig_name)` at simulate.rs line ~578, - // so both commands resolve IUPAC codes identically for the same - // `--seed` and reference. - let ref_seed = derive_seed(seed, contig_name); - let mut ref_rng = rand::rngs::SmallRng::seed_from_u64(ref_seed); - let reference = fasta.load_contig(contig_name, &mut ref_rng)?; - - // Look this contig's variants up from the once-parsed map. An - // absent key or `None` map (no VCF) both yield an empty slice. - let variants: &[crate::vcf::genotype::VariantRecord] = parsed_variants - .as_ref() - .and_then(|p| p.by_contig.get(contig_name)) - .map_or(&[], Vec::as_slice); - - // Build haplotypes with their own deterministic sub-seed, sized - // by the whole-VCF `sample_ploidy` so variant-free contigs land - // the same number of haplotypes as variant-bearing ones. - let hap_seed = derive_seed(seed, &format!("haps@{contig_name}")); - let mut hap_rng = rand::rngs::SmallRng::seed_from_u64(hap_seed); - let haplotypes = build_haplotypes(variants, sample_ploidy, &mut hap_rng); - - // Draw per-haplotype methylation with another deterministic sub-seed. - let meth_seed = derive_seed(seed, &format!("meth@{contig_name}")); - let mut meth_rng = rand::rngs::SmallRng::seed_from_u64(meth_seed); - let methylation = - ContigMethylation::from_haplotypes(&haplotypes, &reference, &model, &mut meth_rng); - - write_contig( - &mut vcf_out, - contig_name, - &reference, - variants, - &methylation, - sample_ploidy, - )?; - - // Append per-contig bedGraph records, if a bedGraph was requested. - // Pass `&haplotypes` so the writer can map every reference CpG - // through each haplotype's local coordinate system before reading - // the methylation bitmap; variant-shifted and variant-destroyed - // CpGs would otherwise be miscounted. - if let Some(bg) = bedgraph_writer.as_mut() { - crate::output::methylation_bedgraph::write_bedgraph_records( - bg, + // 7. Methylate every contig. Each contig is independent: its RNGs are + // seeded purely from `(seed, contig)` with no state carried between + // contigs (the methylation Markov chain runs within a single contig; + // see `meth.rs`), so the loop is a work-stealing `par_iter` whose + // output is identical to the sequential version regardless of thread + // count. Each worker serializes its contig's VCF (and optional + // bedGraph) into its own in-memory buffer; the buffers are written + // out in dict order below so the VCF stays coordinate-sorted. Peak + // memory therefore holds every contig's serialized output at once — + // O(whole reference), not O(one contig) as the old streaming loop + // did. That is fine for hg38-scale references on a multi-GB host + // (a few GB resident); revisit with a bounded ordered-streaming + // writer if a far larger reference must run on a constrained host. + // + // Each worker thread reuses ONE FASTA handle across all of its + // contigs rather than reopening per contig: `Fasta::from_path` + // rebuilds the whole sequence dictionary, so per-contig opens are + // O(contigs²) and dominate the run on references with thousands of + // contigs (e.g. hs38DH's ~3.4k). The handle is only seeked+read by + // `load_contig`, never shared across threads, so reuse is safe and + // deterministic. + let want_bedgraph = self.bedgraph.is_some(); + let ref_path = self.reference.reference.as_path(); + + // One lazily-opened FASTA handle per thread, indexed by + // `rayon::current_thread_index()`. `map_init` cannot be used for this: + // with `with_max_len(1)` below rayon splits the work into one job per + // contig, so `map_init`'s init closure fires once *per contig* rather + // than once per worker — reopening the reference on every contig, the + // exact O(contigs²) cost this handle reuse is meant to avoid. Each + // thread only ever touches its own slot, so every `Mutex` is + // uncontended; it is present solely to make the shared `Vec` `Sync`. + // + // A small workload may run entirely inline on the calling thread rather + // than on the pool, in which case `current_thread_index()` is `None`; + // the extra trailing slot is that caller's handle. Worker indices are + // `0..num_threads`, so they and the fallback slot never collide. + let num_threads = rayon::current_num_threads(); + let fasta_handles: Vec>> = + std::iter::repeat_with(|| Mutex::new(None)).take(num_threads + 1).collect(); + + let contig_outputs: Vec = contig_names + .par_iter() + // One contig per job. Per-contig cost is wildly uneven (chr1 ≫ a + // 50 kb alt) and the large chromosomes sort first, so rayon's + // default contiguous, count-balanced chunking hands one worker the + // whole block of big contigs while the rest drain thousands of + // trivial alts and starve. Forcing single-item jobs lets every + // heavy contig be stolen independently across the pool. + .with_max_len(1) + .map(|contig_name| { + let slot = rayon::current_thread_index().unwrap_or(num_threads); + let mut handle = + fasta_handles[slot].lock().expect("per-worker FASTA handle mutex poisoned"); + // Open this thread's handle on first use, then reuse it across + // every contig the thread is assigned. + if handle.is_none() { + *handle = Some(Fasta::from_path(ref_path)?); + } + let fasta = handle.as_mut().expect("handle initialized above"); + // Look this contig's variants up from the once-parsed map. + // An absent key or `None` map (no VCF) both yield an empty + // slice. + let variants: &[crate::vcf::genotype::VariantRecord] = parsed_variants + .as_ref() + .and_then(|p| p.by_contig.get(contig_name)) + .map_or(&[], Vec::as_slice); + methylate_contig( + fasta, contig_name, - &reference, - &haplotypes, - &methylation, - )?; + seed, + &model, + variants, + sample_ploidy, + want_bedgraph, + ) + }) + // On a multi-contig failure rayon surfaces one failing contig's + // error, not deterministically the first in dict order, so the + // contig named in a fatal error may vary run to run. The job aborts + // either way, so this only affects which contig the message cites. + .collect::>>()?; + + // Write each contig's serialized output in dict order, matching the + // sequential writer exactly: the VCF stays coordinate-sorted and the + // bedGraph contig-ordered. + for (vcf_buf, bedgraph_buf) in &contig_outputs { + vcf_out.write_all(vcf_buf)?; + if let (Some(bg), Some(buf)) = (bedgraph_writer.as_mut(), bedgraph_buf.as_ref()) { + bg.write_all(buf)?; } } @@ -267,6 +312,78 @@ impl Command for Methylate { } } +/// Methylate one contig and return its serialized VCF body (and optional +/// bedGraph body) as in-memory buffers. +/// +/// The contig sequence is loaded through the caller-supplied `fasta` handle +/// (`Fasta::load_contig` seeks the underlying reader, so one handle may be +/// reused across contigs — but a single handle cannot be shared across threads, +/// hence the caller gives each worker its own; see [`Methylate::execute`]). +/// Every RNG is seeded purely from `(seed, contig_name)` with no state carried +/// in from any other contig, so the returned buffers depend only on +/// `(seed, contig_name, model, variants, sample_ploidy)` — never on the handle's +/// prior position or on processing order. Callers may therefore run this across +/// contigs in any order — sequentially or in parallel — and concatenate the +/// buffers in dict order for output identical to a single-threaded run. +fn methylate_contig( + fasta: &mut crate::fasta::Fasta, + contig_name: &str, + seed: u64, + model: &crate::meth::MethylationModel, + variants: &[crate::vcf::genotype::VariantRecord], + sample_ploidy: usize, + want_bedgraph: bool, +) -> Result { + use crate::haplotype::build_haplotypes; + use crate::meth::ContigMethylation; + use crate::seed::derive_seed; + use crate::vcf::methylation::write_contig; + + // Per-contig deterministic seed for reference normalization. Use the bare + // contig name (no prefix) to match `simulate`'s `derive_seed(main_seed, + // contig_name)` at simulate.rs line ~578, so both commands resolve IUPAC + // codes identically for the same `--seed` and reference. + let ref_seed = derive_seed(seed, contig_name); + let mut ref_rng = rand::rngs::SmallRng::seed_from_u64(ref_seed); + let reference = fasta.load_contig(contig_name, &mut ref_rng)?; + + // Build haplotypes with their own deterministic sub-seed, sized by the + // whole-VCF `sample_ploidy` so variant-free contigs land the same number of + // haplotypes as variant-bearing ones. + let hap_seed = derive_seed(seed, &format!("haps@{contig_name}")); + let mut hap_rng = rand::rngs::SmallRng::seed_from_u64(hap_seed); + let haplotypes = build_haplotypes(variants, sample_ploidy, &mut hap_rng); + + // Draw per-haplotype methylation with another deterministic sub-seed. + let meth_seed = derive_seed(seed, &format!("meth@{contig_name}")); + let mut meth_rng = rand::rngs::SmallRng::seed_from_u64(meth_seed); + let methylation = + ContigMethylation::from_haplotypes(&haplotypes, &reference, model, &mut meth_rng); + + let mut vcf_buf = Vec::new(); + write_contig(&mut vcf_buf, contig_name, &reference, variants, &methylation, sample_ploidy)?; + + // Append per-contig bedGraph records, if requested. Pass `&haplotypes` so + // the writer can map every reference CpG through each haplotype's local + // coordinate system before reading the methylation bitmap; variant-shifted + // and variant-destroyed CpGs would otherwise be miscounted. + let bedgraph_buf = if want_bedgraph { + let mut buf = Vec::new(); + crate::output::methylation_bedgraph::write_bedgraph_records( + &mut buf, + contig_name, + &reference, + &haplotypes, + &methylation, + )?; + Some(buf) + } else { + None + }; + + Ok((vcf_buf, bedgraph_buf)) +} + /// Capture the current process's command-line arguments as a space-joined /// string, using lossy UTF-8 for non-Unicode arguments. fn capture_command_line() -> String { diff --git a/src/vcf/methylation.rs b/src/vcf/methylation.rs index ba59877..1c161a8 100644 --- a/src/vcf/methylation.rs +++ b/src/vcf/methylation.rs @@ -173,7 +173,12 @@ fn classify_cpgs_with_haplotypes( var_hap_ranges.push((vi, hap_start, hap_end)); } - // Scan materialized haplotype for CpG dinucleotides. + // Scan materialized haplotype for CpG dinucleotides. `var_cursor` walks + // `var_hap_ranges` (sorted, disjoint) in lock-step with the ascending + // scan so `base_source` is O(1) amortized rather than O(variants) per + // CpG; it is reset per haplotype because `var_hap_ranges` is rebuilt + // above. + let mut var_cursor = 0usize; for h in 0..len - 1 { let c0 = hap_bases[h].to_ascii_uppercase(); let c1 = hap_bases[h + 1].to_ascii_uppercase(); @@ -186,9 +191,10 @@ fn classify_cpgs_with_haplotypes( #[expect(clippy::cast_possible_truncation, reason = "haplotype length fits in u32")] let hpos = h as u32; - // Determine source of each base: variant or reference. - let src_h = base_source(hpos, &var_hap_ranges); - let src_h1 = base_source(hpos + 1, &var_hap_ranges); + // Determine source of each base: variant or reference. Query `hpos` + // before `hpos + 1` so the shared cursor sees non-decreasing inputs. + let src_h = base_source(hpos, &var_hap_ranges, &mut var_cursor); + let src_h1 = base_source(hpos + 1, &var_hap_ranges, &mut var_cursor); let placement = match (src_h, src_h1) { // Both bases come from reference → standalone at the ref @@ -303,17 +309,103 @@ enum BaseSource { } /// Determine whether haplotype coordinate `h` falls inside any variant's -/// alt-allele span. Returns `BaseSource::Variant` for the first matching -/// range, or `BaseSource::Reference` if no variant spans `h`. -fn base_source(h: u32, var_hap_ranges: &[(usize, u32, u32)]) -> BaseSource { - for &(vi, hap_start, hap_end) in var_hap_ranges { - if h >= hap_start && h < hap_end { - return BaseSource::Variant { variant_index: vi, hap_start }; - } +/// alt-allele span, returning `BaseSource::Variant` for the spanning range or +/// `BaseSource::Reference` if none spans `h`. +/// +/// `var_hap_ranges` is sorted by `hap_start` and the ranges are disjoint (one +/// haplotype never carries two overlapping alt spans — enforced upstream by +/// [`check_overlaps_on_shared_haplotypes`]). `cursor` is a monotonic index into +/// `var_hap_ranges` that the caller threads across a strictly **non-decreasing** +/// sequence of `h` queries: it advances past every range ending at or before +/// `h` and never rewinds. This turns what would be an O(ranges) scan per query +/// into O(1) amortized across a full CpG scan, so the per-haplotype CpG +/// classification is O(haplotype length + variants) instead of O(CpGs × +/// variants). +/// +/// The caller must (a) reset `cursor` to `0` before each independent ascending +/// scan and (b) only ever pass `h` values that do not decrease within that +/// scan. Querying `h` then `h + 1` within one loop iteration, with `h` +/// incrementing by one across iterations, satisfies this. +fn base_source(h: u32, var_hap_ranges: &[(usize, u32, u32)], cursor: &mut usize) -> BaseSource { + // Skip ranges that end at or before `h`; they can never contain `h` and, + // because queries are non-decreasing, can never contain any later query + // either, so dropping them permanently is safe. + while *cursor < var_hap_ranges.len() && var_hap_ranges[*cursor].2 <= h { + *cursor += 1; + } + if let Some(&(vi, hap_start, hap_end)) = var_hap_ranges.get(*cursor) + && h >= hap_start + && h < hap_end + { + return BaseSource::Variant { variant_index: vi, hap_start }; } BaseSource::Reference } +#[cfg(test)] +mod base_source_cursor_tests { + //! The monotonic-cursor [`base_source`] must return exactly what a naive + //! linear scan would for the same query, under the real caller access + //! pattern (query `h` then `h + 1` per CpG, `h` ascending by one). This + //! pins the O(CpGs + variants) fast path to the O(CpGs × variants) + //! reference semantics it replaced. + + use super::{BaseSource, base_source}; + use rand::rngs::SmallRng; + use rand::{Rng as _, SeedableRng as _}; + + /// Naive reference implementation: the first (only, since disjoint) range + /// spanning `h`, as `(variant_index, hap_start)`. + fn naive(h: u32, ranges: &[(usize, u32, u32)]) -> Option<(usize, u32)> { + ranges.iter().find(|&&(_, s, e)| h >= s && h < e).map(|&(vi, s, _)| (vi, s)) + } + + fn to_pair(src: BaseSource) -> Option<(usize, u32)> { + match src { + BaseSource::Variant { variant_index, hap_start } => Some((variant_index, hap_start)), + BaseSource::Reference => None, + } + } + + #[test] + fn cursor_matches_linear_scan_under_real_access_pattern() { + let mut rng = SmallRng::seed_from_u64(0xC0FF_EE99); + for _ in 0..500 { + // Build random disjoint, ascending ranges shaped like the + // per-haplotype alt spans `var_hap_ranges` holds: gaps of reference + // between spans, each span at least one base wide. + let n = rng.random_range(0..12usize); + let mut ranges: Vec<(usize, u32, u32)> = Vec::new(); + let mut pos = 0u32; + for vi in 0..n { + let start = pos + rng.random_range(0..5u32); + let end = start + rng.random_range(1..6u32); + ranges.push((vi, start, end)); + pos = end; + } + let max_h = pos + 6; + + // One cursor threaded across the whole ascending scan, querying + // `h` then `h + 1` per step exactly as `classify_cpgs_with_haplotypes` + // does. + let mut cursor = 0usize; + for h in 0..max_h { + assert_eq!( + to_pair(base_source(h, &ranges, &mut cursor)), + naive(h, &ranges), + "query h={h} ranges={ranges:?}", + ); + assert_eq!( + to_pair(base_source(h + 1, &ranges, &mut cursor)), + naive(h + 1, &ranges), + "query h+1={} ranges={ranges:?}", + h + 1, + ); + } + } + } +} + /// Compute a sort key for a [`CpgPlacement`] so that placements are ordered /// by their effective reference position. Standalone records use their own /// `ref_pos`; `OnVariant` records use the variant's reference position as the @@ -691,7 +783,10 @@ fn per_variant_per_hap_cpg_offsets_with_haplotypes( } // Scan for CpG dinucleotides and assign them to the owning variant - // using the same upstream-wins rule as classify_cpgs. + // using the same upstream-wins rule as classify_cpgs. `var_cursor` + // advances with the ascending scan to keep `base_source` O(1) amortized; + // reset per haplotype alongside the rebuilt `var_hap_ranges`. + let mut var_cursor = 0usize; for h in 0..len - 1 { let c0 = hap_bases[h].to_ascii_uppercase(); let c1 = hap_bases[h + 1].to_ascii_uppercase(); @@ -702,8 +797,10 @@ fn per_variant_per_hap_cpg_offsets_with_haplotypes( #[expect(clippy::cast_possible_truncation, reason = "haplotype length fits in u32")] let hpos = h as u32; - let src_h = base_source(hpos, &var_hap_ranges); - let src_h1 = base_source(hpos + 1, &var_hap_ranges); + // Query `hpos` before `hpos + 1` so the shared cursor advances over + // non-decreasing inputs. + let src_h = base_source(hpos, &var_hap_ranges, &mut var_cursor); + let src_h1 = base_source(hpos + 1, &var_hap_ranges, &mut var_cursor); // Determine the owning variant index using the same upstream-wins // rule as classify_cpgs. None means both bases are from reference diff --git a/tests/test_methylate.rs b/tests/test_methylate.rs index 19f6c32..e27528f 100644 --- a/tests/test_methylate.rs +++ b/tests/test_methylate.rs @@ -528,6 +528,52 @@ fn methylate_applies_variants_from_input_vcf() { ); } +#[test] +fn methylate_output_is_identical_across_thread_counts() { + // Several contigs of varied length so the per-contig parallel map spans + // multiple work units; non-repetitive sequence carries CpGs on each. The + // output must be byte-identical whether methylated on one thread or many: + // each contig's RNGs are seeded purely from `(seed, contig)` with no state + // carried between contigs, so processing order cannot affect the result. + let c1 = helpers::non_repetitive_seq(8000); + let c2 = helpers::non_repetitive_seq(5000); + let c3 = helpers::non_repetitive_seq(3000); + let c4 = helpers::non_repetitive_seq(1500); + let env = TestEnv::new(&[("chr1", &c1), ("chr2", &c2), ("chr3", &c3), ("chr4", &c4)]); + + let run = |tag: &str, threads: usize| -> (String, Vec) { + let vcf = env.dir.path().join(format!("{tag}.vcf.gz")); + let bg = env.dir.path().join(format!("{tag}.bedgraph")); + let (ok, _, stderr) = run_methylate_with_threads( + threads, + &[ + "methylate", + "--reference", + env.fasta_path.to_str().unwrap(), + "--output", + vcf.to_str().unwrap(), + "--bedgraph", + bg.to_str().unwrap(), + "--seed", + "42", + ], + ); + assert!(ok, "methylate (threads={threads}) failed: {stderr}"); + (strip_command_line(&decompress_vcf(&vcf)), std::fs::read(&bg).expect("read bedgraph")) + }; + + let (vcf_single, bg_single) = run("single", 1); + let (vcf_multi, bg_multi) = run("multi", 8); + + assert_eq!(vcf_single, vcf_multi, "VCF content differs between 1 and 8 threads"); + assert_eq!(bg_single, bg_multi, "bedGraph differs between 1 and 8 threads"); + // Guard against a vacuous pass on two equal-but-empty outputs. + assert!( + vcf_single.lines().any(|l| l.starts_with("chr1\t") && l.contains("MT:MB")), + "expected CpG methylation rows in the output" + ); +} + // ── Command runner ─────────────────────────────────────────────────────────── /// Run `holodeck methylate` with the given arguments and return @@ -542,6 +588,20 @@ fn run_methylate(args: &[&str]) -> (bool, String, String) { (output.status.success(), stdout, stderr) } +/// Run `holodeck methylate` with `RAYON_NUM_THREADS` pinned, returning +/// `(success, stdout, stderr)`. Used to prove per-contig parallelism does not +/// change the output. +fn run_methylate_with_threads(threads: usize, args: &[&str]) -> (bool, String, String) { + let output = std::process::Command::new(env!("CARGO_BIN_EXE_holodeck")) + .env("RAYON_NUM_THREADS", threads.to_string()) + .args(args) + .output() + .expect("Failed to run holodeck"); + let stdout = String::from_utf8_lossy(&output.stdout).to_string(); + let stderr = String::from_utf8_lossy(&output.stderr).to_string(); + (output.status.success(), stdout, stderr) +} + /// Decompress a BGZF-compressed VCF into a String. The methylate command /// always writes its output as BGZF; tests that need to inspect the text /// use this helper.