feat: split methylate from simulate; methylation truth via VCF MT/MB#11
Merged
Conversation
|
Note Reviews pausedIt looks like this branch is under active development. To avoid overwhelming you with review comments due to an influx of new commits, CodeRabbit has automatically paused this review. You can configure this behavior by changing the Use the following commands to manage reviews:
Use the checkboxes below for quick actions:
📝 WalkthroughWalkthroughAdds a methylation feature: a new ChangesMethylation feature
Sequence Diagram(s)sequenceDiagram
participant CLI as holodeck CLI
participant Methylate as methylate subcommand
participant VCF as MT/MB VCF (BGZF)
participant Simulate as simulate subcommand
participant ReadGen as read generation (generate_read_pair)
participant BAM as Golden BAM writer
CLI->>Methylate: run methylate --reference --methylation-model --output
Methylate->>VCF: write MT/MB records (per-haplotype)
CLI->>Simulate: run simulate --methylation-mode --vcf MT/MB.vcf
Simulate->>VCF: parse MT/MB once, build ContigMethylation
Simulate->>ReadGen: generate_read_pair(..., methylation_config)
ReadGen->>ReadGen: apply_methylation_conversion (fragment-scale)
ReadGen->>BAM: provide methylation annotations
BAM->>BAM: emit XG/XR/YS and optional XM/YM/NM/MD tags
Simulate->>Simulate: record per-mate CpG tallies
Simulate->>File: write cpg-truth bedGraph
Estimated code review effort🎯 5 (Critical) | ⏱️ ~120 minutes Possibly related PRs
✨ Finishing Touches📝 Generate docstrings
|
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 2
🧹 Nitpick comments (1)
README.md (1)
155-155: ⚡ Quick winConsider clarifying "reference orientation" for the
YS:Ztag.The phrase "the pre-conversion read sequence in reference orientation" could be ambiguous. For a read derived from the bottom strand, does "reference orientation" mean:
- The sequence as it would appear in the forward reference strand?
- The sequence before reverse-complementing for BAM storage?
- Something else?
Consider adding a brief clarification, such as: "in reference (forward-strand) orientation" or "oriented to the positive strand" to remove ambiguity.
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@README.md` at line 155, Clarify the ambiguous phrase describing the YS:Z tag by specifying what "reference orientation" means: update the README line for `YS:Z` to say it is the pre-conversion read sequence oriented to the forward (positive) reference strand (e.g., "in reference (forward-strand) orientation" or "oriented to the positive strand") so readers understand how `YS:Z` relates to `SEQ` for bottom-strand-derived reads; keep the note that no bisulfite aligner uses this tag and that it exists to allow downstream evaluators to diff `SEQ` against `YS` base-for-base to recover ground-truth chemistry events.
🤖 Prompt for all review comments with AI agents
Verify each finding against the current code and only fix it if needed.
Inline comments:
In `@src/commands/simulate.rs`:
- Around line 204-210: The validations for methylation_rate and
methylation_conversion_rate should only run when a methylation mode is enabled;
modify the logic around the existing checks (methylation_rate,
methylation_conversion_rate) to first test methylation_mode (e.g., if let
Some(_) = self.methylation_mode or self.methylation_mode.is_some()) and perform
the finite/0..=1.0 validation only in that branch, otherwise reject or warn when
the user supplied those flags without a mode (bail with a clear message asking
for --methylation-mode) so the flags aren’t silently ignored.
- Around line 521-544: The code builds a single methylation table from the
normalized reference (methylation / methylation_config using derive_seed,
main_seed, contig_name and self.methylation_mode) but later samples reads from
per-haplotype sequences (haplotypes[hap_idx]), so CpG sites created/destroyed by
variants will be mis-modeled; either make methylation table construction
haplotype-aware (build a table per haplotype using the haplotype sequence used
for reads) or, if you can't, add an explicit guard that rejects running with
both VCF-driven haplotypes and methylation enabled: detect the combination
(self.methylation_mode.is_some() && vcf/VCF-provided haplotypes in this
command), and return an error/exit with a clear message instead of silently
proceeding; update the codepaths that create methylation/methylation_config (the
block that calls crate::meth::MethylationTable::from_reference and uses
derive_seed/main_seed/contig_name) to implement this check.
---
Nitpick comments:
In `@README.md`:
- Line 155: Clarify the ambiguous phrase describing the YS:Z tag by specifying
what "reference orientation" means: update the README line for `YS:Z` to say it
is the pre-conversion read sequence oriented to the forward (positive) reference
strand (e.g., "in reference (forward-strand) orientation" or "oriented to the
positive strand") so readers understand how `YS:Z` relates to `SEQ` for
bottom-strand-derived reads; keep the note that no bisulfite aligner uses this
tag and that it exists to allow downstream evaluators to diff `SEQ` against `YS`
base-for-base to recover ground-truth chemistry events.
🪄 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: defaults
Review profile: CHILL
Plan: Pro
Run ID: f0879437-524f-43ea-85e6-5ca3af1a6437
⛔ Files ignored due to path filters (1)
Cargo.lockis excluded by!**/*.lock
📒 Files selected for processing (9)
Cargo.tomlREADME.mdsrc/commands/simulate.rssrc/lib.rssrc/meth.rssrc/output/golden_bam.rssrc/read.rstests/helpers/mod.rstests/test_simulate_meth.rs
nh13
force-pushed
the
nh/bisulfite-conversion
branch
from
May 3, 2026 06:08
d1c1a5a to
460d211
Compare
There was a problem hiding this comment.
Actionable comments posted: 3
🤖 Prompt for all review comments with AI agents
Verify each finding against the current code and only fix it if needed.
Inline comments:
In `@src/haplotype.rs`:
- Around line 185-195: hap_position_for currently does a full linear scan over
self.variant_data for every call, causing O(fragments×variants) cost; fix by
precomputing a cumulative delta index and using a binary search instead of
scanning. Add a cached vector (e.g. cumulative_deltas or index_points)
computed/updated whenever variant_data changes that stores for each variant the
ref_end (var.ref_pos+var.ref_len) and the running sum of (alt_bases.len() -
ref_len) as i64, then change hap_position_for to binary_search the largest
ref_end <= ref_pos and read the precomputed cumulative delta to compute hp;
reference the methods/fields variant_data and hap_position_for (and update
construction/update points where variant_data is mutated, e.g. where
extract_fragment or hap_start logic may rely on variant_data) so the cache stays
in sync.
- Around line 127-129: The hap_start is computed using the original ref_start
which is wrong when the pre-loop advances ref_pos past a deletion; change the
computation to call self.hap_position_for using the post-skip coordinate
(ref_pos) instead of ref_start so downstream methylation lookups align with the
actual first surviving reference base—update the hap_start assignment to use
ref_pos and keep references to hap_position_for, ref_start, ref_pos and the
pre-loop behavior in mind when making the change.
In `@src/meth.rs`:
- Around line 96-113: The current cap calculation (let cap =
reference.len().saturating_add(1)) can truncate haplotypes that contain >1 base
of net insertion because extract_fragment(reference, 0, cap) limits emitted
haplotype bases; change the cap so it will not artificially truncate inserted
sequence (e.g. pass a very large cap or usize::MAX) so extract_fragment can walk
until it consumes the whole reference; update the code around cap,
haplotype.extract_fragment, and any related comments so hap_bases and _hap_start
reflect the full materialized haplotype rather than being truncated.
🪄 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: defaults
Review profile: CHILL
Plan: Pro
Run ID: 6785f185-64c8-4955-b973-31e4dd83282e
⛔ Files ignored due to path filters (1)
Cargo.lockis excluded by!**/*.lock
📒 Files selected for processing (11)
Cargo.tomlREADME.mdsrc/commands/simulate.rssrc/fragment.rssrc/haplotype.rssrc/lib.rssrc/meth.rssrc/output/golden_bam.rssrc/read.rstests/helpers/mod.rstests/test_simulate_meth.rs
✅ Files skipped from review due to trivial changes (3)
- src/lib.rs
- Cargo.toml
- README.md
🚧 Files skipped from review as they are similar to previous changes (3)
- src/read.rs
- tests/helpers/mod.rs
- src/output/golden_bam.rs
nh13
force-pushed
the
nh/bisulfite-conversion
branch
2 times, most recently
from
May 4, 2026 05:00
ae0f48e to
4576552
Compare
There was a problem hiding this comment.
Actionable comments posted: 3
🤖 Prompt for all review comments with AI agents
Verify each finding against the current code and only fix it if needed.
Inline comments:
In `@src/commands/simulate.rs`:
- Around line 212-232: The current guard uses float equality on
self.methylation_rate and self.methylation_conversion_rate (via rate_overridden
and conversion_overridden) so explicitly passing “1.0” bypasses validation;
change the logic to check whether the flags were supplied rather than their
numeric value by introducing presence booleans (e.g. track flag presence from
the CLI parser or change the struct fields to Option<f64> and use is_some()) for
methylation_rate and methylation_conversion_rate and then require
self.methylation_mode.is_some() if either presence flag is true; update the
rejection branch that bails with "--methylation-rate and
--methylation-conversion-rate require --methylation-mode" and add regression
tests covering explicitly supplied 1.0 for both methylation_rate and
methylation_conversion_rate to ensure the new presence-based check catches them.
- Around line 128-137: The help text for the methylation probability incorrectly
states methylation is drawn "in the reference" once per genomic position; update
the documentation in simulate.rs to explain methylation is sampled per
haplotype/allele (not just the reference) to match the implementation that uses
ContigMethylation::from_haplotypes(...). Edit the comment block describing the
flag (the paragraph starting "Probability that any CpG-context cytosine...") to
say the methylation state is drawn separately for each haplotype/allele at each
genomic position (deterministic from --seed), so allele-specific/variant-aware
methylation and hemimethylation are correctly communicated to CLI users.
In `@tests/test_simulate_meth.rs`:
- Around line 154-187: The empirical-band test in
test_em_seq_partial_conversion_rate is flaky because it relies on compute_seed
via the temp FASTA path; explicitly pass a fixed seed to run_simulate to
stabilize the RNG (e.g., add the "--seed" and "42" arguments to the run_simulate
args vector used in test_em_seq_partial_conversion_rate) so the derived band
matches the comment's seed; update the test invocation in the
test_em_seq_partial_conversion_rate function accordingly.
🪄 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: defaults
Review profile: CHILL
Plan: Pro
Run ID: 019b0da5-6612-444d-8c5c-c01253dab386
⛔ Files ignored due to path filters (1)
Cargo.lockis excluded by!**/*.lock
📒 Files selected for processing (11)
Cargo.tomlREADME.mdsrc/commands/simulate.rssrc/fragment.rssrc/haplotype.rssrc/lib.rssrc/meth.rssrc/output/golden_bam.rssrc/read.rstests/helpers/mod.rstests/test_simulate_meth.rs
✅ Files skipped from review due to trivial changes (2)
- Cargo.toml
- src/fragment.rs
🚧 Files skipped from review as they are similar to previous changes (6)
- src/lib.rs
- tests/helpers/mod.rs
- src/output/golden_bam.rs
- README.md
- src/meth.rs
- src/haplotype.rs
nh13
force-pushed
the
nh/bisulfite-conversion
branch
from
May 5, 2026 03:19
4576552 to
37592f0
Compare
There was a problem hiding this comment.
♻️ Duplicate comments (3)
src/commands/simulate.rs (2)
128-137:⚠️ Potential issue | 🟡 Minor | ⚡ Quick winDocument the per-haplotype sampling model.
This help text still says methylation is drawn once “in the reference” per genomic position, but the implementation now samples per haplotype × strand CpG via
ContigMethylation::from_haplotypes(...). That hides allele-specific methylation and variant-created/destroyed CpGs from CLI users.✏️ Suggested wording
- /// Probability that any CpG-context cytosine in the reference is - /// methylated. Drawn once per genomic position at simulation start - /// (deterministic from `--seed`). The two C's of each CpG site are - /// drawn independently per strand, so hemimethylation is possible. + /// Probability that any CpG-context cytosine on a haplotype is + /// methylated. Drawn once per haplotype × strand CpG site at + /// simulation start (deterministic from `--seed`), so allele-specific + /// methylation and variant-created/destroyed CpGs are modeled. The two + /// C's of each CpG site are drawn independently per strand, so + /// hemimethylation is possible.🤖 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/simulate.rs` around lines 128 - 137, Update the long help text for the methylation probability flag to document that methylation is sampled per haplotype × strand (not once "in the reference") by the ContigMethylation::from_haplotypes model, so allele-specific methylation and CpGs created/destroyed by variants are captured; change any wording that says "drawn once per genomic position in the reference" to explicitly say "sampled per haplotype and strand (via ContigMethylation::from_haplotypes), deterministic from --seed" and retain the notes about hemimethylation, non-CpG behavior, default value, and examples.
224-244:⚠️ Potential issue | 🟠 Major | ⚡ Quick winReject explicitly supplied default-valued methylation flags too.
rate_overriddenandconversion_overriddenonly detect numeric changes, so--methylation-rate 1.0and--methylation-conversion-rate 1.0still slip through without--methylation-mode. Those invocations silently fall back to the unmethylated path even though this PR advertises that supplying either rate flag without a mode is invalid. Please track flag presence instead of comparing against1.0, and add regressions for the explicit-1.0cases.🤖 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/simulate.rs` around lines 224 - 244, rate_overridden/conversion_overridden currently compare numeric values to 1.0 and therefore miss cases where the user explicitly passed “--methylation-rate 1.0” or “--methylation-conversion-rate 1.0”; change the flag tracking so you detect presence instead of value: make the parsed inputs for methylation_rate and methylation_conversion_rate carry presence (e.g. Option<f64> or a companion bool like methylation_rate_provided/methylation_conversion_rate_provided set by the CLI parser), then replace the comparisons with checks for presence (e.g. let rate_provided = self.methylation_rate.is_some() or self.methylation_rate_provided; let conversion_provided = ... ) and keep the existing guard that if self.methylation_mode.is_none() && (rate_provided || conversion_provided) bail(...). Also add unit/regression tests that call the CLI with explicit “1.0” for each rate flag (and combinations) to verify they are rejected when methylation_mode is absent.tests/test_simulate_meth.rs (1)
154-187:⚠️ Potential issue | 🟠 Major | ⚡ Quick winSeed this empirical-band test explicitly.
The comment says this band was derived from seed
42, but the command omits--seed. Becausecompute_seed()hashes the temp FASTA path, the effective RNG varies between runs and can make this assertion flaky.🔧 Minimal fix
"--methylation-conversion-rate", "0.5", + "--seed", + "42", "--threads", "1",🤖 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 `@tests/test_simulate_meth.rs` around lines 154 - 187, The test test_em_seq_partial_conversion_rate is flaky because the run_simulate call omits an explicit RNG seed (compute_seed() hashes the temp FASTA path), so add a deterministic seed flag (e.g., "--seed", "42") to the run_simulate argument list used in this test to match the empirical band derived with SmallRng::seed_from_u64(42); update the run_simulate invocation in test_em_seq_partial_conversion_rate to include the seed so the RNG output is stable across runs.
🤖 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.
Duplicate comments:
In `@src/commands/simulate.rs`:
- Around line 128-137: Update the long help text for the methylation probability
flag to document that methylation is sampled per haplotype × strand (not once
"in the reference") by the ContigMethylation::from_haplotypes model, so
allele-specific methylation and CpGs created/destroyed by variants are captured;
change any wording that says "drawn once per genomic position in the reference"
to explicitly say "sampled per haplotype and strand (via
ContigMethylation::from_haplotypes), deterministic from --seed" and retain the
notes about hemimethylation, non-CpG behavior, default value, and examples.
- Around line 224-244: rate_overridden/conversion_overridden currently compare
numeric values to 1.0 and therefore miss cases where the user explicitly passed
“--methylation-rate 1.0” or “--methylation-conversion-rate 1.0”; change the flag
tracking so you detect presence instead of value: make the parsed inputs for
methylation_rate and methylation_conversion_rate carry presence (e.g.
Option<f64> or a companion bool like
methylation_rate_provided/methylation_conversion_rate_provided set by the CLI
parser), then replace the comparisons with checks for presence (e.g. let
rate_provided = self.methylation_rate.is_some() or
self.methylation_rate_provided; let conversion_provided = ... ) and keep the
existing guard that if self.methylation_mode.is_none() && (rate_provided ||
conversion_provided) bail(...). Also add unit/regression tests that call the CLI
with explicit “1.0” for each rate flag (and combinations) to verify they are
rejected when methylation_mode is absent.
In `@tests/test_simulate_meth.rs`:
- Around line 154-187: The test test_em_seq_partial_conversion_rate is flaky
because the run_simulate call omits an explicit RNG seed (compute_seed() hashes
the temp FASTA path), so add a deterministic seed flag (e.g., "--seed", "42") to
the run_simulate argument list used in this test to match the empirical band
derived with SmallRng::seed_from_u64(42); update the run_simulate invocation in
test_em_seq_partial_conversion_rate to include the seed so the RNG output is
stable across runs.
ℹ️ Review info
⚙️ Run configuration
Configuration used: defaults
Review profile: CHILL
Plan: Pro
Run ID: 5ad8e9d6-ca2a-493b-b546-0b4d5af557ad
⛔ Files ignored due to path filters (1)
Cargo.lockis excluded by!**/*.lock
📒 Files selected for processing (13)
Cargo.tomlREADME.mdsrc/commands/simulate.rssrc/fragment.rssrc/haplotype.rssrc/lib.rssrc/meth.rssrc/output/cpg_truth.rssrc/output/golden_bam.rssrc/output/mod.rssrc/read.rstests/helpers/mod.rstests/test_simulate_meth.rs
✅ Files skipped from review due to trivial changes (4)
- src/output/mod.rs
- src/lib.rs
- README.md
- Cargo.toml
🚧 Files skipped from review as they are similar to previous changes (5)
- tests/helpers/mod.rs
- src/fragment.rs
- src/read.rs
- src/meth.rs
- src/output/golden_bam.rs
nh13
force-pushed
the
nh/bisulfite-conversion
branch
from
May 5, 2026 05:29
37592f0 to
71de073
Compare
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
nh13
force-pushed
the
nh/bisulfite-conversion
branch
from
May 9, 2026 17:57
05af86a to
d240816
Compare
nh13
changed the title
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
There was a problem hiding this comment.
Actionable comments posted: 2
🤖 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 115-117: Update the README option descriptions for
--methylation-mode, --methylation-rate and --methylation-conversion-rate to
state that methylation is sampled per haplotype × strand (not once per reference
position), and that this enables allele-specific methylation and
variant-created/destroyed CpGs to be modeled; keep the existing default values
and chemistry names (em-seq, bisulfite, taps) but change the --methylation-rate
phrasing from "drawn once per genomic position" to "per-CpG methylation
probability sampled independently for each haplotype and strand" and add a brief
note that this allows allele-specific methylation patterns and variant-induced
CpG changes to be captured.
In `@src/commands/simulate.rs`:
- Around line 242-244: The code only discovers a missing parent directory for
--cpg-truth-bedgraph at write_bedgraph() after the expensive simulation; add an
early check right after option parsing (in the simulate command setup in
src/commands/simulate.rs) that, when self.cpg_truth_bedgraph.is_some(), verifies
the file's parent directory exists and is writable (e.g., obtain
Path::new(&self.cpg_truth_bedgraph.unwrap()).parent(), return an error via bail!
if parent is None or !parent.exists() or !is_writable); do the same guard where
similar outputs are validated (the other blocks referenced around lines ~264-270
and ~409-411) so we fail fast instead of discovering the problem in
write_bedgraph().
🪄 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: defaults
Review profile: CHILL
Plan: Pro
Run ID: 2e6c4956-50d3-436d-918b-001bf54d3e36
⛔ Files ignored due to path filters (1)
Cargo.lockis excluded by!**/*.lock
📒 Files selected for processing (13)
Cargo.tomlREADME.mdsrc/commands/simulate.rssrc/fragment.rssrc/haplotype.rssrc/lib.rssrc/meth.rssrc/output/cpg_truth.rssrc/output/golden_bam.rssrc/output/mod.rssrc/read.rstests/helpers/mod.rstests/test_simulate_meth.rs
✅ Files skipped from review due to trivial changes (2)
- Cargo.toml
- src/output/mod.rs
🚧 Files skipped from review as they are similar to previous changes (9)
- tests/helpers/mod.rs
- src/lib.rs
- src/output/cpg_truth.rs
- src/fragment.rs
- src/output/golden_bam.rs
- src/read.rs
- src/meth.rs
- src/haplotype.rs
- tests/test_simulate_meth.rs
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
There was a problem hiding this comment.
Actionable comments posted: 2
♻️ Duplicate comments (2)
tests/test_simulate_meth.rs (1)
161-187:⚠️ Potential issue | 🟡 Minor | ⚡ Quick winSeed this empirical-band test explicitly.
The comment says the
0.20..0.30band was derived from seed42, but the invocation never passes--seed. That makes the assertion depend on whatever auto-seeding path the CLI currently uses.Suggested fix
"--methylation-conversion-rate", "0.5", + "--seed", + "42", "--threads", "1",🤖 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 `@tests/test_simulate_meth.rs` around lines 161 - 187, The test invocation in tests/test_simulate_meth.rs calls run_simulate(...) without a fixed seed, but the assertion relies on an empirical methylation band derived from seed 42; add a deterministic seed flag to the CLI args (e.g. include "--seed" and "42" in the argument list passed to run_simulate) so the empirical-band assertion remains stable; update the run_simulate call where the arguments array is constructed to include the seed parameter.README.md (1)
115-117:⚠️ Potential issue | 🟡 Minor | ⚡ Quick winDocument methylation sampling as per haplotype × strand, not per reference position.
The README still says the draw happens once per genomic position / “on the reference”, but the implementation samples each haplotype and strand independently. As written, this hides allele-specific methylation and variant-created/destroyed CpGs.
Suggested wording
-| `--methylation-rate` | 1.0 | Per-CpG methylation probability (drawn once per genomic position) | +| `--methylation-rate` | 1.0 | Per-CpG methylation probability sampled independently for each haplotype and strand |-CpG dinucleotides on the reference are independently methylated on each strand (top and bottom) with probability `--methylation-rate` (default `1.0`); the draw happens once per genomic position at simulation start and is deterministic from `--seed`. +CpG dinucleotides on each haplotype are independently methylated on each strand (top and bottom) with probability `--methylation-rate` (default `1.0`); the draw happens once per haplotype × strand CpG site at simulation start and is deterministic from `--seed`. This allows allele-specific methylation and variant-created/destroyed CpGs to be modeled correctly.Also applies to: 151-151
🤖 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 `@README.md` around lines 115 - 117, Update the README description for methylation sampling to state that methylation is sampled independently per haplotype and strand (not once per reference position); specifically revise the lines describing `--methylation-rate` and `--methylation-conversion-rate` to indicate "Per-CpG methylation probability drawn independently for each haplotype × strand" and mention that this allows allele-specific methylation and variant-created/destroyed CpGs; apply the same wording change to the other occurrence noted around line 151.
🤖 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 185-187: The third invariant is incorrect: samtools calmd
recomputes standard SAM tags and does not preserve Bismark-style suppression of
allowed C>T / G>A events, so the README should be updated to avoid recommending
calmd as a validator; replace or reword the `samtools calmd` bullet (the text
mentioning `samtools calmd`, `MD:Z`, and `NM:i`) to state that calmd will
recompute standard NM/MD and can yield false failures on bisulfite-suppressed
BAMs and instead recommend using a bisulfite-aware comparator or a custom
validation script that treats C>T/G>A at CpGs as allowed (or simply remove the
calmd check).
In `@src/methylation_tags.rs`:
- Around line 379-382: The branch in populate_pair_call_tags that only sets
annotation.r2_call_tags when (read2, r2_truth, r2_cigar) are all Some(...)
leaves previous r2_call_tags intact when any of those are None; update the logic
so annotation.r2_call_tags is explicitly cleared (set to None) whenever the
tuple pattern does not match, i.e. when any of read2, r2_truth, or r2_cigar is
None, instead of leaving stale data—keep the compute_one_mate(...) assignment
for the Some case and add an explicit annotation.r2_call_tags = None in the
else/non-matching path to ensure stored tags are always overwritten.
---
Duplicate comments:
In `@README.md`:
- Around line 115-117: Update the README description for methylation sampling to
state that methylation is sampled independently per haplotype and strand (not
once per reference position); specifically revise the lines describing
`--methylation-rate` and `--methylation-conversion-rate` to indicate "Per-CpG
methylation probability drawn independently for each haplotype × strand" and
mention that this allows allele-specific methylation and
variant-created/destroyed CpGs; apply the same wording change to the other
occurrence noted around line 151.
In `@tests/test_simulate_meth.rs`:
- Around line 161-187: The test invocation in tests/test_simulate_meth.rs calls
run_simulate(...) without a fixed seed, but the assertion relies on an empirical
methylation band derived from seed 42; add a deterministic seed flag to the CLI
args (e.g. include "--seed" and "42" in the argument list passed to
run_simulate) so the empirical-band assertion remains stable; update the
run_simulate call where the arguments array is constructed to include the seed
parameter.
🪄 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: defaults
Review profile: CHILL
Plan: Pro
Run ID: c4715dcd-6fd9-45cf-951b-313167a47369
⛔ Files ignored due to path filters (1)
Cargo.lockis excluded by!**/*.lock
📒 Files selected for processing (14)
Cargo.tomlREADME.mdsrc/commands/simulate.rssrc/fragment.rssrc/haplotype.rssrc/lib.rssrc/meth.rssrc/methylation_tags.rssrc/output/cpg_truth.rssrc/output/golden_bam.rssrc/output/mod.rssrc/read.rstests/helpers/mod.rstests/test_simulate_meth.rs
✅ Files skipped from review due to trivial changes (3)
- src/lib.rs
- Cargo.toml
- src/output/mod.rs
🚧 Files skipped from review as they are similar to previous changes (7)
- src/fragment.rs
- tests/helpers/mod.rs
- src/output/cpg_truth.rs
- src/haplotype.rs
- src/meth.rs
- src/commands/simulate.rs
- src/read.rs
nh13
force-pushed
the
nh/bisulfite-conversion
branch
from
May 10, 2026 14:09
d240816 to
eba73ea
Compare
Contributor
Author
|
@coderabbitai resume |
✅ Actions performedReviews resumed. |
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
✅ Actions performedReview triggered.
|
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Contributor
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Real bisulfite/EM-seq conversion is bimodal at the molecule level: most
molecules convert near-completely while a small fraction (fragments that
fail to denature or re-anneal too fast) escape conversion as a coherent
unit. The flat per-cytosine rate could not reproduce this, which is the
exact failure mode a conversion-failure detector must be benchmarked
against.
Draw a per-fragment camp once before the per-base loop: with probability
--methylation-failure-rate the molecule fails and converts its
should-convert cytosines at 1 - conversion_rate (near-zero), otherwise at
conversion_rate. The draw lives in apply_methylation_conversion, which runs
once per fragment, so both mates derive from the same converted buffer and
stay coherent. The golden BAM stamps every record with cf:i:{0|1} -- the
one piece of failure ground truth not recoverable from SEQ/YS alone.
Defaults now reflect realistic chemistry (the feature is unreleased, so no
byte-identical compatibility is owed): conversion-rate 1.0 -> 0.999,
failure-rate -> 0.01. Pass --methylation-failure-rate 0.0 with
--methylation-conversion-rate 1.0 for perfectly deterministic conversion.
Contributor
Author
|
@coderabbitai review |
✅ Action performedReview finished.
|
There was a problem hiding this comment.
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (2)
src/meth.rs (1)
92-109:⚠️ Potential issue | 🟠 Major | 🏗️ Heavy liftReturn typed errors from these public methylation APIs instead of panicking.
MethylationTable::from_haplotypeandapply_methylation_conversionare bothpub, but invalid caller-supplied rates currently tripassert!and abort the process. That turns a recoverable library-input error into a crash path for external callers. Please move this validation onto a typed error surface (or an infallible validated config/newtype) and threadResultthrough the callers/tests instead. As per coding guidelines, "Use anyhow for application errors and thiserror for library error types".Also applies to: 327-353
🤖 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/meth.rs` around lines 92 - 109, Replace the panic-based validation in MethylationTable::from_haplotype (and the same pattern in apply_methylation_conversion) with a typed error return: remove the assert! and validate methylation_rate, returning Result<Self, ErrorType> where ErrorType is a library error defined with thiserror (e.g., MethylationError::InvalidRate { value: f64 }), or alternatively accept a validated newtype; update the function signatures to return Result and propagate that Result through callers and tests using anyhow at application boundaries as needed, ensuring error messages include the invalid value and use the existing function names (MethylationTable::from_haplotype, apply_methylation_conversion) to locate changes.Source: Coding guidelines
src/commands/simulate.rs (1)
258-280:⚠️ Potential issue | 🟠 Major | 🏗️ Heavy liftValidate MT/MB coverage per simulated contig, not just once per VCF.
vcf_has_mt_mb_records(...)only proves that some annotated record exists somewhere in the file.run_simulation()still walks every reference contig, and the first contig without MT/MB falls into the "internal invariant" error path insimulate_contiginstead of failing cleanly up front. That means a subset-contig methylated VCF can still waste a long run and leave partial outputs behind. Please preflight the parsed methylation records against the contigs you intend to simulate, or turn the per-contigNonecase into a normal user-facing validation error before opening writers.🤖 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/simulate.rs` around lines 258 - 280, The current check using crate::vcf::methylation::vcf_has_mt_mb_records only ensures some MT/MB records exist anywhere in the VCF; update run_simulation() to preflight the parsed methylation records against each reference contig you will simulate (use the same contig list used by simulate_contig) and ensure every contig has MT/MB coverage when self.methylation_mode.is_some(); alternatively, modify simulate_contig to return a user-facing validation error (instead of an internal panic) when a contig has no methylation records (i.e., the per-contig None case), and surface that error before opening writers so the run fails fast with a clear message.
🤖 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.
Outside diff comments:
In `@src/commands/simulate.rs`:
- Around line 258-280: The current check using
crate::vcf::methylation::vcf_has_mt_mb_records only ensures some MT/MB records
exist anywhere in the VCF; update run_simulation() to preflight the parsed
methylation records against each reference contig you will simulate (use the
same contig list used by simulate_contig) and ensure every contig has MT/MB
coverage when self.methylation_mode.is_some(); alternatively, modify
simulate_contig to return a user-facing validation error (instead of an internal
panic) when a contig has no methylation records (i.e., the per-contig None
case), and surface that error before opening writers so the run fails fast with
a clear message.
In `@src/meth.rs`:
- Around line 92-109: Replace the panic-based validation in
MethylationTable::from_haplotype (and the same pattern in
apply_methylation_conversion) with a typed error return: remove the assert! and
validate methylation_rate, returning Result<Self, ErrorType> where ErrorType is
a library error defined with thiserror (e.g., MethylationError::InvalidRate {
value: f64 }), or alternatively accept a validated newtype; update the function
signatures to return Result and propagate that Result through callers and tests
using anyhow at application boundaries as needed, ensuring error messages
include the invalid value and use the existing function names
(MethylationTable::from_haplotype, apply_methylation_conversion) to locate
changes.
ℹ️ Review info
⚙️ Run configuration
Configuration used: Organization UI
Review profile: CHILL
Plan: Pro
Run ID: 72444353-4ae4-46fd-bacb-bf44b84d0ecf
📒 Files selected for processing (8)
README.mdsrc/commands/simulate.rssrc/meth.rssrc/methylation_tags.rssrc/output/golden_bam.rssrc/read.rstests/test_methylate.rstests/test_simulate_meth.rs
✅ Files skipped from review due to trivial changes (1)
- README.md
🚧 Files skipped from review as they are similar to previous changes (4)
- src/output/golden_bam.rs
- src/methylation_tags.rs
- src/read.rs
- tests/test_methylate.rs
Cleanup pass from a detailed branch review (opus/sonnet/haiku + docs/test agents). No behavior change to the simulator outputs beyond the bedGraph rate format fix. - bedGraph: round the population-fraction rate to an integer percentage so methylate --bedgraph and simulate --cpg-truth-bedgraph share the MethylDackel `extract` format (a 1/3 site emitted 33.333... before). - Dedupe the two byte-identical reference-CpG scanners into a single pub(crate) meth::find_reference_cpgs. - Hoist the methylation default rates (0.999 / 0.01) to named constants so the seed string can't drift from the effective config defaults. - Document the bitvec dependency rationale (8x memory at chromosome scale) at MethylationTable and in Cargo.toml. Docs: - CLAUDE.md: add the methylate subcommand and the five new methylation modules to the architecture table. - README: document --min-fragment-length, --adapter-r1/r2, --golden-vcf (not yet implemented), and methylate's default METHYLATE sample name. Tests (+6): - conversion-rate out-of-range and NaN rejection - bedGraph non-integer rate rounds to 33 - apply_methylation_conversion no-op at n_genomic == 0 - triploid MT/MB round-trip (ploidy > 2 per-haplotype path) - TAPS+VCF integration: methylated SNP-created CpG converts to T
Contributor
Author
|
@coderabbitai review |
✅ Action performedReview finished.
|
* docs(simulate): document golden BAM header and per-record tags * docs(simulate): document cf:i tag in golden BAM tag count --------- Co-authored-by: Tim Fennell <tfenne@tfenne.com>
* feat(methylate): context-aware Markov methylation model Replace the independent per-CpG Bernoulli draw with a biologically realistic generator. This is a methylate-only change: the binary per-haplotype MT/MB representation is unchanged, so simulate, the golden BAM, chemistry, and the bedGraph writers are untouched. Model: - Each CpG is classified de novo (Gardiner-Garden criteria) into CpG-island interior / shore / open-sea. Islands are hypomethylated, open-sea hypermethylated, shores intermediate. - A two-state (methylated/unmethylated) Markov chain walks each haplotype's CpG list with a keep-or-redraw transition: stationary mean equals the context's target rate, and autocorrelation decays with genomic distance per a per-context correlation length. This yields realistic spatial structure (runs, island/shore/sea level differences) instead of salt-and-pepper noise. - Methylation is symmetric across strands by default; a low sporadic hemimethylation rate drops one strand per CpG. (The old model's independent strand draws overstated hemimethylation at 2p(1-p).) - Per-haplotype independent chains preserve variant-aware CpG handling and give allele-specific methylation. CLI: --methylation-rate is replaced by per-context rate flags (--methylation-rate-island/-shore/-open-sea), per-context correlation lengths, and --hemimethylation-rate. The no-flags default is now biologically realistic rather than fully methylated; set the three rate flags equal for uniform methylation. Adds detector + walk unit tests (stationary mean, autocorrelation, hemi rate, island hypomethylation, determinism) and methylate integration tests; migrates existing call sites and the test helper. * test(methylate): close coverage gaps from adversarial review Multi-agent verification of the context-aware model found no correctness or math defects; these address the Minor coverage/robustness findings: - Add a shore-rate test: island + shore + open-sea with distinct rates, asserting the shore mean falls between island and open-sea. Guards the Shore branch of params_for, which every other test excluded. - Add an allele-specific-methylation test: two haplotypes with identical CpG content must diverge, proving the per-haplotype chains are independent (rate-1.0/0.0 tests can't show this). - Strengthen methylate_default_is_no_longer_fully_methylated: assert the mean per-CpG rate sits near the open-sea target (structural) instead of "at least one CpG < 100" (probabilistic); note it is seed-pinned. - Doc precision: note the inclusive shore boundary; note the realized per-strand fraction m*(1-hemi_rate/2) when hemi_rate > 0. Fix two stale test comments referencing the removed --methylation-rate flag. * refactor(meth): demote new methylation-model surface to pub(crate) holodeck_lib is not consumed outside the crate, so the new model types need no public API surface. Demote MethylationModel, ContextParams, their fields, the DEFAULT_* constants, validate(), and the from_haplotype/ from_haplotypes constructors to pub(crate), matching CpgContext. No behavior change. * docs(meth): sync module overview with the Markov methylation model The top-of-file //! overview still described the old independent per-strand Bernoulli draws (and hemimethylation falling out of them), which contradicted the updated item-level docs. Rewrite it to describe the context-aware symmetric Markov model: CpgContext classification, the two-state chain over ContextParams, symmetric methylation with sporadic hemi_rate, and per-haplotype independence. (CodeRabbit, PR #14.)
Contributor
Author
|
@coderabbitai review |
✅ Action performedReview finished.
|
There was a problem hiding this comment.
Actionable comments posted: 1
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/commands/methylate.rs (1)
162-189:⚠️ Potential issue | 🟠 Major | ⚡ Quick winReject output paths that alias the inputs or each other.
File::create(&self.output)and the optional bedGraphFile::create(...)run before the later VCF parse and before contig loading finishes. If either output path points at--vcfor--reference, this command can clobber its own input; if--output == --bedgraph, both writers target the same file and corrupt each other. Fail fast on any path collision before opening either writer.🤖 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/methylate.rs` around lines 162 - 189, Detect and reject any conflicting paths before opening writers: check that self.output does not equal the VCF input path (self.vcf.vcf) or the reference path (reference field on the command struct), and if self.bedgraph.is_some() also ensure self.bedgraph.as_ref() != Some(&self.output) and != self.vcf.vcf and != reference; if any equality is found return an Err with a clear message instead of calling File::create. Make this check immediately before the current File::create(&self.output) and bedGraph creation block (the code using File::create(&self.output), the optional self.bedgraph creation and before calling crate::vcf::parse_variants_by_contig/write_bedgraph_header) so you fail fast on aliased input/output paths.
🧹 Nitpick comments (1)
src/commands/methylate.rs (1)
240-247: ⚡ Quick winUse the no-variant bedGraph fast path here.
build_haplotypes(..., sample_ploidy, ...)makeshaplotypesnon-empty even on variant-free contigs, so this call never reacheswrite_bedgraph_records'shaplotypes.is_empty()branch and rematerializes every reference-identical haplotype again. Passing an empty slice whenvariants.is_empty()preserves the output but avoids the extra per-haplotype extraction work on the common no-variant path.♻️ Proposed fix
if let Some(bg) = bedgraph_writer.as_mut() { + let bedgraph_haplotypes = if variants.is_empty() { + &[][..] + } else { + haplotypes.as_slice() + }; crate::output::methylation_bedgraph::write_bedgraph_records( bg, contig_name, &reference, - &haplotypes, + bedgraph_haplotypes, &methylation, )?; }🤖 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/methylate.rs` around lines 240 - 247, The call into write_bedgraph_records should use the no-variant fast path: when the local `variants.is_empty()` is true, pass an empty haplotype slice instead of `&haplotypes` to avoid rematerializing reference-identical haplotypes; change the argument to something like `if variants.is_empty() { &[] } else { &haplotypes }` in the block where you call crate::output::methylation_bedgraph::write_bedgraph_records (the surrounding code that previously used `build_haplotypes(..., sample_ploidy, ...)` can remain unchanged).
🤖 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 `@tests/test_methylate.rs`:
- Around line 224-234: The helper parse_bedgraph currently swallows malformed
rows by using filter_map; change it to fail tests on bad input by replacing the
filter_map with a mapping that validates shape and parse success (e.g., split
the line into cols, assert the expected column count and use expect/unwrap with
a clear message when parsing cols[1] or cols[3] fails) so parse_bedgraph
explicitly panics on malformed bedGraph lines instead of skipping them; update
the function named parse_bedgraph to validate cols.len() (and/or check indices
exist) and to unwrap/expect parse results with informative error text.
---
Outside diff comments:
In `@src/commands/methylate.rs`:
- Around line 162-189: Detect and reject any conflicting paths before opening
writers: check that self.output does not equal the VCF input path (self.vcf.vcf)
or the reference path (reference field on the command struct), and if
self.bedgraph.is_some() also ensure self.bedgraph.as_ref() != Some(&self.output)
and != self.vcf.vcf and != reference; if any equality is found return an Err
with a clear message instead of calling File::create. Make this check
immediately before the current File::create(&self.output) and bedGraph creation
block (the code using File::create(&self.output), the optional self.bedgraph
creation and before calling
crate::vcf::parse_variants_by_contig/write_bedgraph_header) so you fail fast on
aliased input/output paths.
---
Nitpick comments:
In `@src/commands/methylate.rs`:
- Around line 240-247: The call into write_bedgraph_records should use the
no-variant fast path: when the local `variants.is_empty()` is true, pass an
empty haplotype slice instead of `&haplotypes` to avoid rematerializing
reference-identical haplotypes; change the argument to something like `if
variants.is_empty() { &[] } else { &haplotypes }` in the block where you call
crate::output::methylation_bedgraph::write_bedgraph_records (the surrounding
code that previously used `build_haplotypes(..., sample_ploidy, ...)` can remain
unchanged).
🪄 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: b98e37d2-ddd9-4823-ab8c-f80f28cdff3f
📒 Files selected for processing (5)
README.mdsrc/commands/methylate.rssrc/meth.rstests/helpers/mod.rstests/test_methylate.rs
🚧 Files skipped from review as they are similar to previous changes (2)
- tests/helpers/mod.rs
- README.md
Comment on lines
+224
to
+234
| /// Parse a population-fraction bedGraph into `(start, rate)` pairs, skipping | ||
| /// the track header line. | ||
| fn parse_bedgraph(contents: &str) -> Vec<(usize, u32)> { | ||
| contents | ||
| .lines() | ||
| .filter(|l| !l.starts_with("track")) | ||
| .filter_map(|l| { | ||
| let cols: Vec<&str> = l.split('\t').collect(); | ||
| Some((cols.get(1)?.parse().ok()?, cols.get(3)?.parse().ok()?)) | ||
| }) | ||
| .collect() |
There was a problem hiding this comment.
Make malformed bedGraph rows fail the test instead of disappearing.
filter_map turns parse/shape errors into skipped lines, so a broken writer can still satisfy the downstream mean/range assertions as long as some rows remain parseable. This helper should assert the expected column count and parse success so format regressions fail loudly.
💚 Proposed fix
fn parse_bedgraph(contents: &str) -> Vec<(usize, u32)> {
contents
.lines()
.filter(|l| !l.starts_with("track"))
- .filter_map(|l| {
- let cols: Vec<&str> = l.split('\t').collect();
- Some((cols.get(1)?.parse().ok()?, cols.get(3)?.parse().ok()?))
+ .map(|l| {
+ let cols: Vec<&str> = l.split('\t').collect();
+ assert_eq!(cols.len(), 6, "unexpected bedGraph row shape: {l}");
+ let start = cols[1].parse().expect("bedGraph start must be an integer");
+ let rate = cols[3].parse().expect("bedGraph rate must be an integer");
+ (start, rate)
})
.collect()
}📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
Suggested change
| /// Parse a population-fraction bedGraph into `(start, rate)` pairs, skipping | |
| /// the track header line. | |
| fn parse_bedgraph(contents: &str) -> Vec<(usize, u32)> { | |
| contents | |
| .lines() | |
| .filter(|l| !l.starts_with("track")) | |
| .filter_map(|l| { | |
| let cols: Vec<&str> = l.split('\t').collect(); | |
| Some((cols.get(1)?.parse().ok()?, cols.get(3)?.parse().ok()?)) | |
| }) | |
| .collect() | |
| /// Parse a population-fraction bedGraph into `(start, rate)` pairs, skipping | |
| /// the track header line. | |
| fn parse_bedgraph(contents: &str) -> Vec<(usize, u32)> { | |
| contents | |
| .lines() | |
| .filter(|l| !l.starts_with("track")) | |
| .map(|l| { | |
| let cols: Vec<&str> = l.split('\t').collect(); | |
| assert_eq!(cols.len(), 6, "unexpected bedGraph row shape: {l}"); | |
| let start = cols[1].parse().expect("bedGraph start must be an integer"); | |
| let rate = cols[3].parse().expect("bedGraph rate must be an integer"); | |
| (start, rate) | |
| }) | |
| .collect() | |
| } |
🤖 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 `@tests/test_methylate.rs` around lines 224 - 234, The helper parse_bedgraph
currently swallows malformed rows by using filter_map; change it to fail tests
on bad input by replacing the filter_map with a mapping that validates shape and
parse success (e.g., split the line into cols, assert the expected column count
and use expect/unwrap with a clear message when parsing cols[1] or cols[3]
fails) so parse_bedgraph explicitly panics on malformed bedGraph lines instead
of skipping them; update the function named parse_bedgraph to validate
cols.len() (and/or check indices exist) and to unwrap/expect parse results with
informative error text.
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
Splits methylation simulation into two composable commands and moves the methylation truth into the input VCF:
holodeck methylate(new): annotates a VCF (or the bare reference) with per-haplotype, per-strand CpG methylation truth in newMT/MBFORMAT fields. Optional--bedgraphwrites a closed-form, MethylDackel-format population-fraction bedGraph.holodeck simulate(refactored): readsMT/MBfrom the input VCF and applies em-seq or TAPS chemistry.--methylation-rateis removed — methylation state is now an inspectable artifact rather than a transient simulation parameter. The full Bismark-compatible golden-BAM tag set (XG/XR/XM/YM/NM/MD/YS) is preserved.The split separates biology (which CpGs are methylated) from chemistry (em-seq vs TAPS conversion at some efficiency), so a genome can be methylated once and re-sequenced under different chemistries.
Architecture
Four commits, organized by layer:
feat(vcf): add methylation classifier and MT/MB VCF I/O— per-haplotype CpG ownership classifier (Standalone vs OnVariant with upstream-wins for adjacent variants), writer + reader forMT/MB-bearing VCFs, all 13 documented edge cases plus closed-loop fuzz round-trip.feat(methylate): add subcommand for producing methylation-annotated VCFs—holodeck methylateend-to-end with Bernoulli draws, optional population-fraction bedGraph.refactor(simulate): load methylation truth from VCF MT/MB; drop --methylation-rate— simulate now sources truth exclusively from the input VCF'sMT/MBfields. Adds a validation matrix governing the(VCF has MT/MB, --methylation-mode set)cells.docs(readme): document the methylate + simulate two-command flow— README rewrite covering the schema, validation matrix, bedGraph comparison, and a complete worked example.Validation matrix (
simulate)MT/MB?--methylation-modeset?run \holodeck methylate` first`MT/MB schema
CpGs straddling a variant boundary are owned by the variant. When two adjacent variants jointly form a CpG, the upstream variant wins (deterministic tiebreaker). Variants must be phased; overlapping
REFspans on the same haplotype are rejected upfront withClassifyError.Test plan
358 tests passing on
cargo ci-fmt && cargo ci-lint && cargo ci-test(baselinemain= 248; PR #11's original commit added the chemistry tests; this refactor adds the classifier + writer + reader + matrix + bedGraph tests and migrates the existing chemistry tests to drive methylation truth viamethylate).ReadErrorvariant tests,header_has_mt_mbprobe tests, plus two closed-loop fuzz round-trips (10 kb random reference, then with phased SNPs).methylate: skeleton smoke test,MT/MBVCF output, end-to-end pipeline (methylate→simulate→ golden BAMYM:Ztag verification), bedGraph,--methylation-raterange validation,--seeddeterminism,--vcfvariants-input path.simulate: four validation-matrix cells,--cpg-truth-bedgraphintegration, byte-identity across same-seed runs in the no-methylation cell, plus all existing chemistry / Bismark-tag tests migrated to use themethylate_to_vcftest helper.Test plan checklist
cargo ci-test)cargo ci-fmt && cargo ci-lintclean (clippy pedantic)methylate→simulate --methylation-mode em-seq --golden-bamproduces a valid Bismark-compatible BAM withXM/YMreflecting the input methylation rateSummary by CodeRabbit
New Features
CLI
Documentation
Tests
Chore