Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 37 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,43 @@ Sections commonly used: Features, Bug fixes, Other changes.

### Features

- **`--clipAdapterType CellRanger4` now matches STAR exactly**, both halves of the
clip. The 5' TSO trim is an overlap alignment against the first 91 bases of the
read, replicating `ClipCR4`'s Opal call (ACGTN alphabet; match +1, mismatch -2,
any-vs-N -2, N-vs-N 0; gap open = gap extend = 2; overlap mode with end
tracking), followed by STAR's acceptance gate
`S<20 || (S==20 && L>26) || (S==21 && L>30)`. The 3' trim is STAR's scored
`polyTail3p` scan, replacing a "trailing run of A >= 8" approximation that was
wrong in both directions (a read ending in 10 A's was trimmed by 10 where STAR
trims 0; `A*15 + C + A*15` was trimmed by 15 where STAR trims 31). The poly-A
transcription is taken from Benjamin Demaille's #148 and moved onto the solo
path, which is where `CellRanger4` is actually reached; #148's `clip_mate`
wiring is not used and that PR is superseded.

The overlap alignment comes from [`hyalite`](https://crates.io/crates/hyalite)
0.2 — a new dependency, pure Rust with no dependencies of its own. Reads are
scanned a batch at a time via `Database::scan_all`, mirroring STAR's
`ClipMate::clipChunk`, which aligns the adapter against a chunk of reads in one
Opal call.

Both halves are gated against STAR's own C++: `tests/data/cr4_opal_oracle.cpp`
links STAR's `opal.cpp` and reproduces `ClipCR4` + `ClipMate::clipChunk`
verbatim, and the committed `cr4_opal_oracle.tsv` is its output over 938 reads
chosen to straddle both decision boundaries. `cr4_tso_matches_star_opal_oracle`
checks every one against the scalar path, a forced-SIMD path, and the
production batch path.

Measured on 10x mouse chr19 with both tools under `--clipAdapterType
CellRanger4`: clip-amount differences against STARsolo drop from 154 to **0**,
and the soft-clipped share moves from 39% (STAR 41%) to 34.9% (STAR 35.0%).
Default SE/PE alignment is untouched, as `CellRanger4` is opt-in: SE 8788/8926
and PE 8390 both-mapped / 0 half-mapped, both unchanged.

Known limitation: `CellRanger4` combined with a non-zero `--clip5pNbases` /
`--clip3pNbases` misplaces most reads by `clip5pNbases`. That is a pre-existing
bug, not introduced here, and is tracked separately; `CellRanger4` on its own
and the fixed clips on their own are both unaffected.

- **STARsolo single-cell quantification (`--soloType`)** — the 10x
Chromium / plate-based count-matrix pipeline, ported from STAR and
verified against real STARsolo (#90).
Expand Down
7 changes: 7 additions & 0 deletions Cargo.lock

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

5 changes: 5 additions & 0 deletions Cargo.toml
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,11 @@ mimalloc = { version = "0.1", default-features = false }
libmimalloc-sys = { version = "0.1.49", features = ["extended"] } # mi_option_set (purge_delay); see main.rs
libdeflater = "1.25.2"
noodles-bgzf = { version = "0.49", features = ["libdeflate"] }
# Deterministic SIMD overlap aligner for the CellRanger4 5' TSO clip (STAR
# links Opal for this). Pure Rust, zero dependencies, bit-identical across
# backends. 0.2 is the first release with `Database::scan_all`, the per-target
# scan the batched clip path needs.
hyalite = "0.2"

[dev-dependencies]
assert_cmd = "2"
Expand Down
18 changes: 17 additions & 1 deletion src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2135,6 +2135,18 @@ fn align_reads_solo<W: AlignmentWriter + ?Sized>(
move |base: u64, batch: Vec<crate::solo::SoloRead>| -> BatchOut<SoloReadProduct> {
let params: &Parameters = &params_arc;
let index = &index;
// CellRanger4 5' TSO clip: resolved for the whole batch in one
// SIMD pass (hyalite `Database::scan_all`), mirroring STAR's
// `ClipMate::clipChunk`, which aligns the adapter against a
// chunk of reads in a single Opal call. Bit-identical to the
// per-read path, ~7x faster on AVX2.
let cr4_tso: Vec<usize> = if cr4_clip {
let seqs: Vec<&[u8]> =
batch.iter().map(|s| s.cdna.sequence.as_slice()).collect();
crate::solo::tso_clip_lens_cr4_batch(&seqs)
} else {
Vec::new()
};
batch
.par_iter()
.enumerate()
Expand All @@ -2155,7 +2167,11 @@ fn align_reads_solo<W: AlignmentWriter + ?Sized>(
// CellRanger4 adapter clipping (TSO 5' + polyA 3') runs before
// the fixed clip5p/clip3p Nbases trimming.
let (cr_seq, cr_qual, cr4_5p, cr4_3p) = if cr4_clip {
crate::solo::clip_adapter_cr4(&read.sequence, &read.quality)
crate::solo::clip_adapter_cr4_with_tso(
&read.sequence,
&read.quality,
cr4_tso[read_idx],
)
} else {
(read.sequence.clone(), read.quality.clone(), 0, 0)
};
Expand Down
Loading
Loading