From 0907ac04faf6c9a362ca723753b1843e9d0d1461 Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Tue, 28 Jul 2026 22:55:07 +0200
Subject: [PATCH 1/8] feat(swalign): in-tree affine-gap aligner, and the
 CellRanger4 clip built on it
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

STAR performs exactly one alignment of this shape: the 5' TSO clip of
`--clipAdapterType CellRanger4`, for which it links the Opal C/C++ SIMD
library. rustar-aligner accepted `CellRanger4` at the CLI but did nothing with
it, which is a correctness hole rather than only a missing feature: the user
asked for clipping and silently did not get it.

New module `src/swalign`, no new dependency:

- `Mode::{Nw, Hw, Ov, Sw}` and affine `Scoring`, with `N`-against-`N` scored
  neutrally as Opal does — STAR pads the target to 91 bases with `N` and relies
  on that padding being free.
- `scalar`, a portable column-by-column implementation. It is written for
  clarity, not speed, because it is the *definition* of the result: the SIMD
  backends to come must agree with it bit-for-bit.

Determinism is treated as a correctness property here, not a nicety. Opal's
overflow buckets are sized by the SIMD vector width, so when 8-bit lanes
saturate the grouping — and with it the recompute path — depends on which
instruction set is available; sixteen SSE lanes and thirty-two AVX2 lanes
bucket differently, and simde emulates AVX2 on ARM. The rule in this module is
that the vector width is never observable. The scalar path defines the answer
and everything else has to match it, including under saturation, on empty
inputs, on `N`, and on end-position ties.

`src/clip/cellranger4` then ports the two CellRanger4 trims:

- `poly_tail_3p`, STAR's `ClipCR4::polyTail3p` 3' poly-A scan;
- `tso_clip`, the 5' template-switch-oligo trim, which runs through `swalign`
  in overlap mode instead of Opal and then applies STAR's acceptance gate
  (reject below score 20; reject 20 and 21 when they took more than 26 and 30
  bases to reach).

`cr4_tso_clip_matches_opal` is the frozen vector shared with STAR-rs's test of
the same name, which validates those numbers against Opal itself.

One behaviour is pinned explicitly because it reads as an off-by-one otherwise:
the poly-A scan does not stop at the tail boundary, so A-rich sequence just
upstream legitimately extends the trim (`ACGT`×5 + 30 A's trims 34, not 30).
That is STAR's behaviour, and the test says so in as many words.

The SIMD backends and the `--clipAdapterType CellRanger4` wiring follow; this
commit is the scalar oracle they will be checked against.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/clip/cellranger4.rs | 170 +++++++++++++++++++++++++++++++++
 src/clip/mod.rs         |   8 +-
 src/lib.rs              |   1 +
 src/swalign/mod.rs      | 148 +++++++++++++++++++++++++++++
 src/swalign/scalar.rs   | 205 ++++++++++++++++++++++++++++++++++++++++
 5 files changed, 529 insertions(+), 3 deletions(-)
 create mode 100644 src/clip/cellranger4.rs
 create mode 100644 src/swalign/mod.rs
 create mode 100644 src/swalign/scalar.rs

diff --git a/src/clip/cellranger4.rs b/src/clip/cellranger4.rs
new file mode 100644
index 0000000..5006c64
--- /dev/null
+++ b/src/clip/cellranger4.rs
@@ -0,0 +1,170 @@
+//! `--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.
+//!
+//! STAR does the 5' alignment with the Opal SIMD library. Here it goes through
+//! [`crate::swalign`], which is required to be bit-identical across
+//! instruction sets, so the clip length cannot depend on the machine.
+
+use crate::swalign::{self, Mode, Scoring};
+
+/// 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
+    }
+}
+
+/// How much of the read STAR aligns the TSO against (`ClipCR4::opalFillOneSeq`).
+const CR4_TARGET_LEN: usize = 91;
+
+/// Number of 5' bases to trim as the 10x TSO.
+///
+/// STAR aligns the TSO against the first 91 bases of the read in overlap mode,
+/// asking for the score and the position in the target where it ends, then
+/// applies an acceptance gate: a score below 20 is rejected outright, and
+/// scores of exactly 20 or 21 are rejected if they took too long to reach
+/// (more than 26 and 30 bases respectively). A weak alignment that happens to
+/// run a long way is what that gate is there to catch.
+///
+/// The read is padded to 91 bases with `N` when it is shorter, which is why
+/// the scoring scheme has to treat `N` against `N` as neutral rather than as a
+/// mismatch: otherwise the padding would drag every score down.
+///
+/// Both arguments are numeric base codes; the return value is a count of 5'
+/// bases to clip, `0` when the alignment is rejected.
+pub fn tso_clip(read: &[u8], tso: &[u8]) -> usize {
+    if tso.is_empty() {
+        return 0;
+    }
+    let take = read.len().min(CR4_TARGET_LEN);
+    let mut target = Vec::with_capacity(CR4_TARGET_LEN);
+    target.extend_from_slice(&read[..take]);
+    target.resize(CR4_TARGET_LEN, 4); // N padding
+
+    let a = swalign::align(tso, &target, Mode::Ov, &Scoring::CLIP_CR4);
+    let clip = a.target_end as i64 + 1; // 1-based end == number of bases covered
+
+    let reject = a.score < 20 || (a.score == 20 && clip > 26) || (a.score == 21 && clip > 30);
+    if reject { 0 } else { clip 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()
+    }
+
+    /// The 10x template switch oligo.
+    const TSO: &str = "AAGCAGTGGTATCAACGCAGAGTACATGGG";
+
+    #[test]
+    fn cr4_tso_clip_matches_opal() {
+        // Frozen vector shared with STAR-rs's `cr4_tso_clip_matches_opal`,
+        // which validates the same numbers against Opal itself. A read that
+        // starts with the TSO is clipped by exactly its length.
+        let read = code(&format!("{TSO}ACGTACGTACGTACGTACGTACGTACGTAC"));
+        assert_eq!(tso_clip(&read, &code(TSO)), 30);
+
+        // A read with no TSO is not clipped at all.
+        let read = code("ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT");
+        assert_eq!(tso_clip(&read, &code(TSO)), 0);
+    }
+
+    #[test]
+    fn tso_clip_is_inert_without_an_adapter() {
+        let read = code(&format!("{TSO}ACGTACGT"));
+        assert_eq!(tso_clip(&read, &[]), 0);
+    }
+
+    #[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)
+        );
+    }
+}
diff --git a/src/clip/mod.rs b/src/clip/mod.rs
index 20901ec..27de1e8 100644
--- a/src/clip/mod.rs
+++ b/src/clip/mod.rs
@@ -17,9 +17,11 @@
 //! `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.
+
+pub mod cellranger4;
 
 use crate::io::fastq::encode_base;
 use crate::params::Parameters;
diff --git a/src/lib.rs b/src/lib.rs
index 6fce173..b73fb01 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -40,6 +40,7 @@ pub mod rng;
 pub mod signal;
 pub mod solo;
 pub mod stats;
+pub mod swalign;
 pub mod wasp;
 
 use log::info;
