Repository navigation
feat(simulate): add --methylation-mode to simulate subcommands - #263
Merged
Merged
Conversation
nh13
marked this pull request as ready for review
April 13, 2026 07:55
Member
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #263 +/- ##
==========================================
+ Coverage 89.24% 89.85% +0.60%
==========================================
Files 120 120
Lines 59642 61359 +1717
==========================================
+ Hits 53228 55133 +1905
+ Misses 6414 6226 -188 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
No actionable comments were generated in the recent review. 🎉 ℹ️ Recent review info⚙️ Run configurationConfiguration used: Organization UI Review profile: CHILL Plan: Pro Run ID: 📒 Files selected for processing (7)
✅ Files skipped from review due to trivial changes (2)
🚧 Files skipped from review as they are similar to previous changes (4)
📝 WalkthroughWalkthroughThe PR makes methylation simulation available across simulate subcommands and changes the crate feature 🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✏️ Tip: You can configure your own custom pre-merge checks in the settings. ✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 5
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/lib/commands/simulate/consensus_reads.rs (1)
359-385:⚠️ Potential issue | 🟠 MajorR2 methylation tags are emitted in the wrong orientation.
generate_consensus_pair()passes a reverse-complemented R2 sequence, butbuild_consensus_record()reuses the unreversed annotation and always callsbuild_mm_ml_tags(..., true, ...)/build_mm_tag_no_ml(..., true, ...). For R2 that makes MM/ML/am/bm point at the wrong bases or disappear entirely. Reverse/reorient the annotation for R2, or only emit methylation tags on R1.Also applies to: 579-598
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@src/lib/commands/simulate/consensus_reads.rs` around lines 359 - 385, generate_consensus_pair() currently reverse-complements the R2 sequence (r2_seq) but passes the original methylation annotation into build_consensus_record(), causing build_mm_ml_tags / build_mm_tag_no_ml to emit MM/ML/am/bm for the wrong bases; to fix, produce a reoriented methylation annotation for R2 before constructing r2_record (or skip emitting methylation for R2), e.g., compute an R2_methylation by reversing the per-base coordinates and flipping strand/orientation flags consistent with reverse_complement, then pass R2_methylation.as_ref() (instead of methylation.as_ref()) into build_consensus_record() when creating r2_record; alternatively, update build_consensus_record() to detect is_first==false and internally call a reorient_methylation helper prior to calling build_mm_ml_tags / build_mm_tag_no_ml (ensure build_mm_ml_tags/build_mm_tag_no_ml receive correctly oriented coords).
🤖 Prompt for all review comments with AI agents
Verify each finding against the current code and only fix it if needed.
Inline comments:
In `@crates/fgumi-consensus/src/methylation.rs`:
- Around line 77-85: The public helper is_cpg_context can panic when pos >=
ref_seq.len(); add an early guard at the start of is_cpg_context (e.g., if pos
>= ref_seq.len() { return false; }) so both top- and bottom-strand branches
never index out of bounds; keep the existing pos+1 < ref_seq.len() and pos > 0
checks but ensure the initial length check prevents panics for invalid
coordinates.
In `@src/lib/commands/simulate/common.rs`:
- Around line 331-339: The sampling range incorrectly omits the last valid start
and treats exact-fit (length == self.total_length) as empty; compute max_start =
self.total_length.saturating_sub(length) and if length > self.total_length
return None, otherwise pick genome_pos from 0..=max_start (inclusive) so that
when max_start == 0 you still get start 0; update the call using the sampler
(e.g., change rng.random_range(0..self.total_length.saturating_sub(length)) in
sample_sequence to an inclusive range or equivalent logic that uses max_start
and handles exact-fit correctly.
In `@src/lib/commands/simulate/consensus_reads.rs`:
- Around line 112-115: The new flag methylation_depth_mean is not validated and
passing it into create_depth_distribution(methylation_depth_mean,
methylation_depth_mean / 2.0) can panic for zero, negative, or non-finite
values; add the same upfront validation used for the rate parameters to reject
values that are <= 0.0 or not finite before constructing the log-normal
(validate methylation_depth_mean in the argument parsing or immediately after
parsing), return an error or exit with a clear message if invalid, and apply the
same check for the other occurrence noted around the create_depth_distribution
call at the later block (lines referenced 188–191).
In `@src/lib/commands/simulate/grouped_reads.rs`:
- Around line 142-147: The grouped-read code loads family templates from
ReferenceGenome (ReferenceGenome::load) but still writes BAM headers and record
positions in the synthetic --ref-name/--ref-length coordinate space; update the
write path that builds the BAM header and record positions so it maps each
template's locus/contig from the loaded ReferenceGenome into the BAM ref-name
and ref-length used in the header, and set each Record's reference id and POS
from the template's actual locus; likewise, when the fallback sampler replaces a
template or sequence, ensure it also updates the Record.sequence, Record.pos and
reference id to match the sampled locus (this affects the places that
construct/write BAM headers and records around the grouped-MI output and the
fallback-sampler branches referenced in the diff).
In `@src/lib/commands/simulate/mapped_reads.rs`:
- Around line 134-139: The BAM records are being assigned a synthetic
contig/position that can mismatch the actual FASTA locus used for read bases;
when loading a reference via ReferenceGenome::load (ref_genome /
self.reference.reference) ensure the sampled locus (contig/chromosome name and
start coordinate) is preserved and written into the mapped record metadata
instead of always using the synthetic --ref-name/--ref-length. Concretely: when
you sample sequence from the ReferenceGenome, capture the actual contig name and
reference start/length and set the BAM record's reference ID/name and position
to those values (update the code paths mentioned around 145-164 and 305-310 that
construct mapped records) and only fall back to --ref-name/--ref-length when no
reference was provided; keep ReferenceGenome::load usage and propagate the
sampled locus through to the record creation so coordinates and bases remain
consistent.
---
Outside diff comments:
In `@src/lib/commands/simulate/consensus_reads.rs`:
- Around line 359-385: generate_consensus_pair() currently reverse-complements
the R2 sequence (r2_seq) but passes the original methylation annotation into
build_consensus_record(), causing build_mm_ml_tags / build_mm_tag_no_ml to emit
MM/ML/am/bm for the wrong bases; to fix, produce a reoriented methylation
annotation for R2 before constructing r2_record (or skip emitting methylation
for R2), e.g., compute an R2_methylation by reversing the per-base coordinates
and flipping strand/orientation flags consistent with reverse_complement, then
pass R2_methylation.as_ref() (instead of methylation.as_ref()) into
build_consensus_record() when creating r2_record; alternatively, update
build_consensus_record() to detect is_first==false and internally call a
reorient_methylation helper prior to calling build_mm_ml_tags /
build_mm_tag_no_ml (ensure build_mm_ml_tags/build_mm_tag_no_ml receive correctly
oriented coords).
🪄 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: 7f6e3d5a-88f5-486e-95f7-15632b46d387
📒 Files selected for processing (7)
Cargo.tomlcrates/fgumi-consensus/src/methylation.rssrc/lib/commands/simulate/common.rssrc/lib/commands/simulate/consensus_reads.rssrc/lib/commands/simulate/fastq_reads.rssrc/lib/commands/simulate/grouped_reads.rssrc/lib/commands/simulate/mapped_reads.rs
Add methylation-aware simulation to fastq-reads, mapped-reads, grouped-reads, and consensus-reads subcommands. For base-level commands (fastq-reads, mapped-reads, grouped-reads): - New --methylation-mode <em-seq|taps> enables C→T / G→A conversion - --cpg-methylation-rate and --conversion-rate control stochastic behavior - Requires --reference for realistic template sequences - Per-read independent conversion produces realistic cu/ct count ratios - Strand-aware: top strand converts C→T, bottom strand converts G→A For consensus-reads: - Generates MM/ML/cu/ct tags directly (no base conversion) - Duplex mode adds per-strand au/at/bu/bt/am/bm tags - --methylation-depth-mean controls count sampling depth Shared infrastructure: - MethylationConfig struct bundles mode + rates to reduce parameter sprawl - MethylationArgs CLI args with validation (rate bounds checking) - is_cpg_context moved to fgumi-consensus::methylation for reuse - ReferenceGenome gains sequence_at_genome_pos for multi-chromosome support - simulate feature now depends on simplex for methylation types
Remove the synthetic reference mode from all simulate subcommands. ReferenceArgs now has a single required --reference PathBuf field. fastq-reads, mapped-reads, grouped-reads, and consensus-reads all require --reference; ref_name and ref_length are derived from the loaded ReferenceGenome instead of CLI flags.
nh13
force-pushed
the
nh/simulate-methylation
branch
from
April 13, 2026 21:16
7b17b46 to
0b88de9
Compare
There was a problem hiding this comment.
Actionable comments posted: 10
♻️ Duplicate comments (2)
src/lib/commands/simulate/grouped_reads.rs (1)
349-353:⚠️ Potential issue | 🟠 MajorKeep fallback templates and reported loci in sync.
If the original
(chrom_idx, local_pos)cannot provideinsert_sizebases, this falls back to a different sampled locus but throws away the new coordinates. The BAM records and truth TSV then advertise a position that does not match the emitted sequence.🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@src/lib/commands/simulate/grouped_reads.rs` around lines 349 - 353, The code currently uses ref_genome.sequence_at(chrom_idx, local_pos, insert_size).or_else(|| ref_genome.sample_sequence(insert_size).map(|(_, _, seq)| seq)) to get shared_template but discards the coordinates when falling back; update this so that when sample_sequence is used you also capture and propagate the returned (chrom_idx, local_pos) (or a new vars like sampled_chrom, sampled_pos) and replace the original chrom_idx/local_pos used for BAM record and truth TSV emission. In practice change the selection to return both coordinates and sequence (e.g., Option<(chrom_idx, local_pos, seq)>) and ensure all downstream uses (shared_template, BAM write, truth TSV creation) read coordinates from that single source so reported loci match the emitted template sequence.src/lib/commands/simulate/mapped_reads.rs (1)
324-328:⚠️ Potential issue | 🟠 MajorKeep fallback templates and coordinates together.
If
sequence_at()misses, the fallback samples a different locus but only keeps the bases. Records and truth still use the originalchrom_idx/local_pos, so coordinates can disagree with the emitted sequence.🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@src/lib/commands/simulate/mapped_reads.rs` around lines 324 - 328, The fallback currently discards the sampled locus and only keeps the bases, causing coordinates (chrom_idx/local_pos) to mismatch the emitted sequence; change the logic around ref_genome.sequence_at and ref_genome.sample_sequence so the fallback preserves both coordinates and sequence (e.g., return a tuple with (chrom_idx, local_pos, sequence) or a small struct) instead of mapping to just seq, and update uses of shared_template and any record/truth emission to read coordinates from that returned value so sequence and coordinates always come from the same source (refer to sequence_at, sample_sequence, shared_template, chrom_idx, local_pos, insert_size).
🧹 Nitpick comments (1)
tests/integration/test_e2e_regression.rs (1)
439-445: Assert on methylation tags here.Comparing
simplex.bamto itself only checks that the BAM is readable. Missing or wrongMM/ML/cu/cttags — and any TAPS-specific regression — still pass.🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@tests/integration/test_e2e_regression.rs` around lines 439 - 445, The test currently just calls assert_bams_identical(&simplex, &simplex, ...) which only checks BAM readability; replace or augment that with an assertion that the simplex BAM actually contains methylation/TAPS tags: open the simplex BAM and iterate records (or a representative sample) and assert each record (or sample) has non-empty MM and ML tags and the TAPS-specific cu and ct tags present and well-formed; either add a helper like assert_bam_contains_methylation_tags(simplex) or extend the existing assert_bams_identical call to validate tags 'MM', 'ML', 'cu', and 'ct' are present and valid (e.g., ML is a byte array of expected length, MM is a string, cu/ct exist) rather than comparing the file to itself.
🤖 Prompt for all review comments with AI agents
Verify each finding against the current code and only fix it if needed.
Inline comments:
In `@docs/simulate-cli.md`:
- Around line 41-42: The CLI reference is missing documentation for the new
methylation flags; add entries for --methylation-mode, --cpg-methylation-rate,
and --conversion-rate wherever flag tables are defined (e.g., near the
`--reference` entry and at the other referenced sections) by describing accepted
values and defaults: --methylation-mode (choices like "none"/"global"/"file" or
the actual modes implemented), --cpg-methylation-rate (float between 0–1,
default value), and --conversion-rate (float between 0–1, default value); ensure
the flags are listed in the main simulate CLI reference tables and include short
descriptions and default values to match how other flags (e.g., `--reference`)
are documented so the feature is discoverable.
- Line 17: The docs still describe the old truth-file schema (1-based `position`
and `A/B` strand) but the code in src/lib/commands/simulate/mapped_reads.rs and
src/lib/commands/simulate/grouped_reads.rs emits raw 0-based positions and
'+'/'-' strand values; update docs/simulate-cli.md to change the schema
description and example tables to state positions are 0-based (raw coordinates)
and strands use '+' or '-' (not A/B), and replace any example rows or column
descriptions that show 1-based positions or A/B with representative 0-based
numbers and +/- values; also apply the same change to the other occurrence noted
at line 148 so both places consistently reflect the new schema.
In `@src/lib/commands/simulate/consensus_reads.rs`:
- Around line 337-341: The code currently calls
params.ref_genome.sample_sequence(...).expect(...) inside
generate_consensus_pair (or the rayon worker), which will panic on bad user
input; change this to propagate a proper error instead of panicking: make
generate_consensus_pair return a Result (or otherwise fallible) and replace the
expect with match/? handling of sample_sequence's Option/Result so you return a
CLI-friendly error; alternatively validate params.ref_genome (e.g., ensure at
least one contig ≥ read_length) once up front before spawning rayon workers and
error out with a clear message if sampling is impossible (references:
generate_consensus_pair, params.ref_genome.sample_sequence, read_length).
- Around line 208-219: GenerationParams never receives the CLI/resolved
conversion_rate and the simulation still uses a hard-coded total_depth / 5
window; add a conversion_rate field to GenerationParams, thread the resolved
conversion_rate into GenerationParams where params are constructed, and replace
all uses of the fixed "total_depth / 5" split (notably the logic in the simulate
consensus code around the current windows at the locations referenced, e.g., the
blocks using total_depth/5 and the split logic at 375-384 and 511-533) to use
params.conversion_rate when determining the converted/unconverted read counts so
the --conversion-rate flag actually controls CU/CT and MM/ML output.
- Around line 389-427: The methylation and duplex-methylation annotations must
be reversed based on each record's orientation, not read number; update the
calls to build_consensus_record so that the methylation arg for R1 is
methylation.as_ref().map(|m| if r1_is_reverse { m.reverse() } else { m.clone()
}) and for R2 uses the record's reverse flag (e.g., r2_is_reverse =
!r1_is_reverse) to decide whether to reverse, and similarly if params.duplex is
true pass the duplex tuple reversed for the record that is reverse-complemented
instead of always passing the same tuple for both reads; adjust usage of
r1_is_reverse, r2_is_reverse, r1_methylation and r2_methylation when
constructing the two build_consensus_record calls.
- Around line 599-607: The MM/ML and am/bm tag builders currently hard-code
is_top_strand = true, which ignores the is_reverse flag and causes reverse
records to use the wrong base track; update each call that sets is_top_strand
(the MM/ML and am/bm tag helper invocations in the block handling
methylation/methylation_mode around where methylation: Option<&MethylationData>,
methylation_mode: MethylationMode and is_reverse: bool are in scope) to derive
the strand from is_reverse (e.g., is_top_strand = !is_reverse or equivalent) so
the helpers switch between C (top strand) and G (bottom strand) correctly for
reverse-complemented records; make this change consistently for all occurrences
in the 642–671 range.
In `@src/lib/commands/simulate/fastq_reads.rs`:
- Around line 389-392: sample_sequence() can return None (oversized inserts or
exhausted retries), so replace the expect on sample_sequence with proper None
handling: check the Option returned by reference.sample_sequence(template_len,
&mut rng) and return a user-facing error (Err) instead of panicking; include
contextual info (template_len and that sampling failed on the given reference)
in the error message and propagate it from the current function so Rayon workers
don't panic. Ensure you reference the same variables (chrom_idx, pos, template)
after successful sampling and keep usage of reference.name(chrom_idx) unchanged.
- Around line 408-410: The code always chooses strand/family with
rng.random_bool regardless of the CLI flag, so ensure the --duplex flag governs
this: change the assignment for (strand, family_idx) to check the duplex flag
(e.g., if duplex { if rng.random_bool(0.5) { ("A",0) } else { ("B",1) } } else {
("A",0) }) so non-duplex runs always produce strand "A" and family_idx 0 while
duplex runs keep the 50/50 random behavior; update the expression where (strand,
family_idx) is set (currently using rng.random_bool) to this conditional.
In `@src/lib/commands/simulate/grouped_reads.rs`:
- Around line 176-198: Reject zero positions before building the position table:
if position_dist.num_positions (computed as num_positions) is 0 while
self.common.num_molecules > 0, return an error (anyhow::bail!) explaining that
--num-positions must be >0 (or reduce molecule count) instead of continuing to
call ref_genome.sample_positions and later doing mol_id % num_positions or
indexing position_table[pos_idx]; add this check immediately before creating
pos_rng / calling ref_genome.sample_positions to fail fast.
In `@src/lib/commands/simulate/mapped_reads.rs`:
- Around line 169-187: Check for num_positions == 0 right after computing let
num_positions =
self.position_dist.num_positions.unwrap_or(self.common.num_molecules); and
return an error if zero (e.g., using anyhow::bail!) before any use of
num_positions; update the block around usable_bases/pos_rng/sample_positions to
bail with a clear message when num_positions == 0 so ref_genome.sample_positions
and later code that uses mol_id % num_positions never see zero. Use the existing
symbols num_positions, self.position_dist.num_positions,
self.common.num_molecules, ref_genome.sample_positions and create_rng to locate
where to insert this check.
---
Duplicate comments:
In `@src/lib/commands/simulate/grouped_reads.rs`:
- Around line 349-353: The code currently uses ref_genome.sequence_at(chrom_idx,
local_pos, insert_size).or_else(||
ref_genome.sample_sequence(insert_size).map(|(_, _, seq)| seq)) to get
shared_template but discards the coordinates when falling back; update this so
that when sample_sequence is used you also capture and propagate the returned
(chrom_idx, local_pos) (or a new vars like sampled_chrom, sampled_pos) and
replace the original chrom_idx/local_pos used for BAM record and truth TSV
emission. In practice change the selection to return both coordinates and
sequence (e.g., Option<(chrom_idx, local_pos, seq)>) and ensure all downstream
uses (shared_template, BAM write, truth TSV creation) read coordinates from that
single source so reported loci match the emitted template sequence.
In `@src/lib/commands/simulate/mapped_reads.rs`:
- Around line 324-328: The fallback currently discards the sampled locus and
only keeps the bases, causing coordinates (chrom_idx/local_pos) to mismatch the
emitted sequence; change the logic around ref_genome.sequence_at and
ref_genome.sample_sequence so the fallback preserves both coordinates and
sequence (e.g., return a tuple with (chrom_idx, local_pos, sequence) or a small
struct) instead of mapping to just seq, and update uses of shared_template and
any record/truth emission to read coordinates from that returned value so
sequence and coordinates always come from the same source (refer to sequence_at,
sample_sequence, shared_template, chrom_idx, local_pos, insert_size).
---
Nitpick comments:
In `@tests/integration/test_e2e_regression.rs`:
- Around line 439-445: The test currently just calls
assert_bams_identical(&simplex, &simplex, ...) which only checks BAM
readability; replace or augment that with an assertion that the simplex BAM
actually contains methylation/TAPS tags: open the simplex BAM and iterate
records (or a representative sample) and assert each record (or sample) has
non-empty MM and ML tags and the TAPS-specific cu and ct tags present and
well-formed; either add a helper like
assert_bam_contains_methylation_tags(simplex) or extend the existing
assert_bams_identical call to validate tags 'MM', 'ML', 'cu', and 'ct' are
present and valid (e.g., ML is a byte array of expected length, MM is a string,
cu/ct exist) rather than comparing the file to itself.
🪄 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: f2761f07-8ae4-4f52-bd1b-09b815b45c67
📒 Files selected for processing (9)
Cargo.tomlcrates/fgumi-consensus/src/methylation.rsdocs/simulate-cli.mdsrc/lib/commands/simulate/common.rssrc/lib/commands/simulate/consensus_reads.rssrc/lib/commands/simulate/fastq_reads.rssrc/lib/commands/simulate/grouped_reads.rssrc/lib/commands/simulate/mapped_reads.rstests/integration/test_e2e_regression.rs
🚧 Files skipped from review as they are similar to previous changes (1)
- Cargo.toml
nh13
force-pushed
the
nh/simulate-methylation
branch
from
April 13, 2026 21:47
0b88de9 to
16b491d
Compare
Member
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
There was a problem hiding this comment.
Actionable comments posted: 2
♻️ Duplicate comments (3)
src/lib/commands/simulate/fastq_reads.rs (1)
233-240:⚠️ Potential issue | 🟠 MajorReturn a user error here instead of panicking.
The preflight only checks
read_length, but sampling usesinsert_size - umi_length. Oversized inserts or exhausted retries can still hit thisexpect, which will panic a Rayon worker instead of surfacing a normal CLI failure.Also applies to: 402-405
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@src/lib/commands/simulate/fastq_reads.rs` around lines 233 - 240, The code currently panics via an expect when a contig is too short (e.g., the check using reference.max_contig_length() vs self.common.read_length and the sampling code that uses insert_size - umi_length), which can crash a Rayon worker; replace those expect/panic sites with returning a user-facing error (use anyhow::bail or return Err) that includes the offending values (requested read_length or effective insert length and longest contig) so the CLI surfaces a clean failure; update the checks around reference.max_contig_length(), the sampling logic that computes insert_size - umi_length, and the other duplicate site (the block at ~402-405) to return a descriptive error instead of panicking.src/lib/commands/simulate/grouped_reads.rs (1)
354-356:⚠️ Potential issue | 🟠 MajorDon't change the template without changing the locus.
When
sequence_atmisses, this can pull bases from another locus — or even generate a fresh random template per read pair — whilechrom_idx/local_posstay fixed. That makes the BAM/truth coordinates disagree with the reference context, and families can stop sharing one underlying template.Also applies to: 490-495
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@src/lib/commands/simulate/grouped_reads.rs` around lines 354 - 356, The code calls ref_genome.sequence_at(...) and falls back to ref_genome.sample_sequence(...) but keeps chrom_idx and local_pos unchanged, causing templates to come from a different locus while coordinates stay the same; update the logic in the generation of shared_template (and the same pattern around lines with sequence_at/sample_sequence) so that when sequence_at returns None you call sample_sequence and also set chrom_idx and local_pos (or an equivalent locus identifier) to the sampled locus before using the returned seq, ensuring the BAM/truth coordinates and shared_template are consistent and that sample_sequence is invoked once per family (not per read pair); reference symbols: sequence_at, sample_sequence, shared_template, chrom_idx, local_pos, insert_size.src/lib/commands/simulate/mapped_reads.rs (1)
329-331:⚠️ Potential issue | 🟠 MajorFallback sampling still breaks coordinate fidelity.
On a
sequence_atmiss, this can substitute bases from another locus, and if that also fails it invents a random template, butchrom_idx/local_posand truth output never change. The emitted alignments can therefore disagree with the reference bases implied by their coordinates.Also applies to: 341-346
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@src/lib/commands/simulate/mapped_reads.rs` around lines 329 - 331, The fallback sampling in shared_template (call to ref_genome.sequence_at(...) .or_else(|| ref_genome.sample_sequence(...))) can produce a sequence from a different locus (or random) without updating chrom_idx/local_pos or the truth output, breaking coordinate fidelity; change the logic so that when sequence_at returns None you either (A) do not fallback to sample_sequence and instead treat the fragment as unmappable/skip emitting it, or (B) if you must sample, capture the sampled locus returned by ref_genome.sample_sequence and update chrom_idx and local_pos and any truth/alignment fields accordingly so coordinates match the sequence; modify the code around shared_template, ref_genome.sequence_at, ref_genome.sample_sequence, chrom_idx, local_pos, and the truth output generation to implement one of these two behaviors consistently (also apply same fix to the similar block at lines 341-346).
🤖 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/lib/commands/simulate/mapped_reads.rs`:
- Around line 273-275: The loop writing pairs currently always emits r1 then r2
which can violate template-coordinate order when generate_molecule_reads
produces an is_top_strand == false pair (R1 is rightmost); change the write
logic so you compute each read's template/start coordinate (e.g.,
record.reference_start/pos) and write the record with the smaller coordinate
first: if r1.start <= r2.start write r1 then r2, otherwise write r2 then r1;
keep using writer.write_alignment_record and preserve read pairing/metadata
(read_name, umi_str, is_top_strand) but swap the write order based on the
coordinate comparison to guarantee template-coordinate-sorted output.
In `@tests/integration/test_e2e_regression.rs`:
- Around line 439-445: The test currently calls assert_bams_identical(&simplex,
&simplex, ...) which only verifies BAM readability; update the test to validate
methylation emission by either comparing &simplex to a known-good reference BAM
(not itself) or by adding explicit checks for methylation tags/values in the
simplex file: open simplex and assert presence and non-empty values for tags
like "cu"/"ct" or the MM/ML tagpair for at least one read (use the same helper
functions used elsewhere in tests or add a small helper that scans reads for
MM/ML and fails if missing or empty).
---
Duplicate comments:
In `@src/lib/commands/simulate/fastq_reads.rs`:
- Around line 233-240: The code currently panics via an expect when a contig is
too short (e.g., the check using reference.max_contig_length() vs
self.common.read_length and the sampling code that uses insert_size -
umi_length), which can crash a Rayon worker; replace those expect/panic sites
with returning a user-facing error (use anyhow::bail or return Err) that
includes the offending values (requested read_length or effective insert length
and longest contig) so the CLI surfaces a clean failure; update the checks
around reference.max_contig_length(), the sampling logic that computes
insert_size - umi_length, and the other duplicate site (the block at ~402-405)
to return a descriptive error instead of panicking.
In `@src/lib/commands/simulate/grouped_reads.rs`:
- Around line 354-356: The code calls ref_genome.sequence_at(...) and falls back
to ref_genome.sample_sequence(...) but keeps chrom_idx and local_pos unchanged,
causing templates to come from a different locus while coordinates stay the
same; update the logic in the generation of shared_template (and the same
pattern around lines with sequence_at/sample_sequence) so that when sequence_at
returns None you call sample_sequence and also set chrom_idx and local_pos (or
an equivalent locus identifier) to the sampled locus before using the returned
seq, ensuring the BAM/truth coordinates and shared_template are consistent and
that sample_sequence is invoked once per family (not per read pair); reference
symbols: sequence_at, sample_sequence, shared_template, chrom_idx, local_pos,
insert_size.
In `@src/lib/commands/simulate/mapped_reads.rs`:
- Around line 329-331: The fallback sampling in shared_template (call to
ref_genome.sequence_at(...) .or_else(|| ref_genome.sample_sequence(...))) can
produce a sequence from a different locus (or random) without updating
chrom_idx/local_pos or the truth output, breaking coordinate fidelity; change
the logic so that when sequence_at returns None you either (A) do not fallback
to sample_sequence and instead treat the fragment as unmappable/skip emitting
it, or (B) if you must sample, capture the sampled locus returned by
ref_genome.sample_sequence and update chrom_idx and local_pos and any
truth/alignment fields accordingly so coordinates match the sequence; modify the
code around shared_template, ref_genome.sequence_at, ref_genome.sample_sequence,
chrom_idx, local_pos, and the truth output generation to implement one of these
two behaviors consistently (also apply same fix to the similar block at lines
341-346).
🪄 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: a5a3e554-e813-4f6e-843a-a2d709207954
📒 Files selected for processing (7)
docs/simulate-cli.mdsrc/lib/commands/simulate/common.rssrc/lib/commands/simulate/consensus_reads.rssrc/lib/commands/simulate/fastq_reads.rssrc/lib/commands/simulate/grouped_reads.rssrc/lib/commands/simulate/mapped_reads.rstests/integration/test_e2e_regression.rs
✅ Files skipped from review due to trivial changes (1)
- docs/simulate-cli.md
🚧 Files skipped from review as they are similar to previous changes (2)
- src/lib/commands/simulate/consensus_reads.rs
- src/lib/commands/simulate/common.rs
…ReferenceGenome Add build_bam_header, num_chromosomes, chromosome_length, and sample_positions methods to ReferenceGenome. sample_positions pre-samples random loci genome-wide, skipping N-containing windows, for efficient position-table lookup in upcoming callers. Include tests for all new methods covering correctness, determinism, and multi-chromosome coverage.
Remove the compute_position function and its test. Replace call sites in
mapped_reads and grouped_reads with todo!("position table lookup")
stubs, pending the position-table approach that will use the newly added
sample_positions method. Remove the now-unused num_positions and
ref_length fields from GenerationParams in both callers.
…nd coin flip Replace the todo!() stubs in mapped-reads with real reference coordinate support. The BAM header now includes all reference contigs from the loaded FASTA. Positions are pre-sampled from the reference genome using a dedicated RNG, and molecule seeds use a separate RNG to avoid correlation. Each molecule gets a 50/50 strand coin flip producing both F1R2 and R1F2 orientations. The truth TSV now reports real chromosome names, positions, and strand (+/-).
…ce coordinates Load reference genome and build BAM header from contigs. Sample template sequences from the reference instead of generating random ones. Add 50/50 strand coin flip per molecule for F1R2/R1F2 orientation. Produce mapped records with reference ID, position, and mapping quality. Update truth TSV with chrom, pos, and strand columns.
Remove Option wrappers around reference genome since --reference is now required. Remove random sequence fallback — all templates are sampled from the reference. Add 50/50 strand coin flip for both duplex and non-duplex modes so molecules are randomly assigned F1R2 or R1F2 orientation. Simplify truth TSV writing (always has chrom/pos).
All simulate subcommands now require --reference. Update the e2e regression test helpers to create a test reference FASTA and pass it to simulate grouped-reads and fastq-reads invocations.
Mark --reference as required for all simulate subcommands (except correct-reads). Remove --ref-name and --ref-length documentation. Update consensus-reads from "unmapped" to "mapped" BAM. Update all examples to include --reference. Note strand coin flip behavior.
Verify that simulated reads with real reference coordinates can be consumed by downstream methylation-aware tools: grouped-reads with --methylation-mode em-seq -> simplex with --methylation-mode em-seq --ref. This exercises the full pipeline to confirm contigs match.
nh13
force-pushed
the
nh/simulate-methylation
branch
from
April 13, 2026 23:12
16b491d to
41cfa35
Compare
Member
Author
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
Merged
Merged
This branch was previously deployed
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
1 participant
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
--methylation-mode <em-seq|taps>tofastq-reads,mapped-reads,grouped-reads, andconsensus-readssimulate subcommands--cpg-methylation-rateand--conversion-rate, requiring--referencefor realistic templatesconsensus-readsgenerates MM/ML/cu/ct tags directly (plus per-strand au/at/bu/bt/am/bm for duplex)MethylationConfigstruct andMethylationArgsCLI args reduce parameter sprawl across subcommandsis_cpg_contextmoved tofgumi-consensus::methylationfor reuse across cratesReferenceGenomegainssequence_at_genome_posfor correct multi-chromosome template samplingsimulatefeature now depends onsimplexfor access to methylation typesCloses #262
Test plan
cargo ci-test)cargo ci-fmt)cargo ci-lint)simulate grouped-reads --methylation-mode em-seq→simplex --methylation-mode em-seq→filter --methylation-mode em-seq--methylation-mode tapssimulate consensus-reads --methylation-mode em-seq --duplexproduces au/at/bu/bt/am/bm tags