clip: CellRanger4 3' poly-A trim - #148
Open
BenjaminDEMAILLE wants to merge 7 commits into
Open
Conversation
…uilt on it
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>
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>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
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>
… case 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>
… 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>
Collaborator
|
This is too much to be handled within rustar. This is why Star imports opal. I have already written a tool/library that does deterministic sw with full simd and rustar integration in rust as an opal replacement. This keeps the complexity of a whole different aligner in its own real of tests and verification. I mentioned this in the scverse chat as well as in the PR handling CR4 and that I would be making that tool. James |
Per review on scverse#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>
BenjaminDEMAILLE
changed the title
Collaborator
|
@BenjaminDEMAILLE you need to be the one taking responsibility and reading things, not claude. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
CellRanger4 3' poly-A trim.
What changed
src/clip/cellranger4.rs, ported from STAR'sClipCR4.cppandClipMate_clipChunk.cpp:poly_tail_3pscores in from the 3' end,+1perAand-2per non-A, keeps the longest prefix of that walk still clearing a 70% density threshold, gives up once the score falls more than 27 behind, and returns nothing unless the kept score reached 20. Wired intoclip_mate_cellranger4and reached by--clipAdapterType CellRanger4.The 5' TSO clip is not in this PR. It needs a Smith-Waterman aligner, which per @Psy-Fer's review belongs in a separate library rather than inside rustar;
src/swalignwas removed from this branch.--clipAdapterType CellRanger4therefore performs the poly-A trim and logs a warning that the TSO clip is not yet applied, rather than silently doing half the job.Why
--clipAdapterType CellRanger4parsed and validated, then did nothing:src/clip/mod.rsdeclared the mode out of scope, so a user who asked for it got no clipping and the run reported success.One behaviour worth pointing at
The poly-A scan does not stop at the tail boundary, so A-rich sequence upstream legitimately extends the trim:
ACGT×5 followed by 30As trims 34, not 30. That is STAR's behaviour and reads as an off-by-one;poly_tail_scan_does_not_stop_at_the_tail_boundarypins it.Verification
cr4_polya_trim_matches_star— the frozen vectorpoly_tail_scan_does_not_stop_at_the_tail_boundary— the behaviour abovetest_clip_adapter_type_cellranger4— end to end through the CLI-D warnings,cargo fmt --check, MSRV 1.89Output-neutral unless
--clipAdapterType CellRanger4is passed.🤖 Generated with Claude Code