diff --git a/src/swalign/mod.rs b/src/swalign/mod.rs
new file mode 100644
index 0000000..99dd0bb
--- /dev/null
+++ b/src/swalign/mod.rs
@@ -0,0 +1,148 @@
+//! Deterministic Smith-Waterman-family alignment with affine gaps.
+//!
+//! STAR performs exactly one alignment of this shape: the 5' TSO clip of
+//! `--clipAdapterType CellRanger4`, for which it links the Opal C/C++ SIMD
+//! library (`OPAL_MODE_OV` + `OPAL_SEARCH_SCORE_END`). This module provides
+//! the same capability in-tree, with no new dependency.
+//!
+//! # Determinism is a correctness property here
+//!
+//! Opal's overflow handling groups database sequences into buckets sized by the
+//! SIMD vector width, so when 8-bit lanes saturate, the grouping — and with it
+//! the recompute path — depends on which instruction set is available. Sixteen
+//! SSE lanes and thirty-two AVX2 lanes bucket differently, and simde emulates
+//! AVX2 on ARM, so the same input can take a different path on a different
+//! machine.
+//!
+//! That is not acceptable for an aligner whose output is supposed to be
+//! reproducible. Here the rule is: **the vector width is never observable.**
+//! [`scalar`] defines the result; every other backend must agree with it
+//! bit-for-bit, including under saturation, on empty inputs, on `N`, and on
+//! end-position ties. The differential test in [`tests`] is what enforces it.
+//!
+//! # Coordinates and conventions
+//!
+//! Sequences are numeric base codes, `0..=3` for ACGT and `4` for `N`, the same
+//! encoding the rest of the aligner uses. Scores are `i32`; the caller supplies
+//! the scoring scheme.
+
+pub mod scalar;
+
+/// How the ends of the two sequences are treated.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Mode {
+    /// Global: both sequences must be consumed end to end (Needleman-Wunsch).
+    Nw,
+    /// Semi-global on the target: the query must be consumed entirely, the
+    /// target may extend past it on both sides.
+    Hw,
+    /// Overlap: gaps at the start of either sequence and at the end of either
+    /// sequence are free. This is the mode STAR's TSO clip uses.
+    Ov,
+    /// Local (Smith-Waterman): the best-scoring subalignment.
+    Sw,
+}
+
+/// Affine-gap scoring.
+///
+/// `gap_open` and `gap_extend` are the penalties *subtracted*, so both are
+/// given as positive numbers. STAR's ClipCR4 uses `match_score = 1`,
+/// `mismatch = -2`, `gap_open = gap_extend = 2`, and scores `N` against `N` as
+/// zero rather than as a mismatch.
+#[derive(Debug, Clone, Copy)]
+pub struct Scoring {
+    /// Added when two bases are equal.
+    pub match_score: i32,
+    /// Added when two bases differ (negative).
+    pub mismatch: i32,
+    /// Subtracted to open a gap.
+    pub gap_open: i32,
+    /// Subtracted for each base a gap is extended by.
+    pub gap_extend: i32,
+    /// Added when both positions are `N`. Opal treats this pairing as neutral
+    /// rather than as a mismatch, and STAR relies on that when it pads the
+    /// target with `N`.
+    pub n_vs_n: i32,
+}
+
+impl Scoring {
+    /// The scoring STAR uses for the CellRanger4 TSO clip
+    /// (`ClipCR4.cpp`: match +1, mismatch -2, gaps 2, `N`/`N` neutral).
+    pub const CLIP_CR4: Self = Self {
+        match_score: 1,
+        mismatch: -2,
+        gap_open: 2,
+        gap_extend: 2,
+        n_vs_n: 0,
+    };
+
+    /// Score one aligned pair of base codes.
+    #[inline]
+    pub fn pair(&self, q: u8, t: u8) -> i32 {
+        if q == 4 && t == 4 {
+            self.n_vs_n
+        } else if q == t {
+            self.match_score
+        } else {
+            self.mismatch
+        }
+    }
+}
+
+/// The outcome of one alignment.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub struct Alignment {
+    /// Best score under the chosen mode.
+    pub score: i32,
+    /// Zero-based position in the target where that score is reached.
+    ///
+    /// Ties are broken towards the **earlier** column, except that a score
+    /// reached in the final column wins outright: that is what
+    /// `OPAL_SEARCH_SCORE_END` means, and STAR's clip length depends on it.
+    pub target_end: usize,
+}
+
+/// Align `query` against `target` under `mode`.
+///
+/// Dispatches to the fastest backend available for this machine. Every backend
+/// is required to return exactly what [`scalar::align`] would, so the choice is
+/// invisible in the output.
+pub fn align(query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
+    // Only the scalar backend exists so far. The SIMD backends land next, and
+    // the differential test is already written so they cannot land silently
+    // wrong.
+    scalar::align(query, target, mode, scoring)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn empty_inputs_score_zero_in_every_mode() {
+        for mode in [Mode::Nw, Mode::Hw, Mode::Ov, Mode::Sw] {
+            let a = align(&[], &[], mode, &Scoring::CLIP_CR4);
+            assert_eq!(a.score, 0, "{mode:?} on two empty sequences");
+            let a = align(&[0, 1, 2], &[], mode, &Scoring::CLIP_CR4);
+            assert_eq!(a.score, 0, "{mode:?} on an empty target");
+        }
+    }
+
+    #[test]
+    fn n_against_n_is_neutral_not_a_mismatch() {
+        let s = Scoring::CLIP_CR4;
+        assert_eq!(s.pair(4, 4), 0);
+        assert_eq!(s.pair(4, 0), -2);
+        assert_eq!(s.pair(0, 4), -2);
+        assert_eq!(s.pair(0, 0), 1);
+        assert_eq!(s.pair(0, 1), -2);
+    }
+
+    #[test]
+    fn identical_sequences_score_one_per_base() {
+        let seq = [0u8, 1, 2, 3, 0, 1, 2, 3];
+        let a = align(&seq, &seq, Mode::Ov, &Scoring::CLIP_CR4);
+        assert_eq!(a.score, seq.len() as i32);
+        assert_eq!(a.target_end, seq.len() - 1);
+    }
+}
diff --git a/src/swalign/scalar.rs b/src/swalign/scalar.rs
new file mode 100644
index 0000000..53930e6
--- /dev/null
+++ b/src/swalign/scalar.rs
@@ -0,0 +1,205 @@
+//! Portable reference implementation.
+//!
+//! This is the definition of the result. It is written for clarity rather than
+//! speed: the SIMD backends must agree with it bit-for-bit, so it needs to be
+//! obviously correct more than it needs to be fast.
+
+use super::{Alignment, Mode, Scoring};
+
+/// Sentinel for "unreachable". Far below any real score, and far enough from
+/// `i32::MIN` that subtracting a gap penalty cannot wrap.
+const NEG: i32 = i32::MIN / 4;
+
+/// Align `query` against `target`, filling the DP column by column.
+///
+/// Only two columns are ever live, so the working set is `O(|query|)` rather
+/// than the full matrix. That is also what makes the striped SIMD form a
+/// drop-in replacement later: it computes the same columns in the same order.
+pub fn align(query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
+    let ql = query.len();
+    let tl = target.len();
+    if ql == 0 || tl == 0 {
+        return Alignment {
+            score: 0,
+            target_end: 0,
+        };
+    }
+
+    // Whether a gap before the start of each sequence is free.
+    let (free_query_start, free_target_start) = match mode {
+        Mode::Nw => (false, false),
+        Mode::Hw => (false, true),
+        Mode::Ov | Mode::Sw => (true, true),
+    };
+
+    // Column 0 of the previous iteration: H is the score of aligning the first
+    // `r+1` query bases against nothing.
+    let mut prev_h = vec![0i32; ql];
+    let mut prev_e = vec![NEG; ql];
+    if !free_query_start {
+        for (r, h) in prev_h.iter_mut().enumerate() {
+            *h = -(scoring.gap_open + scoring.gap_extend * r as i32);
+        }
+    }
+
+    // Best score seen on the final query row, and the column that first
+    // achieved it. Ties go to the earlier column.
+    let mut best_last_row = NEG;
+    let mut best_last_row_col = 0usize;
+    // Best score anywhere in the most recent column.
+    let mut last_col_max = NEG;
+    // Best score anywhere in the matrix, for local mode.
+    let mut best_anywhere = 0i32;
+    let mut best_anywhere_col = 0usize;
+
+    for (c, &tc) in target.iter().enumerate() {
+        // Top of the column: aligning the first `c+1` target bases against
+        // nothing.
+        let mut up_h = if free_target_start {
+            0
+        } else {
+            -(scoring.gap_open + scoring.gap_extend * c as i32)
+        };
+        // The diagonal predecessor, i.e. the cell up and to the left. Column 0
+        // has no predecessor and scores 0 either way, so it falls out of the
+        // free-start case.
+        let mut diag_h = if free_target_start || c == 0 {
+            0
+        } else {
+            -(scoring.gap_open + scoring.gap_extend * (c as i32 - 1))
+        };
+        let mut up_f = NEG;
+        let mut col_max = NEG;
+        let mut h = NEG;
+
+        for r in 0..ql {
+            // Gap in the query (moving right): open from H, or extend E.
+            let e = (prev_h[r] - scoring.gap_open).max(prev_e[r] - scoring.gap_extend);
+            // Gap in the target (moving down): open from H, or extend F.
+            let f = (up_h - scoring.gap_open).max(up_f - scoring.gap_extend);
+            // Match or mismatch on the diagonal.
+            let d = diag_h + scoring.pair(query[r], tc);
+
+            h = e.max(f).max(d);
+            if mode == Mode::Sw {
+                // Local alignment never carries a negative prefix forward.
+                h = h.max(0);
+                if h > best_anywhere {
+                    best_anywhere = h;
+                    best_anywhere_col = c;
+                }
+            }
+            if h > col_max {
+                col_max = h;
+            }
+
+            up_f = f;
+            up_h = h;
+            diag_h = prev_h[r];
+            prev_e[r] = e;
+            prev_h[r] = h;
+        }
+
+        // `h` now holds the last query row for this column.
+        if h > best_last_row {
+            best_last_row = h;
+            best_last_row_col = c;
+        }
+        last_col_max = col_max;
+    }
+
+    match mode {
+        Mode::Sw => Alignment {
+            score: best_anywhere,
+            target_end: best_anywhere_col,
+        },
+        Mode::Nw => Alignment {
+            // Global: the corner cell, which is the last row of the last
+            // column.
+            score: prev_h[ql - 1],
+            target_end: tl - 1,
+        },
+        Mode::Hw => Alignment {
+            // Semi-global on the target: the query must be consumed, so the
+            // answer lives on the last query row.
+            score: best_last_row,
+            target_end: best_last_row_col,
+        },
+        Mode::Ov => {
+            // Overlap: either the query ran out (best on the last row) or the
+            // target ran out (best in the last column). A score reached in the
+            // final column wins the tie, which is what `OPAL_SEARCH_SCORE_END`
+            // specifies and what STAR's clip length depends on.
+            let score = last_col_max.max(best_last_row);
+            let target_end = if last_col_max >= best_last_row {
+                tl - 1
+            } else {
+                best_last_row_col
+            };
+            Alignment { score, target_end }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::swalign::Scoring;
+
+    const S: Scoring = Scoring::CLIP_CR4;
+
+    #[test]
+    fn global_mode_pays_for_every_unaligned_base() {
+        // Query is a prefix of the target; NW must pay to skip the tail.
+        let q = [0u8, 1, 2];
+        let t = [0u8, 1, 2, 3, 3];
+        let nw = align(&q, &t, Mode::Nw, &S);
+        let ov = align(&q, &t, Mode::Ov, &S);
+        assert_eq!(ov.score, 3, "overlap: the tail is free");
+        assert!(
+            nw.score < ov.score,
+            "global should pay for the tail, got {} vs {}",
+            nw.score,
+            ov.score
+        );
+    }
+
+    #[test]
+    fn local_mode_ignores_flanking_mismatch() {
+        // A clean 4-base core buried in mismatching flanks.
+        let q = [3u8, 3, 0, 1, 2, 3, 3, 3];
+        let t = [1u8, 1, 0, 1, 2, 3, 1, 1];
+        let sw = align(&q, &t, Mode::Sw, &S);
+        assert!(
+            sw.score >= 4,
+            "local should find the shared core, got {}",
+            sw.score
+        );
+    }
+
+    #[test]
+    fn local_score_is_never_negative() {
+        let q = [0u8, 0, 0, 0];
+        let t = [3u8, 3, 3, 3];
+        assert_eq!(align(&q, &t, Mode::Sw, &S).score, 0);
+    }
+
+    #[test]
+    fn overlap_prefers_the_final_column_on_a_tie() {
+        // A query that matches equally well at two positions: the one that
+        // runs to the end of the target must win, because that is the tie-break
+        // STAR's clip length is built on.
+        let q = [0u8, 1];
+        let t = [0u8, 1, 4, 4, 0, 1];
+        let a = align(&q, &t, Mode::Ov, &S);
+        assert_eq!(a.target_end, t.len() - 1);
+    }
+
+    #[test]
+    fn affine_gaps_cost_less_than_repeated_opens() {
+        // One 3-base gap should beat three separate 1-base gaps.
+        let one_long = S.gap_open + S.gap_extend * 3;
+        let three_short = 3 * (S.gap_open + S.gap_extend);
+        assert!(one_long < three_short);
+    }
+}

From bf61301977da1c958c43e512ab8f2a8a3e6af8bb Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Tue, 28 Jul 2026 22:57:58 +0200
Subject: [PATCH 2/8] feat(clip): --clipAdapterType CellRanger4 now actually
 clips
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The value parsed and validated but did nothing: the clip module declared the
mode out of scope, so a user who asked for CellRanger4 clipping silently got
none. That is a correctness hole rather than a missing feature, since the
run appears to succeed.

`ClipParams` gains the mode selector and the 5' adapter, and `clip_mate`
branches to `clip_mate_cellranger4`, which replaces the Hamming rules with
STAR's two CR4 trims: the TSO overlap alignment at the 5' end and the poly-A
scan at the 3'. The fixed `--clip{5,3}pNbases` still apply first, as in the
Hamming path.

Adds `--clip5pAdapterSeq` and `--clip5pAdapterMMp`. The 5' adapter is only
configured for the first mate under STAR, so mate 2 simply has no TSO and its
5' trim reduces to the fixed clip — covered by a test rather than left implicit.

The wiring tests check the whole shape: a read with a TSO at the 5' end, a
poly-A tail at the 3' and mappable sequence between comes back with exactly
that middle section surviving.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/clip/mod.rs   | 118 +++++++++++++++++++++++++++++++++++++++++++++-
 src/params/mod.rs |  10 ++++
 2 files changed, 127 insertions(+), 1 deletion(-)

diff --git a/src/clip/mod.rs b/src/clip/mod.rs
index 27de1e8..973fdfb 100644
--- a/src/clip/mod.rs
+++ b/src/clip/mod.rs
@@ -19,7 +19,8 @@
 //!
 //! `--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.
+//! [`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;
 
@@ -50,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`
@@ -67,7 +74,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(),
@@ -122,6 +136,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;
@@ -157,6 +175,98 @@ 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, then the TSO overlap alignment on what is left.
+    let mut c5 = p.five.n.min(len);
+    if !p.five_adapter.is_empty() {
+        c5 += cellranger4::tso_clip(&read[c5..], &p.five_adapter).min(len - c5);
+    }
+    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_clips_the_tso_and_the_polya_tail() {
+        // TSO at the 5' end, a clean poly-A tail at the 3', mappable sequence
+        // in between.
+        let body = "CGTCGTCGTCGTCGTCGTCGTCGTCGTCGT";
+        let read = code(&format!("{TSO}{body}{}", "A".repeat(30)));
+        let (c5, c3) = clip_mate(&read, &cr4_params(TSO));
+        assert_eq!(c5, 30, "the TSO should be clipped from the 5' end");
+        assert_eq!(c3, 30, "the poly-A tail should be clipped from the 3' end");
+        // What survives is exactly the body.
+        assert_eq!(&read[c5..read.len() - c3], code(body).as_slice());
+    }
+
+    #[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_clips_first() {
+        let body = "CGTCGTCGTCGTCGTCGTCGTCGTCGTCGT";
+        let read = code(&format!("{TSO}{body}{}", "A".repeat(30)));
+        let mut p = cr4_params(TSO);
+        p.five.n = 5;
+        let (c5, _) = clip_mate(&read, &p);
+        // The fixed 5 bases come off, then the rest of the TSO is still found.
+        assert_eq!(c5, 30);
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -185,6 +295,8 @@ mod tests {
     #[test]
     fn fixed_5p_3p() {
         let p = ClipParams {
+            cellranger4: false,
+            five_adapter: Vec::new(),
             five: ClipEnd {
                 n: 3,
                 ..Default::default()
@@ -219,6 +331,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()
@@ -242,6 +356,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()
diff --git a/src/params/mod.rs b/src/params/mod.rs
index 0536a85..7f5cf4a 100644
--- a/src/params/mod.rs
+++ b/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,

From 079d07a22af8f7eb6279cfbab57e914b0274e307 Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Tue, 28 Jul 2026 22:58:08 +0200
Subject: [PATCH 3/8] docs: changelog for CellRanger4 clipping

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 CHANGELOG.md | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index b6e4893..56bc7c0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,18 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Features
 
+- **`--clipAdapterType CellRanger4`** — the 10x Chromium v4 clipping rules: a
+  5' template-switch-oligo trim and a 3' poly-A trim, with `--clip5pAdapterSeq`
+  and `--clip5pAdapterMMp`.
+
+  The 5' trim is an overlap alignment, for which STAR links the Opal C/C++ SIMD
+  library. New in-tree module `swalign` provides it with no new dependency:
+  affine-gap alignment in `NW`/`HW`/`OV`/`SW` modes, with a portable scalar
+  path that defines the result. Determinism is treated as a correctness
+  property — Opal's overflow buckets are sized by the SIMD vector width, so its
+  recompute path depends on the available instruction set; here the vector
+  width is never observable.
+
 - **STARsolo single-cell quantification (`--soloType`)** — the 10x
   Chromium / plate-based count-matrix pipeline, ported from STAR and
   verified against real STARsolo (#90).

From 1efcca5f149fcfbbe55aaa82396f1054ff93a06c Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Tue, 28 Jul 2026 23:08:33 +0200
Subject: [PATCH 4/8] test(swalign): differential harness that every backend
 has to pass
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Adds `Backend` and `differential_check`. The contract a backend signs is
narrow: given the same inputs it returns exactly what the scalar reference
returns — not the same score with a different end position, the same
`Alignment`.

The harness sweeps 3612 cases: every mode against every pair of lengths from
{0, 1, 2, 7, 8, 9, 15, 16, 17, 30, 31, 33, 64, 91, 128}, at four `N` densities
from none to all. Those lengths sit either side of the lane counts a vectorised
kernel cares about, and 91 is the size STAR's TSO clip actually uses.

Then the case that matters most: long exact matches scoring one per base, which
overflow an 8-bit lane long before they trouble the scalar's `i32`. A backend
that escalates lane width incorrectly fails there and nowhere else, so it is
checked explicitly rather than left to the random draw.

Sequences come from the in-tree splitmix64, so a failure reproduces exactly and
the error names the case, the mode and both sequences.

`Backend::detect()` currently returns `Scalar`, and the harness proves it
consistent with itself. That is tautological today and deliberately so: it
proves the harness runs, reaches the saturation cases and can name a failure,
so the SIMD backends inherit machinery that is known to work rather than
machinery written at the same time as the code it judges.

The backends themselves (SSE2 baseline, runtime-detected AVX2, aarch64 NEON)
are next. `differential_check` is the gate they have to clear before
`detect()` will offer them.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/swalign/backend.rs | 163 +++++++++++++++++++++++++++++++++++++++++
 src/swalign/mod.rs     |   8 +-
 2 files changed, 167 insertions(+), 4 deletions(-)
 create mode 100644 src/swalign/backend.rs

diff --git a/src/swalign/backend.rs b/src/swalign/backend.rs
new file mode 100644
index 0000000..9e9209f
--- /dev/null
+++ b/src/swalign/backend.rs
@@ -0,0 +1,163 @@
+//! Backend selection, and the differential harness that keeps backends honest.
+//!
+//! The contract every backend signs: given the same inputs it returns exactly
+//! what [`super::scalar::align`] returns. Not "within rounding", not "the same
+//! score with a different end position" — the same [`Alignment`].
+//!
+//! [`differential_check`] is what enforces that. It is written here rather than
+//! in a test module so a backend can be checked from a test, a benchmark or a
+//! debug session without duplicating the generator.
+
+use super::{Alignment, Mode, Scoring, scalar};
+
+/// Which implementation actually ran.
+///
+/// Exposed so tests can assert that a machine with the hardware really used it,
+/// rather than silently falling back and reporting a pass that proves nothing.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Backend {
+    /// The portable reference. Always available, and the definition of the
+    /// result.
+    Scalar,
+}
+
+impl Backend {
+    /// The best backend this machine can run.
+    pub fn detect() -> Self {
+        // SIMD backends land here as they are written. Each one has to pass
+        // `differential_check` on this machine before being offered.
+        Backend::Scalar
+    }
+
+    /// Run this backend.
+    pub fn align(self, query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
+        match self {
+            Backend::Scalar => scalar::align(query, target, mode, scoring),
+        }
+    }
+}
+
+/// A deterministic sequence generator for the differential harness.
+///
+/// Seeded from the in-tree splitmix64 so the cases are identical on every
+/// machine and every run: a backend that fails does so reproducibly, on a case
+/// the report can name.
+struct SeqGen(crate::rng::SplitMix64);
+
+impl SeqGen {
+    fn new(seed: u64) -> Self {
+        Self(crate::rng::SplitMix64::seed(seed))
+    }
+
+    /// A sequence of `len` base codes. `n_rate` out of 16 bases are `N`, so the
+    /// generator covers the `N`-heavy inputs STAR's padding produces as well as
+    /// clean ones.
+    fn seq(&mut self, len: usize, n_rate: u64) -> Vec<u8> {
+        (0..len)
+            .map(|_| {
+                let r = self.0.next_u64();
+                if r % 16 < n_rate {
+                    4
+                } else {
+                    (r >> 8) as u8 % 4
+                }
+            })
+            .collect()
+    }
+}
+
+/// Compare `backend` against the scalar reference over a spread of inputs.
+///
+/// Returns `Err` with a description of the first disagreement, naming the case
+/// so it can be reproduced. `Ok(n)` reports how many cases were checked.
+///
+/// The spread is deliberately awkward: empty and length-1 sequences, queries
+/// longer than targets, all-`N` inputs, and long runs that push scores far
+/// enough to saturate a narrow lane. Saturation is where a SIMD backend is most
+/// likely to diverge, so it is not left to chance.
+pub fn differential_check(backend: Backend, seed: u64) -> Result<usize, String> {
+    let scoring = Scoring::CLIP_CR4;
+    let mut rng = SeqGen::new(seed);
+    let mut checked = 0usize;
+
+    // Lengths chosen around the boundaries a vectorised kernel cares about:
+    // zero, one, just under and just over a typical lane count, and the 91
+    // STAR actually uses.
+    const LENS: &[usize] = &[0, 1, 2, 7, 8, 9, 15, 16, 17, 30, 31, 33, 64, 91, 128];
+
+    for &ql in LENS {
+        for &tl in LENS {
+            for &n_rate in &[0u64, 1, 8, 16] {
+                for mode in [Mode::Nw, Mode::Hw, Mode::Ov, Mode::Sw] {
+                    let q = rng.seq(ql, n_rate);
+                    let t = rng.seq(tl, n_rate);
+                    let want = scalar::align(&q, &t, mode, &scoring);
+                    let got = backend.align(&q, &t, mode, &scoring);
+                    if got != want {
+                        return Err(format!(
+                            "{backend:?} disagrees with scalar on {mode:?}, \
+                             |q|={ql} |t|={tl} n_rate={n_rate}/16: \
+                             scalar {want:?}, backend {got:?}\n  query  {q:?}\n  target {t:?}"
+                        ));
+                    }
+                    checked += 1;
+                }
+            }
+        }
+    }
+
+    // Saturation: a long exact match scores one per base, which overflows an
+    // 8-bit lane well before it overflows the scalar's i32. A backend that
+    // escalates lane width incorrectly fails here and nowhere else.
+    for &len in &[200usize, 400, 1000] {
+        let q = rng.seq(len, 0);
+        for mode in [Mode::Nw, Mode::Hw, Mode::Ov, Mode::Sw] {
+            let want = scalar::align(&q, &q, mode, &scoring);
+            let got = backend.align(&q, &q, mode, &scoring);
+            if got != want {
+                return Err(format!(
+                    "{backend:?} disagrees with scalar under saturation, \
+                     {mode:?}, len={len}: scalar {want:?}, backend {got:?}"
+                ));
+            }
+            checked += 1;
+        }
+    }
+
+    Ok(checked)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn scalar_is_consistent_with_itself() {
+        // Tautological for the scalar backend, but it proves the harness runs,
+        // covers every mode and length pair, and reaches the saturation cases.
+        // A backend added later inherits a harness that is known to work.
+        let n = differential_check(Backend::Scalar, 0x5EED).expect("scalar vs scalar");
+        assert!(n > 900, "harness covered only {n} cases");
+    }
+
+    #[test]
+    fn the_detected_backend_agrees_with_scalar_on_this_machine() {
+        // The test that matters once SIMD backends exist: whatever this
+        // machine selects must match the reference on this machine.
+        let backend = Backend::detect();
+        if let Err(e) = differential_check(backend, 0x00C0_FFEE) {
+            panic!("{e}");
+        }
+    }
+
+    #[test]
+    fn generator_is_reproducible() {
+        // The harness is only useful if a failure can be reproduced, which
+        // needs the sequences to be identical run to run.
+        let a = SeqGen::new(7).seq(64, 4);
+        let b = SeqGen::new(7).seq(64, 4);
+        assert_eq!(a, b);
+        assert!(a.contains(&4), "n_rate 4/16 should produce Ns");
+        assert!(a.iter().any(|&b| b < 4), "and also real bases");
+    }
+}
diff --git a/src/swalign/mod.rs b/src/swalign/mod.rs
index 99dd0bb..04412a7 100644
--- a/src/swalign/mod.rs
+++ b/src/swalign/mod.rs
@@ -26,8 +26,11 @@
 //! encoding the rest of the aligner uses. Scores are `i32`; the caller supplies
 //! the scoring scheme.
 
+pub mod backend;
 pub mod scalar;
 
+pub use backend::Backend;
+
 /// How the ends of the two sequences are treated.
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
 pub enum Mode {
@@ -108,10 +111,7 @@ pub struct Alignment {
 /// is required to return exactly what [`scalar::align`] would, so the choice is
 /// invisible in the output.
 pub fn align(query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
-    // Only the scalar backend exists so far. The SIMD backends land next, and
-    // the differential test is already written so they cannot land silently
-    // wrong.
-    scalar::align(query, target, mode, scoring)
+    Backend::detect().align(query, target, mode, scoring)
 }
 
 #[cfg(test)]

From bfec07e12f5fa01549c714f999b01d8f64234b0d Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Wed, 29 Jul 2026 00:50:35 +0200
Subject: [PATCH 5/8] feat(swalign): NEON backend, checked against the scalar
 path on every case
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Adds an aarch64 NEON backend and selects it via `Backend::detect()`.

**Anti-diagonals, not Farrar stripes.** The striped layout is faster, but its
lazy-F correction and stripe-relative indexing make the end position awkward to
extract — and the end position is exactly what STAR's clip length is built on.
A backend that is quicker but disagrees about where an alignment ends is
worthless here.

Cells on one anti-diagonal are mutually independent: `H(r,c)` depends on
`E(r,c-1)` and `F(r-1,c)` on `d-1`, and `H(r-1,c-1)` on `d-2`. So a whole
anti-diagonal computes in parallel with no correction pass and no reordering,
which makes agreement with the scalar path a property of the layout rather than
something to test for and hope about. The cost is strided access, which is
irrelevant at the 30×91 matrix STAR's TSO clip uses.

The harness earned its place immediately. The first version disagreed on local
alignment: `|q|=7 |t|=2`, same score, different end. The scalar sweeps
column-major and keeps the first cell to reach a new maximum, so among equal
scores the smallest `(c, r)` wins; visiting anti-diagonals changes the arrival
order. The tie-break key is now compared explicitly instead of being inferred
from arrival. That is precisely the class of bug the differential check exists
to catch, and it would have been invisible in a benchmark.

`the_detected_backend_agrees_with_scalar_on_this_machine` also now asserts that
aarch64 really selected NEON. A pass proves nothing if the machine quietly fell
back to the reference.

x86 backends (SSE2 baseline, runtime-detected AVX2) still to come; they will use
the same layout and the same harness, which is now known to catch a real
disagreement rather than merely being present.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/swalign/backend.rs |  21 +++-
 src/swalign/mod.rs     |   2 +
 src/swalign/neon.rs    | 238 +++++++++++++++++++++++++++++++++++++++++
 3 files changed, 259 insertions(+), 2 deletions(-)
 create mode 100644 src/swalign/neon.rs

diff --git a/src/swalign/backend.rs b/src/swalign/backend.rs
index 9e9209f..b6bcd1d 100644
--- a/src/swalign/backend.rs
+++ b/src/swalign/backend.rs
@@ -19,13 +19,18 @@ pub enum Backend {
     /// The portable reference. Always available, and the definition of the
     /// result.
     Scalar,
+    /// aarch64 NEON, computing one anti-diagonal at a time.
+    #[cfg(target_arch = "aarch64")]
+    Neon,
 }
 
 impl Backend {
     /// The best backend this machine can run.
     pub fn detect() -> Self {
-        // SIMD backends land here as they are written. Each one has to pass
-        // `differential_check` on this machine before being offered.
+        #[cfg(target_arch = "aarch64")]
+        if super::neon::is_available() {
+            return Backend::Neon;
+        }
         Backend::Scalar
     }
 
@@ -33,6 +38,8 @@ impl Backend {
     pub fn align(self, query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
         match self {
             Backend::Scalar => scalar::align(query, target, mode, scoring),
+            #[cfg(target_arch = "aarch64")]
+            Backend::Neon => super::neon::align(query, target, mode, scoring),
         }
     }
 }
@@ -145,6 +152,16 @@ mod tests {
         // The test that matters once SIMD backends exist: whatever this
         // machine selects must match the reference on this machine.
         let backend = Backend::detect();
+
+        // A pass proves nothing if the machine quietly fell back to the
+        // reference, so assert that the hardware backend really was selected.
+        #[cfg(target_arch = "aarch64")]
+        assert_eq!(
+            backend,
+            Backend::Neon,
+            "aarch64 must select NEON; a silent fallback would make this test vacuous"
+        );
+
         if let Err(e) = differential_check(backend, 0x00C0_FFEE) {
             panic!("{e}");
         }
diff --git a/src/swalign/mod.rs b/src/swalign/mod.rs
index 04412a7..ed305e9 100644
--- a/src/swalign/mod.rs
+++ b/src/swalign/mod.rs
@@ -27,6 +27,8 @@
 //! the scoring scheme.
 
 pub mod backend;
+#[cfg(target_arch = "aarch64")]
+pub mod neon;
 pub mod scalar;
 
 pub use backend::Backend;
diff --git a/src/swalign/neon.rs b/src/swalign/neon.rs
new file mode 100644
index 0000000..1121b95
--- /dev/null
+++ b/src/swalign/neon.rs
@@ -0,0 +1,238 @@
+//! aarch64 NEON backend.
+//!
+//! # Why anti-diagonals rather than Farrar stripes
+//!
+//! The striped layout is faster, but its lazy-F correction and its
+//! stripe-relative indexing make the *end position* awkward to extract, and the
+//! end position is exactly what STAR's clip length is built on. A backend that
+//! is quicker but disagrees with [`super::scalar`] about where an alignment
+//! ends is worthless here.
+//!
+//! Cells on one anti-diagonal `d = r + c` are mutually independent: `H(r,c)`
+//! depends on `E(r,c-1)` and `F(r-1,c)`, both on `d-1`, and on `H(r-1,c-1)` on
+//! `d-2`. So a whole anti-diagonal can be computed in parallel with no
+//! correction pass and no reordering, which makes agreement with the scalar
+//! path a property of the layout rather than something to be tested for and
+//! hoped about.
+//!
+//! The cost is strided access. For the 30×91 matrix STAR's TSO clip actually
+//! uses, that is irrelevant.
+
+use std::arch::aarch64::{vaddq_s32, vdupq_n_s32, vld1q_s32, vmaxq_s32, vst1q_s32, vsubq_s32};
+
+use super::{Alignment, Mode, Scoring};
+
+/// Matches `scalar::NEG`: far below any real score, far enough from `i32::MIN`
+/// that subtracting a gap penalty cannot wrap.
+const NEG: i32 = i32::MIN / 4;
+
+/// Lanes per NEON vector at 32-bit width.
+const LANES: usize = 4;
+
+/// Align `query` against `target` using NEON.
+///
+/// # Safety
+///
+/// The caller must have established that NEON is available. On aarch64 it is
+/// architecturally guaranteed, so [`is_available`] is a constant.
+pub fn align(query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
+    let ql = query.len();
+    let tl = target.len();
+    if ql == 0 || tl == 0 {
+        return Alignment {
+            score: 0,
+            target_end: 0,
+        };
+    }
+
+    let (free_query_start, free_target_start) = match mode {
+        Mode::Nw => (false, false),
+        Mode::Hw => (false, true),
+        Mode::Ov | Mode::Sw => (true, true),
+    };
+
+    // Score of aligning the first `r+1` query bases against nothing, i.e. the
+    // cell one column to the left of column 0.
+    let left_h = |r: usize| -> i32 {
+        if free_query_start {
+            0
+        } else {
+            -(scoring.gap_open + scoring.gap_extend * r as i32)
+        }
+    };
+    // Score of aligning the first `c+1` target bases against nothing, i.e. the
+    // cell one row above row 0.
+    let top_h = |c: i64| -> i32 {
+        // A column before the first scores 0 either way, so it folds into the
+        // free-start case.
+        if free_target_start || c < 0 {
+            0
+        } else {
+            -(scoring.gap_open + scoring.gap_extend * c as i32)
+        }
+    };
+
+    // Rows of the two previous anti-diagonals, indexed by `r`.
+    let mut h1 = vec![NEG; ql]; // H on d-1
+    let mut e1 = vec![NEG; ql]; // E on d-1
+    let mut f1 = vec![NEG; ql]; // F on d-1
+    let mut h2 = vec![NEG; ql]; // H on d-2
+    let mut h0 = vec![NEG; ql];
+    let mut e0 = vec![NEG; ql];
+    let mut f0 = vec![NEG; ql];
+
+    let mut best_last_row = NEG;
+    let mut best_last_row_col = 0usize;
+    let mut last_col_max = NEG;
+    let mut best_anywhere = 0i32;
+    let mut best_anywhere_col = 0usize;
+    // The scalar sweeps column-major and keeps the first cell to reach a new
+    // maximum, so among equal scores the smallest `(c, r)` wins. Visiting
+    // anti-diagonals changes the order, so the key has to be compared
+    // explicitly rather than inferred from arrival.
+    let mut best_anywhere_row = 0usize;
+
+    let go = scoring.gap_open;
+    let ge = scoring.gap_extend;
+
+    for d in 0..(ql + tl - 1) {
+        let r_lo = d.saturating_sub(tl - 1);
+        let r_hi = d.min(ql - 1);
+
+        // SAFETY: every load and store below is a plain lane-wise op on stack
+        // scalars gathered by index; no pointer arithmetic escapes the slices,
+        // which are all length `ql` and indexed within `r_lo..=r_hi`.
+        unsafe {
+            let vgo = vdupq_n_s32(go);
+            let vge = vdupq_n_s32(ge);
+
+            let mut row = r_lo;
+            while row <= r_hi {
+                let lanes = LANES.min(r_hi - row + 1);
+
+                // Gather the four predecessor terms for lanes r..r+n.
+                let mut e_prev_h = [NEG; LANES]; // H(r, c-1)
+                let mut e_prev_e = [NEG; LANES]; // E(r, c-1)
+                let mut f_prev_h = [NEG; LANES]; // H(r-1, c)
+                let mut f_prev_f = [NEG; LANES]; // F(r-1, c)
+                let mut diag = [NEG; LANES]; // H(r-1, c-1)
+                let mut sub = [0i32; LANES]; // score(q[r], t[c])
+
+                for k in 0..lanes {
+                    let rr = row + k;
+                    let col = d - rr;
+                    // (rr, c-1) lives on d-1 at row rr.
+                    if col == 0 {
+                        e_prev_h[k] = left_h(rr);
+                        e_prev_e[k] = NEG;
+                    } else {
+                        e_prev_h[k] = h1[rr];
+                        e_prev_e[k] = e1[rr];
+                    }
+                    // (rr-1, c) lives on d-1 at row rr-1.
+                    if rr == 0 {
+                        f_prev_h[k] = top_h(col as i64);
+                        f_prev_f[k] = NEG;
+                    } else {
+                        f_prev_h[k] = h1[rr - 1];
+                        f_prev_f[k] = f1[rr - 1];
+                    }
+                    // (rr-1, c-1) lives on d-2 at row rr-1.
+                    diag[k] = if rr == 0 {
+                        top_h(col as i64 - 1)
+                    } else if col == 0 {
+                        left_h(rr - 1)
+                    } else {
+                        h2[rr - 1]
+                    };
+                    sub[k] = scoring.pair(query[rr], target[col]);
+                }
+
+                let ins = vmaxq_s32(
+                    vsubq_s32(vld1q_s32(e_prev_h.as_ptr()), vgo),
+                    vsubq_s32(vld1q_s32(e_prev_e.as_ptr()), vge),
+                );
+                let del = vmaxq_s32(
+                    vsubq_s32(vld1q_s32(f_prev_h.as_ptr()), vgo),
+                    vsubq_s32(vld1q_s32(f_prev_f.as_ptr()), vge),
+                );
+                let diagv = vaddq_s32(vld1q_s32(diag.as_ptr()), vld1q_s32(sub.as_ptr()));
+                let mut cur = vmaxq_s32(vmaxq_s32(ins, del), diagv);
+                if mode == Mode::Sw {
+                    cur = vmaxq_s32(cur, vdupq_n_s32(0));
+                }
+
+                let mut hs = [0i32; LANES];
+                let mut es = [0i32; LANES];
+                let mut fs = [0i32; LANES];
+                vst1q_s32(hs.as_mut_ptr(), cur);
+                vst1q_s32(es.as_mut_ptr(), ins);
+                vst1q_s32(fs.as_mut_ptr(), del);
+
+                for k in 0..lanes {
+                    let rr = row + k;
+                    h0[rr] = hs[k];
+                    e0[rr] = es[k];
+                    f0[rr] = fs[k];
+                }
+                row += lanes;
+            }
+        }
+
+        // Bookkeeping, in increasing column order so ties go to the earlier
+        // column exactly as the scalar's column loop does.
+        for (rr, &cell) in h0.iter().enumerate().take(r_hi + 1).skip(r_lo) {
+            let col = d - rr;
+            if mode == Mode::Sw
+                && (cell > best_anywhere
+                    || (cell == best_anywhere
+                        && (col, rr) < (best_anywhere_col, best_anywhere_row)))
+            {
+                best_anywhere = cell;
+                best_anywhere_col = col;
+                best_anywhere_row = rr;
+            }
+            if rr == ql - 1 && cell > best_last_row {
+                best_last_row = cell;
+                best_last_row_col = col;
+            }
+            if col == tl - 1 {
+                last_col_max = last_col_max.max(cell);
+            }
+        }
+
+        std::mem::swap(&mut h2, &mut h1);
+        std::mem::swap(&mut h1, &mut h0);
+        std::mem::swap(&mut e1, &mut e0);
+        std::mem::swap(&mut f1, &mut f0);
+    }
+
+    match mode {
+        Mode::Sw => Alignment {
+            score: best_anywhere,
+            target_end: best_anywhere_col,
+        },
+        Mode::Nw => Alignment {
+            score: h1[ql - 1],
+            target_end: tl - 1,
+        },
+        Mode::Hw => Alignment {
+            score: best_last_row,
+            target_end: best_last_row_col,
+        },
+        Mode::Ov => {
+            let score = last_col_max.max(best_last_row);
+            let target_end = if last_col_max >= best_last_row {
+                tl - 1
+            } else {
+                best_last_row_col
+            };
+            Alignment { score, target_end }
+        }
+    }
+}
+
+/// NEON is architecturally guaranteed on aarch64.
+pub const fn is_available() -> bool {
+    true
+}

From 8632649f11d94ad9e48b5e55d5711896a2408ccb Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Wed, 29 Jul 2026 00:57:40 +0200
Subject: [PATCH 6/8] test(swalign): check every backend this machine can run,
 not just the detected one

`detect()` returns one backend, so testing only that one leaves the others
unexercised. On a machine with AVX2 the SSE2 baseline would never run under
test, despite being what older hardware executes.

`every_available_backend_agrees_with_scalar` enumerates every backend the
machine can actually run and puts each through the differential check. Today
that is scalar and, on aarch64, NEON; it gains the x86 backends without further
edits when they land.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/swalign/backend.rs | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/src/swalign/backend.rs b/src/swalign/backend.rs
index b6bcd1d..48ddf49 100644
--- a/src/swalign/backend.rs
+++ b/src/swalign/backend.rs
@@ -138,6 +138,30 @@ pub fn differential_check(backend: Backend, seed: u64) -> Result<usize, String>
 mod tests {
     use super::*;
 
+    /// Every backend this machine can run, not just the one `detect()` picks.
+    ///
+    /// On x86 that matters: a machine with AVX2 would otherwise never exercise
+    /// SSE2, and the baseline is what runs on everything older.
+    fn available_backends() -> Vec<Backend> {
+        let mut v = vec![Backend::Scalar];
+        #[cfg(target_arch = "aarch64")]
+        if super::super::neon::is_available() {
+            v.push(Backend::Neon);
+        }
+        v
+    }
+
+    #[test]
+    fn every_available_backend_agrees_with_scalar() {
+        // Not just the detected one: on a machine with AVX2, SSE2 would
+        // otherwise go unchecked despite being what older hardware runs.
+        for backend in available_backends() {
+            if let Err(e) = differential_check(backend, 0x0BAD_5EED) {
+                panic!("{e}");
+            }
+        }
+    }
+
     #[test]
     fn scalar_is_consistent_with_itself() {
         // Tautological for the scalar backend, but it proves the harness runs,

From 6dd63aed30a660cce97b07b9edb8d35b2744c29e Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Wed, 29 Jul 2026 01:40:56 +0200
Subject: [PATCH 7/8] refactor(clip): drop the in-tree SW engine, keep the
 poly-A trim
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Per review on #148: a full deterministic-SIMD Smith-Waterman is a second
aligner and belongs in its own crate, with its own tests and verification,
rather than inside rustar. That is where this started before I chose in-tree to
avoid a dependency, and I had the trade the wrong way round.

Removes `src/swalign` entirely — scalar, NEON, differential harness — and the
5' TSO trim that depended on it.

What remains needs no alignment engine and is rustar's own:

- `poly_tail_3p`, the 3' poly-A trim, a plain scan;
- the `clip_mate_cellranger4` plumbing and `ClipParams` fields;
- `--clip5pAdapterSeq` and `--clip5pAdapterMMp`.

`--clipAdapterType CellRanger4` stays accepted, as on main, because an existing
integration test depends on that and rejecting it would break a contract this
branch has no business changing. On main it does nothing at all; here it now
performs the poly-A trim, which is strictly closer to correct.

Since it is no longer a complete implementation of the mode, `clip_params_from`
warns once when a TSO is configured that the 5' trim is not applied. Half a
mode running silently is the failure I wanted to avoid; a warning is the honest
version of it until the SW crate exists.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 CHANGELOG.md            |  16 +--
 src/clip/cellranger4.rs |  66 +----------
 src/clip/mod.rs         |  43 ++++----
 src/lib.rs              |   1 -
 src/swalign/backend.rs  | 204 ----------------------------------
 src/swalign/mod.rs      | 150 -------------------------
 src/swalign/neon.rs     | 238 ----------------------------------------
 src/swalign/scalar.rs   | 205 ----------------------------------
 8 files changed, 29 insertions(+), 894 deletions(-)
 delete mode 100644 src/swalign/backend.rs
 delete mode 100644 src/swalign/mod.rs
 delete mode 100644 src/swalign/neon.rs
 delete mode 100644 src/swalign/scalar.rs

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 56bc7c0..f5129e1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,17 +13,11 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Features
 
-- **`--clipAdapterType CellRanger4`** — the 10x Chromium v4 clipping rules: a
-  5' template-switch-oligo trim and a 3' poly-A trim, with `--clip5pAdapterSeq`
-  and `--clip5pAdapterMMp`.
-
-  The 5' trim is an overlap alignment, for which STAR links the Opal C/C++ SIMD
-  library. New in-tree module `swalign` provides it with no new dependency:
-  affine-gap alignment in `NW`/`HW`/`OV`/`SW` modes, with a portable scalar
-  path that defines the result. Determinism is treated as a correctness
-  property — Opal's overflow buckets are sized by the SIMD vector width, so its
-  recompute path depends on the available instruction set; here the vector
-  width is never observable.
+- **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
diff --git a/src/clip/cellranger4.rs b/src/clip/cellranger4.rs
index 5006c64..e0d7a45 100644
--- a/src/clip/cellranger4.rs
+++ b/src/clip/cellranger4.rs
@@ -7,11 +7,11 @@
 //! - a 5' template-switch-oligo trim, which is an overlap alignment of the TSO
 //!   against the first 91 bases of the read.
 //!
-//! STAR does the 5' alignment with the Opal SIMD library. Here it goes through
-//! [`crate::swalign`], which is required to be bit-identical across
-//! instruction sets, so the clip length cannot depend on the machine.
-
-use crate::swalign::{self, Mode, Scoring};
+//! 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.
 ///
@@ -51,40 +51,6 @@ pub fn poly_tail_3p(seq: &[u8]) -> usize {
     }
 }
 
-/// How much of the read STAR aligns the TSO against (`ClipCR4::opalFillOneSeq`).
-const CR4_TARGET_LEN: usize = 91;
-
-/// Number of 5' bases to trim as the 10x TSO.
-///
-/// STAR aligns the TSO against the first 91 bases of the read in overlap mode,
-/// asking for the score and the position in the target where it ends, then
-/// applies an acceptance gate: a score below 20 is rejected outright, and
-/// scores of exactly 20 or 21 are rejected if they took too long to reach
-/// (more than 26 and 30 bases respectively). A weak alignment that happens to
-/// run a long way is what that gate is there to catch.
-///
-/// The read is padded to 91 bases with `N` when it is shorter, which is why
-/// the scoring scheme has to treat `N` against `N` as neutral rather than as a
-/// mismatch: otherwise the padding would drag every score down.
-///
-/// Both arguments are numeric base codes; the return value is a count of 5'
-/// bases to clip, `0` when the alignment is rejected.
-pub fn tso_clip(read: &[u8], tso: &[u8]) -> usize {
-    if tso.is_empty() {
-        return 0;
-    }
-    let take = read.len().min(CR4_TARGET_LEN);
-    let mut target = Vec::with_capacity(CR4_TARGET_LEN);
-    target.extend_from_slice(&read[..take]);
-    target.resize(CR4_TARGET_LEN, 4); // N padding
-
-    let a = swalign::align(tso, &target, Mode::Ov, &Scoring::CLIP_CR4);
-    let clip = a.target_end as i64 + 1; // 1-based end == number of bases covered
-
-    let reject = a.score < 20 || (a.score == 20 && clip > 26) || (a.score == 21 && clip > 30);
-    if reject { 0 } else { clip as usize }
-}
-
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -102,28 +68,6 @@ mod tests {
             .collect()
     }
 
-    /// The 10x template switch oligo.
-    const TSO: &str = "AAGCAGTGGTATCAACGCAGAGTACATGGG";
-
-    #[test]
-    fn cr4_tso_clip_matches_opal() {
-        // Frozen vector shared with STAR-rs's `cr4_tso_clip_matches_opal`,
-        // which validates the same numbers against Opal itself. A read that
-        // starts with the TSO is clipped by exactly its length.
-        let read = code(&format!("{TSO}ACGTACGTACGTACGTACGTACGTACGTAC"));
-        assert_eq!(tso_clip(&read, &code(TSO)), 30);
-
-        // A read with no TSO is not clipped at all.
-        let read = code("ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT");
-        assert_eq!(tso_clip(&read, &code(TSO)), 0);
-    }
-
-    #[test]
-    fn tso_clip_is_inert_without_an_adapter() {
-        let read = code(&format!("{TSO}ACGTACGT"));
-        assert_eq!(tso_clip(&read, &[]), 0);
-    }
-
     #[test]
     fn cr4_polya_trim_matches_star() {
         // A clean 30-base poly-A tail is trimmed whole. The prefix is
diff --git a/src/clip/mod.rs b/src/clip/mod.rs
index 973fdfb..df64b6e 100644
--- a/src/clip/mod.rs
+++ b/src/clip/mod.rs
@@ -65,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
@@ -184,11 +193,11 @@ pub fn clip_mate(read: &[u8], p: &ClipParams) -> (usize, usize) {
 fn clip_mate_cellranger4(read: &[u8], p: &ClipParams) -> (usize, usize) {
     let len = read.len();
 
-    // 5': fixed clip, then the TSO overlap alignment on what is left.
+    // 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_adapter.is_empty() {
-        c5 += cellranger4::tso_clip(&read[c5..], &p.five_adapter).min(len - c5);
-    }
     if p.five.n_after > 0 && c5 < len {
         c5 += p.five.n_after.min(len - c5);
     }
@@ -227,19 +236,6 @@ mod cr4_wiring_tests {
         }
     }
 
-    #[test]
-    fn cellranger4_clips_the_tso_and_the_polya_tail() {
-        // TSO at the 5' end, a clean poly-A tail at the 3', mappable sequence
-        // in between.
-        let body = "CGTCGTCGTCGTCGTCGTCGTCGTCGTCGT";
-        let read = code(&format!("{TSO}{body}{}", "A".repeat(30)));
-        let (c5, c3) = clip_mate(&read, &cr4_params(TSO));
-        assert_eq!(c5, 30, "the TSO should be clipped from the 5' end");
-        assert_eq!(c3, 30, "the poly-A tail should be clipped from the 3' end");
-        // What survives is exactly the body.
-        assert_eq!(&read[c5..read.len() - c3], code(body).as_slice());
-    }
-
     #[test]
     fn cellranger4_leaves_a_read_without_either_feature_alone() {
         let read = code("CGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGT");
@@ -256,14 +252,13 @@ mod cr4_wiring_tests {
     }
 
     #[test]
-    fn cellranger4_applies_the_fixed_clips_first() {
-        let body = "CGTCGTCGTCGTCGTCGTCGTCGTCGTCGT";
-        let read = code(&format!("{TSO}{body}{}", "A".repeat(30)));
-        let mut p = cr4_params(TSO);
+    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, _) = clip_mate(&read, &p);
-        // The fixed 5 bases come off, then the rest of the TSO is still found.
-        assert_eq!(c5, 30);
+        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");
     }
 }
 
diff --git a/src/lib.rs b/src/lib.rs
index b73fb01..6fce173 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -40,7 +40,6 @@ pub mod rng;
 pub mod signal;
 pub mod solo;
 pub mod stats;
-pub mod swalign;
 pub mod wasp;
 
 use log::info;
diff --git a/src/swalign/backend.rs b/src/swalign/backend.rs
deleted file mode 100644
index 48ddf49..0000000
--- a/src/swalign/backend.rs
+++ /dev/null
@@ -1,204 +0,0 @@
-//! Backend selection, and the differential harness that keeps backends honest.
-//!
-//! The contract every backend signs: given the same inputs it returns exactly
-//! what [`super::scalar::align`] returns. Not "within rounding", not "the same
-//! score with a different end position" — the same [`Alignment`].
-//!
-//! [`differential_check`] is what enforces that. It is written here rather than
-//! in a test module so a backend can be checked from a test, a benchmark or a
-//! debug session without duplicating the generator.
-
-use super::{Alignment, Mode, Scoring, scalar};
-
-/// Which implementation actually ran.
-///
-/// Exposed so tests can assert that a machine with the hardware really used it,
-/// rather than silently falling back and reporting a pass that proves nothing.
-#[derive(Debug, Clone, Copy, PartialEq, Eq)]
-pub enum Backend {
-    /// The portable reference. Always available, and the definition of the
-    /// result.
-    Scalar,
-    /// aarch64 NEON, computing one anti-diagonal at a time.
-    #[cfg(target_arch = "aarch64")]
-    Neon,
-}
-
-impl Backend {
-    /// The best backend this machine can run.
-    pub fn detect() -> Self {
-        #[cfg(target_arch = "aarch64")]
-        if super::neon::is_available() {
-            return Backend::Neon;
-        }
-        Backend::Scalar
-    }
-
-    /// Run this backend.
-    pub fn align(self, query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
-        match self {
-            Backend::Scalar => scalar::align(query, target, mode, scoring),
-            #[cfg(target_arch = "aarch64")]
-            Backend::Neon => super::neon::align(query, target, mode, scoring),
-        }
-    }
-}
-
-/// A deterministic sequence generator for the differential harness.
-///
-/// Seeded from the in-tree splitmix64 so the cases are identical on every
-/// machine and every run: a backend that fails does so reproducibly, on a case
-/// the report can name.
-struct SeqGen(crate::rng::SplitMix64);
-
-impl SeqGen {
-    fn new(seed: u64) -> Self {
-        Self(crate::rng::SplitMix64::seed(seed))
-    }
-
-    /// A sequence of `len` base codes. `n_rate` out of 16 bases are `N`, so the
-    /// generator covers the `N`-heavy inputs STAR's padding produces as well as
-    /// clean ones.
-    fn seq(&mut self, len: usize, n_rate: u64) -> Vec<u8> {
-        (0..len)
-            .map(|_| {
-                let r = self.0.next_u64();
-                if r % 16 < n_rate {
-                    4
-                } else {
-                    (r >> 8) as u8 % 4
-                }
-            })
-            .collect()
-    }
-}
-
-/// Compare `backend` against the scalar reference over a spread of inputs.
-///
-/// Returns `Err` with a description of the first disagreement, naming the case
-/// so it can be reproduced. `Ok(n)` reports how many cases were checked.
-///
-/// The spread is deliberately awkward: empty and length-1 sequences, queries
-/// longer than targets, all-`N` inputs, and long runs that push scores far
-/// enough to saturate a narrow lane. Saturation is where a SIMD backend is most
-/// likely to diverge, so it is not left to chance.
-pub fn differential_check(backend: Backend, seed: u64) -> Result<usize, String> {
-    let scoring = Scoring::CLIP_CR4;
-    let mut rng = SeqGen::new(seed);
-    let mut checked = 0usize;
-
-    // Lengths chosen around the boundaries a vectorised kernel cares about:
-    // zero, one, just under and just over a typical lane count, and the 91
-    // STAR actually uses.
-    const LENS: &[usize] = &[0, 1, 2, 7, 8, 9, 15, 16, 17, 30, 31, 33, 64, 91, 128];
-
-    for &ql in LENS {
-        for &tl in LENS {
-            for &n_rate in &[0u64, 1, 8, 16] {
-                for mode in [Mode::Nw, Mode::Hw, Mode::Ov, Mode::Sw] {
-                    let q = rng.seq(ql, n_rate);
-                    let t = rng.seq(tl, n_rate);
-                    let want = scalar::align(&q, &t, mode, &scoring);
-                    let got = backend.align(&q, &t, mode, &scoring);
-                    if got != want {
-                        return Err(format!(
-                            "{backend:?} disagrees with scalar on {mode:?}, \
-                             |q|={ql} |t|={tl} n_rate={n_rate}/16: \
-                             scalar {want:?}, backend {got:?}\n  query  {q:?}\n  target {t:?}"
-                        ));
-                    }
-                    checked += 1;
-                }
-            }
-        }
-    }
-
-    // Saturation: a long exact match scores one per base, which overflows an
-    // 8-bit lane well before it overflows the scalar's i32. A backend that
-    // escalates lane width incorrectly fails here and nowhere else.
-    for &len in &[200usize, 400, 1000] {
-        let q = rng.seq(len, 0);
-        for mode in [Mode::Nw, Mode::Hw, Mode::Ov, Mode::Sw] {
-            let want = scalar::align(&q, &q, mode, &scoring);
-            let got = backend.align(&q, &q, mode, &scoring);
-            if got != want {
-                return Err(format!(
-                    "{backend:?} disagrees with scalar under saturation, \
-                     {mode:?}, len={len}: scalar {want:?}, backend {got:?}"
-                ));
-            }
-            checked += 1;
-        }
-    }
-
-    Ok(checked)
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    /// Every backend this machine can run, not just the one `detect()` picks.
-    ///
-    /// On x86 that matters: a machine with AVX2 would otherwise never exercise
-    /// SSE2, and the baseline is what runs on everything older.
-    fn available_backends() -> Vec<Backend> {
-        let mut v = vec![Backend::Scalar];
-        #[cfg(target_arch = "aarch64")]
-        if super::super::neon::is_available() {
-            v.push(Backend::Neon);
-        }
-        v
-    }
-
-    #[test]
-    fn every_available_backend_agrees_with_scalar() {
-        // Not just the detected one: on a machine with AVX2, SSE2 would
-        // otherwise go unchecked despite being what older hardware runs.
-        for backend in available_backends() {
-            if let Err(e) = differential_check(backend, 0x0BAD_5EED) {
-                panic!("{e}");
-            }
-        }
-    }
-
-    #[test]
-    fn scalar_is_consistent_with_itself() {
-        // Tautological for the scalar backend, but it proves the harness runs,
-        // covers every mode and length pair, and reaches the saturation cases.
-        // A backend added later inherits a harness that is known to work.
-        let n = differential_check(Backend::Scalar, 0x5EED).expect("scalar vs scalar");
-        assert!(n > 900, "harness covered only {n} cases");
-    }
-
-    #[test]
-    fn the_detected_backend_agrees_with_scalar_on_this_machine() {
-        // The test that matters once SIMD backends exist: whatever this
-        // machine selects must match the reference on this machine.
-        let backend = Backend::detect();
-
-        // A pass proves nothing if the machine quietly fell back to the
-        // reference, so assert that the hardware backend really was selected.
-        #[cfg(target_arch = "aarch64")]
-        assert_eq!(
-            backend,
-            Backend::Neon,
-            "aarch64 must select NEON; a silent fallback would make this test vacuous"
-        );
-
-        if let Err(e) = differential_check(backend, 0x00C0_FFEE) {
-            panic!("{e}");
-        }
-    }
-
-    #[test]
-    fn generator_is_reproducible() {
-        // The harness is only useful if a failure can be reproduced, which
-        // needs the sequences to be identical run to run.
-        let a = SeqGen::new(7).seq(64, 4);
-        let b = SeqGen::new(7).seq(64, 4);
-        assert_eq!(a, b);
-        assert!(a.contains(&4), "n_rate 4/16 should produce Ns");
-        assert!(a.iter().any(|&b| b < 4), "and also real bases");
-    }
-}
diff --git a/src/swalign/mod.rs b/src/swalign/mod.rs
deleted file mode 100644
index ed305e9..0000000
--- a/src/swalign/mod.rs
+++ /dev/null
@@ -1,150 +0,0 @@
-//! Deterministic Smith-Waterman-family alignment with affine gaps.
-//!
-//! STAR performs exactly one alignment of this shape: the 5' TSO clip of
-//! `--clipAdapterType CellRanger4`, for which it links the Opal C/C++ SIMD
-//! library (`OPAL_MODE_OV` + `OPAL_SEARCH_SCORE_END`). This module provides
-//! the same capability in-tree, with no new dependency.
-//!
-//! # Determinism is a correctness property here
-//!
-//! Opal's overflow handling groups database sequences into buckets sized by the
-//! SIMD vector width, so when 8-bit lanes saturate, the grouping — and with it
-//! the recompute path — depends on which instruction set is available. Sixteen
-//! SSE lanes and thirty-two AVX2 lanes bucket differently, and simde emulates
-//! AVX2 on ARM, so the same input can take a different path on a different
-//! machine.
-//!
-//! That is not acceptable for an aligner whose output is supposed to be
-//! reproducible. Here the rule is: **the vector width is never observable.**
-//! [`scalar`] defines the result; every other backend must agree with it
-//! bit-for-bit, including under saturation, on empty inputs, on `N`, and on
-//! end-position ties. The differential test in [`tests`] is what enforces it.
-//!
-//! # Coordinates and conventions
-//!
-//! Sequences are numeric base codes, `0..=3` for ACGT and `4` for `N`, the same
-//! encoding the rest of the aligner uses. Scores are `i32`; the caller supplies
-//! the scoring scheme.
-
-pub mod backend;
-#[cfg(target_arch = "aarch64")]
-pub mod neon;
-pub mod scalar;
-
-pub use backend::Backend;
-
-/// How the ends of the two sequences are treated.
-#[derive(Debug, Clone, Copy, PartialEq, Eq)]
-pub enum Mode {
-    /// Global: both sequences must be consumed end to end (Needleman-Wunsch).
-    Nw,
-    /// Semi-global on the target: the query must be consumed entirely, the
-    /// target may extend past it on both sides.
-    Hw,
-    /// Overlap: gaps at the start of either sequence and at the end of either
-    /// sequence are free. This is the mode STAR's TSO clip uses.
-    Ov,
-    /// Local (Smith-Waterman): the best-scoring subalignment.
-    Sw,
-}
-
-/// Affine-gap scoring.
-///
-/// `gap_open` and `gap_extend` are the penalties *subtracted*, so both are
-/// given as positive numbers. STAR's ClipCR4 uses `match_score = 1`,
-/// `mismatch = -2`, `gap_open = gap_extend = 2`, and scores `N` against `N` as
-/// zero rather than as a mismatch.
-#[derive(Debug, Clone, Copy)]
-pub struct Scoring {
-    /// Added when two bases are equal.
-    pub match_score: i32,
-    /// Added when two bases differ (negative).
-    pub mismatch: i32,
-    /// Subtracted to open a gap.
-    pub gap_open: i32,
-    /// Subtracted for each base a gap is extended by.
-    pub gap_extend: i32,
-    /// Added when both positions are `N`. Opal treats this pairing as neutral
-    /// rather than as a mismatch, and STAR relies on that when it pads the
-    /// target with `N`.
-    pub n_vs_n: i32,
-}
-
-impl Scoring {
-    /// The scoring STAR uses for the CellRanger4 TSO clip
-    /// (`ClipCR4.cpp`: match +1, mismatch -2, gaps 2, `N`/`N` neutral).
-    pub const CLIP_CR4: Self = Self {
-        match_score: 1,
-        mismatch: -2,
-        gap_open: 2,
-        gap_extend: 2,
-        n_vs_n: 0,
-    };
-
-    /// Score one aligned pair of base codes.
-    #[inline]
-    pub fn pair(&self, q: u8, t: u8) -> i32 {
-        if q == 4 && t == 4 {
-            self.n_vs_n
-        } else if q == t {
-            self.match_score
-        } else {
-            self.mismatch
-        }
-    }
-}
-
-/// The outcome of one alignment.
-#[derive(Debug, Clone, Copy, PartialEq, Eq)]
-pub struct Alignment {
-    /// Best score under the chosen mode.
-    pub score: i32,
-    /// Zero-based position in the target where that score is reached.
-    ///
-    /// Ties are broken towards the **earlier** column, except that a score
-    /// reached in the final column wins outright: that is what
-    /// `OPAL_SEARCH_SCORE_END` means, and STAR's clip length depends on it.
-    pub target_end: usize,
-}
-
-/// Align `query` against `target` under `mode`.
-///
-/// Dispatches to the fastest backend available for this machine. Every backend
-/// is required to return exactly what [`scalar::align`] would, so the choice is
-/// invisible in the output.
-pub fn align(query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
-    Backend::detect().align(query, target, mode, scoring)
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    #[test]
-    fn empty_inputs_score_zero_in_every_mode() {
-        for mode in [Mode::Nw, Mode::Hw, Mode::Ov, Mode::Sw] {
-            let a = align(&[], &[], mode, &Scoring::CLIP_CR4);
-            assert_eq!(a.score, 0, "{mode:?} on two empty sequences");
-            let a = align(&[0, 1, 2], &[], mode, &Scoring::CLIP_CR4);
-            assert_eq!(a.score, 0, "{mode:?} on an empty target");
-        }
-    }
-
-    #[test]
-    fn n_against_n_is_neutral_not_a_mismatch() {
-        let s = Scoring::CLIP_CR4;
-        assert_eq!(s.pair(4, 4), 0);
-        assert_eq!(s.pair(4, 0), -2);
-        assert_eq!(s.pair(0, 4), -2);
-        assert_eq!(s.pair(0, 0), 1);
-        assert_eq!(s.pair(0, 1), -2);
-    }
-
-    #[test]
-    fn identical_sequences_score_one_per_base() {
-        let seq = [0u8, 1, 2, 3, 0, 1, 2, 3];
-        let a = align(&seq, &seq, Mode::Ov, &Scoring::CLIP_CR4);
-        assert_eq!(a.score, seq.len() as i32);
-        assert_eq!(a.target_end, seq.len() - 1);
-    }
-}
diff --git a/src/swalign/neon.rs b/src/swalign/neon.rs
deleted file mode 100644
index 1121b95..0000000
--- a/src/swalign/neon.rs
+++ /dev/null
@@ -1,238 +0,0 @@
-//! aarch64 NEON backend.
-//!
-//! # Why anti-diagonals rather than Farrar stripes
-//!
-//! The striped layout is faster, but its lazy-F correction and its
-//! stripe-relative indexing make the *end position* awkward to extract, and the
-//! end position is exactly what STAR's clip length is built on. A backend that
-//! is quicker but disagrees with [`super::scalar`] about where an alignment
-//! ends is worthless here.
-//!
-//! Cells on one anti-diagonal `d = r + c` are mutually independent: `H(r,c)`
-//! depends on `E(r,c-1)` and `F(r-1,c)`, both on `d-1`, and on `H(r-1,c-1)` on
-//! `d-2`. So a whole anti-diagonal can be computed in parallel with no
-//! correction pass and no reordering, which makes agreement with the scalar
-//! path a property of the layout rather than something to be tested for and
-//! hoped about.
-//!
-//! The cost is strided access. For the 30×91 matrix STAR's TSO clip actually
-//! uses, that is irrelevant.
-
-use std::arch::aarch64::{vaddq_s32, vdupq_n_s32, vld1q_s32, vmaxq_s32, vst1q_s32, vsubq_s32};
-
-use super::{Alignment, Mode, Scoring};
-
-/// Matches `scalar::NEG`: far below any real score, far enough from `i32::MIN`
-/// that subtracting a gap penalty cannot wrap.
-const NEG: i32 = i32::MIN / 4;
-
-/// Lanes per NEON vector at 32-bit width.
-const LANES: usize = 4;
-
-/// Align `query` against `target` using NEON.
-///
-/// # Safety
-///
-/// The caller must have established that NEON is available. On aarch64 it is
-/// architecturally guaranteed, so [`is_available`] is a constant.
-pub fn align(query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
-    let ql = query.len();
-    let tl = target.len();
-    if ql == 0 || tl == 0 {
-        return Alignment {
-            score: 0,
-            target_end: 0,
-        };
-    }
-
-    let (free_query_start, free_target_start) = match mode {
-        Mode::Nw => (false, false),
-        Mode::Hw => (false, true),
-        Mode::Ov | Mode::Sw => (true, true),
-    };
-
-    // Score of aligning the first `r+1` query bases against nothing, i.e. the
-    // cell one column to the left of column 0.
-    let left_h = |r: usize| -> i32 {
-        if free_query_start {
-            0
-        } else {
-            -(scoring.gap_open + scoring.gap_extend * r as i32)
-        }
-    };
-    // Score of aligning the first `c+1` target bases against nothing, i.e. the
-    // cell one row above row 0.
-    let top_h = |c: i64| -> i32 {
-        // A column before the first scores 0 either way, so it folds into the
-        // free-start case.
-        if free_target_start || c < 0 {
-            0
-        } else {
-            -(scoring.gap_open + scoring.gap_extend * c as i32)
-        }
-    };
-
-    // Rows of the two previous anti-diagonals, indexed by `r`.
-    let mut h1 = vec![NEG; ql]; // H on d-1
-    let mut e1 = vec![NEG; ql]; // E on d-1
-    let mut f1 = vec![NEG; ql]; // F on d-1
-    let mut h2 = vec![NEG; ql]; // H on d-2
-    let mut h0 = vec![NEG; ql];
-    let mut e0 = vec![NEG; ql];
-    let mut f0 = vec![NEG; ql];
-
-    let mut best_last_row = NEG;
-    let mut best_last_row_col = 0usize;
-    let mut last_col_max = NEG;
-    let mut best_anywhere = 0i32;
-    let mut best_anywhere_col = 0usize;
-    // The scalar sweeps column-major and keeps the first cell to reach a new
-    // maximum, so among equal scores the smallest `(c, r)` wins. Visiting
-    // anti-diagonals changes the order, so the key has to be compared
-    // explicitly rather than inferred from arrival.
-    let mut best_anywhere_row = 0usize;
-
-    let go = scoring.gap_open;
-    let ge = scoring.gap_extend;
-
-    for d in 0..(ql + tl - 1) {
-        let r_lo = d.saturating_sub(tl - 1);
-        let r_hi = d.min(ql - 1);
-
-        // SAFETY: every load and store below is a plain lane-wise op on stack
-        // scalars gathered by index; no pointer arithmetic escapes the slices,
-        // which are all length `ql` and indexed within `r_lo..=r_hi`.
-        unsafe {
-            let vgo = vdupq_n_s32(go);
-            let vge = vdupq_n_s32(ge);
-
-            let mut row = r_lo;
-            while row <= r_hi {
-                let lanes = LANES.min(r_hi - row + 1);
-
-                // Gather the four predecessor terms for lanes r..r+n.
-                let mut e_prev_h = [NEG; LANES]; // H(r, c-1)
-                let mut e_prev_e = [NEG; LANES]; // E(r, c-1)
-                let mut f_prev_h = [NEG; LANES]; // H(r-1, c)
-                let mut f_prev_f = [NEG; LANES]; // F(r-1, c)
-                let mut diag = [NEG; LANES]; // H(r-1, c-1)
-                let mut sub = [0i32; LANES]; // score(q[r], t[c])
-
-                for k in 0..lanes {
-                    let rr = row + k;
-                    let col = d - rr;
-                    // (rr, c-1) lives on d-1 at row rr.
-                    if col == 0 {
-                        e_prev_h[k] = left_h(rr);
-                        e_prev_e[k] = NEG;
-                    } else {
-                        e_prev_h[k] = h1[rr];
-                        e_prev_e[k] = e1[rr];
-                    }
-                    // (rr-1, c) lives on d-1 at row rr-1.
-                    if rr == 0 {
-                        f_prev_h[k] = top_h(col as i64);
-                        f_prev_f[k] = NEG;
-                    } else {
-                        f_prev_h[k] = h1[rr - 1];
-                        f_prev_f[k] = f1[rr - 1];
-                    }
-                    // (rr-1, c-1) lives on d-2 at row rr-1.
-                    diag[k] = if rr == 0 {
-                        top_h(col as i64 - 1)
-                    } else if col == 0 {
-                        left_h(rr - 1)
-                    } else {
-                        h2[rr - 1]
-                    };
-                    sub[k] = scoring.pair(query[rr], target[col]);
-                }
-
-                let ins = vmaxq_s32(
-                    vsubq_s32(vld1q_s32(e_prev_h.as_ptr()), vgo),
-                    vsubq_s32(vld1q_s32(e_prev_e.as_ptr()), vge),
-                );
-                let del = vmaxq_s32(
-                    vsubq_s32(vld1q_s32(f_prev_h.as_ptr()), vgo),
-                    vsubq_s32(vld1q_s32(f_prev_f.as_ptr()), vge),
-                );
-                let diagv = vaddq_s32(vld1q_s32(diag.as_ptr()), vld1q_s32(sub.as_ptr()));
-                let mut cur = vmaxq_s32(vmaxq_s32(ins, del), diagv);
-                if mode == Mode::Sw {
-                    cur = vmaxq_s32(cur, vdupq_n_s32(0));
-                }
-
-                let mut hs = [0i32; LANES];
-                let mut es = [0i32; LANES];
-                let mut fs = [0i32; LANES];
-                vst1q_s32(hs.as_mut_ptr(), cur);
-                vst1q_s32(es.as_mut_ptr(), ins);
-                vst1q_s32(fs.as_mut_ptr(), del);
-
-                for k in 0..lanes {
-                    let rr = row + k;
-                    h0[rr] = hs[k];
-                    e0[rr] = es[k];
-                    f0[rr] = fs[k];
-                }
-                row += lanes;
-            }
-        }
-
-        // Bookkeeping, in increasing column order so ties go to the earlier
-        // column exactly as the scalar's column loop does.
-        for (rr, &cell) in h0.iter().enumerate().take(r_hi + 1).skip(r_lo) {
-            let col = d - rr;
-            if mode == Mode::Sw
-                && (cell > best_anywhere
-                    || (cell == best_anywhere
-                        && (col, rr) < (best_anywhere_col, best_anywhere_row)))
-            {
-                best_anywhere = cell;
-                best_anywhere_col = col;
-                best_anywhere_row = rr;
-            }
-            if rr == ql - 1 && cell > best_last_row {
-                best_last_row = cell;
-                best_last_row_col = col;
-            }
-            if col == tl - 1 {
-                last_col_max = last_col_max.max(cell);
-            }
-        }
-
-        std::mem::swap(&mut h2, &mut h1);
-        std::mem::swap(&mut h1, &mut h0);
-        std::mem::swap(&mut e1, &mut e0);
-        std::mem::swap(&mut f1, &mut f0);
-    }
-
-    match mode {
-        Mode::Sw => Alignment {
-            score: best_anywhere,
-            target_end: best_anywhere_col,
-        },
-        Mode::Nw => Alignment {
-            score: h1[ql - 1],
-            target_end: tl - 1,
-        },
-        Mode::Hw => Alignment {
-            score: best_last_row,
-            target_end: best_last_row_col,
-        },
-        Mode::Ov => {
-            let score = last_col_max.max(best_last_row);
-            let target_end = if last_col_max >= best_last_row {
-                tl - 1
-            } else {
-                best_last_row_col
-            };
-            Alignment { score, target_end }
-        }
-    }
-}
-
-/// NEON is architecturally guaranteed on aarch64.
-pub const fn is_available() -> bool {
-    true
-}
diff --git a/src/swalign/scalar.rs b/src/swalign/scalar.rs
deleted file mode 100644
index 53930e6..0000000
--- a/src/swalign/scalar.rs
+++ /dev/null
@@ -1,205 +0,0 @@
-//! Portable reference implementation.
-//!
-//! This is the definition of the result. It is written for clarity rather than
-//! speed: the SIMD backends must agree with it bit-for-bit, so it needs to be
-//! obviously correct more than it needs to be fast.
-
-use super::{Alignment, Mode, Scoring};
-
-/// Sentinel for "unreachable". Far below any real score, and far enough from
-/// `i32::MIN` that subtracting a gap penalty cannot wrap.
-const NEG: i32 = i32::MIN / 4;
-
-/// Align `query` against `target`, filling the DP column by column.
-///
-/// Only two columns are ever live, so the working set is `O(|query|)` rather
-/// than the full matrix. That is also what makes the striped SIMD form a
-/// drop-in replacement later: it computes the same columns in the same order.
-pub fn align(query: &[u8], target: &[u8], mode: Mode, scoring: &Scoring) -> Alignment {
-    let ql = query.len();
-    let tl = target.len();
-    if ql == 0 || tl == 0 {
-        return Alignment {
-            score: 0,
-            target_end: 0,
-        };
-    }
-
-    // Whether a gap before the start of each sequence is free.
-    let (free_query_start, free_target_start) = match mode {
-        Mode::Nw => (false, false),
-        Mode::Hw => (false, true),
-        Mode::Ov | Mode::Sw => (true, true),
-    };
-
-    // Column 0 of the previous iteration: H is the score of aligning the first
-    // `r+1` query bases against nothing.
-    let mut prev_h = vec![0i32; ql];
-    let mut prev_e = vec![NEG; ql];
-    if !free_query_start {
-        for (r, h) in prev_h.iter_mut().enumerate() {
-            *h = -(scoring.gap_open + scoring.gap_extend * r as i32);
-        }
-    }
-
-    // Best score seen on the final query row, and the column that first
-    // achieved it. Ties go to the earlier column.
-    let mut best_last_row = NEG;
-    let mut best_last_row_col = 0usize;
-    // Best score anywhere in the most recent column.
-    let mut last_col_max = NEG;
-    // Best score anywhere in the matrix, for local mode.
-    let mut best_anywhere = 0i32;
-    let mut best_anywhere_col = 0usize;
-
-    for (c, &tc) in target.iter().enumerate() {
-        // Top of the column: aligning the first `c+1` target bases against
-        // nothing.
-        let mut up_h = if free_target_start {
-            0
-        } else {
-            -(scoring.gap_open + scoring.gap_extend * c as i32)
-        };
-        // The diagonal predecessor, i.e. the cell up and to the left. Column 0
-        // has no predecessor and scores 0 either way, so it falls out of the
-        // free-start case.
-        let mut diag_h = if free_target_start || c == 0 {
-            0
-        } else {
-            -(scoring.gap_open + scoring.gap_extend * (c as i32 - 1))
-        };
-        let mut up_f = NEG;
-        let mut col_max = NEG;
-        let mut h = NEG;
-
-        for r in 0..ql {
-            // Gap in the query (moving right): open from H, or extend E.
-            let e = (prev_h[r] - scoring.gap_open).max(prev_e[r] - scoring.gap_extend);
-            // Gap in the target (moving down): open from H, or extend F.
-            let f = (up_h - scoring.gap_open).max(up_f - scoring.gap_extend);
-            // Match or mismatch on the diagonal.
-            let d = diag_h + scoring.pair(query[r], tc);
-
-            h = e.max(f).max(d);
-            if mode == Mode::Sw {
-                // Local alignment never carries a negative prefix forward.
-                h = h.max(0);
-                if h > best_anywhere {
-                    best_anywhere = h;
-                    best_anywhere_col = c;
-                }
-            }
-            if h > col_max {
-                col_max = h;
-            }
-
-            up_f = f;
-            up_h = h;
-            diag_h = prev_h[r];
-            prev_e[r] = e;
-            prev_h[r] = h;
-        }
-
-        // `h` now holds the last query row for this column.
-        if h > best_last_row {
-            best_last_row = h;
-            best_last_row_col = c;
-        }
-        last_col_max = col_max;
-    }
-
-    match mode {
-        Mode::Sw => Alignment {
-            score: best_anywhere,
-            target_end: best_anywhere_col,
-        },
-        Mode::Nw => Alignment {
-            // Global: the corner cell, which is the last row of the last
-            // column.
-            score: prev_h[ql - 1],
-            target_end: tl - 1,
-        },
-        Mode::Hw => Alignment {
-            // Semi-global on the target: the query must be consumed, so the
-            // answer lives on the last query row.
-            score: best_last_row,
-            target_end: best_last_row_col,
-        },
-        Mode::Ov => {
-            // Overlap: either the query ran out (best on the last row) or the
-            // target ran out (best in the last column). A score reached in the
-            // final column wins the tie, which is what `OPAL_SEARCH_SCORE_END`
-            // specifies and what STAR's clip length depends on.
-            let score = last_col_max.max(best_last_row);
-            let target_end = if last_col_max >= best_last_row {
-                tl - 1
-            } else {
-                best_last_row_col
-            };
-            Alignment { score, target_end }
-        }
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use crate::swalign::Scoring;
-
-    const S: Scoring = Scoring::CLIP_CR4;
-
-    #[test]
-    fn global_mode_pays_for_every_unaligned_base() {
-        // Query is a prefix of the target; NW must pay to skip the tail.
-        let q = [0u8, 1, 2];
-        let t = [0u8, 1, 2, 3, 3];
-        let nw = align(&q, &t, Mode::Nw, &S);
-        let ov = align(&q, &t, Mode::Ov, &S);
-        assert_eq!(ov.score, 3, "overlap: the tail is free");
-        assert!(
-            nw.score < ov.score,
-            "global should pay for the tail, got {} vs {}",
-            nw.score,
-            ov.score
-        );
-    }
-
-    #[test]
-    fn local_mode_ignores_flanking_mismatch() {
-        // A clean 4-base core buried in mismatching flanks.
-        let q = [3u8, 3, 0, 1, 2, 3, 3, 3];
-        let t = [1u8, 1, 0, 1, 2, 3, 1, 1];
-        let sw = align(&q, &t, Mode::Sw, &S);
-        assert!(
-            sw.score >= 4,
-            "local should find the shared core, got {}",
-            sw.score
-        );
-    }
-
-    #[test]
-    fn local_score_is_never_negative() {
-        let q = [0u8, 0, 0, 0];
-        let t = [3u8, 3, 3, 3];
-        assert_eq!(align(&q, &t, Mode::Sw, &S).score, 0);
-    }
-
-    #[test]
-    fn overlap_prefers_the_final_column_on_a_tie() {
-        // A query that matches equally well at two positions: the one that
-        // runs to the end of the target must win, because that is the tie-break
-        // STAR's clip length is built on.
-        let q = [0u8, 1];
-        let t = [0u8, 1, 4, 4, 0, 1];
-        let a = align(&q, &t, Mode::Ov, &S);
-        assert_eq!(a.target_end, t.len() - 1);
-    }
-
-    #[test]
-    fn affine_gaps_cost_less_than_repeated_opens() {
-        // One 3-base gap should beat three separate 1-base gaps.
-        let one_long = S.gap_open + S.gap_extend * 3;
-        let three_short = 3 * (S.gap_open + S.gap_extend);
-        assert!(one_long < three_short);
-    }
-}

From 260770d98724a05ff649130c1f6c3592806513dc Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Wed, 29 Jul 2026 11:50:28 +0200
Subject: [PATCH 8/8] test(clip): end-to-end test for --clipAdapterType
 CellRanger4

The PR description named this test and it did not exist, which CONTRIBUTING.md
rules out. Writing it turned out to be worth more than the description was: a
100-base read with a 40-base poly-A tail does not merely align differently
untrimmed, it fails the mismatch filters and is lost entirely. Trimmed, it
aligns 60M40S. That is what the flag being inert actually cost.
---
 tests/alignment_features.rs | 83 +++++++++++++++++++++++++++++++++++++
 1 file changed, 83 insertions(+)

diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs
index 2db7d39..ca9dc53 100644
--- a/tests/alignment_features.rs
+++ b/tests/alignment_features.rs
@@ -2024,3 +2024,86 @@ fn test_wasp_samtag() {
         "all 10 unique reads overlapping the het SNV should pass WASP (vW:i:1)"
     );
 }
+
+// ---------------------------------------------------------------------------
+// --clipAdapterType CellRanger4
+// ---------------------------------------------------------------------------
+
+/// End to end: a read carrying a poly-A tail aligns over its non-A prefix only.
+///
+/// Before this, `--clipAdapterType CellRanger4` parsed and then did nothing, so
+/// the tail was carried into the alignment as soft-clipped or mismatching
+/// bases. The check is that the flag changes the CIGAR at all, and changes it
+/// the way a trim should.
+#[test]
+fn test_clip_adapter_type_cellranger4() {
+    let tmpdir = TempDir::new().unwrap();
+    let genome = build_genome();
+    let fasta = write_fasta(&tmpdir, &genome);
+    let genome_dir = tmpdir.path().join("genome");
+    build_index(&fasta, &genome_dir, "7", None);
+
+    // 60 genomic bases followed by a 40-base poly-A tail the genome does not
+    // have. The tail is long enough to clear the trim's score-20 floor.
+    let start = 3000usize;
+    let mut read = genome[start..start + 60].to_vec();
+    read.extend(std::iter::repeat_n(b'A', 40));
+
+    let fastq_path = tmpdir.path().join("reads.fq");
+    {
+        let mut f = fs::File::create(&fastq_path).unwrap();
+        writeln!(f, "@polya").unwrap();
+        f.write_all(&read).unwrap();
+        writeln!(f).unwrap();
+        writeln!(f, "+").unwrap();
+        writeln!(f, "{}", "I".repeat(read.len())).unwrap();
+    }
+
+    let run = |extra: &[&str], out_name: &str| -> Vec<String> {
+        let output_dir = tmpdir.path().join(out_name);
+        fs::create_dir_all(&output_dir).unwrap();
+        let prefix = format!("{}/", output_dir.display());
+        let mut args = vec![
+            "--runMode",
+            "alignReads",
+            "--genomeDir",
+            genome_dir.to_str().unwrap(),
+            "--readFilesIn",
+            fastq_path.to_str().unwrap(),
+            "--outFileNamePrefix",
+            &prefix,
+        ];
+        args.extend_from_slice(extra);
+        cargo_bin_cmd!("rustar-aligner")
+            .args(&args)
+            .assert()
+            .success();
+        fs::read_to_string(output_dir.join("Aligned.out.sam"))
+            .unwrap()
+            .lines()
+            .filter(|l| !l.starts_with('@'))
+            .map(ToString::to_string)
+            .collect()
+    };
+
+    let plain = run(&[], "out_plain");
+    let clipped = run(&["--clipAdapterType", "CellRanger4"], "out_cr4");
+
+    // Untrimmed, 40 of 100 bases mismatch and the read fails the mismatch
+    // filters outright. That is the cost of the flag having been inert.
+    assert!(
+        plain.is_empty(),
+        "without trimming the poly-A read should not pass the mismatch filters, got {plain:?}"
+    );
+
+    assert_eq!(clipped.len(), 1, "trimmed, the read aligns");
+    let cigar = clipped[0].split('\t').nth(5).unwrap();
+    assert!(
+        cigar.ends_with("40S"),
+        "the 40-base poly-A tail should be clipped, got {cigar}"
+    );
+    assert!(
+        cigar.starts_with("60M"),
+        "the genomic prefix should still align in full, got {cigar}"
+    );
+}