clip: CellRanger4 3' poly-A trim

CHANGELOG.md

@@ -13,6 +13,12 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Features
 
+- **CellRanger4 3' poly-A trim** (`clip::cellranger4::poly_tail_3p`) and the
+  `--clipAdapterType CellRanger4` clip plumbing, plus `--clip5pAdapterSeq` and
+  `--clip5pAdapterMMp`. The 5' TSO trim is an overlap alignment and waits on a
+  dedicated deterministic-SIMD crate, so `CellRanger4` is still rejected at
+  parse time rather than half-applied.
+
 - **STARsolo single-cell quantification (`--soloType`)** — the 10x
   Chromium / plate-based count-matrix pipeline, ported from STAR and
   verified against real STARsolo (#90).

src/clip/cellranger4.rs

@@ -0,0 +1,114 @@
+//! `--clipAdapterType CellRanger4`: the 10x Chromium v4 clipping rules.
+//!
+//! Two independent trims, ported from STAR's `ClipCR4.cpp` and
+//! `ClipMate_clipChunk.cpp`:
+//!
+//! - a 3' poly-A tail trim, scored base by base from the 3' end;
+//! - a 5' template-switch-oligo trim, which is an overlap alignment of the TSO
+//!   against the first 91 bases of the read.
+//!
+//! Only the poly-A trim is implemented here. The 5' TSO trim is an overlap
+//! alignment, for which STAR links the Opal SIMD library; rustar will take
+//! that from a dedicated deterministic-SIMD crate rather than carrying a
+//! second aligner in-tree. Until then `--clipAdapterType CellRanger4` is
+//! rejected rather than silently doing half the job.
+
+/// Number of 3' bases to trim as a CellRanger4 poly-A tail.
+///
+/// STAR `ClipCR4::polyTail3p`. Walks in from the 3' end scoring `+1` per `A`
+/// and `-2` per non-`A`, and remembers the longest prefix of that walk whose
+/// running score still clears a 70% density threshold (`score * 10 >= ib * 7`).
+/// It gives up once the score has fallen more than 27 behind the position, and
+/// returns nothing unless the remembered score reached 20.
+///
+/// `seq` is numeric base codes, so `A == 0`.
+pub fn poly_tail_3p(seq: &[u8]) -> usize {
+    let seq_len = seq.len();
+    if seq_len < 20 {
+        return 0;
+    }
+    let mut best_len: i64 = seq_len as i64 - 1;
+    let mut score: i64 = 0;
+    let mut best_score: i64 = 0;
+    for ib in 1..=seq_len as i64 {
+        if seq[seq_len - ib as usize] == 0 {
+            score += 1;
+            if score * 10 >= ib * 7 {
+                best_len = ib;
+                best_score = score;
+            }
+        } else {
+            score -= 2;
+            if ib - score > 27 {
+                break;
+            }
+        }
+    }
+    if best_score < 20 {
+        0
+    } else {
+        best_len as usize
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// ACGT text to base codes.
+    fn code(s: &str) -> Vec<u8> {
+        s.bytes()
+            .map(|b| match b {
+                b'A' => 0,
+                b'C' => 1,
+                b'G' => 2,
+                b'T' => 3,
+                _ => 4,
+            })
+            .collect()
+    }
+
+    #[test]
+    fn cr4_polya_trim_matches_star() {
+        // A clean 30-base poly-A tail is trimmed whole. The prefix is
+        // deliberately A-free: the scan does not stop at the tail boundary, so
+        // an `A` just upstream of it would legitimately extend the trim.
+        let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{}", "A".repeat(30)));
+        assert_eq!(poly_tail_3p(&read), 30);
+
+        // No tail: nothing to trim.
+        let read = code("ACGTACGTACGTACGTACGTACGTACGTACGT");
+        assert_eq!(poly_tail_3p(&read), 0);
+
+        // A tail shorter than the score-20 floor is not trimmed, however clean.
+        let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{}", "A".repeat(10)));
+        assert_eq!(poly_tail_3p(&read), 0);
+    }
+
+    #[test]
+    fn poly_tail_needs_twenty_bases_of_read() {
+        assert_eq!(poly_tail_3p(&code("AAAAAAAAAAAAAAAAAAA")), 0); // 19
+    }
+
+    #[test]
+    fn poly_tail_scan_does_not_stop_at_the_tail_boundary() {
+        // STAR keeps scoring past the run of A's, so A-rich sequence just
+        // upstream extends the trim. Worth pinning: it looks like an off-by-one
+        // otherwise.
+        let read = code(&format!("ACGTACGTACGTACGTACGT{}", "A".repeat(30)));
+        assert_eq!(poly_tail_3p(&read), 34);
+    }
+
+    #[test]
+    fn poly_tail_tolerates_a_single_interruption() {
+        // 70% density is the threshold, so one non-A inside a long tail is
+        // survivable.
+        let tail = format!("{}C{}", "A".repeat(15), "A".repeat(15));
+        let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{tail}"));
+        assert!(
+            poly_tail_3p(&read) >= 15,
+            "one mismatch should not abandon the tail, got {}",
+            poly_tail_3p(&read)
+        );
+    }
+}

src/clip/mod.rs

@@ -17,9 +17,12 @@
 //! `clip5pAfterAdapterNbases`), then 3' on the 5'-clipped read (`clip3pNbases`,
 //! then the 3' adapter Hamming scan, then `clip3pAfterAdapterNbases`).
 //!
-//! Only `--clipAdapterType Hamming` (STAR's default) is supported; a 5' Hamming
-//! adapter is not a thing STAR itself supports either (only `CellRanger4` mode
-//! clips a 5' adapter, the 10x TSO) — that mode is out of scope here.
+//! `--clipAdapterType Hamming` (STAR's default) uses the 3' Hamming scan above.
+//! `CellRanger4` instead trims a 3' poly-A tail and a 5' TSO; see
+//! [`cellranger4`]. A 5' *Hamming* adapter is not a thing STAR supports either:
+//! the 5' end only ever carries an adapter under CellRanger4.
+
+pub mod cellranger4;
 
 use crate::io::fastq::encode_base;
 use crate::params::Parameters;
@@ -48,6 +51,12 @@ pub struct ClipParams {
     pub five: ClipEnd,
     /// 3' end (STAR `ClipMate` type 1).
     pub three: ClipEnd,
+    /// `--clipAdapterType CellRanger4`: replace the Hamming rules with the 10x
+    /// poly-A / TSO trims. See [`cellranger4`].
+    pub cellranger4: bool,
+    /// `--clip5pAdapterSeq` as base codes: the 10x TSO, clipped from the 5' end
+    /// under CellRanger4. Empty when none is configured.
+    pub five_adapter: Vec<u8>,
 }
 
 /// Build [`ClipParams`] for `mate` (0 or 1) from the run's `--clip{5,3}pNbases`
@@ -56,6 +65,15 @@ pub struct ClipParams {
 /// per-mate (`clip5p(mate)`/`clip3p(mate)`); the adapter / mmp / after-adapter
 /// clips are single-valued and apply to both mates. Cheap to build per batch.
 pub fn clip_params_from(params: &Parameters, mate: usize) -> ClipParams {
+    // The 5' TSO trim of CellRanger4 needs an overlap alignment, which waits on
+    // a dedicated deterministic-SIMD crate. Say so once, loudly, rather than
+    // leaving the user to infer from the output that half the mode ran.
+    if params.clip_adapter_type == "CellRanger4" && params.clip5p_adapter_seq != "-" && mate == 0 {
+        log::warn!(
+            "--clipAdapterType CellRanger4: the 3' poly-A trim is applied, but the 5' \
+             adapter (TSO) trim is not yet implemented and --clip5pAdapterSeq is ignored"
+        );
+    }
     // Encode the adapter to base codes (A=0..T=3) ONCE here. The read reaching
     // clip_mate is already numeric (io::fastq encodes at read time), so the 3'
     // Hamming scan compares numeric-vs-numeric — re-encoding the read there turned
@@ -65,7 +83,14 @@ pub fn clip_params_from(params: &Parameters, mate: usize) -> ClipParams {
     } else {
         params.clip3p_adapter_seq.bytes().map(encode_base).collect()
     };
+    let five_adapter = if params.clip5p_adapter_seq == "-" {
+        Vec::new()
+    } else {
+        params.clip5p_adapter_seq.bytes().map(encode_base).collect()
+    };
     ClipParams {
+        cellranger4: params.clip_adapter_type == "CellRanger4",
+        five_adapter,
         five: ClipEnd {
             n: params.clip5p(mate),
             adapter: Vec::new(),
@@ -120,6 +145,10 @@ fn local_search(x: &[u8], y: &[u8], p_mm: f64) -> usize {
 pub fn clip_mate(read: &[u8], p: &ClipParams) -> (usize, usize) {
     let len = read.len();
 
+    if p.cellranger4 {
+        return clip_mate_cellranger4(read, p);
+    }
+
     // ---- 5' end (STAR ClipMate type 0) ----
     let five_active = p.five.n > 0;
     let mut c5 = 0;
@@ -155,6 +184,84 @@ pub fn clip_mate(read: &[u8], p: &ClipParams) -> (usize, usize) {
     (c5, c3)
 }
 
+/// `--clipAdapterType CellRanger4`, which replaces the Hamming rules entirely
+/// (STAR `ClipCR4`): a 5' TSO trim and a 3' poly-A trim.
+///
+/// The fixed `--clip{5,3}pNbases` still apply first, as in the Hamming path.
+/// The 5' TSO is only configured for the first mate, matching STAR, so mate 2
+/// simply has no 5' adapter and the 5' trim reduces to the fixed clip.
+fn clip_mate_cellranger4(read: &[u8], p: &ClipParams) -> (usize, usize) {
+    let len = read.len();
+
+    // 5': fixed clip only. The TSO trim is an overlap alignment and waits on a
+    // dedicated deterministic-SIMD crate. A configured TSO is therefore not
+    // applied, and the caller is warned once rather than left to infer it from
+    // the output.
+    let mut c5 = p.five.n.min(len);
+    if p.five.n_after > 0 && c5 < len {
+        c5 += p.five.n_after.min(len - c5);
+    }
+
+    // 3': fixed clip, then the poly-A scan on the 5'-clipped read.
+    let s = &read[c5..];
+    let sl = s.len();
+    let mut c3 = p.three.n.min(sl);
+    let remaining = sl - c3;
+    if remaining > 0 {
+        c3 += cellranger4::poly_tail_3p(&s[..remaining]).min(remaining);
+    }
+    if p.three.n_after > 0 && c3 < sl {
+        c3 += p.three.n_after.min(sl - c3);
+    }
+
+    (c5, c3)
+}
+
+#[cfg(test)]
+mod cr4_wiring_tests {
+    use super::*;
+
+    fn code(s: &str) -> Vec<u8> {
+        s.bytes().map(encode_base).collect()
+    }
+
+    const TSO: &str = "AAGCAGTGGTATCAACGCAGAGTACATGGG";
+
+    fn cr4_params(tso: &str) -> ClipParams {
+        ClipParams {
+            cellranger4: true,
+            five_adapter: code(tso),
+            five: ClipEnd::default(),
+            three: ClipEnd::default(),
+        }
+    }
+
+    #[test]
+    fn cellranger4_leaves_a_read_without_either_feature_alone() {
+        let read = code("CGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGT");
+        assert_eq!(clip_mate(&read, &cr4_params(TSO)), (0, 0));
+    }
+
+    #[test]
+    fn cellranger4_without_a_tso_still_trims_polya() {
+        // Mate 2 has no 5' adapter under STAR, so only the 3' trim applies.
+        let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{}", "A".repeat(30)));
+        let (c5, c3) = clip_mate(&read, &cr4_params("-"));
+        assert_eq!(c5, 0);
+        assert_eq!(c3, 30);
+    }
+
+    #[test]
+    fn cellranger4_applies_the_fixed_clip_before_the_polya_trim() {
+        let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{}", "A".repeat(30)));
+        let mut p = cr4_params("-");
+        p.five.n = 5;
+        let (c5, c3) = clip_mate(&read, &p);
+        assert_eq!(c5, 5, "the fixed 5' clip still applies");
+        assert_eq!(c3, 30, "and the poly-A tail is trimmed from what remains");
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -183,6 +290,8 @@ mod tests {
     #[test]
     fn fixed_5p_3p() {
         let p = ClipParams {
+            cellranger4: false,
+            five_adapter: Vec::new(),
             five: ClipEnd {
                 n: 3,
                 ..Default::default()
@@ -217,6 +326,8 @@ mod tests {
     fn after_adapter_alone_is_noop() {
         // STAR's inactive-end short-circuit: n_after with no fixed clip and no adapter clips nothing.
         let p = ClipParams {
+            cellranger4: false,
+            five_adapter: Vec::new(),
             five: ClipEnd {
                 n_after: 4,
                 ..Default::default()
@@ -240,6 +351,8 @@ mod tests {
     #[test]
     fn no_adapter_configured_only_fixed_clips() {
         let p = ClipParams {
+            cellranger4: false,
+            five_adapter: Vec::new(),
             five: ClipEnd {
                 n: 2,
                 ..Default::default()

src/params/mod.rs

@@ -555,6 +555,16 @@ pub struct Parameters {
     #[arg(long = "clipAdapterType", default_value = "Hamming")]
     pub clip_adapter_type: String,
 
+    /// 5' adapter sequence to clip, one per mate. Only used by
+    /// `--clipAdapterType CellRanger4`, where it is the 10x template switch
+    /// oligo. `-` (the default) means none.
+    #[arg(long = "clip5pAdapterSeq", default_value = "-")]
+    pub clip5p_adapter_seq: String,
+
+    /// Max mismatch fraction for the 5' adapter.
+    #[arg(long = "clip5pAdapterMMp", default_value_t = 0.1)]
+    pub clip5p_adapter_mmp: f64,
+
     /// 3' adapter sequence to clip (Hamming scan), `-` = none
     #[arg(long = "clip3pAdapterSeq", default_value = "-")]
     pub clip3p_adapter_seq: String,