Roadmap: remaining STAR-rs port work, split by theme

[BenjaminDEMAILLE]

Following the consolidation PRs (#133, #134, #135, #137, #138, #139), here is a proposed roadmap for the remaining STAR-rs port work, organised by theme rather than by flag, which is the shape the consolidations converged on.

Status: all ten themes now have a PR open. The table below carries the current state; the sections after it are the original analysis and remain accurate except where a measurement overturned them (P9, noted inline).

Everything below branches from 4961e3d4. Per your note on #118 I will keep PRs on the same base commit, stacking only where two PRs rewrite the same functions (stitch.rs, src/solo/, src/io/).

Themes

#	Branch	Theme	Size
P1	bd/aligner-core-faithfulness	Annotated-junction stitching (sjA fast path + boundary snapping), alignEndsType, in-recursion genomic-length penalty	~+1100 / -200
P2	bd/starlong (#116)	STARlong chaining-DP stitcher + rustar-aligner-long binary	~+900 / -70
P3	bd/perf-seed-batch	Batched prefetching MMP search, allocation-free stitch recursion (output-neutral)	~+750 / -250
P4	bd/solo-barcode-chemistry	Chemistry presets, soloCBtype String, adapter-anchored complex geometry, cbMinP posterior threshold	~+1100 / -150
P5	bd/solo-transcript3p-cellfilter	Transcript3p, CellReads.stats, STAR-faithful EmptyDrops_CR, --runMode soloCellFiltering	~+2200 / -350
P6	bd/cli-output-parity	SAM/SJ/read-input knobs and the STAR limit* surface	~+1200 / -150
P7	bd/index-libsais (#109, retargeted)	SuperTranscriptome genome type, transform back-output, RAM-aware SA builder	~+1400 / -100
P8	bd/chimeric-multimap	--chimMultimapNmax detection, chimFilter / score-drop gates	~+750 / -130
P9	bd/perf-parallel-output	Multithreaded BGZF, parallel BAM record encoding, binned coordinate sort (output-neutral)	~+500 / -250
P10	bd/clip-cellranger4	--clipAdapterType CellRanger4, via an in-tree deterministic SIMD Smith-Waterman engine	~+2200 / -40

Where P1 comes from

P1 is not a flag PR. It closes the open items in your own STAR-RS-COMPARISON.md §7.2: #4b (annotated-boundary snapping, which §7.4 calls "the best remaining lever"), #7 (genomic-length penalty inside the recursion rather than a post-pass), #8 (alignEndsType / extend-to-end), #9 (same_structure dedup guard, stitch mismatch cap, alignInsertionFlush Right).

§7.2 item #2 (the gLength tie-break key) stays dropped: §7.5 records that it regresses.

Today stitch.rs:1389-1407 calls db.is_annotated(...) for strand 0/1/2 and adds sjdb_score on a hit, but it never snaps exon coordinates onto the annotated donor/acceptor, and WindowAlignment has no sj_a field so the sjAB fast path cannot exist. Enabling work: SpliceJunctionDb gains a find(stored_start, stored_end) next to the existing boolean is_annotated, and WindowAlignment gains sj_a.

This does not change the index format. PreparedJunction (src/junction/sjdb_insert.rs:258-273) already carries motif, shift_left, shift_right and strand; sjdbInfo.txt already writes and reads all four; sort_and_dedup (:368) already sorts by (stored_start, stored_end), which is exactly the ordering STAR's binarySearch2 needs; and decode_gsj_hit (:157) already computes the junction index at :174 and then discards it. So find() is a search over data that is already on disk and already loaded, and sj_a is a value that is already computed. No version bump, no rebuild of existing indices.

One latent bug surfaces along the way: the align-time --sjdbGTFfile path (src/index/io.rs:72-79) goes through from_raw_junctions, which throws away motif and shift, while the index-loaded path keeps them. The two paths should build identical tables; P1 makes them do so.

P1 also adds a yeast_100_gtf tier to test/run_tests.sh and test/ci.sh. §7.5 documents that rustar recognised 0 annotated junctions against STAR's 71 for months precisely because no integration tier passed a GTF at align time. Without an annotated tier, #4b is untestable in CI.

Two performance PRs

Both are strictly output-neutral, with an empty diff against the previous binary on yeast 10k SE and PE as the acceptance criterion.

P9 is the bigger wall-clock win and touches nothing in the aligner core. make_bgzf_writer (src/io/bam.rs:48) builds a single-threaded bgzf::io::Writer, even though noodles-bgzf 0.49 is already declared with features = ["libdeflate"] and libdeflater 1.25.2 is already a dependency (used only in src/solo/count.rs). So this is wiring on dependencies that already pass the 5-platform matrix, not a new dependency. On the STAR-rs side the equivalent single-threaded BGZF cost ~39 s of an ~81 s run on 8M GIAB pairs, that is, it alone made the aligner ~2x slower than STAR while alignment was already at parity.

P3 is the seed-search side: rustar already has SuffixArray::prefetch, what is missing is interleaving 32 independent MMP searches so memory-level parallelism goes from 1 to 32. It is deliberately sequenced last, since an A/B measurement taken while P1 and P2 are still moving the transcript set is worthless.

STAR C++ bugs

Where STAR 2.7.11b has undefined or plainly wrong behaviour, I propose fixing it by default, with a documented entry in a new docs-old/dev/divergences.md and a locking test asserting the correct result (never asserting STAR's UB). Each fix is a separately titled commit inside its theme PR, not its own PR.

• --soloUMIfiltering MultiGeneUMI_All is a no-op in STAR (its consumption test checks only the MultiGeneUMI flag). rustar currently aliases it to MultiGeneUMI (src/solo/count.rs:125), which is neither STAR's behaviour nor the documented one. Separate issue to follow. Lands in P5.
• --genomeType SuperTranscriptome minus-strand sequence reads uninitialised memory in STAR. Cannot land anywhere else, the feature is born in P7.
• EmptyDrops_CR Simple-Good-Turing with fewer than 5 distinct frequencies. Lands in P5 alongside the faithful libc++ mt19937 / discrete_distribution port.
• CB_UMI_Complex homopolymer UMI (umiL=0 quirk). Lands in P4.

Two are deliberately not ported: --outMultimapperOrder Random is already resolved here by #111/#139, and the transcriptome-BAM primary-flag choice is already thread-invariant here via the per-read seed, so there is nothing to gain from adopting the STAR-rs j == 0 convention.

Proposed order

P1, P6, P7, P9, P2, P8, P4, P5, P3, P10, with at most 3 or 4 open at a time and stacked PRs opened only after their base merges.

Happy to reshuffle, split further, or drop any of these. If you would rather I bundle differently, say so before I open the first one.

Current state

#	PR	State
P1	#146	open — 835 → 933/1000 reads matching STAR on the new annotated-junction tier
P2	#116	open — rebased
P3	—	not started; deliberately last, since an A/B taken while #146/#116 still move the transcript set is worthless
P4	#150	open — solo differential against real STARsolo still owed
P5	#152	open — MultiGeneUMI_All + libc++ RNG; Transcript3p and CellReads.stats not yet done
P6	#147	open — parameter surface 143/200 → 174/200, machine-checked
P7	#109	open — libsais selected by --limitGenomeGenerateRAM; SuperTranscriptome not yet done
P8	#149	open — chimericDetectionMult + chimFilter
P9	#151	open — binned sort only. The multithreaded-BGZF half was measured ~37% slower and dropped (comment below); #153 supersedes it with an optional htslib backend, ~17% faster
P10	#148	open — narrowed to the poly-A trim; the SW engine belongs in a separate library per review on that PR

[BenjaminDEMAILLE]

Negative result on the P9 BGZF item. Recorded so it is not repeated.

The claim in the issue body that single-threaded BGZF is the biggest wall-clock lever is wrong for this codebase. Swapping to noodles_bgzf::MultithreadedWriter is a regression.

Measured on yeast, --outSAMtype BAM Unsorted --outBAMcompression 6, same binary otherwise:

reads	BAM out	threads	main	MultithreadedWriter
120k	7 MB	8	0.60 s	0.84 s
120k	7 MB	16	0.91 s	1.20 s
720k	44 MB	8	2.00 s	2.74 s

Consistently ~37% slower, with no cross-over as output grows.

Cause: noodles-bgzf is declared here with its libdeflate feature, and both Writer and MultithreadedWriter go through the same deflate::encode. The single-threaded path was never the naive-deflate bottleneck the STAR-rs comparison describes; the multithreaded writer only adds per-block channel and ordering overhead.

Decoded records were identical across 1/4/16 threads and against main, so the change was correct — just slower. The branch was dropped rather than opened as a PR.

Not ruled out: a run large enough that per-block compression outweighs coordination, or a writer batching many blocks per hand-off. Either needs measuring at GIAB scale rather than on yeast.

The other parts of that theme are unaffected: binned coordinate sorting landed as #151, and parallel record encoding is untouched.

[BenjaminDEMAILLE]

Follow-up to the BGZF negative result: rust-htslib. Implemented as #153.

htslib's hts_set_threads is a genuinely parallel deflate rather than a layer over one compressor, which is what the noodles route failed to deliver. Measured at 16 threads, five interleaved A/B pairs: noodles median 2.51 s, htslib median 2.08 s, ~17% faster, non-overlapping. Decoded records identical.

It cannot be a mandatory dependency, and that decides the design:

• rust-htslib ships pre-built bindings for Mac and Linux only; its README calls the Windows bindgen path untested.
• Its CI runs ubuntu-22.04, ubuntu-22.04-arm and macOS. No Windows job.
• This project's CI matrix requires windows-x86_64.

So #153 declares it under cfg(not(windows)) behind an off-by-default htslib-bam feature, with default-features = false dropping bzip2-sys and lzma-sys. noodles stays the default and the only path on Windows; the published crate and the CI matrix are unchanged.

The trade that remains open regardless of gating: it pulls htslib (C) into a project otherwise pure Rust apart from libdeflater. That is a dependency-story change on a crates.io crate. If it is unwanted, the measurement stands on its own as a record of where the BAM write cost sits.