feat(simulate): model terminal soft-clip artifacts - #19
Conversation
|
Important Review skippedDraft detected. Please check the settings in the CodeRabbit UI or the ⚙️ Run configurationConfiguration used: Organization UI Review profile: CHILL Plan: Pro Run ID: You can disable this status message by setting the Use the checkbox below for a quick retry:
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 |
Add an optional terminal soft-clip model so simulated reads no longer align essentially end-to-end. Real libraries soft-clip a few bases off one or both read ends on a meaningful fraction of reads (end-repair fill-in, damaged or non-templated ends); Twist EM-seq shows ~8% of reads clipped, strongly 5'-biased and modal at 6-10 bp. Four new `simulate` flags control it: `--clip-5p-rate` and `--clip-3p-rate` (both default 0.0, disabled) give independent per-end probabilities; `--clip-length-mean` (8) and `--clip-length-max` (20) shape a shared truncated-geometric clip-length distribution. Separate 5'/3' rates capture the asymmetry real libraries show; the length distribution is shared because observed lengths do not differ meaningfully between ends. Clipped bases are replaced with random sequence so a downstream aligner re-derives the clip, and the soft-clip is recorded in the golden-BAM CIGAR (advancing the alignment start when a forward read's 5' end is clipped) so ground truth stays exact. The model is protocol-agnostic (independent of --methylation-mode); a disabled or omitted model draws no randomness, leaving existing seeded runs byte-identical. The CIGAR builder now takes explicit leading/trailing soft-clip lengths instead of (adapter_bases, strand), folding adapter read-through and terminal clips into the correct record ends per strand.
1cc8f47 to
8ce562a
Compare
|
@coderabbitai pause |
Rate Limit Exceeded
|
2 participants
Closes #18.
Adds an optional terminal soft-clip artifact model to
simulate. By default holodeck's reads align essentially end-to-end (aside from VCF indels and 3' adapter read-through), but real libraries soft-clip a few bases off one or both read ends on a meaningful fraction of reads. Measured against real Twist EM-seq (meth-twist-emseq-5M,bwa-mem3 mem --meth, ~2M primary reads): ~8.7% of reads carry a soft clip, strongly 5'-biased (8.2% vs 1.8%), modal at 6–10 bp, with near-random clipped-base composition. The mechanism (end-repair fill-in, damaged/non-templated ends) is protocol-agnostic, so this is useful for both methylation and non-methylation simulation.Behavior
Four new flags, all off by default:
--clip-5p-rate0.0--clip-3p-rate0.0--clip-length-mean8--clip-length-max20Separate 5'/3' rates capture the asymmetry real libraries show; the length distribution is shared because observed lengths don't differ meaningfully between ends. This resolves open question 1 in #18 (separate rates) and 2 (no per-mate R1/R2 split for now). Defaults remain placeholders pending a full-length untrimmed lane (open question 3) — they only take effect once a rate is set.
When an end clips, the terminal bases are replaced with random sequence so a downstream aligner re-derives the clip, and the soft-clip is recorded in the golden-BAM CIGAR (advancing the alignment start when a forward read's 5' end is clipped). Ground truth therefore stays exact. A disabled or omitted model draws no randomness, so existing seeded runs are byte-identical — verified by a test.
Implementation
src/clip.rs:TerminalClipConfigwith per-end Bernoulli draws and a truncated-geometric length sampler.cigar_from_ref_positionsnow takes explicit(lead_clip, trail_clip)lengths instead of(adapter_bases, negative_strand);build_matefolds adapter read-through and terminal clips into the correct record ends per strand and trims the reference-position slice so POS and CIGAR stay consistent.simulatewith validation.Testing
clip.rs(disabled-is-inert/no-RNG-draw, certain-clip, rate tracking, mean length, never-consumes-whole-alignment).read.rsintegration: 5'-forward clip advances POS + leading clip; 3' is trailing; disabled config is byte-identical toNone.test_simulate.rs: golden BAM shows the expected clip fraction, baseline (no flags) is clip-free, clips are terminal-only with clean M interiors, and lengths respect the cap.cargo fmt,cargo clippy --all-features --all-targets -- -D warnings, and the full test suite all pass.Smoke run (
--clip-5p-rate 0.5 --clip-3p-rate 0.1): 55.9% of records soft-clipped, matching the analytic 55%, with both leading and trailing clips of varied lengths.