perf(align): resolve the genome storage once per extension and junction scan - #257
Open
BenjaminDEMAILLE wants to merge 2 commits into
Open
perf(align): resolve the genome storage once per extension and junction scan#257BenjaminDEMAILLE wants to merge 2 commits into
BenjaminDEMAILLE wants to merge 2 commits into
Conversation
…igns
`find_best_junction_position` is the single hottest function in the
aligner: on a 50k-pair yeast PE run it accounts for ~60% of alignment
time. The scan itself is not wasteful, but it is repeated.
`stitchWindowAligns`' include/exclude recursion reaches the same
(exon A end, seed B) pair through many different branch paths, and each
path re-runs the identical scan. The scan is a pure function of its
arguments, and within one window `read_seq`, the genome, `is_reverse`
and `n_genome` are all fixed, so six coordinates identify a scan
completely: the exon A read end and genome end, the read and genome
gaps, the previous exon length and the next seed length.
Add `JunctionScanCache`, a per-window `FxHashMap` on that key, and a
`find_best_junction_position_cached` wrapper that consults it. The
uncached function is untouched, so a hit returns exactly what a fresh
scan would have. The cache is created per window in `stitch_seeds_core`
and threaded down the recursion, which is what keeps the key complete:
it never outlives the read, genome and strand it was filled for. An
empty `HashMap` does not allocate, so windows that stitch nothing pay
nothing.
Measured on 50k yeast read pairs (Apple M4 Max, quiet machine,
best-of-6 interleaved rounds, `--outSAMtype None`):
threads before after change
1 19.51s 16.76s -14.1%
8 2.55s 2.18s -14.5%
Output is byte-identical: `Aligned.out.sam` (records and header alike,
modulo the `@PG` CL line naming the binary) and `SJ.out.tab` compare
equal against the pre-change binary on the same input. 592 tests pass,
0 clippy warnings.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…on scan `Genome::get_base` matches on the `GenomeSeq` discriminant, bounds-checks and, for a memory-mapped genome's reverse-complement half, recomputes the mirrored index and complements the byte. That is fine per call, but the alignment inner loops call it once per base: the junction-position scan reads two bases per candidate position across three loops, and `extend_alignment` reads one per extended base. Together those two functions are ~55% of alignment time after the scan memo. Add `GenomeSeq::view()`, returning a `Copy` `SeqView` that resolves the variant once. The view is a slice plus one integer, so the loops keep it in registers and each base costs a bounds check and a load. Hoist it out of the junction scan (including the sliding motif window) and out of both `extend_alignment` loops. Out-of-range reads return the `OUT_OF_RANGE` sentinel instead of `None`. Every call site converted here already treated "not one of A/C/G/T" the same way a `None` was treated, so the branch structure is preserved exactly; `score.rs` already had this sentinel locally for the motif window and it moves next to the view it belongs to. `SeqView::base` duplicates the reverse-complement arithmetic in `GenomeSeq::base`, so a unit test asserts the two agree on every index in `0..2n` plus the first out-of-range one, for both storage variants. Measured on 50k yeast read pairs (Apple M4 Max, quiet machine, best-of-6 interleaved rounds, `--outSAMtype None`, 8 threads), on top of the junction-scan memo: 2.18s to 2.12s, and 2.60s to 2.12s against the pre-memo baseline (-18%). Output is byte-identical: `Aligned.out.sam` and `SJ.out.tab` compare equal against the pre-change binary on the same input. 593 tests pass, 0 clippy warnings. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
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
1 participant
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.
Stacked on #256, which cannot be used as a base branch here (it lives on a fork), so this PR is opened against
mainand its diff contains #256's commit as well. Review only the second commit (perf(align): resolve the genome storage once...); merge #256 first and this one rebases down to that single commit.What
Genome::get_basematches on theGenomeSeqdiscriminant, bounds-checks, and for a memory-mapped genome's reverse-complement half recomputes the mirrored index and complements the byte. Fine per call, but the alignment inner loops call it once per base:extend_alignmentreads one per extended baseAfter #256 those two functions are still ~55% of alignment time (31% + 24%).
This PR adds
GenomeSeq::view(), returning aCopySeqViewthat resolves the variant once. The view is a slice plus one integer, so the loops keep it in registers and each base costs a bounds check and a load. It is hoisted out of the junction scan (motif window included) and out of bothextend_alignmentloops.Sentinel instead of Option
Out-of-range reads return
OUT_OF_RANGErather thanNone. Every call site converted here already treated "not one of A/C/G/T" the same way it treated aNone, so the branch structure is preserved exactly.score.rsalready had this sentinel locally for the motif window; it moves next to the view it belongs to.SeqView::baseduplicates the reverse-complement arithmetic inGenomeSeq::base, so a unit test asserts the two agree on every index in0..2nplus the first out-of-range one, for both storage variants. That test is the guard against the two definitions drifting apart.Benchmark
50k yeast read pairs, Apple M4 Max, quiet machine, best-of-6 interleaved rounds,
--outSAMtype None, 8 threads:Smaller than #256 on its own, and reported as such: most of the remaining cost in those loops is the scanning work itself, not the dispatch.
Correctness
Output is byte-identical.
Aligned.out.samandSJ.out.tabcompare equal against the pre-change binary on the same input.--all-targets --release)cargo fmt --checkclean🤖 Generated with Claude Code