feat(simulate): model terminal soft-clip artifacts#19
Draft
nh13 wants to merge 1 commit into
Draft
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.
nh13
force-pushed
the
nh/terminal-clip-model
branch
from
June 15, 2026 19:45
1cc8f47 to
8ce562a
Compare
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
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.
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.