From 677f8ba7fe41d3e2667782b58bf0ccd55c5005f6 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 15:20:52 +0200 Subject: [PATCH 1/5] perf(score): slide the splice-motif window across the junction scan `detect_splice_motif` and `find_best_junction_position` were the two largest self-time entries in a sampled 2M-read run, together about 45% of the on-CPU samples. The scan moves the junction one base per iteration, so the four bases that decide the motif (the intron's first two and last two) overlap the previous position in two of four places. Carrying them in a `MotifWindow` and sliding turns four genome reads per position into two. Whether the motif branch runs at all is decided once before the loop rather than per iteration, since `del` does not change across the scan. An out-of-range position becomes a sentinel that no motif arm matches, which is what `get_base` returning `None` already meant. Output-neutral: SAM byte-identical to the parent commit on 200k real reads. Interleaved A/B at 16 threads on 2M reads, machine at 87-95% CPU idle: 22.94s median before, 22.29s after, the same direction in all four rounds. About 2.8%. Co-Authored-By: Claude Opus 5 (1M context) --- src/align/score.rs | 113 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 97 insertions(+), 16 deletions(-) diff --git a/src/align/score.rs b/src/align/score.rs index 7a34f55..0e89881 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -353,6 +353,18 @@ impl AlignmentScorer { let mut best_motif = SpliceMotif::NonCanonical; let mut best_motif_score = self.score_gap_noncan; + // `del` does not change across the scan, so whether the motif branch + // runs is decided once rather than per iteration. + let motif_in_range = + del >= self.align_intron_min as i64 && del <= self.align_intron_max as i64; + + // The four motif bases sit at `donor`, `donor+1`, `donor+del-2` and + // `donor+del-1`, and `donor` moves exactly one base per iteration, so + // consecutive iterations share two of the four. Carrying the window + // instead of re-fetching it turns four genome reads per position into + // two. `window` is `None` until the first motif iteration primes it. + let mut window: Option = None; + loop { let ri = r_a_end_inc as i64 + jr1 as i64; if ri >= 0 && (ri as usize) < read_seq.len() { @@ -378,7 +390,7 @@ impl AlignmentScorer { } // Check splice motif at this junction position - if del >= self.align_intron_min as i64 && del <= self.align_intron_max as i64 { + if motif_in_range { // Donor position in SA space: one past the last donor-exon base let donor_sa = (g_a_end_inc as i64 + jr1 as i64 + 1) as u64; // Convert to forward genome coordinates for motif detection @@ -387,7 +399,14 @@ impl AlignmentScorer { } else { donor_sa }; - let motif = self.detect_splice_motif(donor_fwd, del as u32, genome); + let w = match window.as_mut() { + Some(w) => { + w.slide_to(donor_fwd, del as u64, genome); + &*w + } + None => window.insert(MotifWindow::at(donor_fwd, del as u64, genome)), + }; + let motif = w.motif(); let motif_score = self.score_splice_junction(motif); let score2 = score1 + motif_score; @@ -535,20 +554,82 @@ impl AlignmentScorer { /// `donor_pos` is the 0-based position of the intron's first base on the /// forward strand; `intron_len` is the intron length in bases. pub fn detect_splice_motif(donor_pos: u64, intron_len: u32, genome: &Genome) -> SpliceMotif { - let d1 = genome.get_base(donor_pos); - let d2 = genome.get_base(donor_pos + 1); - let a1 = genome.get_base(donor_pos + intron_len as u64 - 2); - let a2 = genome.get_base(donor_pos + intron_len as u64 - 1); - - // Base encoding: A=0, C=1, G=2, T=3. - match (d1, d2, a1, a2) { - (Some(2), Some(3), Some(0), Some(2)) => SpliceMotif::GtAg, - (Some(2), Some(1), Some(0), Some(2)) => SpliceMotif::GcAg, - (Some(0), Some(3), Some(0), Some(1)) => SpliceMotif::AtAc, - (Some(1), Some(3), Some(0), Some(1)) => SpliceMotif::CtAc, - (Some(1), Some(3), Some(2), Some(1)) => SpliceMotif::CtGc, - (Some(2), Some(3), Some(0), Some(3)) => SpliceMotif::GtAt, - _ => SpliceMotif::NonCanonical, + MotifWindow::at(donor_pos, intron_len as u64, genome).motif() +} + +/// A position off the end of the genome. No motif arm matches it, so it falls +/// through to `NonCanonical` exactly as `get_base` returning `None` did. +const OUT_OF_RANGE: u8 = u8::MAX; + +#[inline] +fn base_or_out_of_range(genome: &Genome, pos: u64) -> u8 { + genome.get_base(pos).unwrap_or(OUT_OF_RANGE) +} + +/// The four bases that decide a splice motif: the intron's first two and last +/// two, on the forward strand. +/// +/// Kept as a struct so the junction scan can slide it. Successive junction +/// positions differ by one base, so three of the four positions overlap the +/// previous window and only two bases have to be read from the genome. +#[derive(Clone, Copy)] +struct MotifWindow { + donor: u64, + d1: u8, + d2: u8, + a1: u8, + a2: u8, +} + +impl MotifWindow { + #[inline] + fn at(donor: u64, intron_len: u64, genome: &Genome) -> Self { + Self { + donor, + d1: base_or_out_of_range(genome, donor), + d2: base_or_out_of_range(genome, donor + 1), + a1: base_or_out_of_range(genome, donor + intron_len - 2), + a2: base_or_out_of_range(genome, donor + intron_len - 1), + } + } + + /// Move the window to `donor`, reusing what overlaps. + /// + /// A step of exactly one base either way shares two of the four: moving + /// right, the old `d2`/`a2` become the new `d1`/`a1`; moving left, the old + /// `d1`/`a1` become the new `d2`/`a2`. Any other step is rare enough that + /// re-reading all four is the simpler answer. + #[inline] + fn slide_to(&mut self, donor: u64, intron_len: u64, genome: &Genome) { + if donor == self.donor + 1 { + self.d1 = self.d2; + self.a1 = self.a2; + self.d2 = base_or_out_of_range(genome, donor + 1); + self.a2 = base_or_out_of_range(genome, donor + intron_len - 1); + self.donor = donor; + } else if donor + 1 == self.donor { + self.d2 = self.d1; + self.a2 = self.a1; + self.d1 = base_or_out_of_range(genome, donor); + self.a1 = base_or_out_of_range(genome, donor + intron_len - 2); + self.donor = donor; + } else if donor != self.donor { + *self = Self::at(donor, intron_len, genome); + } + } + + /// Base encoding: A=0, C=1, G=2, T=3. + #[inline] + fn motif(&self) -> SpliceMotif { + match (self.d1, self.d2, self.a1, self.a2) { + (2, 3, 0, 2) => SpliceMotif::GtAg, + (2, 1, 0, 2) => SpliceMotif::GcAg, + (0, 3, 0, 1) => SpliceMotif::AtAc, + (1, 3, 0, 1) => SpliceMotif::CtAc, + (1, 3, 2, 1) => SpliceMotif::CtGc, + (2, 3, 0, 3) => SpliceMotif::GtAt, + _ => SpliceMotif::NonCanonical, + } } } From 56c1829e1e650cc3e6189a594faa3a339739ae9b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 20:00:04 +0200 Subject: [PATCH 2/5] perf(score): look the splice motif up in a table instead of matching The four-base motif decision compiled to a chain of comparisons whose outcome is data-dependent and close to unpredictable, and the junction scan pays it at every position. A 256-entry table indexed by the four bases packed two bits each replaces it with one load. `examples/jrbench.rs` measures the scan in isolation: 8.5 ns per iteration before, 5.2 ns after, a 38% cut, with identical results. That example exists because the earlier attempts on this path were guesses about whether the loop was memory-bound, and two of them were wrong. It answers the question by construction, holding the iteration count and the instruction mix fixed and changing only how far the acceptor sits from the donor: del ns/iter 64 8.42 1024 8.47 16384 8.52 262144 8.41 4194304 8.62 67108864 8.56 Flat from 64 bases to 67 million, so the acceptor distance costs nothing and the loop is not stalled on that stream. Which is why skipping the acceptor read behind a donor-pair test measured *slower*: it traded a prefetchable access for an unpredictable branch. The branch was the cost all along, and this removes it. `the_motif_table_agrees_with_the_original_match_on_every_input` checks the table against the match it replaces over every combination of `A`/`C`/`G`/`T`, `N`, the chromosome-boundary byte and the out-of-range sentinel, so the inputs that no match arm covered are covered explicitly. Output-neutral: SAM byte-identical on 200k real reads. Co-Authored-By: Claude Opus 5 (1M context) --- examples/jrbench.rs | 97 +++++++++++++++++++++++++++++++++++++++++++++ src/align/score.rs | 82 ++++++++++++++++++++++++++++++++++---- 2 files changed, 171 insertions(+), 8 deletions(-) create mode 100644 examples/jrbench.rs diff --git a/examples/jrbench.rs b/examples/jrbench.rs new file mode 100644 index 0000000..d987c9b --- /dev/null +++ b/examples/jrbench.rs @@ -0,0 +1,97 @@ +//! Is the junction scan stalled on memory, and is it the acceptor stream? +//! +//! The stack sampler gives self time, not cache misses, so this asks the +//! question by construction instead: run the same scan over the same genome, +//! changing only how far the acceptor sits from the donor. Everything else +//! (iteration count, branch pattern, instruction mix) is held fixed, because +//! `del` enters the loop only as an address offset once it is inside the +//! intron-length range. +//! +//! If time per iteration is flat across `del`, the loop is not memory-bound on +//! the acceptor and prefetching it would be wasted work. If it climbs with +//! `del`, that stream is the cost. +//! +//! Run: cargo run --release --example jrbench + +use rustar_aligner::align::score::AlignmentScorer; +use rustar_aligner::genome::{Genome, GenomeSeq}; +use std::time::Instant; + +/// Deterministic pseudo-random bases, so the run is reproducible and the motif +/// hit rate is the same at every `del`. +fn synthetic_genome(n: usize) -> Genome { + let mut state = 0x2545_F491_4F6C_DD1Du64; + let mut seq = Vec::with_capacity(n); + for _ in 0..n { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + seq.push((state >> 33) as u8 & 3); + } + Genome { + transform_blocks: None, + sequence: GenomeSeq::Owned(seq), + n_genome: n as u64, + n_genome_real: n as u64, + n_chr_real: 1, + chr_name: vec!["chr1".to_string()], + chr_length: vec![n as u64], + chr_start: vec![0, n as u64], + } +} + +fn main() { + // Large enough that the donor and acceptor streams cannot both stay in + // cache when they are far apart. + const N: usize = 256 << 20; // 256 MB of bases, one byte each + const READ_LEN: usize = 150; + const SCANS: usize = 20_000; + + let genome = synthetic_genome(N); + let mut scorer = AlignmentScorer::from_params_minimal(); + scorer.align_intron_min = 20; + scorer.align_intron_max = u32::MAX; + + let mut read = vec![0u8; READ_LEN]; + let mut state = 0x9E37_79B9_7F4A_7C15u64; + for b in read.iter_mut() { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + *b = (state >> 33) as u8 & 3; + } + + println!("{:>12} {:>10} {:>12}", "del", "wall", "ns/iter"); + for del in [64u64, 1_024, 16_384, 262_144, 4_194_304, 67_108_864] { + // Spread the scans over the genome so no single window stays resident. + let stride = (N as u64 - del - 4096) / SCANS as u64; + let iters_per_scan = 100usize; // r_gap 0 + next_seed_len 100 + let t0 = Instant::now(); + let mut sink = 0i64; + for k in 0..SCANS { + let g_a_end = 2048 + k as u64 * stride; + let (jr, _motif, score, _l, _r) = scorer.find_best_junction_position( + &read, + 75, + g_a_end, + 0, + del as i64, + &genome, + false, + N as u64, + 75, + iters_per_scan, + ); + sink += jr as i64 + score as i64; + } + let dt = t0.elapsed(); + let iters = (SCANS * iters_per_scan) as f64; + println!( + "{:>12} {:>9.3}s {:>11.2} (sink {})", + del, + dt.as_secs_f64(), + dt.as_secs_f64() * 1e9 / iters, + sink + ); + } +} diff --git a/src/align/score.rs b/src/align/score.rs index 0e89881..408a2d2 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -619,20 +619,45 @@ impl MotifWindow { } /// Base encoding: A=0, C=1, G=2, T=3. + /// + /// A table lookup rather than a match on the four bases. The match compiled + /// to a chain of compares whose outcome is data-dependent and close to + /// unpredictable, which the scan pays at every junction position; the table + /// is one load at a computed index. `|` binds tighter than `>=` in Rust, so + /// the guard tests the union of the four bases and rejects anything that is + /// not `A`, `C`, `G` or `T` — `N`, the chromosome-boundary byte, and the + /// out-of-range sentinel all take that path, exactly as no match arm + /// covered them. #[inline] fn motif(&self) -> SpliceMotif { - match (self.d1, self.d2, self.a1, self.a2) { - (2, 3, 0, 2) => SpliceMotif::GtAg, - (2, 1, 0, 2) => SpliceMotif::GcAg, - (0, 3, 0, 1) => SpliceMotif::AtAc, - (1, 3, 0, 1) => SpliceMotif::CtAc, - (1, 3, 2, 1) => SpliceMotif::CtGc, - (2, 3, 0, 3) => SpliceMotif::GtAt, - _ => SpliceMotif::NonCanonical, + let (d1, d2, a1, a2) = (self.d1, self.d2, self.a1, self.a2); + if d1 | d2 | a1 | a2 >= 4 { + return SpliceMotif::NonCanonical; } + MOTIF_TABLE[motif_index(d1, d2, a1, a2)] } } +/// Every four-base combination, packed `d1 d2 a1 a2` at two bits each. +static MOTIF_TABLE: [SpliceMotif; 256] = build_motif_table(); + +/// Pack the four bases into the table index, two bits each. +const fn motif_index(d1: u8, d2: u8, a1: u8, a2: u8) -> usize { + ((d1 << 6) | (d2 << 4) | (a1 << 2) | a2) as usize +} + +const fn build_motif_table() -> [SpliceMotif; 256] { + // A=0, C=1, G=2, T=3. + let mut t = [SpliceMotif::NonCanonical; 256]; + t[motif_index(2, 3, 0, 2)] = SpliceMotif::GtAg; + t[motif_index(2, 1, 0, 2)] = SpliceMotif::GcAg; + t[motif_index(0, 3, 0, 1)] = SpliceMotif::AtAc; + t[motif_index(1, 3, 0, 1)] = SpliceMotif::CtAc; + t[motif_index(1, 3, 2, 1)] = SpliceMotif::CtGc; + t[motif_index(2, 3, 0, 3)] = SpliceMotif::GtAt; + t +} + /// Splice junction motif types #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] pub enum SpliceMotif { @@ -705,6 +730,47 @@ mod tests { use super::*; + /// The table has to agree with the match it replaced on every input the + /// scan can present, not merely on the six motifs. That includes `N` (4), + /// the chromosome-boundary byte (5) and the out-of-range sentinel, all of + /// which no match arm covered and which must therefore be `NonCanonical`. + #[test] + fn the_motif_table_agrees_with_the_original_match_on_every_input() { + fn original(d1: u8, d2: u8, a1: u8, a2: u8) -> SpliceMotif { + match (d1, d2, a1, a2) { + (2, 3, 0, 2) => SpliceMotif::GtAg, + (2, 1, 0, 2) => SpliceMotif::GcAg, + (0, 3, 0, 1) => SpliceMotif::AtAc, + (1, 3, 0, 1) => SpliceMotif::CtAc, + (1, 3, 2, 1) => SpliceMotif::CtGc, + (2, 3, 0, 3) => SpliceMotif::GtAt, + _ => SpliceMotif::NonCanonical, + } + } + + let values = [0u8, 1, 2, 3, 4, 5, OUT_OF_RANGE]; + for &d1 in &values { + for &d2 in &values { + for &a1 in &values { + for &a2 in &values { + let w = MotifWindow { + donor: 0, + d1, + d2, + a1, + a2, + }; + assert_eq!( + w.motif(), + original(d1, d2, a1, a2), + "disagreement at ({d1}, {d2}, {a1}, {a2})" + ); + } + } + } + } + } + fn make_test_genome(seq: &[u8]) -> Genome { // Create simple genome with one chromosome let n_genome = ((seq.len() as u64 + 1) / 64 + 1) * 64; // Pad to 64-byte boundary From f3760e08493ac0f2062af131d40621b926ed9e31 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 20:12:05 +0200 Subject: [PATCH 3/5] docs(changelog): record the junction-scan motif work Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..ddf5041 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -131,3 +131,11 @@ Sections commonly used: Features, Bug fixes, Other changes. needs random access to the SA in RAM. Initial release of Rust rewrite of STAR. +### Other changes + +- The splice-motif check in the junction scan carries a sliding window + and looks the motif up in a table instead of re-reading four genome + bases and matching on them at every position. Output is unchanged + (SAM byte-identical on 200k reads); the scan itself goes from 8.5 ns + to 5.2 ns per iteration, measured by `examples/jrbench.rs`. + From ad21f29f81e8d9e07ea838cd0c3fe417b686b342 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 20:58:26 +0200 Subject: [PATCH 4/5] fix(test): make the unmapped-reason test drive the shipped logic The test at `read_align.rs` built its own `HashMap<&str, i32>` and re-implemented the `has_mismatch` / `has_other` decision inside the test body, then asserted on its own copy. It would have passed unchanged if the shipped code had been deleted. It now drives `FilterReasons`, which is the same set of counters the aligner uses, so a change to that decision fails the test. `FilterReasons` replaces the `HashMap<&str, i32>` the filter loop used: a fixed `[u32; 9]` indexed by constant. That was a heap allocation per read and a string hash per filtered transcript, for something whose only consumers ask whether a count is non-zero (two debug logs and the `unmapped_reason` decision), so nothing here can reach the output except through that decision. Its `Debug` prints only the non-zero reasons, in a fixed order, which is more readable than the map's arbitrary one. `SpliceJunctionDb` moves to `FxHashMap` in the same spirit, matching what `cluster_seeds` already does for its integer-keyed maps. Only `get` and `insert` are called on it and it is never iterated. Honest about the performance: **neither change is measurable end to end.** Six interleaved rounds on 2M reads gave medians 20.70 s against 20.48 s with the direction mixed, which is inside the run-to-run spread. They are here because they remove work and because the test was not testing anything, not because they show up on a clock. Output-neutral: SAM byte-identical on 200k real reads. Co-Authored-By: Claude Opus 5 (1M context) --- src/align/read_align.rs | 142 +++++++++++++++++++++++++++++++--------- src/junction/mod.rs | 18 +++-- 2 files changed, 123 insertions(+), 37 deletions(-) diff --git a/src/align/read_align.rs b/src/align/read_align.rs index b4d2f20..eaffd37 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -12,6 +12,81 @@ use crate::params::{IntronMotifFilter, IntronStrandFilter, MultimapperOrder, Par use crate::stats::UnmappedReason; use std::hash::{DefaultHasher, Hash, Hasher}; +/// Why a transcript was dropped by the quality filters. +/// +/// A fixed set of counters rather than a `HashMap<&str, i32>`: this is built +/// once per read and incremented per filtered transcript, so the map cost a +/// heap allocation per read and a string hash per increment. The counts feed +/// two debug logs and the `unmapped_reason` decision below, all of which ask +/// only whether a count is non-zero, so nothing here can reach the output +/// except through that decision. +#[derive(Default, Clone, Copy)] +struct FilterReasons { + counts: [u32; FilterReasons::N], +} + +impl FilterReasons { + const N: usize = 9; + const NAMES: [&'static str; Self::N] = [ + "score_min", + "score_min_relative", + "mismatch_max", + "mismatch_rate", + "match_min", + "match_min_relative", + "noncanonical_junction", + "noncanonical_unannotated_junction", + "inconsistent_strand", + ]; + + const SCORE_MIN: usize = 0; + const SCORE_MIN_RELATIVE: usize = 1; + const MISMATCH_MAX: usize = 2; + const MISMATCH_RATE: usize = 3; + const MATCH_MIN: usize = 4; + const MATCH_MIN_RELATIVE: usize = 5; + const NONCANONICAL_JUNCTION: usize = 6; + const NONCANONICAL_UNANNOTATED_JUNCTION: usize = 7; + const INCONSISTENT_STRAND: usize = 8; + + #[inline] + fn add(&mut self, reason: usize) { + self.counts[reason] += 1; + } + + #[inline] + fn has(&self, reason: usize) -> bool { + self.counts[reason] > 0 + } + + /// Any reason other than the two mismatch ones. + #[inline] + fn has_non_mismatch(&self) -> bool { + self.counts + .iter() + .enumerate() + .any(|(i, &c)| c > 0 && i != Self::MISMATCH_MAX && i != Self::MISMATCH_RATE) + } +} + +impl std::fmt::Debug for FilterReasons { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let mut first = true; + write!(f, "{{")?; + for (i, &c) in self.counts.iter().enumerate() { + if c == 0 { + continue; + } + if !first { + write!(f, ", ")?; + } + first = false; + write!(f, "{}: {}", Self::NAMES[i], c)?; + } + write!(f, "}}") + } +} + /// Derive a deterministic per-read RNG seed from `run_rng_seed` + the read name. /// /// STAR seeds `std::mt19937` once per chunk/thread (`runRNGseed*(iChunk+1)`), @@ -409,24 +484,24 @@ pub fn align_read( // Log filtering statistics let pre_filter_count = transcripts.len(); - let mut filter_reasons = std::collections::HashMap::new(); + let mut filter_reasons = FilterReasons::default(); transcripts.retain(|t| { // Absolute score threshold if t.score < params.out_filter_score_min { - *filter_reasons.entry("score_min").or_insert(0) += 1; + filter_reasons.add(FilterReasons::SCORE_MIN); return false; } // Relative score threshold: STAR casts to intScore (i32) if t.score < (params.out_filter_score_min_over_lread * lread_m1) as i32 { - *filter_reasons.entry("score_min_relative").or_insert(0) += 1; + filter_reasons.add(FilterReasons::SCORE_MIN_RELATIVE); return false; } // Absolute mismatch count if t.n_mismatch > params.out_filter_mismatch_nmax { - *filter_reasons.entry("mismatch_max").or_insert(0) += 1; + filter_reasons.add(FilterReasons::MISMATCH_MAX); log::debug!( "Filtered {}: {} mismatches > {} max (read_len={}, score={})", read_name, @@ -441,7 +516,7 @@ pub fn align_read( // Relative mismatch count (mismatches / read_length) let mismatch_rate = t.n_mismatch as f64 / read_length; if mismatch_rate > params.out_filter_mismatch_nover_lmax { - *filter_reasons.entry("mismatch_rate").or_insert(0) += 1; + filter_reasons.add(FilterReasons::MISMATCH_RATE); log::debug!( "Filtered {}: {:.1}% mismatch rate > {:.1}% max ({}/{} bases, score={})", read_name, @@ -457,13 +532,13 @@ pub fn align_read( // Absolute matched bases let n_matched = t.n_matched(); if n_matched < params.out_filter_match_nmin as usize { - *filter_reasons.entry("match_min").or_insert(0) += 1; + filter_reasons.add(FilterReasons::MATCH_MIN); return false; } // Relative matched bases: STAR casts to uint (u32) if (n_matched as f64) < params.out_filter_match_nmin_over_lread * lread_m1 { - *filter_reasons.entry("match_min_relative").or_insert(0) += 1; + filter_reasons.add(FilterReasons::MATCH_MIN_RELATIVE); return false; } @@ -475,7 +550,7 @@ pub fn align_read( IntronMotifFilter::RemoveNoncanonical => { // Reject if any junction is non-canonical if t.junction_motifs.contains(&SpliceMotif::NonCanonical) { - *filter_reasons.entry("noncanonical_junction").or_insert(0) += 1; + filter_reasons.add(FilterReasons::NONCANONICAL_JUNCTION); return false; } } @@ -486,9 +561,7 @@ pub fn align_read( .zip(t.junction_annotated.iter()) .any(|(m, annotated)| *m == SpliceMotif::NonCanonical && !annotated) { - *filter_reasons - .entry("noncanonical_unannotated_junction") - .or_insert(0) += 1; + filter_reasons.add(FilterReasons::NONCANONICAL_UNANNOTATED_JUNCTION); return false; } } @@ -512,7 +585,7 @@ pub fn align_read( } } if has_plus && has_minus { - *filter_reasons.entry("inconsistent_strand").or_insert(0) += 1; + filter_reasons.add(FilterReasons::INCONSISTENT_STRAND); return false; } } @@ -639,11 +712,9 @@ pub fn align_read( // if the quality filter emptied the set, the best transcript failed → classify // short/mismatch; otherwise, if too many loci survive, it's multi; else mapped. let unmapped_reason = if transcripts.is_empty() { - let has_mismatch = filter_reasons.contains_key("mismatch_max") - || filter_reasons.contains_key("mismatch_rate"); - let has_other = filter_reasons - .keys() - .any(|k| *k != "mismatch_max" && *k != "mismatch_rate"); + let has_mismatch = filter_reasons.has(FilterReasons::MISMATCH_MAX) + || filter_reasons.has(FilterReasons::MISMATCH_RATE); + let has_other = filter_reasons.has_non_mismatch(); if has_mismatch && !has_other { Some(UnmappedReason::TooManyMismatches) } else { @@ -1708,14 +1779,17 @@ mod tests { fn test_unmapped_reason_mismatch_classification() { // Verify the filter_reasons logic used to derive unmapped_reason. // mismatch-only → TooManyMismatches; anything else (or mixed) → TooShort. - let classify = |reasons: &[&str]| -> UnmappedReason { - let filter_reasons: std::collections::HashMap<&str, i32> = - reasons.iter().map(|k| (*k, 1)).collect(); - let has_mismatch = filter_reasons.contains_key("mismatch_max") - || filter_reasons.contains_key("mismatch_rate"); - let has_other = filter_reasons - .keys() - .any(|k| *k != "mismatch_max" && *k != "mismatch_rate"); + // Drives `FilterReasons` itself rather than re-implementing the + // decision: the previous version of this test kept its own copy of the + // logic, so it would have passed even if the shipped code had changed. + let classify = |reasons: &[usize]| -> UnmappedReason { + let mut filter_reasons = FilterReasons::default(); + for &r in reasons { + filter_reasons.add(r); + } + let has_mismatch = filter_reasons.has(FilterReasons::MISMATCH_MAX) + || filter_reasons.has(FilterReasons::MISMATCH_RATE); + let has_other = filter_reasons.has_non_mismatch(); if has_mismatch && !has_other { UnmappedReason::TooManyMismatches } else { @@ -1724,25 +1798,31 @@ mod tests { }; assert_eq!( - classify(&["mismatch_max"]), + classify(&[FilterReasons::MISMATCH_MAX]), UnmappedReason::TooManyMismatches ); assert_eq!( - classify(&["mismatch_rate"]), + classify(&[FilterReasons::MISMATCH_RATE]), UnmappedReason::TooManyMismatches ); assert_eq!( - classify(&["mismatch_max", "mismatch_rate"]), + classify(&[FilterReasons::MISMATCH_MAX, FilterReasons::MISMATCH_RATE]), UnmappedReason::TooManyMismatches ); // Mixed with score filter → TooShort assert_eq!( - classify(&["mismatch_max", "score_min"]), + classify(&[FilterReasons::MISMATCH_MAX, FilterReasons::SCORE_MIN]), UnmappedReason::TooShort ); // Score-only → TooShort - assert_eq!(classify(&["score_min"]), UnmappedReason::TooShort); - assert_eq!(classify(&["match_min"]), UnmappedReason::TooShort); + assert_eq!( + classify(&[FilterReasons::SCORE_MIN]), + UnmappedReason::TooShort + ); + assert_eq!( + classify(&[FilterReasons::MATCH_MIN]), + UnmappedReason::TooShort + ); } #[test] diff --git a/src/junction/mod.rs b/src/junction/mod.rs index 13776ce..040af01 100644 --- a/src/junction/mod.rs +++ b/src/junction/mod.rs @@ -17,7 +17,7 @@ use crate::params::Parameters; use crate::error::Error; use crate::genome::Genome; -use std::collections::HashMap; +use rustc_hash::FxHashMap; use std::path::Path; /// Key for junction lookup: (chr_idx, intron_start, intron_end, strand). @@ -54,14 +54,19 @@ pub struct NovelJunctionKey { #[derive(Clone)] pub struct SpliceJunctionDb { /// Map: (chr_idx, intron_start, intron_end, strand) → annotated - junctions: HashMap, + /// + /// `FxHashMap`, not the default `HashMap`: this is queried once per + /// candidate junction while stitching, and SipHash of a 25-byte key showed + /// up in the profile. Only `get` and `insert` are ever called on it and it + /// is never iterated, so the hasher cannot reach the output. + junctions: FxHashMap, } impl SpliceJunctionDb { /// Create empty database (for no-GTF mode) pub fn empty() -> Self { Self { - junctions: HashMap::new(), + junctions: FxHashMap::default(), } } @@ -95,7 +100,8 @@ impl SpliceJunctionDb { /// `TranscriptomeIndex` and the `sjdb_insert` pipeline without /// re-parsing the file. pub fn from_raw_junctions(raw: &[(usize, u64, u64, u8)]) -> Self { - let mut junctions = HashMap::with_capacity(raw.len()); + let mut junctions = + FxHashMap::with_capacity_and_hasher(raw.len(), rustc_hash::FxBuildHasher); for &(chr_idx, intron_start, intron_end, strand) in raw { let key = JunctionKey { chr_idx, @@ -411,7 +417,6 @@ mod tests { #[test] fn test_db_keyed_in_genome_absolute_zero_based_multi_chr() { use crate::junction::gtf::{GtfRecord, extract_junctions_configured}; - use std::collections::HashMap; // Two-chromosome toy genome so chr_start[1] != 0. let genome = Genome { @@ -426,7 +431,8 @@ mod tests { }; let make_exon = |seqname: &str, start: u64, end: u64, transcript: &str| -> GtfRecord { - let mut attrs: HashMap = HashMap::new(); + let mut attrs: std::collections::HashMap = + std::collections::HashMap::new(); attrs.insert("gene_id".to_string(), "G".to_string()); attrs.insert("transcript_id".to_string(), transcript.to_string()); GtfRecord { From 2e20e0d0306398aa9e48a69861e9850588994990 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 20:58:26 +0200 Subject: [PATCH 5/5] docs(changelog): record the filter-counter change Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index ddf5041..bb92cb8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -133,6 +133,10 @@ Sections commonly used: Features, Bug fixes, Other changes. Initial release of Rust rewrite of STAR. ### Other changes +- The per-read transcript-filter bookkeeping is a fixed set of counters + rather than a `HashMap<&str, i32>`, removing a heap allocation per read + and a string hash per filtered transcript. Output is unchanged. + - The splice-motif check in the junction scan carries a sliding window and looks the motif up in a table instead of re-reading four genome bases and matching on them at every position. Output is unchanged