feat(simulate): model terminal soft-clip artifacts

CHANGELOG.md

@@ -7,6 +7,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+- Optional terminal soft-clip artifact model for `simulate`. 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, for
+  example, shows ~8% of reads clipped, strongly 5'-biased and modal at 6–10 bp.
+  Four new flags control it: `--clip-5p-rate` and `--clip-3p-rate` (both
+  default `0.0`, i.e. disabled) give independent per-end probabilities, and
+  `--clip-length-mean` (default `8`) / `--clip-length-max` (default `20`) shape
+  a shared truncated-geometric clip-length distribution. Clipped bases are
+  replaced with random sequence 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.
+
 ## [0.3.0] - 2026-06-13
 
 ### Added

CLAUDE.md

@@ -40,6 +40,7 @@ Holodeck is a single-crate Rust project with a binary (`holodeck`) and library (
 | `vcf/methylation.rs` | MT/MB FORMAT field read/write; CpG classification across haplotypes |
 | `output/cpg_truth.rs` | Coverage-weighted per-CpG truth bedGraph (from simulated reads) |
 | `output/methylation_bedgraph.rs` | Closed-form population-fraction methylation bedGraph |
+| `clip.rs` | Terminal soft-clip artifact model (per-end 5'/3' clip sampling) |
 | `fragment.rs` | Fragment extraction, reverse complement, adapter padding |
 | `read.rs` | Read pair generation combining fragments + error model + naming |
 | `error_model/mod.rs` | ErrorModel trait + apply_errors free function |

README.md

@@ -125,6 +125,10 @@ holodeck simulate -r ref.fa -v methylated.vcf.gz -o output_taps \
 | `--min-error-rate` | 0.001 | Error rate at start of reads |
 | `--max-error-rate` | 0.01 | Error rate at end of reads |
 | `--max-n-frac` | 0.02 | Reject reads with >this fraction of bases from ambiguous reference positions (see [Ambiguous reference bases](#ambiguous-reference-bases)) |
+| `--clip-5p-rate` | 0.0 | Probability a read's 5' end carries a terminal soft-clip artifact (see [Terminal soft-clip artifacts](#terminal-soft-clip-artifacts)) |
+| `--clip-3p-rate` | 0.0 | Probability a read's 3' end carries a terminal soft-clip artifact, independent of 3' adapter read-through |
+| `--clip-length-mean` | 8 | Mean injected clip length (truncated geometric); only used when a clip rate is non-zero |
+| `--clip-length-max` | 20 | Maximum injected clip length; only used when a clip rate is non-zero |
 | `--methylation-mode` | none | Methylation chemistry: `em-seq` (or `bisulfite`) or `taps`. Presence enables methylation simulation |
 | `--methylation-conversion-rate` | 0.999 | Chemistry efficiency for normally-converting molecules (probability that the converting class of C converts to T) |
 | `--methylation-failure-rate` | 0.01 | Fraction of molecules that are whole-molecule conversion failures (convert at `1 − conversion-rate`). Requires `--methylation-mode` |
@@ -155,6 +159,14 @@ Real references contain a mix of `A`/`C`/`G`/`T`, large stretches of `N` (assemb
 
 **Known limitation:** requested `--coverage` is computed from raw contig/BED lengths, not from the non-ambiguous territory. For a reference like hs38DH (~5% N), rejection is noise and coverage lands where you'd expect. For simulations targeted at heavily-N contigs (or with `--max-n-frac 0.0` in N-dense regions), effective coverage will be slightly below the requested value; a warning is logged if the resampling budget is exhausted.
 
+### Terminal soft-clip artifacts
+
+By default holodeck's reads align essentially end-to-end: aside from VCF-derived indels and 3' adapter read-through, the simulated CIGAR is one long match block. Real libraries are messier — end-repair fill-in, damaged or non-templated read ends, and similar effects make aligners soft-clip a few bases off one or both ends of a meaningful fraction of reads. Twist EM-seq, for example, shows ~8% of reads carrying a soft clip, strongly 5'-biased and modal at 6–10 bp. The mechanism is protocol-agnostic, so the model is independent of `--methylation-mode` and useful for both methylation and non-methylation simulation.
+
+`--clip-5p-rate` and `--clip-3p-rate` turn this on (both default to `0.0`, i.e. disabled). For each read, the 5' and 3' ends independently receive a clip with their configured probabilities; a clipped end's length is drawn from a truncated geometric distribution with mean `--clip-length-mean`, capped at `--clip-length-max`. The 5' and 3' rates are separate because real libraries are markedly asymmetric (the 5' end clips several times more often), while the length distribution is shared across both ends.
+
+When an end is clipped, holodeck replaces those terminal bases with random sequence — so a downstream aligner re-derives the soft-clip rather than forcing a match — and records the soft-clip in the golden-BAM CIGAR (and advances the alignment start when a forward read's 5' end is clipped). Ground truth therefore stays exact: the aligned block still maps to the correct reference span. A disabled (or omitted) model draws no randomness, so existing seeded runs remain byte-identical.
+
 ### Golden BAM
 
 Pass `--golden-bam` to write a perfect-truth BAM alongside the FASTQ output. Every record carries the alignment the simulator generated it from — contig, position, strand, CIGAR — at MAPQ 60, in unsorted (generation) order. Sort downstream if a coordinate-sorted BAM is needed.

src/clip.rs

@@ -0,0 +1,224 @@
+//! Terminal soft-clip artifact model.
+//!
+//! Real sequencing libraries frequently produce reads whose extreme 5' or 3'
+//! ends diverge from the reference enough that an aligner soft-clips them
+//! (end-repair fill-in, damaged or otherwise non-templated read ends). Twist
+//! EM-seq, for example, shows ~8% of reads carrying a soft clip, strongly
+//! 5'-biased and modal at 6-10 bp.
+//!
+//! Holodeck's simulated reads otherwise align essentially end-to-end, so this
+//! module injects that artifact: for a configurable fraction of reads it
+//! marks the first/last `K` bases of the *aligned* portion as soft-clipped.
+//! The caller corrupts those bases with random sequence (so a downstream
+//! aligner re-derives the clip) and records the soft-clip in the golden BAM
+//! CIGAR, keeping ground truth exact.
+//!
+//! The mechanism is protocol-agnostic — it is independent of methylation
+//! chemistry and useful for any library type.
+
+use rand::Rng;
+
+/// Configuration for the terminal soft-clip artifact model.
+///
+/// The 5' and 3' ends clip independently, each with its own probability,
+/// because real libraries are markedly asymmetric (the 5' end clips several
+/// times more often than the 3' end). The clip-length distribution is shared
+/// across both ends, as observed lengths do not differ meaningfully between
+/// them.
+#[derive(Debug, Clone, Copy)]
+pub struct TerminalClipConfig {
+    /// Probability in `[0.0, 1.0]` that a read's 5' end receives a clip.
+    pub rate_5p: f64,
+    /// Probability in `[0.0, 1.0]` that a read's 3' end receives a clip.
+    pub rate_3p: f64,
+    /// Mean clip length in bases. Lengths are drawn from a truncated
+    /// geometric distribution with this mean.
+    pub length_mean: usize,
+    /// Maximum clip length in bases; sampled lengths are clamped to this.
+    pub length_max: usize,
+}
+
+impl TerminalClipConfig {
+    /// Whether either end can produce a clip. When this is `false`,
+    /// [`Self::sample_clips`] consumes no randomness and returns `(0, 0)`, so
+    /// the simulator's output is byte-identical to having no clip model at
+    /// all.
+    #[must_use]
+    pub fn is_enabled(&self) -> bool {
+        self.rate_5p > 0.0 || self.rate_3p > 0.0
+    }
+
+    /// Sample `(clip_5p, clip_3p)` lengths for one read whose aligned
+    /// (genomic) portion is `aligned_len` bases long.
+    ///
+    /// Each end independently clips with its configured probability; a chosen
+    /// clip's length is drawn from a truncated geometric with mean
+    /// `length_mean`, clamped to `[1, length_max]`. The combined clip is then
+    /// capped so at least one aligned base survives (an alignment needs a
+    /// match block), favouring the 5' end when both are drawn and the budget
+    /// is tight.
+    ///
+    /// Returns `(0, 0)` without drawing any randomness when the model is
+    /// disabled or `aligned_len < 2`.
+    pub fn sample_clips(&self, aligned_len: usize, rng: &mut impl Rng) -> (usize, usize) {
+        // No clip is representable without leaving an aligned base behind, and
+        // a disabled model must not perturb the RNG stream.
+        if !self.is_enabled() || aligned_len < 2 {
+            return (0, 0);
+        }
+
+        let mut clip_5p = self.sample_end(self.rate_5p, rng);
+        let mut clip_3p = self.sample_end(self.rate_3p, rng);
+
+        // Cap so at least one aligned base remains. Trim the 3' end first
+        // (the rarer, less characteristic artifact), then the 5' end.
+        let budget = aligned_len - 1;
+        if clip_5p + clip_3p > budget {
+            clip_3p = clip_3p.min(budget.saturating_sub(clip_5p));
+            clip_5p = clip_5p.min(budget);
+        }
+
+        (clip_5p, clip_3p)
+    }
+
+    /// Draw a clip length for a single end: `0` if the end does not clip,
+    /// otherwise a truncated-geometric length in `[1, length_max]`. Consumes
+    /// no randomness when `rate <= 0.0`.
+    fn sample_end(&self, rate: f64, rng: &mut impl Rng) -> usize {
+        if rate <= 0.0 {
+            return 0;
+        }
+        if rng.random::<f64>() >= rate {
+            return 0;
+        }
+        self.sample_length(rng)
+    }
+
+    /// Sample a clip length from a truncated geometric distribution with mean
+    /// [`Self::length_mean`], clamped to `[1, length_max]`.
+    fn sample_length(&self, rng: &mut impl Rng) -> usize {
+        let max = self.length_max.max(1);
+        let mean = self.length_mean.max(1);
+        if mean <= 1 {
+            // Degenerate distribution: always the shortest representable clip.
+            return 1;
+        }
+        #[expect(clippy::cast_precision_loss, reason = "small clip means")]
+        let p = 1.0 / mean as f64;
+        let u: f64 = rng.random(); // [0, 1)
+        // Inverse-CDF of the geometric (trials >= 1): k = floor(ln(1-u)/ln(1-p)) + 1.
+        let k = ((1.0 - u).ln() / (1.0 - p).ln()).floor();
+        #[expect(
+            clippy::cast_possible_truncation,
+            clippy::cast_sign_loss,
+            reason = "k is finite, non-negative after floor"
+        )]
+        let k = k as usize + 1;
+        k.clamp(1, max)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use rand::SeedableRng;
+    use rand::rngs::SmallRng;
+
+    fn cfg(rate_5p: f64, rate_3p: f64) -> TerminalClipConfig {
+        TerminalClipConfig { rate_5p, rate_3p, length_mean: 8, length_max: 20 }
+    }
+
+    #[test]
+    fn disabled_config_is_not_enabled() {
+        assert!(!cfg(0.0, 0.0).is_enabled());
+        assert!(cfg(0.1, 0.0).is_enabled());
+        assert!(cfg(0.0, 0.1).is_enabled());
+    }
+
+    #[test]
+    fn disabled_config_returns_zero_and_consumes_no_randomness() {
+        // A disabled model must leave the RNG stream untouched so existing
+        // simulations stay byte-identical.
+        let mut rng_a = SmallRng::seed_from_u64(7);
+        let clips = cfg(0.0, 0.0).sample_clips(100, &mut rng_a);
+        assert_eq!(clips, (0, 0));
+        let after: u64 = rng_a.random();
+
+        let mut rng_b = SmallRng::seed_from_u64(7);
+        let direct: u64 = rng_b.random();
+        assert_eq!(after, direct, "disabled clip model must not advance the RNG");
+    }
+
+    #[test]
+    fn short_alignment_cannot_clip() {
+        let mut rng = SmallRng::seed_from_u64(1);
+        assert_eq!(cfg(1.0, 1.0).sample_clips(0, &mut rng), (0, 0));
+        assert_eq!(cfg(1.0, 1.0).sample_clips(1, &mut rng), (0, 0));
+    }
+
+    #[test]
+    fn certain_5p_clip_always_clips_5p_only() {
+        let mut rng = SmallRng::seed_from_u64(42);
+        for _ in 0..200 {
+            let (c5, c3) = cfg(1.0, 0.0).sample_clips(150, &mut rng);
+            assert!(c5 >= 1, "rate_5p=1.0 must always clip the 5' end");
+            assert_eq!(c3, 0, "rate_3p=0.0 must never clip the 3' end");
+            assert!(c5 <= 20, "clip must respect length_max");
+        }
+    }
+
+    #[test]
+    fn certain_3p_clip_always_clips_3p_only() {
+        let mut rng = SmallRng::seed_from_u64(42);
+        for _ in 0..200 {
+            let (c5, c3) = cfg(0.0, 1.0).sample_clips(150, &mut rng);
+            assert_eq!(c5, 0);
+            assert!((1..=20).contains(&c3));
+        }
+    }
+
+    #[test]
+    fn clip_never_consumes_whole_alignment() {
+        let mut rng = SmallRng::seed_from_u64(99);
+        // length_max larger than the alignment forces the cap to engage.
+        let config =
+            TerminalClipConfig { rate_5p: 1.0, rate_3p: 1.0, length_mean: 8, length_max: 50 };
+        for aligned in 2..30 {
+            for _ in 0..50 {
+                let (c5, c3) = config.sample_clips(aligned, &mut rng);
+                assert!(c5 + c3 <= aligned - 1, "must leave >=1 aligned base (aligned={aligned})");
+            }
+        }
+    }
+
+    #[test]
+    fn observed_rate_tracks_configured_rate() {
+        let mut rng = SmallRng::seed_from_u64(2024);
+        let mut clipped = 0;
+        let n = 20_000;
+        for _ in 0..n {
+            let (c5, _) = cfg(0.3, 0.0).sample_clips(150, &mut rng);
+            if c5 > 0 {
+                clipped += 1;
+            }
+        }
+        let observed = f64::from(clipped) / f64::from(n);
+        assert!((observed - 0.3).abs() < 0.02, "expected ~0.30 clip rate, got {observed}");
+    }
+
+    #[test]
+    fn mean_length_is_in_the_right_ballpark() {
+        let mut rng = SmallRng::seed_from_u64(2025);
+        let config =
+            TerminalClipConfig { rate_5p: 1.0, rate_3p: 0.0, length_mean: 8, length_max: 1000 };
+        let mut sum = 0usize;
+        let n = 50_000;
+        for _ in 0..n {
+            let (c5, _) = config.sample_clips(2000, &mut rng);
+            sum += c5;
+        }
+        #[expect(clippy::cast_precision_loss, reason = "test arithmetic")]
+        let mean = sum as f64 / f64::from(n);
+        assert!((mean - 8.0).abs() < 1.0, "expected mean ~8, got {mean}");
+    }
+}

src/commands/simulate.rs

@@ -123,6 +123,34 @@ pub struct Simulate {
     #[arg(long, default_value_t = 0.02, value_name = "FLOAT")]
     pub max_n_frac: f64,
 
+    /// Probability in `[0.0, 1.0]` that a read's 5' end carries a terminal
+    /// soft-clip artifact (end-repair fill-in, damaged or non-templated read
+    /// ends). The clipped bases are replaced with random sequence so a
+    /// downstream aligner soft-clips them, and the golden BAM records the clip
+    /// in the CIGAR so ground truth stays exact. Real Twist EM-seq shows ~8%
+    /// of reads 5'-clipped. Defaults to `0.0` (disabled). Protocol-agnostic:
+    /// independent of `--methylation-mode`.
+    #[arg(long, default_value_t = 0.0, value_name = "FLOAT")]
+    pub clip_5p_rate: f64,
+
+    /// Probability in `[0.0, 1.0]` that a read's 3' end carries a terminal
+    /// soft-clip artifact, independent of (and in addition to) 3' adapter
+    /// read-through. Real Twist EM-seq shows ~2% of reads 3'-clipped. Defaults
+    /// to `0.0` (disabled).
+    #[arg(long, default_value_t = 0.0, value_name = "FLOAT")]
+    pub clip_3p_rate: f64,
+
+    /// Mean length in bases of an injected terminal soft-clip. Lengths are
+    /// drawn from a truncated geometric distribution centered here and capped
+    /// by `--clip-length-max`. Only takes effect when a clip rate is non-zero.
+    #[arg(long, default_value_t = 8, value_name = "INT")]
+    pub clip_length_mean: usize,
+
+    /// Maximum length in bases of an injected terminal soft-clip. Only takes
+    /// effect when a clip rate is non-zero.
+    #[arg(long, default_value_t = 20, value_name = "INT")]
+    pub clip_length_max: usize,
+
     /// Enable methylation chemistry simulation. `em-seq` (or `bisulfite`)
     /// converts unmethylated cytosines to thymine and preserves methylated
     /// ones (matches both bisulfite and em-seq protocols). `taps` is the
@@ -223,6 +251,22 @@ impl Simulate {
         if !self.max_n_frac.is_finite() || !(0.0..=1.0).contains(&self.max_n_frac) {
             bail!("--max-n-frac must be in [0.0, 1.0]");
         }
+        if !self.clip_5p_rate.is_finite() || !(0.0..=1.0).contains(&self.clip_5p_rate) {
+            bail!("--clip-5p-rate must be in [0.0, 1.0]");
+        }
+        if !self.clip_3p_rate.is_finite() || !(0.0..=1.0).contains(&self.clip_3p_rate) {
+            bail!("--clip-3p-rate must be in [0.0, 1.0]");
+        }
+        // Validate clip-length parameters only when the model is enabled, so a
+        // disabled model (both rates zero) never rejects on unused values.
+        if self.clip_5p_rate > 0.0 || self.clip_3p_rate > 0.0 {
+            if self.clip_length_max == 0 {
+                bail!("--clip-length-max must be > 0 when a terminal clip rate is set");
+            }
+            if self.clip_length_mean > self.clip_length_max {
+                bail!("--clip-length-mean must be <= --clip-length-max");
+            }
+        }
         if let Some(rate) = self.methylation_conversion_rate
             && (!rate.is_finite() || !(0.0..=1.0).contains(&rate))
         {
@@ -349,6 +393,18 @@ impl Simulate {
 
         let error_model =
             IlluminaErrorModel::new(self.read_length, self.min_error_rate, self.max_error_rate);
+
+        // Terminal soft-clip artifact model. Passed as `Some` only when at
+        // least one end can clip; a disabled model is equivalent to `None`
+        // (it draws no randomness) but keeping it `None` makes intent clear.
+        let clip_config = (self.clip_5p_rate > 0.0 || self.clip_3p_rate > 0.0).then_some(
+            crate::clip::TerminalClipConfig {
+                rate_5p: self.clip_5p_rate,
+                rate_3p: self.clip_3p_rate,
+                length_mean: self.clip_length_mean,
+                length_max: self.clip_length_max,
+            },
+        );
         let frag_dist = Normal::new(self.fragment_mean as f64, self.fragment_stddev as f64)
             .map_err(|e| anyhow::anyhow!("Invalid fragment distribution parameters: {e}"))?;
 
@@ -457,6 +513,7 @@ impl Simulate {
                 targets.as_ref(),
                 total_reads,
                 &error_model,
+                clip_config.as_ref(),
                 &frag_dist,
                 adapter_r1.as_bytes(),
                 adapter_r2.as_bytes(),
@@ -599,6 +656,7 @@ impl Simulate {
         targets: Option<&TargetRegions>,
         total_reads: u64,
         error_model: &IlluminaErrorModel,
+        clip_config: Option<&crate::clip::TerminalClipConfig>,
         frag_dist: &Normal<f64>,
         adapter_r1: &[u8],
         adapter_r2: &[u8],
@@ -803,6 +861,7 @@ impl Simulate {
                 // YS:Z tag — skip the per-mate clone when no golden BAM is
                 // requested.
                 self.golden_bam,
+                clip_config,
                 rng,
             ) else {
                 // Too many ambiguity-resolved bases in this read pair; resample.
@@ -926,6 +985,10 @@ mod tests {
             min_error_rate: 0.001,
             max_error_rate: 0.01,
             max_n_frac: 0.02,
+            clip_5p_rate: 0.0,
+            clip_3p_rate: 0.0,
+            clip_length_mean: 8,
+            clip_length_max: 20,
             methylation_mode: None,
             methylation_conversion_rate: None,
             methylation_failure_rate: None,