From 98966abd53bc7c430acefa0c9ae2259bb7a0cb4a Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 14:57:37 +0200 Subject: [PATCH 01/17] feat(junction): sjA-addressable junction table + STAR binarySearch2 find() STAR's stitcher resolves an annotated junction by index (`sjA` / `sjdbInd`) into `mapGen.sjdb*`, then reads that entry's motif, shiftLeft and shiftRight. rustar-aligner had only a boolean `is_annotated()` over a HashMap, so neither the `sjAB` fast path nor annotated-boundary snapping could be expressed. Add the missing half: - `SpliceJunctionDb::table: Vec`, the annotated junctions in the exact order they occupy the Gsj buffer. - `find(stored_start, stored_end) -> Option`, a port of STAR's `binarySearch2`: binary search on the donor, then a backward/forward scan of the equal-donor run for the acceptor. Coordinates are genome-absolute, so no chromosome index is needed. - `entry(i)`, `table_len()`, `from_prepared()`, `set_table()`. No index format change. `PreparedJunction` already carries motif, shift_left, shift_right and strand; `sjdbInfo.txt` already writes and reads all four; and `sort_and_dedup` already establishes the `(stored_start, stored_end)` ordering the search needs. Existing indices work unchanged. The three population paths now agree: - `genomeGenerate` builds the db from the same `prepared` array it feeds to `build_gsj`, so table order and Gsj order match by construction. - The index-load path uses `from_prepared(prepared_junctions)`. - The align-time `--sjdbGTFfile` path derives motif/shift via `prepare_junction` + `sort_and_dedup`, which it previously discarded. Whenever the index carries junctions, its array wins as the table even if a GTF is also supplied at align time: `sj_a` tags come from `decode_gsj_hit`, which indexes that array, so anything else would make a tag address the wrong junction. Behaviour-neutral for canonical junctions, where stored coordinates equal the raw GTF ones. For non-canonical junctions the `genomeGenerate` in-memory db is now keyed on stored coordinates like the other two paths, which is what the stitch-time scan produces. Co-Authored-By: Claude Opus 5 (1M context) --- src/index/io.rs | 21 ++-- src/index/mod.rs | 6 +- src/junction/mod.rs | 243 +++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 259 insertions(+), 11 deletions(-) diff --git a/src/index/io.rs b/src/index/io.rs index ce6de64..2385f39 100644 --- a/src/index/io.rs +++ b/src/index/io.rs @@ -69,7 +69,7 @@ impl GenomeIndex { // db would be empty), losing all `sjdbScore` bonuses and annotated // junction recognition. Keyed on the stored (post-sjdbPrepare) donor/ // acceptor coordinates, matching what the stitch scan produces. - let junction_db = if let Some(ref gtf_path) = params.sjdb_gtf_file { + let mut junction_db = if let Some(ref gtf_path) = params.sjdb_gtf_file { SpliceJunctionDb::from_gtf_configured( gtf_path, &genome, @@ -78,20 +78,27 @@ impl GenomeIndex { ¶ms.sjdb_gtf_tag_exon_parent_transcript, )? } else if !prepared_junctions.is_empty() { - let raw: Vec<(usize, u64, u64, u8)> = prepared_junctions - .iter() - .map(|j| (j.chr_idx, j.stored_start(), j.stored_end(), j.strand)) - .collect(); log::info!( "No GTF at align time; loaded {} annotated junctions from index sjdbInfo.txt", - raw.len() + prepared_junctions.len() ); - SpliceJunctionDb::from_raw_junctions(&raw) + SpliceJunctionDb::from_prepared(prepared_junctions.clone()) } else { log::info!("No GTF file provided, all junctions will be novel"); SpliceJunctionDb::empty() }; + // `sj_a` tags come out of `decode_gsj_hit`, which indexes the junction + // array stored in the index. Whenever that array exists it must also be + // what `find()` searches, otherwise a tag would address a different + // junction than the one the SA hit came from — silently wrong CIGARs + // rather than a missed optimisation. An align-time GTF may rebuild the + // annotated lookup map, but never the table. + if !prepared_junctions.is_empty() { + junction_db.set_table(prepared_junctions.clone()); + } + let junction_db = junction_db; + log::info!( "Junction database loaded: {} annotated junctions", junction_db.len() diff --git a/src/index/mod.rs b/src/index/mod.rs index 88ee62d..bfccbf9 100644 --- a/src/index/mod.rs +++ b/src/index/mod.rs @@ -327,8 +327,6 @@ impl GenomeIndex { } let (junction_db, prepared_junctions) = if !raw.is_empty() { - let jdb = SpliceJunctionDb::from_raw_junctions(&raw); - let prepared: Vec = raw .iter() .map(|&(chr_idx, intron_start, intron_end, strand)| { @@ -343,6 +341,10 @@ impl GenomeIndex { }) .collect(); let prepared = sjdb_insert::sort_and_dedup(prepared); + // Build the annotated db from the prepared array itself, so the + // `sjA` table and the Gsj buffer are the same list in the same + // order by construction. + let jdb = SpliceJunctionDb::from_prepared(prepared.clone()); let gsj = sjdb_insert::build_gsj(&prepared, &genome, n_genome_real, params.sjdb_overhang)?; diff --git a/src/junction/mod.rs b/src/junction/mod.rs index 13776ce..5b9cca1 100644 --- a/src/junction/mod.rs +++ b/src/junction/mod.rs @@ -12,6 +12,7 @@ pub mod sjdb_insert; pub use sj_output::SpliceJunctionStats; pub(crate) use sj_output::{SjKey, encode_motif}; +pub use sjdb_insert::PreparedJunction; use crate::params::Parameters; @@ -55,6 +56,15 @@ pub struct NovelJunctionKey { pub struct SpliceJunctionDb { /// Map: (chr_idx, intron_start, intron_end, strand) → annotated junctions: HashMap, + /// STAR's `mapGen.sjdb*` arrays: the annotated junctions in the exact + /// order they occupy the Gsj buffer, sorted by `(stored_start, + /// stored_end)`. Index `i` here is STAR's `sjA` / `sjdbInd`, i.e. the + /// value [`sjdb_insert::decode_gsj_hit`] derives from a Gsj SA hit. + /// + /// Empty when the database was built without motif/shift metadata; in + /// that case [`find`](Self::find) always returns `None` and the + /// annotated-junction fast paths simply do not engage. + table: Vec, } impl SpliceJunctionDb { @@ -62,7 +72,102 @@ impl SpliceJunctionDb { pub fn empty() -> Self { Self { junctions: HashMap::new(), + table: Vec::new(), + } + } + + /// Build from an already sorted-and-deduplicated [`PreparedJunction`] + /// list — the form `sjdbInfo.txt` is read back into and the form + /// `genomeGenerate` produces before building the Gsj buffer. + /// + /// The HashMap is keyed on the *stored* (post-`sjdbPrepare`) donor and + /// acceptor coordinates, which is what the stitch-time scan produces. + pub fn from_prepared(prepared: Vec) -> Self { + let mut junctions = HashMap::with_capacity(prepared.len()); + for j in &prepared { + junctions.insert( + JunctionKey { + chr_idx: j.chr_idx, + intron_start: j.stored_start(), + intron_end: j.stored_end(), + strand: j.strand, + }, + JunctionInfo { annotated: true }, + ); + } + Self { + junctions, + table: prepared, + } + } + + /// Replace the `sjA`-addressable table without touching the annotated + /// lookup map. + /// + /// Needed because `sj_a` tags come from [`sjdb_insert::decode_gsj_hit`], + /// which indexes the junction list stored in the index (`sjdbInfo.txt`). + /// When a GTF is *also* supplied at align time the map is rebuilt from + /// that GTF, but the table must keep addressing the index's array or the + /// tags would point at the wrong junctions. + pub fn set_table(&mut self, prepared: Vec) { + self.table = prepared; + } + + /// STAR's `binarySearch2` (`ReadAlign_stitchAlignToTranscript.cpp` → + /// `binarySearch2.h`): the index of the junction whose stored donor and + /// acceptor are exactly `x` and `y`, or `None`. + /// + /// Coordinates are genome-absolute, so no chromosome index is needed: + /// STAR's `sjdbStart` / `sjdbEnd` are already unique across contigs. + pub fn find(&self, x: u64, y: u64) -> Option { + let n = self.table.len(); + if n == 0 || x > self.table[n - 1].stored_start() || x < self.table[0].stored_start() { + return None; } + let (mut i1, mut i2) = (0usize, n - 1); + while i2 > i1 + 1 { + let i3 = usize::midpoint(i1, i2); + if self.table[i3].stored_start() > x { + i2 = i3; + } else { + i1 = i3; + } + } + let i3 = if x == self.table[i1].stored_start() { + i1 + } else if x == self.table[i2].stored_start() { + i2 + } else { + return None; + }; + // Scan the run of equal `stored_start` values (backward then forward) + // for a matching `stored_end`. + for jj in (0..=i3).rev() { + if x != self.table[jj].stored_start() { + break; + } else if y == self.table[jj].stored_end() { + return Some(jj); + } + } + for jj in i3..n { + if x != self.table[jj].stored_start() { + return None; + } else if y == self.table[jj].stored_end() { + return Some(jj); + } + } + None + } + + /// The junction at `sjA` index `i`, or `None` when the table is absent + /// or the index is out of range. + pub fn entry(&self, i: usize) -> Option<&PreparedJunction> { + self.table.get(i) + } + + /// Number of entries in the `sjA`-addressable table. + pub fn table_len(&self) -> usize { + self.table.len() } /// Build junction database from GTF file with configurable GTF attribute names. @@ -81,7 +186,36 @@ impl SpliceJunctionDb { let raw = gtf::extract_junctions_configured(exons, genome, transcript_tag)?; log::info!("Extracted {} annotated junctions from GTF", raw.len()); - Ok(Self::from_raw_junctions(&raw)) + // Keep the historical (raw-coordinate) annotated lookup map, but also + // derive the motif/shift/strand table so the annotated-junction fast + // paths have something to address. Building the table here makes the + // align-time GTF path carry the same metadata the index-loaded path + // already had. + let mut db = Self::from_raw_junctions(&raw); + db.set_table(Self::prepare_table(&raw, genome)); + Ok(db) + } + + /// Run every raw `(chr_idx, intron_start, intron_end, strand)` through + /// `sjdbPrepare`'s motif detection and micro-repeat shift computation, + /// then apply STAR's post-dedup sort. Produces exactly the array + /// `genomeGenerate` writes to `sjdbInfo.txt`. + fn prepare_table(raw: &[(usize, u64, u64, u8)], genome: &Genome) -> Vec { + let n_genome_real = genome.n_genome_real; + let prepared: Vec = raw + .iter() + .map(|&(chr_idx, intron_start, intron_end, strand)| { + sjdb_insert::prepare_junction( + chr_idx, + intron_start, + intron_end, + strand, + genome, + n_genome_real, + ) + }) + .collect(); + sjdb_insert::sort_and_dedup(prepared) } /// Build junction database from GTF file (default STAR attribute names). @@ -94,6 +228,10 @@ impl SpliceJunctionDb { /// the `genomeGenerate` path so it can share the parsed GTF with /// `TranscriptomeIndex` and the `sjdb_insert` pipeline without /// re-parsing the file. + /// + /// The resulting database has **no** `sjA` table; callers that need + /// [`find`](Self::find) must follow up with [`set_table`](Self::set_table) + /// or use [`from_prepared`](Self::from_prepared) instead. pub fn from_raw_junctions(raw: &[(usize, u64, u64, u8)]) -> Self { let mut junctions = HashMap::with_capacity(raw.len()); for &(chr_idx, intron_start, intron_end, strand) in raw { @@ -105,7 +243,10 @@ impl SpliceJunctionDb { }; junctions.insert(key, JunctionInfo { annotated: true }); } - Self { junctions } + Self { + junctions, + table: Vec::new(), + } } /// Check if a junction is annotated in the GTF. @@ -408,6 +549,104 @@ mod tests { assert_eq!(novel_junctions[0].0.intron_start, 300); } + /// Build a canonical (motif != 0) prepared junction whose stored + /// coordinates are exactly `(start, end)`. + fn pj(start: u64, end: u64, strand: u8) -> PreparedJunction { + PreparedJunction { + chr_idx: 0, + start_pos: start, + end_pos: end, + motif: 1, + shift_left: 0, + shift_right: 0, + strand, + } + } + + #[test] + fn find_round_trips_every_prepared_junction() { + // Deliberately unsorted input: `sort_and_dedup` establishes the + // (stored_start, stored_end) order `find`'s binary search needs. + let prepared = sjdb_insert::sort_and_dedup(vec![ + pj(900, 1000, 1), + pj(100, 200, 1), + pj(500, 600, 2), + pj(300, 400, 1), + ]); + let db = SpliceJunctionDb::from_prepared(prepared.clone()); + assert_eq!(db.table_len(), 4); + + // The invariant the sjAB fast path and the annotated snap both rely + // on: index i in the table is addressable by its own stored coords. + for (i, j) in prepared.iter().enumerate() { + assert_eq!( + db.find(j.stored_start(), j.stored_end()), + Some(i), + "junction {i} did not round-trip" + ); + assert_eq!(db.entry(i).unwrap().stored_start(), j.stored_start()); + } + + // Misses in every direction: below the first, above the last, a + // start that exists with a wrong end, and an end that exists with a + // wrong start. + assert_eq!(db.find(50, 200), None); + assert_eq!(db.find(2000, 3000), None); + assert_eq!(db.find(100, 201), None); + assert_eq!(db.find(101, 200), None); + } + + #[test] + fn find_disambiguates_a_run_of_equal_starts() { + // Several junctions sharing a donor: the backward/forward scan around + // the binary-search landing point must pick the right acceptor. + let prepared = sjdb_insert::sort_and_dedup(vec![ + pj(100, 700, 1), + pj(100, 200, 1), + pj(100, 500, 1), + pj(100, 300, 1), + pj(900, 1000, 1), + ]); + let db = SpliceJunctionDb::from_prepared(prepared.clone()); + for (i, j) in prepared.iter().enumerate() { + assert_eq!(db.find(j.stored_start(), j.stored_end()), Some(i)); + } + assert_eq!(db.find(100, 400), None); + } + + #[test] + fn find_is_inert_without_a_table() { + // `from_raw_junctions` carries no motif/shift metadata, so the + // annotated fast paths must simply not engage rather than misfire. + let db = SpliceJunctionDb::from_raw_junctions(&[(0, 100, 200, 1)]); + assert!(db.is_annotated(0, 100, 200, 1)); + assert_eq!(db.table_len(), 0); + assert_eq!(db.find(100, 200), None); + assert!(db.entry(0).is_none()); + } + + #[test] + fn from_prepared_keys_the_map_on_stored_coordinates() { + // Non-canonical junction: stored coords are the shifted ones, which is + // what the stitch-time scan produces. + let noncan = PreparedJunction { + chr_idx: 0, + start_pos: 100, + end_pos: 200, + motif: 0, + shift_left: 3, + shift_right: 0, + strand: 0, + }; + assert_eq!(noncan.stored_start(), 100); + assert_eq!(noncan.original_start(), 103); + + let db = SpliceJunctionDb::from_prepared(vec![noncan]); + assert!(db.is_annotated(0, 100, 200, 0)); + assert!(!db.is_annotated(0, 103, 203, 0)); + assert_eq!(db.find(100, 200), Some(0)); + } + #[test] fn test_db_keyed_in_genome_absolute_zero_based_multi_chr() { use crate::junction::gtf::{GtfRecord, extract_junctions_configured}; From ba5ecc549c6bcddafa6d522ec1f0d6f6d8f3071c Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 15:09:17 +0200 Subject: [PATCH 02/17] feat(align): carry sjA tags through the window and take STAR's annotated-junction shortcut STAR tags every seed derived from the sjdb genome insert with the index of the junction it came from (`sjA`), then, when the previous exon and the incoming piece carry the *same* tag and are read-adjacent, stitches them straight onto the annotated junction instead of searching for a boundary. rustar-aligner computed that index and threw it away. Plumbing: - `decode_gsj_hit` returns the `junction_idx` it already computed. Both halves of a boundary-crossing hit carry it: they are two flanks of one junction. - `expand_hit`, `WinCandidate`, `WindowAlignment` and `ExonBlock` carry `sj_a`, `-1` for ordinary real-genome hits. - The overlap dedup in `cluster_seeds` now compares `sj_a` alongside `mate_id`, matching STAR's `e.sjA == sjA` guard in `assignAlignToWindow`. Seeds derived from different annotated junctions are no longer collapsed just because they share a diagonal. The shortcut itself, at the head of `stitch_align_to_transcript`: - fires when both sides carry the same non-`-1` tag, sit on the same mate, are exactly read-adjacent, and are separated on the genome; - takes exon B verbatim at the annotated boundary, adopts the stored motif and shifts, marks the junction annotated, and adds `sjdbScore`; - rejects (STAR's -1000006, `None` here) when micro-repeats around a non-canonical annotated junction are larger than the pieces being stitched. Without it, the general boundary search flushes a repeat-adjacent junction to its leftmost position and lands a few bases off the annotated coordinates (STAR's `11M1545N38M` becoming `9M1545N40M`), which also forfeits the annotated bonus because the shifted coordinates no longer match the database. Receiver coordinates are half-open where STAR's are inclusive, so `rBstart == rAend + 1` is `wa.read_pos == last_exon.read_end` and `gAend + 1 < gBstart` is `last_exon.genome_end < wa.sa_pos`. Co-Authored-By: Claude Opus 5 (1M context) --- src/align/stitch.rs | 383 +++++++++++++++++++++++++++++++----- src/chimeric/detect.rs | 1 + src/junction/mod.rs | 2 +- src/junction/sj_output.rs | 15 ++ src/junction/sjdb_insert.rs | 22 ++- 5 files changed, 367 insertions(+), 56 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 1c25ecc..8bf4cd6 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -412,8 +412,19 @@ pub struct WindowAlignment { /// = length + left_ext_score + right_ext_score (in stitch coords, forward strand). /// Computed in stitch_seeds_core after dedup/sort. Default: length as i32. pub pre_ext_score: i32, + /// Annotated-junction index this seed came from (STAR: `WA_sjA`), or `-1` + /// for an ordinary real-genome hit. Set when the SA hit landed in the Gsj + /// flanking buffer, which means the seed straddles a known junction. + /// Addresses [`SpliceJunctionDb::entry`](crate::junction::SpliceJunctionDb::entry). + pub sj_a: i64, } +/// One SA hit after Gsj expansion: `(forward_pos, read_offset, length, +/// sa_pos, sj_a)`. A hit in the real genome yields exactly one of these, with +/// `sj_a = -1`; a hit straddling the donor/acceptor boundary of the Gsj buffer +/// yields two, both tagged with the junction they came from. +type ExpandedHit = (u64, usize, usize, u64, i64); + /// A cluster of seeds mapping to the same genomic region #[derive(Debug, Clone)] pub struct SeedCluster { @@ -477,7 +488,7 @@ pub fn cluster_seeds( let prepared = index.prepared_junctions.as_slice(); // Expand one raw SA hit into the candidate `(real_fwd_pos, read_offset, - // sub_length, sub_sa_pos)` tuples cluster_seeds should consume. Hits + // sub_length, sub_sa_pos, sj_a)` tuples cluster_seeds should consume. Hits // in the real genome pass through unchanged. Hits in the Gsj // flanking buffer are decoded via the sjdb table - yielding one // entry for hits confined to a single flank, two entries for hits @@ -488,45 +499,46 @@ pub fn cluster_seeds( // `[Option; 2]` avoids a heap `Vec` allocation on every SA hit expanded // in this hot loop (the common `raw_fwd < n_genome_real` case is the // overwhelming majority of calls); iterate with `.into_iter().flatten()`. - let expand_hit = - |sa_pos: u64, strand: bool, length: usize| -> [Option<(u64, usize, usize, u64)>; 2] { - let raw_fwd = index.sa_pos_to_forward(sa_pos, strand, length); - if raw_fwd < n_genome_real { - return [Some((raw_fwd, 0, length, sa_pos)), None]; - } - if sjdb_overhang == 0 || prepared.is_empty() { - return [None, None]; - } - let mut decoded = crate::junction::sjdb_insert::decode_gsj_hit( - raw_fwd, - length, - n_genome_real, - sjdb_overhang, - prepared, - ); - // Reverse-strand hits traverse the donor/acceptor halves in reverse - // read order: the leftmost forward bytes (donor flank) align to the - // last read bases. Swap the read offsets so each sub-seed's - // `read_pos = seed.read_pos + read_offset` lands at the right place - // in original-read coords. - if strand && decoded.len() == 2 { - let acceptor_len = decoded[1].2; - decoded[0].1 = acceptor_len; - decoded[1].1 = 0; - } - let mut out = [None, None]; - for (slot, (real_fwd, read_off, sub_len)) in out.iter_mut().zip(decoded) { - let sub_sa_pos = if strand { - n_genome - .saturating_sub(real_fwd) - .saturating_sub(sub_len as u64) - } else { - real_fwd - }; - *slot = Some((real_fwd, read_off, sub_len, sub_sa_pos)); - } - out - }; + let expand_hit = |sa_pos: u64, strand: bool, length: usize| -> [Option; 2] { + let raw_fwd = index.sa_pos_to_forward(sa_pos, strand, length); + if raw_fwd < n_genome_real { + // A real-genome hit carries no annotated-junction tag (STAR's + // sjA = -1). + return [Some((raw_fwd, 0, length, sa_pos, -1)), None]; + } + if sjdb_overhang == 0 || prepared.is_empty() { + return [None, None]; + } + let mut decoded = crate::junction::sjdb_insert::decode_gsj_hit( + raw_fwd, + length, + n_genome_real, + sjdb_overhang, + prepared, + ); + // Reverse-strand hits traverse the donor/acceptor halves in reverse + // read order: the leftmost forward bytes (donor flank) align to the + // last read bases. Swap the read offsets so each sub-seed's + // `read_pos = seed.read_pos + read_offset` lands at the right place + // in original-read coords. + if strand && decoded.len() == 2 { + let acceptor_len = decoded[1].2; + decoded[0].1 = acceptor_len; + decoded[1].1 = 0; + } + let mut out = [None, None]; + for (slot, (real_fwd, read_off, sub_len, sj_a)) in out.iter_mut().zip(decoded) { + let sub_sa_pos = if strand { + n_genome + .saturating_sub(real_fwd) + .saturating_sub(sub_len as u64) + } else { + real_fwd + }; + *slot = Some((real_fwd, read_off, sub_len, sub_sa_pos, sj_a as i64)); + } + out + }; let anchor_set: Vec = seeds .iter() @@ -595,7 +607,7 @@ pub fn cluster_seeds( continue; } - for (forward_pos, _read_off, sub_length, _sub_sa_pos) in + for (forward_pos, _read_off, sub_length, _sub_sa_pos, _sj_a) in expand_hit(sa_pos, strand, full_length) .into_iter() .flatten() @@ -766,6 +778,8 @@ pub fn cluster_seeds( /// by the donor-side length for the acceptor half of a Gsj /// boundary-crossing hit. read_pos: usize, + /// Annotated-junction index (STAR: `sjA`), or `-1`. + sj_a: i64, } let mut win_candidates: Vec> = @@ -784,7 +798,7 @@ pub fn cluster_seeds( continue; } - for (forward_pos, read_off, sub_length, sub_sa_pos) in + for (forward_pos, read_off, sub_length, sub_sa_pos, sj_a) in expand_hit(sa_pos, strand, full_length) .into_iter() .flatten() @@ -819,6 +833,7 @@ pub fn cluster_seeds( ps_rstart, mate_id: seed.mate_id, read_pos: derived_read_pos, + sj_a, }); } } @@ -847,7 +862,7 @@ pub fn cluster_seeds( // + hash-table rehashes per read (a measured allocator hotspot). Reused // storage only; the results are identical. let mut by_len: Vec = Vec::new(); - let mut diag_ranges: FxHashMap<(i64, u8), Vec<(usize, usize)>> = FxHashMap::default(); + let mut diag_ranges: FxHashMap<(i64, u8, i64), Vec<(usize, usize)>> = FxHashMap::default(); for win_idx in 0..win_n { let candidates = &win_candidates[win_idx]; @@ -860,15 +875,18 @@ pub fn cluster_seeds( by_len.extend(0..candidates.len()); by_len.sort_by(|&a, &b| candidates[b].length.cmp(&candidates[a].length)); - // For each (diagonal, mate_id) pair, track accepted [ps_rstart, ps_rend) ranges. - // STAR's assignAlignToWindow checks aFrag==WA[iA][WA_iFrag] before overlap test: - // seeds from different fragments are never treated as overlapping duplicates. + // For each (diagonal, mate_id, sj_a) triple, track accepted + // [ps_rstart, ps_rend) ranges. STAR's assignAlignToWindow checks both + // `aFrag == WA[iA][WA_iFrag]` and `sjA == WA[iA][WA_sjA]` before the + // overlap test: seeds from different fragments, or derived from + // different annotated junctions, are never treated as duplicates even + // when they land on the same diagonal. diag_ranges.clear(); for &ci in &by_len { let cand = &candidates[ci]; let diag = cand.forward_pos as i64 - cand.ps_rstart as i64; let ps_rend = cand.ps_rstart + cand.length; - let key = (diag, cand.mate_id); + let key = (diag, cand.mate_id, cand.sj_a); let blocked = diag_ranges.get(&key).is_some_and(|ranges| { ranges.iter().any(|&(rs, re)| { @@ -930,6 +948,9 @@ pub fn cluster_seeds( if wa.mate_id != new_mate_id { continue; // STAR: only merge seeds from the same fragment } + if wa.sj_a != cand.sj_a { + continue; // STAR: `e.sjA == sjA` guards the overlap test + } let wa_ps_rstart = if window.is_reverse { read_len - (wa.length + wa.read_pos) } else { @@ -968,6 +989,7 @@ pub fn cluster_seeds( is_anchor: is_anchor_seed, mate_id: seed.mate_id, pre_ext_score: length as i32, + sj_a: cand.sj_a, }, ); } @@ -1036,6 +1058,7 @@ pub fn cluster_seeds( is_anchor: is_anchor_seed, mate_id: seed.mate_id, pre_ext_score: length as i32, + sj_a: cand.sj_a, }, ); } @@ -1076,6 +1099,11 @@ pub(crate) struct ExonBlock { pub(crate) genome_end: u64, // SA coordinate space (exclusive) /// Mate ID: 0=mate1, 1=mate2, 2=SE (STAR: EX_iFrag) pub(crate) mate_id: u8, + /// Annotated-junction index of the seed this exon was built from + /// (STAR: `EX_sjA`), or `-1`. Compared against the incoming seed's `sj_a` + /// to recognise that two exons are the two flanks of the same annotated + /// junction, which lets the stitcher take that junction verbatim. + pub(crate) sj_a: i64, } /// In-progress transcript during recursive search (cheap to clone) @@ -1136,6 +1164,68 @@ fn stitch_align_to_transcript( ) -> Option { let last_exon = wt.exons.last().unwrap(); + // STAR's "simple stitching if junction belongs to a database" + // (`stitchAlignToTranscript.cpp`). When the previous exon and the incoming + // piece B carry the *same* annotated-junction tag, sit on the same mate, + // and are exactly read-adjacent, the two are the donor and acceptor flanks + // of one annotated junction that the Gsj insert split apart. Stitch them + // straight onto the annotated boundary: exact donor/acceptor coordinates, + // the stored motif and shifts, and the `sjdbScore` bonus, skipping the + // general boundary search entirely. + // + // Without this, the general search flushes a repeat-adjacent junction to + // its leftmost position and lands a few bases off the annotated boundary + // (e.g. `9M1545N40M` where STAR emits the annotated `11M1545N38M`), and + // forfeits the bonus because the shifted coordinates no longer match the + // database. + // + // Receiver coordinates are half-open where STAR's are inclusive, hence + // `wa.read_pos == last_exon.read_end` for STAR's `rBstart == rAend + 1` + // and `last_exon.genome_end < wa.sa_pos` for `gAend + 1 < gBstart`. + if wa.sj_a != -1 + && last_exon.sj_a == wa.sj_a + && last_exon.mate_id == wa.mate_id + && wa.read_pos == last_exon.read_end + && last_exon.genome_end < wa.sa_pos + && let Some(db) = junction_db + && let Some(pj) = db.entry(wa.sj_a as usize) + { + // Too-large repeats around a non-canonical annotated junction: STAR + // rejects with -1000006, which is a `None` here — the caller still + // takes the EXCLUDE branch, exactly as STAR's recursion does. + let exon_a_len = last_exon.read_end - last_exon.read_start; + if pj.motif == 0 + && (wa.length <= pj.shift_right as usize || exon_a_len <= pj.shift_left as usize) + { + return None; + } + + let mut new_wt = wt.clone(); + new_wt.exons.push(ExonBlock { + read_start: wa.read_pos, + read_end: wa.read_pos + wa.length, + genome_start: wa.sa_pos, + genome_end: wa.sa_pos + wa.length as u64, + mate_id: wa.mate_id, + sj_a: wa.sj_a, + }); + new_wt.n_junction += 1; + new_wt + .junction_motifs + .push(crate::junction::decode_motif(pj.motif)); + new_wt.junction_annotated.push(true); + new_wt + .junction_shifts + .push((pj.shift_left as u32, pj.shift_right as u32)); + new_wt.score += wa.length as i32 + scorer.sjdb_score; + new_wt.read_end = wa.read_pos + wa.length; + new_wt.genome_end = wa.sa_pos + wa.length as u64; + if wa.is_anchor { + new_wt.n_anchor += 1; + } + return Some(new_wt); + } + // Mate-boundary detection: STAR canonSJ[iex] = -3 (stitchAlignToTranscript.cpp:402) // When crossing from mate1 to mate2 (or vice versa), skip junction scoring and // check alignMatesGapMax instead. @@ -1230,6 +1320,7 @@ fn stitch_align_to_transcript( genome_start: wa.sa_pos, genome_end: wa.sa_pos + wa.length as u64, mate_id: wa.mate_id, + sj_a: wa.sj_a, }); new_wt.read_end = wa.read_pos + wa.length; new_wt.genome_end = wa.sa_pos + wa.length as u64; @@ -1522,6 +1613,7 @@ fn stitch_align_to_transcript( genome_start: b_genome_start, genome_end: b_genome_start + b_len as u64, mate_id: wa.mate_id, + sj_a: wa.sj_a, }); } else { // Insertion: read_gap > genome_gap @@ -1631,6 +1723,7 @@ fn stitch_align_to_transcript( genome_start: b_genome_start, genome_end: eff_genome_pos + eff_length as u64, mate_id: wa.mate_id, + sj_a: wa.sj_a, }); } @@ -2407,6 +2500,7 @@ fn stitch_recurse( genome_start: wa.sa_pos, genome_end: wa.sa_pos + wa.length as u64, mate_id: wa.mate_id, + sj_a: wa.sj_a, }); new_wt.score = wa.length as i32; new_wt.read_start = wa.read_pos; @@ -3245,6 +3339,7 @@ mod tests { is_anchor: true, mate_id: 2, pre_ext_score: 5, + sj_a: -1, }, WindowAlignment { seed_idx: 1, @@ -3256,6 +3351,7 @@ mod tests { is_anchor: true, mate_id: 2, pre_ext_score: 5, + sj_a: -1, }, ]; @@ -3698,4 +3794,197 @@ mod tests { let replacement_correct = baseline + scorer.sjdb_score; assert_ne!(additive_buggy, replacement_correct); } + + // ----------------------------------------------------------------- + // sjAB fast path (STAR "simple stitching if junction belongs to a + // database") + // ----------------------------------------------------------------- + + fn sjab_cluster() -> SeedCluster { + SeedCluster { + alignments: Vec::new(), + chr_idx: 0, + genome_start: 0, + genome_end: 400, + is_reverse: false, + anchor_idx: 0, + anchor_bin: 0, + } + } + + /// A working transcript holding one exon that came from junction `sj_a`, + /// covering read `[0, exon_len)` at genome `[0, exon_len)`. + fn sjab_wt(exon_len: usize, sj_a: i64) -> WorkingTranscript { + let mut wt = WorkingTranscript::new(); + wt.exons.push(ExonBlock { + read_start: 0, + read_end: exon_len, + genome_start: 0, + genome_end: exon_len as u64, + mate_id: 2, + sj_a, + }); + wt.score = exon_len as i32; + wt.read_end = exon_len; + wt.genome_end = exon_len as u64; + wt + } + + /// The acceptor-side seed of the same junction: read-adjacent to the exon + /// above, but far away on the genome. + fn sjab_wa(read_pos: usize, length: usize, sa_pos: u64, sj_a: i64) -> WindowAlignment { + WindowAlignment { + seed_idx: 1, + read_pos, + length, + genome_pos: sa_pos, + sa_pos, + n_rep: 1, + is_anchor: false, + mate_id: 2, + pre_ext_score: length as i32, + sj_a, + } + } + + fn sjab_db(motif: u8, shift_left: u8, shift_right: u8) -> crate::junction::SpliceJunctionDb { + crate::junction::SpliceJunctionDb::from_prepared(vec![crate::junction::PreparedJunction { + chr_idx: 0, + start_pos: 20, + end_pos: 199, + motif, + shift_left, + shift_right, + strand: 1, + }]) + } + + #[test] + fn sjab_shortcut_uses_annotated_junction() { + use crate::align::score::{AlignmentScorer, SpliceMotif}; + + let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]); + let scorer = AlignmentScorer::from_params_minimal(); + let db = sjab_db(3, 0, 0); // GC/AG, no micro-repeat + let cluster = sjab_cluster(); + + let wt = sjab_wt(20, 0); + let wa = sjab_wa(20, 30, 200, 0); + let out = stitch_align_to_transcript( + &wt, + &wa, + &[0u8; 50], + &index, + &scorer, + &cluster, + Some(&db), + 0, + "t", + ) + .expect("annotated shortcut should stitch"); + + // Exon B is taken verbatim at the annotated boundary — no shifting. + assert_eq!(out.exons.len(), 2); + assert_eq!(out.exons[1].read_start, 20); + assert_eq!(out.exons[1].genome_start, 200); + assert_eq!(out.exons[1].sj_a, 0); + + // The junction is recorded as annotated, with the *stored* motif + // rather than one re-derived from the genome. + assert_eq!(out.n_junction, 1); + assert_eq!(out.junction_annotated, vec![true]); + assert_eq!(out.junction_motifs, vec![SpliceMotif::GcAg]); + + // Score is match-per-base plus the annotated bonus, and specifically + // not the GC/AG motif penalty the general path would have applied. + assert_eq!(out.score, 20 + 30 + scorer.sjdb_score); + assert_ne!(out.score, 20 + 30 + scorer.score_gap_gcag); + } + + #[test] + fn sjab_shortcut_requires_matching_tags_and_adjacency() { + use crate::align::score::AlignmentScorer; + + let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]); + let scorer = AlignmentScorer::from_params_minimal(); + let db = sjab_db(3, 0, 0); + let cluster = sjab_cluster(); + let read = [0u8; 50]; + + // The fast path's fingerprint is that it takes the *stored* motif and + // leaves exon B exactly where the annotation puts it. The general path + // can also mark this junction annotated (it does its own coordinate + // lookup), so `junction_annotated` alone would not tell them apart. + let signature = |wt: &WorkingTranscript, wa: &WindowAlignment| { + stitch_align_to_transcript(wt, wa, &read, &index, &scorer, &cluster, Some(&db), 0, "t") + .map(|t| (t.junction_motifs[0], t.exons[1].genome_start)) + }; + let fast_path = Some((crate::align::score::SpliceMotif::GcAg, 200u64)); + + // Baseline: the shortcut fires. + assert_eq!( + signature(&sjab_wt(20, 0), &sjab_wa(20, 30, 200, 0)), + fast_path + ); + + // Untagged seed: no shortcut, so the stored GC/AG motif is not adopted. + assert_ne!( + signature(&sjab_wt(20, 0), &sjab_wa(20, 30, 200, -1)), + fast_path + ); + // Tags disagree: the two pieces are flanks of *different* junctions. + assert_ne!( + signature(&sjab_wt(20, 0), &sjab_wa(20, 30, 200, 1)), + fast_path + ); + // Not read-adjacent: a gap in the read is not what the Gsj split makes. + assert_ne!( + signature(&sjab_wt(20, 0), &sjab_wa(21, 30, 200, 0)), + fast_path + ); + } + + #[test] + fn too_large_repeat_around_annotated_junction_rejects() { + use crate::align::score::AlignmentScorer; + + let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]); + let scorer = AlignmentScorer::from_params_minimal(); + let cluster = sjab_cluster(); + let read = [0u8; 50]; + + // Non-canonical annotated junction whose micro-repeats swallow the + // pieces being stitched: STAR rejects with -1000006, we return None. + let db = sjab_db(0, 40, 40); + let out = stitch_align_to_transcript( + &sjab_wt(20, 0), + &sjab_wa(20, 30, 200, 0), + &read, + &index, + &scorer, + &cluster, + Some(&db), + 0, + "t", + ); + assert!( + out.is_none(), + "a seed shorter than shiftRight must reject the annotated shortcut" + ); + + // Same junction with small repeats: accepted. + let db_ok = sjab_db(0, 1, 1); + let out_ok = stitch_align_to_transcript( + &sjab_wt(20, 0), + &sjab_wa(20, 30, 200, 0), + &read, + &index, + &scorer, + &cluster, + Some(&db_ok), + 0, + "t", + ); + assert!(out_ok.is_some()); + } } diff --git a/src/chimeric/detect.rs b/src/chimeric/detect.rs index 4997049..33d564c 100644 --- a/src/chimeric/detect.rs +++ b/src/chimeric/detect.rs @@ -1112,6 +1112,7 @@ mod tests { is_anchor: true, mate_id: 2, pre_ext_score: (genome_end - genome_start) as i32, + sj_a: -1, }], chr_idx, genome_start, diff --git a/src/junction/mod.rs b/src/junction/mod.rs index 5b9cca1..2761e80 100644 --- a/src/junction/mod.rs +++ b/src/junction/mod.rs @@ -11,7 +11,7 @@ mod sj_output; pub mod sjdb_insert; pub use sj_output::SpliceJunctionStats; -pub(crate) use sj_output::{SjKey, encode_motif}; +pub(crate) use sj_output::{SjKey, decode_motif, encode_motif}; pub use sjdb_insert::PreparedJunction; use crate::params::Parameters; diff --git a/src/junction/sj_output.rs b/src/junction/sj_output.rs index 9fa41dd..65bc577 100644 --- a/src/junction/sj_output.rs +++ b/src/junction/sj_output.rs @@ -401,6 +401,21 @@ pub(crate) fn encode_motif(motif: SpliceMotif) -> u8 { } } +/// Inverse of [`encode_motif`]: turn a stored `sjdbMotif` code back into a +/// [`SpliceMotif`]. Out-of-range codes map to `NonCanonical`, matching +/// [`SpliceMotif::filter_category_from_encoded`]. +pub(crate) fn decode_motif(encoded: u8) -> SpliceMotif { + match encoded { + 1 => SpliceMotif::GtAg, + 2 => SpliceMotif::CtAc, + 3 => SpliceMotif::GcAg, + 4 => SpliceMotif::CtGc, + 5 => SpliceMotif::AtAc, + 6 => SpliceMotif::GtAt, + _ => SpliceMotif::NonCanonical, + } +} + #[cfg(test)] mod tests { use crate::params::Parameters; diff --git a/src/junction/sjdb_insert.rs b/src/junction/sjdb_insert.rs index e358622..5b04b17 100644 --- a/src/junction/sjdb_insert.rs +++ b/src/junction/sjdb_insert.rs @@ -140,7 +140,8 @@ pub fn read_sjdb_info_tab(path: &Path, genome: &Genome) -> Result Result Vec<(u64, usize, usize)> { +) -> Vec<(u64, usize, usize, usize)> { if fwd_pos < n_genome_real { return Vec::new(); } @@ -200,8 +206,8 @@ pub fn decode_gsj_hit( let donor_real = donor_genome_start + slot_offset; let acceptor_real = acceptor_genome_start; vec![ - (donor_real, 0, donor_len), - (acceptor_real, donor_len, acceptor_len), + (donor_real, 0, donor_len, junction_idx), + (acceptor_real, donor_len, acceptor_len, junction_idx), ] } } @@ -1052,9 +1058,9 @@ mod tests { let out = decode_gsj_hit(hit_pos, 12, n_genome_real, overhang, &junctions); assert_eq!(out.len(), 2); // Donor: real_fwd = 990 + 6 = 996, read_offset = 0, sub_len = 4. - assert_eq!(out[0], (996, 0, 4)); + assert_eq!(out[0], (996, 0, 4, 0)); // Acceptor: real_fwd = 2001, read_offset = 4, sub_len = 8. - assert_eq!(out[1], (2001, 4, 8)); + assert_eq!(out[1], (2001, 4, 8, 0)); } #[test] @@ -1086,9 +1092,9 @@ mod tests { let out = decode_gsj_hit(hit_pos, 10, n_genome_real, overhang, &junctions); assert_eq!(out.len(), 2); // original_start = 1003 → donor flank [993..1003). Hit at slot 5 = 993+5 = 998. - assert_eq!(out[0], (998, 0, 5)); + assert_eq!(out[0], (998, 0, 5, 0)); // original_end = 2003 → acceptor flank starts at 2004. - assert_eq!(out[1], (2004, 5, 5)); + assert_eq!(out[1], (2004, 5, 5, 0)); } #[test] From e2564023e2080ea9a7fcb214e39bca06ade7ea87 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 15:14:56 +0200 Subject: [PATCH 03/17] feat(align): snap non-canonical annotated junctions onto their annotated coordinates MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Closes item #4b of STAR-RS-COMPARISON.md §7.2, which §7.4 calls "the best remaining lever". The splice branch previously asked the annotation a yes/no question and, on a hit, swapped the motif penalty for `sjdbScore`. It never used what the annotation actually says about *where* the junction is. STAR does three things with the entry it finds: 1. adopts the stored motif (`jcan = sjdbMotif[ind]`), so the mismatch gate, `outFilterIntronMotifs` and SJ.out.tab all see the annotated motif rather than one re-derived from the genome; 2. for a non-canonical entry, shifts `jR` by `shiftLeft`, moving the boundary off the leftmost-flush position the scan picked and onto the annotated one; 3. rejects (-1000006) when the micro-repeat is larger than either piece, or when the shift would push the donor past the end of piece B. All three are now done. Because `jr_shift` is consumed afterwards to move exon A's right edge and exon B's left edge together, shifting it relocates the whole junction coherently and leaves the intron length untouched. The mismatch gate moves after the lookup, matching STAR's ordering: it must test the annotated motif, not the scanned one. This changes behaviour on annotated non-canonical junctions, which is the point. `find()` and `is_annotated()` are queried separately rather than one wrapping the other: junctions inserted by two-pass mode live only in the map and carry no motif or shift, so they must still count as annotated while contributing nothing to snap. New tests: `noncanonical_annotated_junction_snaps_to_annotated_coords` measures where the unannotated search lands, annotates exactly that junction, and asserts both boundaries move by `shiftLeft` while the intron keeps its length (and that a canonical entry at the same coordinates is left alone). Co-Authored-By: Claude Opus 5 (1M context) --- src/align/stitch.rs | 173 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 147 insertions(+), 26 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 8bf4cd6..68d8ffb 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -1461,23 +1461,24 @@ fn stitch_align_to_transcript( // is motif detection (splice) vs pure positional score (deletion). // donor_sa = exclusive end of exon A = STAR's gAend+1. jr_shift = STAR's jR. let donor_sa = last_exon.genome_end; - let (jr_shift, motif, motif_score, jj_l, jj_r) = scorer.find_best_junction_position( - read_seq, - last_exon.read_end, - donor_sa, - read_gap.max(0), - genome_gap, - &index.genome, - cluster.is_reverse, - index.genome.n_genome, - last_exon.read_end - last_exon.read_start, - eff_length, - ); + let (jr_shift, mut motif, motif_score, mut jj_l, mut jj_r) = scorer + .find_best_junction_position( + read_seq, + last_exon.read_end, + donor_sa, + read_gap.max(0), + genome_gap, + &index.genome, + cluster.is_reverse, + index.genome.n_genome, + last_exon.read_end - last_exon.read_start, + eff_length, + ); // Clamp shift: jr_shift = STAR's jR. Lower bound: can't consume entire exon A. // Upper bound: scan already limited to < shared+eff_length but clamp for safety. let prev_match_len = (last_exon.read_end - last_exon.read_start) as i32; - let jr_shift = jr_shift + let mut jr_shift = jr_shift .max(-prev_match_len) .min((eff_length + shared) as i32); @@ -1560,21 +1561,61 @@ fn stitch_align_to_transcript( // --- Type-specific scoring and tracking --- if is_splice { - // Check stitch mismatch limit + // Look the junction up in the annotation *before* gating on the + // motif, because an annotated junction supplies its own motif and + // that is what STAR's mismatch gate sees + // (`stitchAlignToTranscript.cpp`: `sjdbInd` is resolved, `jcan` is + // overwritten from `sjdbMotif`, and only then is the gate applied). + let junc_donor_sa = (donor_sa as i64 + jr_shift as i64) as u64; + let donor_fwd = + index.sa_pos_to_forward(junc_donor_sa, cluster.is_reverse, del as usize); + let acceptor_fwd = donor_fwd + del as u64 - 1; + + // The metadata-bearing lookup. Junctions inserted by two-pass mode + // carry no motif or shift, so they are only visible through the + // annotated map below; the two lookups are therefore separate. + let sj_entry = junction_db + .and_then(|db| db.find(donor_fwd, acceptor_fwd).and_then(|i| db.entry(i))); + let is_annotated = sj_entry.is_some() + || junction_db.is_some_and(|db| { + db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 0) + || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 1) + || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 2) + }); + + if let Some(pj) = sj_entry { + // The annotated entry's motif wins over the one scanned from + // the genome (STAR: `jcan = sjdbMotif[ind]`). This is what the + // mismatch gate, `outFilterIntronMotifs` and SJ.out.tab see. + motif = crate::junction::decode_motif(pj.motif); + + // Non-canonical annotated junction: snap the boundary onto the + // annotated coordinates instead of leaving it wherever the + // leftmost-flush scan put it (STAR shifts `jR` by `shiftLeft`). + // `jr_shift` is consumed below to move exon A's right edge and + // exon B's left edge together, so shifting it here relocates + // both coherently and leaves the intron length unchanged. + if pj.motif == 0 { + let sl = pj.shift_left as i32; + if (eff_length as i32) <= sl || prev_match_len <= sl { + return None; // STAR -1000006 + } + jr_shift += sl; + if (last_exon.read_end as i64 + jr_shift as i64) + >= (eff_read_pos + eff_length) as i64 + { + return None; // STAR -1000006 + } + jj_l = pj.shift_left as u32; + jj_r = pj.shift_right as u32; + } + } + + // Check stitch mismatch limit (now against the annotated motif). if !scorer.stitch_mismatch_allowed(&motif, gap_mm) { return None; } - let is_annotated = junction_db.is_some_and(|db| { - let junc_donor_sa = (donor_sa as i64 + jr_shift as i64) as u64; - let donor_fwd = - index.sa_pos_to_forward(junc_donor_sa, cluster.is_reverse, del as usize); - let acceptor_fwd = donor_fwd + del as u64 - 1; - db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 0) - || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 1) - || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 2) - }); - if is_annotated { d_score += scorer.sjdb_score; } else { @@ -3847,11 +3888,16 @@ mod tests { } } + /// A one-entry table whose stored coordinates are deliberately nowhere near + /// the seeds being stitched. The `sjA` path addresses entries by *index*, + /// trusting the tag the Gsj hit carried, so it still fires; the general + /// path looks entries up by *coordinate* and finds nothing. That is what + /// makes the two paths distinguishable here. fn sjab_db(motif: u8, shift_left: u8, shift_right: u8) -> crate::junction::SpliceJunctionDb { crate::junction::SpliceJunctionDb::from_prepared(vec![crate::junction::PreparedJunction { chr_idx: 0, - start_pos: 20, - end_pos: 199, + start_pos: 50_000, + end_pos: 50_180, motif, shift_left, shift_right, @@ -3944,6 +3990,81 @@ mod tests { ); } + #[test] + fn noncanonical_annotated_junction_snaps_to_annotated_coords() { + use crate::align::score::AlignmentScorer; + + let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]); + let scorer = AlignmentScorer::from_params_minimal(); + let cluster = sjab_cluster(); + let read = [0u8; 50]; + // No `sj_a` tag, so this exercises the general coordinate-lookup path. + let wt = sjab_wt(20, -1); + let wa = sjab_wa(20, 30, 200, -1); + + let run = |db: Option<&crate::junction::SpliceJunctionDb>| { + stitch_align_to_transcript(&wt, &wa, &read, &index, &scorer, &cluster, db, 0, "t") + }; + + // Calibrate: where does the unannotated boundary search land? + let base = run(None).expect("baseline stitch should succeed"); + let donor_end = base.exons[0].genome_end; + let acceptor_start = base.exons[1].genome_start; + let intron_len = acceptor_start - donor_end; + + // Annotate exactly that junction as non-canonical with a 2-base + // micro-repeat to its left. + let shift_left = 2u8; + let annotated = crate::junction::SpliceJunctionDb::from_prepared(vec![ + crate::junction::PreparedJunction { + chr_idx: 0, + start_pos: donor_end, + end_pos: acceptor_start - 1, + motif: 0, + shift_left, + shift_right: 0, + strand: 0, + }, + ]); + let snapped = run(Some(&annotated)).expect("annotated stitch should succeed"); + + // Both boundaries move right by shift_left; the intron keeps its length. + assert_eq!(snapped.exons[0].genome_end, donor_end + shift_left as u64); + assert_eq!( + snapped.exons[1].genome_start, + acceptor_start + shift_left as u64 + ); + assert_eq!( + snapped.exons[1].genome_start - snapped.exons[0].genome_end, + intron_len, + "snapping must relocate the junction, not resize the intron" + ); + + // It is recorded as annotated, with the stored shifts. + assert_eq!(snapped.junction_annotated, vec![true]); + assert_eq!(snapped.junction_shifts, vec![(shift_left as u32, 0)]); + + // A canonical annotated junction at the same coordinates, carrying the + // same micro-repeat, is *not* snapped: STAR only shifts non-canonical + // entries. Canonical entries store `start_pos + shift_left`, so the + // pre-shift positions are backed off to land on the same stored key. + let canonical = crate::junction::SpliceJunctionDb::from_prepared(vec![ + crate::junction::PreparedJunction { + chr_idx: 0, + start_pos: donor_end - shift_left as u64, + end_pos: acceptor_start - 1 - shift_left as u64, + motif: 1, + shift_left, + shift_right: 0, + strand: 1, + }, + ]); + let unsnapped = run(Some(&canonical)).expect("canonical annotated stitch should succeed"); + assert_eq!(unsnapped.exons[0].genome_end, donor_end); + assert_eq!(unsnapped.exons[1].genome_start, acceptor_start); + assert_eq!(unsnapped.junction_annotated, vec![true]); + } + #[test] fn too_large_repeat_around_annotated_junction_rejects() { use crate::align::score::AlignmentScorer; From 1d9a3face3b7a5b4e66f50f00dddea7c43013d07 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 22:07:21 +0200 Subject: [PATCH 04/17] feat(params): six more STAR align knobs, alignInsertionFlush Right, read-length mismatch cap Adds the flags this theme needs that #145 did not bring: `alignEndsProtrude`, `alignSoftClipAtReferenceEnds`, `alignInsertionFlush`, `alignTranscriptsPerReadNmax`, `outFilterMismatchNoverReadLmax`, `seedNoneLociPerWindow`, `seedSplitMin`. Off-menu values are rejected loudly rather than silently ignored. `--alignEndsType` and its `ext[mate][end]` matrix come from #145; this branch builds on that rather than duplicating it. The scorer gains only the three fields #145 does not have: `flush_right`, `soft_clip_at_reference_ends` and `p_mm_max_read`. Two behaviours land here: - `--alignInsertionFlush Right` now works. The insertion-placement scan accepts ties as well as strict improvements, then walks the insertion further right for as long as the read keeps matching, rejecting when it runs out of read on the B side (STAR -1000009). `None`, the default, is bit-identical to before. - `--outFilterMismatchNoverReadLmax` is enforced. STAR's `outFilterMismatchNmaxTotal` is the tightest of three caps: absolute, a fraction of the mapped length, and a fraction of the read length. Only the first two were applied, so the flag had nothing to bind on. Now in `AlignmentScorer::mismatch_nmax_total`. Note for a follow-up: the pre-existing `--clipAdapterType` validation sits inside the `solo_enabled()` branch of `try_parse_from`, so it only fires for STARsolo runs. Co-Authored-By: Claude Opus 5 (1M context) --- src/align/score.rs | 74 ++++++++++++++++++++++++++ src/align/stitch.rs | 54 +++++++++++++++---- src/params/mod.rs | 124 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 241 insertions(+), 11 deletions(-) diff --git a/src/align/score.rs b/src/align/score.rs index 7a34f55..b906b4d 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -52,6 +52,9 @@ pub struct AlignmentScorer { /// Read-end extension policy (alignEndsType). `ext[iMate][iEnd]==true` forces /// full end-to-end extension (no terminal soft-clip) of that mate/end. pub align_ends_type: crate::params::AlignEndsType, + pub flush_right: bool, + pub soft_clip_at_reference_ends: bool, + pub p_mm_max_read: f64, } impl AlignmentScorer { @@ -80,6 +83,9 @@ impl AlignmentScorer { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, } } @@ -120,6 +126,9 @@ impl AlignmentScorer { out_filter_score_min_over_lread: params.out_filter_score_min_over_lread, // Parsed+validated in Parameters::validate; default to Local if unset. align_ends_type: params.align_ends_type.parse().unwrap_or_default(), + flush_right: params.align_insertion_flush == "Right", + soft_clip_at_reference_ends: params.align_soft_clip_at_reference_ends != "No", + p_mm_max_read: params.out_filter_mismatch_nover_read_lmax, } } @@ -132,6 +141,17 @@ impl AlignmentScorer { ((genomic_span as f64).log2() * self.score_genomic_length_log2_scale - 0.5).ceil() as i32 } + /// STAR's `outFilterMismatchNmaxTotal`: the tightest of the absolute cap, + /// a fraction of the mapped length, and a fraction of the read length. + /// + /// `mapped_len` is the alignment's span in the read; `read_len` is the + /// whole read (both mates summed, for a pair). + pub fn mismatch_nmax_total(&self, mapped_len: usize, read_len: usize) -> u32 { + let by_mapped = (self.p_mm_max * mapped_len as f64) as u32; + let by_read = (self.p_mm_max_read * read_len as f64) as u32; + self.n_mm_max.min(by_mapped).min(by_read) + } + /// Annotated junctions score `sjdb_score`; unannotated junctions score `motif_score`. pub fn score_annotated_junction(&self, motif_score: i32, annotated: bool) -> i32 { if annotated { @@ -693,6 +713,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; // Intron from position 2, length 12 (spans positions 2-13 inclusive) @@ -738,6 +761,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; let motif = scorer.detect_splice_motif(2, 12, &genome); @@ -782,6 +808,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; let motif = scorer.detect_splice_motif(2, 12, &genome); @@ -824,6 +853,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; let motif = scorer.detect_splice_motif(2, 12, &genome); @@ -859,6 +891,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; let (score, gap_type) = scorer.score_gap(0, 5, 0, &genome); @@ -892,6 +927,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; // Small gap (< align_intron_min) is deletion @@ -933,6 +971,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; // Gap starting at position 2 (GT), length 26 (>= 21) is splice junction @@ -972,6 +1013,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; let annotated_score = scorer.score_annotated_junction(0, true); @@ -1012,6 +1056,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; // CT-AC motif: (1,3,0,1) — reverse complement of GT-AG @@ -1120,6 +1167,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; // Gap of exactly 589824 starting at position 100 should be splice junction @@ -1159,6 +1209,24 @@ mod tests { assert_eq!(SpliceMotif::GtAt.filter_category(), 3); } + #[test] + fn mismatch_nmax_total_takes_the_tightest_of_three_caps() { + let mut s = AlignmentScorer::from_params_minimal(); + s.n_mm_max = 10; + s.p_mm_max = 0.3; // fraction of MAPPED length + s.p_mm_max_read = 1.0; // fraction of READ length + + // Mapped fraction binds: 0.3 * 20 = 6. + assert_eq!(s.mismatch_nmax_total(20, 100), 6); + // Absolute cap binds: 0.3 * 100 = 30, capped at 10. + assert_eq!(s.mismatch_nmax_total(100, 100), 10); + + // The read-length cap is the term that was previously missing: with a + // 100-base alignment it now binds before either of the other two. + s.p_mm_max_read = 0.02; + assert_eq!(s.mismatch_nmax_total(100, 100), 2); + } + #[test] fn test_filter_category_from_encoded() { assert_eq!(SpliceMotif::filter_category_from_encoded(0), 0); // non-canonical @@ -1200,6 +1268,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; // Gap of 1001 (> 1000 max) should be deletion, not splice junction @@ -1252,6 +1323,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, } } diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 68d8ffb..c1ee89d 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -1663,14 +1663,15 @@ fn stitch_align_to_transcript( let mut jr = 0i32; // number of shared bases going to A side + let genome_offset: u64 = if cluster.is_reverse { + index.genome.n_genome + } else { + 0 + }; + if shared > 0 { // jR scanning to find optimal insertion placement (STAR lines 265-291) let shared_usize = shared as usize; - let genome_offset: u64 = if cluster.is_reverse { - index.genome.n_genome - } else { - 0 - }; let mut score1 = 0i32; let mut max_score1 = 0i32; @@ -1699,10 +1700,11 @@ fn stitch_align_to_transcript( } } - // STAR default (alignInsertionFlush=None): strict > only. - // First maximum wins = leftmost insertion in current coordinate space. - // For flushRight mode (not yet implemented): Score1 >= maxScore1. - if score1 > max_score1 { + // `alignInsertionFlush None` (STAR's default): strict `>`, so + // the first maximum wins and the insertion sits at its leftmost + // position. `Right` accepts ties too, walking the insertion to + // the rightmost equally-scoring position. + if score1 > max_score1 || (score1 == max_score1 && scorer.flush_right) { max_score1 = score1; jr = jr1 as i32; } @@ -1739,6 +1741,27 @@ fn stitch_align_to_transcript( } // shared == 0: simple insertion, no scanning needed, jr stays 0 + // `alignInsertionFlush Right`: walk the insertion further right for as + // long as the read keeps matching the genome, so a repeat-ambiguous + // insertion lands at the rightmost position it can. Running out of + // read on the B side means there is nowhere left to put the insertion + // (STAR -1000009). + if scorer.flush_right { + let jr_limit = + (eff_read_pos + eff_length).saturating_sub(last_exon.read_end + ins) as i32; + while jr < jr_limit { + let r_idx = last_exon.read_end + jr as usize; + let g_pos = last_exon.genome_end + jr as u64 + genome_offset; + match index.genome.get_base(g_pos) { + Some(gb) if gb < 4 && read_seq.get(r_idx) == Some(&gb) => jr += 1, + _ => break, + } + } + if jr >= jr_limit { + return None; + } + } + let ins_score = scorer.score_ins_open + scorer.score_ins_base * ins as i32; d_score += ins_score; new_wt.n_gap += 1; @@ -1768,10 +1791,13 @@ fn stitch_align_to_transcript( }); } - // Mismatch limit check + // Mismatch limit check. STAR's `outFilterMismatchNmaxTotal` is the minimum + // of three caps: the absolute one, a fraction of the *mapped* length, and + // a fraction of the *read* length. The third was missing here, which made + // `--outFilterMismatchNoverReadLmax` unenforceable. let total_mm = new_wt.n_mismatch + gap_mm; let total_len = new_wt.read_end.max(eff_read_pos + eff_length) - new_wt.read_start; - let mm_limit = ((scorer.p_mm_max * total_len as f64) as u32).min(scorer.n_mm_max); + let mm_limit = scorer.mismatch_nmax_total(total_len, read_seq.len()); if total_mm > mm_limit { return None; } @@ -3681,6 +3707,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; // Left overhang (prev.length) = 3, below min of 5 @@ -3725,6 +3754,9 @@ mod tests { align_spliced_mate_map_lmin_over_lmate: 0.66, out_filter_score_min_over_lread: 0.66, align_ends_type: crate::params::AlignEndsType::default(), + flush_right: false, + soft_clip_at_reference_ends: true, + p_mm_max_read: 1.0, }; // Both overhangs >= 5 diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..7b38133 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -55,6 +55,8 @@ impl std::str::FromStr for RunMode { } } +// --------------------------------------------------------------------------- + impl std::fmt::Display for RunMode { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { match self { @@ -713,6 +715,12 @@ pub struct Parameters { #[arg(long = "outFilterMismatchNoverLmax", default_value_t = 0.3)] pub out_filter_mismatch_nover_lmax: f64, + /// Max mismatches per pair as a fraction of *read* length (both mates + /// summed), as opposed to `outFilterMismatchNoverLmax`, which is a + /// fraction of the *mapped* length. + #[arg(long = "outFilterMismatchNoverReadLmax", default_value_t = 1.0)] + pub out_filter_mismatch_nover_read_lmax: f64, + /// Min alignment score (absolute) #[arg(long = "outFilterScoreMin", default_value_t = 0)] pub out_filter_score_min: i32, @@ -818,6 +826,22 @@ pub struct Parameters { default_values_t = vec![0, -1, 0, 0], allow_hyphen_values = true)] pub align_sj_stitch_mismatch_nmax: Vec, + /// Allow protrusion of one mate's alignment past the other's start: + /// `ConcordantPair` or ` ConcordantPair`. + #[arg(long = "alignEndsProtrude", num_args = 1..=2, + default_values_t = vec!["0".to_string(), "ConcordantPair".to_string()])] + pub align_ends_protrude: Vec, + + /// Allow the soft-clipping of alignments past the ends of chromosomes + /// (`Yes`, default) or forbid it (`No`). + #[arg(long = "alignSoftClipAtReferenceEnds", default_value = "Yes")] + pub align_soft_clip_at_reference_ends: String, + + /// Which end of a repeat-ambiguous insertion the insertion is flushed to: + /// `None` (default, no flushing) or `Right`. + #[arg(long = "alignInsertionFlush", default_value = "None")] + pub align_insertion_flush: String, + /// Splice junction penalty (canonical) #[arg(long = "scoreGap", default_value_t = 0)] pub score_gap: i32, @@ -887,6 +911,14 @@ pub struct Parameters { #[arg(long = "seedPerWindowNmax", default_value_t = 50)] pub seed_per_window_nmax: usize, + /// Max number of one-seed loci per window. + #[arg(long = "seedNoneLociPerWindow", default_value_t = 10)] + pub seed_none_loci_per_window: usize, + + /// Min length of the seed sequences split by the seed search. + #[arg(long = "seedSplitMin", default_value_t = 12)] + pub seed_split_min: usize, + /// Max distance between seed search start positions (defines Nstart = readLen/seedSearchStartLmax + 1) #[arg(long = "seedSearchStartLmax", default_value_t = 50)] pub seed_search_start_lmax: usize, @@ -915,6 +947,12 @@ pub struct Parameters { #[arg(long = "alignTranscriptsPerWindowNmax", default_value_t = 100)] pub align_transcripts_per_window_nmax: usize, + /// Max number of different alignments per read to consider. STAR stops + /// collecting once the running total plus one window's worth would exceed + /// this. + #[arg(long = "alignTranscriptsPerReadNmax", default_value_t = 10000)] + pub align_transcripts_per_read_nmax: usize, + // ── Splice junction database ──────────────────────────────────────── /// GTF file for splice junction annotations #[arg(long = "sjdbGTFfile")] @@ -1551,6 +1589,45 @@ impl Parameters { )); } + // Validate --alignSoftClipAtReferenceEnds. + if !matches!( + params.align_soft_clip_at_reference_ends.as_str(), + "Yes" | "No" + ) { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "unknown --alignSoftClipAtReferenceEnds '{}'; expected Yes or No", + params.align_soft_clip_at_reference_ends + ), + )); + } + // Validate --alignInsertionFlush. + if !matches!(params.align_insertion_flush.as_str(), "None" | "Right") { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "unknown --alignInsertionFlush '{}'; expected None or Right", + params.align_insertion_flush + ), + )); + } + // Validate --alignEndsProtrude: an optional base count followed by + // the protrusion type. + { + let p = ¶ms.align_ends_protrude; + let kind = p.last().map_or("ConcordantPair", String::as_str); + let n_ok = p.len() < 2 || p[0].parse::().is_ok(); + if !n_ok || !matches!(kind, "ConcordantPair") { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "invalid --alignEndsProtrude '{}'; expected [] ConcordantPair", + p.join(" ") + ), + )); + } + } // ── STARsolo validation ───────────────────────────────────────── if params.run_mode == RunMode::AlignReads && params.solo_enabled() { // CB_UMI_Complex needs one CB position + whitelist per segment. @@ -1881,6 +1958,53 @@ mod tests { Parameters::try_parse_from(&full) } + #[test] + fn new_align_knobs_parse_and_reject() { + let p = try_parse(&["--readFilesIn", "r.fq"]).unwrap(); + assert_eq!(p.align_ends_type, "Local"); + assert_eq!(p.align_soft_clip_at_reference_ends, "Yes"); + assert_eq!(p.align_insertion_flush, "None"); + assert_eq!(p.align_transcripts_per_read_nmax, 10000); + assert!((p.out_filter_mismatch_nover_read_lmax - 1.0).abs() < f64::EPSILON); + assert_eq!(p.seed_none_loci_per_window, 10); + assert_eq!(p.seed_split_min, 12); + + assert!(try_parse(&["--readFilesIn", "r.fq", "--alignEndsType", "EndToEnd"]).is_ok()); + assert!(try_parse(&["--readFilesIn", "r.fq", "--alignInsertionFlush", "Right"]).is_ok()); + assert!( + try_parse(&[ + "--readFilesIn", + "r.fq", + "--alignSoftClipAtReferenceEnds", + "No" + ]) + .is_ok() + ); + assert!( + try_parse(&[ + "--readFilesIn", + "r.fq", + "--alignEndsProtrude", + "5", + "ConcordantPair" + ]) + .is_ok() + ); + + // Off-menu values are rejected loudly rather than silently ignored. + assert!(try_parse(&["--readFilesIn", "r.fq", "--alignEndsType", "Global"]).is_err()); + assert!(try_parse(&["--readFilesIn", "r.fq", "--alignInsertionFlush", "Left"]).is_err()); + assert!( + try_parse(&[ + "--readFilesIn", + "r.fq", + "--alignSoftClipAtReferenceEnds", + "Maybe" + ]) + .is_err() + ); + } + #[test] fn defaults() { let p = try_parse(&["--readFilesIn", "reads.fq"]).unwrap(); From 6bfd996b4b2f246827735af9ed636fd545fff3df Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 22:14:43 +0200 Subject: [PATCH 05/17] feat(align): alignSoftClipAtReferenceEnds An alignment may no longer be soft-clipped past the end of its chromosome when `--alignSoftClipAtReferenceEnds No` is given. `Yes`, the default, leaves behaviour unchanged. Sits directly after the end-to-end boundary check from #145: both are "this extension may not run off the reference" rules and share the extension results. Co-Authored-By: Claude Opus 5 (1M context) --- src/align/stitch.rs | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index c1ee89d..402ec2a 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -2041,6 +2041,28 @@ pub(crate) fn finalize_transcript( return None; } + // `--alignSoftClipAtReferenceEnds No`: an alignment may not be + // soft-clipped past the end of its chromosome, so discard any that would + // be. `Yes`, the default, leaves behaviour unchanged. + if !scorer.soft_clip_at_reference_ends { + let chr = cluster.chr_idx; + let chr_start = index.genome.chr_start[chr]; + let chr_end = chr_start + index.genome.chr_length[chr]; + let span_start = wt + .genome_start + .saturating_sub(left_extend.extend_len as u64); + let span_end = wt.genome_end + right_extend.extend_len as u64; + let clipped_left = alignment_start - left_extend.extend_len; + let clipped_right = read_seq + .len() + .saturating_sub(alignment_end + right_extend.extend_len); + if (clipped_left > 0 && span_start <= chr_start) + || (clipped_right > 0 && span_end >= chr_end) + { + return None; + } + } + // STAR finalization check: exon lengths including repeat lengths (shiftSJ) // For non-annotated junctions: exon_len >= alignSJoverhangMin + shiftSJ[side] // For annotated junctions: exon_len >= alignSJDBoverhangMin From 0a771fbc7a50779b0cd01ce2a95fa0be8f6a494c Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 15:41:25 +0200 Subject: [PATCH 06/17] feat(align): apply the genomic-length penalty inside the recursion, not after it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Closes item #7 of STAR-RS-COMPARISON.md §7.2. STAR applies `scoreGenomicLengthLog2scale` in the base case of `stitchWindowAligns`, before the dedup and eviction that decide which transcripts survive the window. rustar-aligner applied it only in `finalize_transcript`, which runs after those decisions, so the recursion ranked candidates on unpenalised scores: a sprawling alignment could displace a compact one of equal raw score, and the penalty arrived too late to matter. `WorkingTranscript` gains `rank_score` = `score` + genomic-length penalty, computed in the base case and used by the three places that compare siblings: the two `same_structure` dominance tests and the `min_by_key` eviction. `rank_score` is kept separate from `score` rather than folded into it. The PE mate-split path derives two per-mate transcripts from one combined score, and splitting an already-penalised score across mates has no clean semantics. `finalize_transcript` still applies the penalty to the final score exactly as before, so nothing is double-counted. This changes which alignment wins among equal-raw-score candidates, which is the point. §7.5's lesson from the dropped item #2 applies: if the differential measures negative, this commit is the one to revert. Co-Authored-By: Claude Opus 5 (1M context) --- src/align/stitch.rs | 37 +++++++++++++++++++++++++++++++++---- 1 file changed, 33 insertions(+), 4 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 402ec2a..31ac476 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -1120,6 +1120,21 @@ pub(crate) struct WorkingTranscript { /// STAR's shiftSJ[isj][0] and shiftSJ[isj][1]. pub(crate) junction_shifts: Vec<(u32, u32)>, pub(crate) n_anchor: u32, + /// Score used to *rank* this transcript against its siblings inside the + /// recursion: `score` plus the genomic-length penalty. + /// + /// STAR applies `scoreGenomicLengthLog2scale` in the base case of + /// `stitchWindowAligns`, before the dedup and eviction that pick which + /// transcripts survive, so a compact alignment can beat a sprawling one of + /// equal raw score. rustar-aligner applied it only at finalization, i.e. + /// after those decisions had already been made. + /// + /// Kept separate from `score` rather than folded into it: the PE mate-split + /// path builds two per-mate transcripts from one combined `score`, and + /// splitting an already-penalised score across mates has no clean meaning. + /// `finalize_transcript` still applies the penalty to the final score, so + /// nothing is counted twice. + pub(crate) rank_score: i32, // Tight bounds for extension at finalization pub(crate) read_start: usize, pub(crate) read_end: usize, @@ -1139,6 +1154,7 @@ impl WorkingTranscript { junction_annotated: Vec::new(), junction_shifts: Vec::new(), n_anchor: 0, + rank_score: 0, read_start: 0, read_end: 0, genome_start: 0, @@ -2520,6 +2536,13 @@ fn stitch_recurse( } } + // STAR applies the genomic-length penalty here, in the base case, + // *before* dedup and eviction. Ranking on the raw score instead + // lets a sprawling alignment displace a compact one of equal + // score, which STAR would never do. + wt.rank_score = wt.score + + scorer.genomic_length_penalty(wt.genome_end.saturating_sub(wt.genome_start)); + // Dedup via blocks_overlap: drop if subset of existing higher-score transcript. // Use same_structure guard: only dedup transcripts with same number of exon // blocks. A non-spliced path should never be killed by a spliced one here @@ -2543,11 +2566,11 @@ fn stitch_recurse( // Only dedup transcripts with same number of exon blocks (junctions). let same_structure = wt.exons.len() == existing.exons.len(); - if same_structure && overlap >= wt_len && existing.score >= wt.score { + if same_structure && overlap >= wt_len && existing.rank_score >= wt.rank_score { dominated = true; break; } - if same_structure && overlap >= ex_len && wt.score >= existing.score { + if same_structure && overlap >= ex_len && wt.rank_score >= existing.rank_score { remove_indices.push(idx); } } @@ -2562,8 +2585,8 @@ fn stitch_recurse( } else if let Some(worst_idx) = transcripts .iter() .enumerate() - .min_by_key(|(_, t)| t.score) - .filter(|(_, t)| t.score < wt.score) + .min_by_key(|(_, t)| t.rank_score) + .filter(|(_, t)| t.rank_score < wt.rank_score) .map(|(i, _)| i) { // STAR-faithful eviction: keep the N best transcripts. @@ -2812,6 +2835,9 @@ pub(crate) fn split_combined_wt( junction_annotated: m1_ja, junction_shifts: m1_js, n_anchor: 0, + // Per-mate transcripts are never ranked against each other inside + // the recursion; `finalize_transcript` applies the penalty. + rank_score: m1_score, read_start: m1_read_start, read_end: m1_read_end, genome_start: m1_genome_start, @@ -2827,6 +2853,9 @@ pub(crate) fn split_combined_wt( junction_annotated: m2_ja, junction_shifts: m2_js, n_anchor: 0, + // Per-mate transcripts are never ranked against each other inside + // the recursion; `finalize_transcript` applies the penalty. + rank_score: m2_score, read_start: m2_read_start, read_end: m2_read_end, genome_start: m2_genome_start, From 2fbe698c32a7ad4ed22c9c5230a01b5e60e9b21b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 15:49:15 +0200 Subject: [PATCH 07/17] feat(test): annotated-junction differential tier, plus alignTranscriptsPerReadNmax MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds `test/simulate_junction_reads.py` and a `yeast_junction_gtf` tier. The existing tiers align real RNA-seq reads, of which only a small fraction cross an intron — yeast has a few hundred — so the annotated-junction code paths are barely exercised and a divergence in them is easy to miss. The simulator draws reads from spliced transcript sequences centred on annotated junctions, half of them reverse-complemented, from a fixed seed so the FASTQ is reproducible. Every read crosses an annotated junction. Also wires `--alignTranscriptsPerReadNmax`: STAR breaks out of the window loop once one more window's worth of transcripts could overrun the per-read cap, so the cap is never exceeded rather than merely noticed afterwards. Measured with this tier, 1000 reads, yeast R64-1-1 with the Ensembl 110 GTF, against native STAR 2.7.11b, records compared on FLAG/RNAME/POS/CIGAR: origin/main 835 / 1000 identical this branch 933 / 1000 identical gained 98, regressed 0 junctions (total / annotated) STAR 155 / 155 origin/main 149 / 131 this branch 149 / 131 The gain is entirely in CIGARs: the junction set is unchanged, but 98 reads now place their exon boundaries where STAR places them. Against an unannotated index both binaries score 938/1000, gained 0, lost 0, confirming the in-recursion length penalty is neutral there and the whole improvement comes from annotated-boundary snapping. Output is byte-identical at `--runThreadN` 1, 4 and 16. Adds `docs-old/dev/divergences.md`, recording the two places this branch knowingly differs from STAR: the extend-to-end rejection signal, and the missing per-junction strand on non-canonical annotated junctions. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 25 +++++ docs-old/dev/divergences.md | 78 ++++++++++++++++ src/align/read_align.rs | 10 ++ test/run_tests.sh | 10 ++ test/simulate_junction_reads.py | 156 ++++++++++++++++++++++++++++++++ 5 files changed, 279 insertions(+) create mode 100644 docs-old/dev/divergences.md create mode 100644 test/simulate_junction_reads.py diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..d189e09 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,31 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- **Annotated-junction stitching, `--alignEndsType`, and the in-recursion + genomic-length penalty** — closes items #4b, #7, #8 and part of #9 of + `STAR-RS-COMPARISON.md` §7.2. + + - `SpliceJunctionDb` gains STAR's `sjdb*` table and a `binarySearch2`-style + `find()`, so an annotated junction can be addressed by index (`sjA`) and + not merely tested for existence. No index format change: `sjdbInfo.txt` + already carried motif, shiftLeft, shiftRight and strand. + - Seeds derived from the sjdb genome insert carry their `sjA` tag through the + window into the stitcher, which takes STAR's annotated-junction shortcut + when two pieces turn out to be the two flanks of one junction. + - Non-canonical annotated junctions are snapped onto their annotated + coordinates instead of being left wherever the leftmost-flush scan put + them, and the stored motif overrides the one re-derived from the genome. + - `--alignEndsType`, `--alignSoftClipAtReferenceEnds`, + `--alignInsertionFlush Right` and `--outFilterMismatchNoverReadLmax` are + implemented; `--alignEndsProtrude`, `--alignTranscriptsPerReadNmax`, + `--seedNoneLociPerWindow` and `--seedSplitMin` are accepted. + - The genomic-length penalty is applied inside the stitch recursion, where it + can affect which transcripts survive, rather than only at finalisation. + + On 1000 junction-spanning yeast reads against native STAR 2.7.11b, records + identical in FLAG/RNAME/POS/CIGAR rise from 835 to 933 out of 1000, with no + read regressing. Output against an unannotated index is unchanged. + - **STARsolo single-cell quantification (`--soloType`)** — the 10x Chromium / plate-based count-matrix pipeline, ported from STAR and verified against real STARsolo (#90). diff --git a/docs-old/dev/divergences.md b/docs-old/dev/divergences.md new file mode 100644 index 0000000..a951cdb --- /dev/null +++ b/docs-old/dev/divergences.md @@ -0,0 +1,78 @@ +# Deliberate divergences from STAR 2.7.11b + +rustar-aligner aims to be a faithful port: given the same inputs and flags, it should +produce the same bytes as STAR 2.7.11b. This file records the places where it +deliberately does not, and why. + +Every entry here is a case where STAR's behaviour is undefined, plainly wrong, or +impossible to reproduce, and where reproducing it would mean shipping a known +defect. Each is covered by a test that asserts the **correct** result, never +STAR's. + +Divergences that arise from bugs on our side are not listed here. They are bugs, +and they get fixed. + +--- + +## Format + +Each entry gives: what STAR does, what rustar-aligner does instead, why, and the +test that locks the behaviour in. + +--- + +## D-01 · `--alignEndsType` failure is signalled by rejection, not by a sentinel score + +**STAR.** When an end cannot be extended to the end of the read (the chromosome +runs out first), `extendAlign` returns a score of `-999999999` and lets that +value propagate; downstream score filters then discard the alignment. + +**rustar-aligner.** `ExtendResult::hit_boundary` says so explicitly, and +`finalize_transcript` drops the transcript. + +**Why.** `finalize_transcript` clamps the final score with `.max(0)`. A +`-999999999` sentinel would be clamped to `0` and the alignment would *survive* +as a zero-score record rather than being rejected. The sentinel only works in +STAR because STAR does not clamp. The observable outcome is the same, without +depending on arithmetic this code path does not perform. + +**Test.** `extend_to_end_reports_boundary_when_the_chromosome_runs_out` +(`src/align/stitch.rs`). + +--- + +## D-02 · Non-canonical annotated junctions carry no strand + +**STAR.** Tracks a per-junction strand (`sjStr`) alongside the motif, and for an +annotated junction takes that strand from `sjdbStrand`. SJ.out.tab column 4 then +reports the annotated strand even when the motif is non-canonical. + +**rustar-aligner.** Has no per-junction strand on the working transcript, so +column 4 reports `0` (undefined) for a non-canonical annotated junction where +STAR would report `1` or `2`. Canonical junctions are unaffected: their strand +is implied by the motif, which both derive identically. + +**Why.** Not a deliberate improvement, just not yet ported. Recorded here so the +difference is not mistaken for a bug in the annotated-junction path, which is +otherwise faithful. Adding `junction_strands` to `WorkingTranscript` and +`Transcript` would close it. + +**Test.** None yet; this entry is the marker. + +--- + +## Divergences considered and rejected + +Cases where STAR-rs, the other Rust port, made a different choice that +rustar-aligner deliberately does not adopt. + +### Transcriptome-BAM primary flag + +STAR picks the primary alignment for the transcriptome BAM with a per-thread +RNG, so the output depends on `--runThreadN`. STAR-rs resolves this by always +taking the first alignment (`j == 0`). + +rustar-aligner already picks by a per-read seed, which is thread-count invariant +and therefore fixes the same defect. The two ports resolve it differently and +both diverge from STAR; there is no fidelity argument for switching to `j == 0`, +so rustar-aligner keeps its own rule. diff --git a/src/align/read_align.rs b/src/align/read_align.rs index b4d2f20..f959cd0 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -299,6 +299,16 @@ pub fn align_read( }; for (ci, cluster) in clusters.iter().enumerate() { + // STAR's `alignTranscriptsPerReadNmax` headroom break + // (`ReadAlign_stitchPieces.cpp`): stop collecting windows once one more + // window's worth of transcripts could overrun the per-read cap. The + // test is on the headroom, not on the running total, so the cap is + // never exceeded rather than merely noticed afterwards. + if transcripts.len() + params.align_transcripts_per_window_nmax + >= params.align_transcripts_per_read_nmax + { + break; + } let debug_name = if debug_read { read_name } else { "" }; let cluster_transcripts = stitch_seeds_with_jdb_debug( cluster, diff --git a/test/run_tests.sh b/test/run_tests.sh index fbfa7ba..5b75599 100755 --- a/test/run_tests.sh +++ b/test/run_tests.sh @@ -44,6 +44,15 @@ TEST_CASES=( "yeast_1k_paired:yeast:ERR12389696_sub_1_1000.fastq.gz,ERR12389696_sub_2_1000.fastq.gz:paired:--outSAMtype SAM" "yeast_1k_twopass:yeast:ERR12389696_sub_1_1000.fastq.gz:single:--twopassMode Basic --outSAMtype SAM" "yeast_10k:yeast:ERR12389696_sub_1_10k.fastq.gz:single:--outSAMtype SAM" + # Annotated-junction tier. The tiers above use real RNA-seq reads, of which + # only a small fraction cross an intron -- yeast has a few hundred -- so the + # annotated-junction code paths (the sjA shortcut, annotated-boundary + # snapping) are barely exercised and a divergence in them is easy to miss. + # These reads are drawn from spliced transcripts and every one of them + # crosses an annotated junction. Generate with: + # test/simulate_junction_reads.py \ + # test/data/small/yeast/reads/sim_junction_1k.fastq 1000 100 + "yeast_junction_gtf:yeast_gtf:sim_junction_1k.fastq:single:--outSAMtype SAM" ) # ============================================================================== @@ -350,6 +359,7 @@ run_test_case() { # Parse test specification IFS=':' read -r name dataset reads mode extra_args <<< "$test_spec" + local test_dir="$RESULTS_DIR/${TIMESTAMP}_${name}" local star_dir="$test_dir/star" local rustar_aligner_dir="$test_dir/rustar-aligner-aligner" diff --git a/test/simulate_junction_reads.py b/test/simulate_junction_reads.py new file mode 100644 index 0000000..8e0b21b --- /dev/null +++ b/test/simulate_junction_reads.py @@ -0,0 +1,156 @@ +#!/usr/bin/env python3 +"""Simulate junction-spanning reads from a genome + GTF. + +The existing test tiers all use an unannotated index, which makes the +annotated-junction code paths (`sjA` shortcut, annotated-boundary snapping) +invisible: they simply never run. This generates reads that are guaranteed to +cross annotated exon-exon junctions, so those paths are exercised and any +divergence from STAR shows up as a junction count or CIGAR difference. + +Reads are drawn from spliced transcript sequences, centred on junctions, with a +fixed seed so the FASTQ is reproducible byte-for-byte. + +Usage: + simulate_junction_reads.py [n_reads] [read_len] +""" + +import gzip +import random +import sys +from collections import defaultdict + + +def read_fasta(path): + """Return {contig_name: sequence}. Names are truncated at first whitespace.""" + opener = gzip.open if path.endswith(".gz") else open + seqs, name, chunks = {}, None, [] + with opener(path, "rt") as fh: + for line in fh: + if line.startswith(">"): + if name is not None: + seqs[name] = "".join(chunks) + name = line[1:].split()[0] + chunks = [] + else: + chunks.append(line.strip()) + if name is not None: + seqs[name] = "".join(chunks) + return seqs + + +def read_exons(path): + """Return {transcript_id: (contig, strand, [(start, end), ...])}, 1-based inclusive.""" + opener = gzip.open if path.endswith(".gz") else open + exons = defaultdict(list) + meta = {} + with opener(path, "rt") as fh: + for line in fh: + if line.startswith("#"): + continue + f = line.rstrip("\n").split("\t") + if len(f) < 9 or f[2] != "exon": + continue + attrs = f[8] + key = 'transcript_id "' + i = attrs.find(key) + if i < 0: + continue + j = attrs.find('"', i + len(key)) + tid = attrs[i + len(key) : j] + exons[tid].append((int(f[3]), int(f[4]))) + meta[tid] = (f[0], f[6]) + out = {} + for tid, blocks in exons.items(): + blocks.sort() + contig, strand = meta[tid] + out[tid] = (contig, strand, blocks) + return out + + +COMPLEMENT = str.maketrans("ACGTNacgtn", "TGCANtgcan") + + +def revcomp(s): + return s.translate(COMPLEMENT)[::-1] + + +def main(): + if len(sys.argv) < 4: + sys.exit(__doc__) + genome_path, gtf_path, out_path = sys.argv[1:4] + n_reads = int(sys.argv[4]) if len(sys.argv) > 4 else 1000 + read_len = int(sys.argv[5]) if len(sys.argv) > 5 else 100 + + genome = read_fasta(genome_path) + transcripts = read_exons(gtf_path) + + # Collect every (transcript, junction) pair with enough flank on both sides + # to host a read that actually crosses the junction. + candidates = [] + for tid, (contig, strand, blocks) in transcripts.items(): + if contig not in genome or len(blocks) < 2: + continue + # Offset of each junction within the spliced transcript. + offset = 0 + for k in range(len(blocks) - 1): + offset += blocks[k][1] - blocks[k][0] + 1 + candidates.append((tid, offset)) + + if not candidates: + sys.exit("no multi-exon transcripts found; nothing to simulate") + + rng = random.Random(20260728) + rng.shuffle(candidates) + + spliced_cache = {} + + def spliced(tid): + if tid not in spliced_cache: + contig, strand, blocks = transcripts[tid] + seq = "".join(genome[contig][s - 1 : e] for s, e in blocks) + spliced_cache[tid] = (seq, strand) + return spliced_cache[tid] + + # Yeast has only a few hundred introns, so one read per junction is not + # enough to be interesting. Cycle the junction list, varying the overhang + # each pass, until the requested count is reached or a pass places nothing. + written = 0 + with open(out_path, "w") as out: + while written < n_reads: + placed_this_pass = 0 + for tid, junction_offset in candidates: + if written >= n_reads: + break + seq, _strand = spliced(tid) + if len(seq) < read_len: + continue + # Put the junction in the middle third of the read so both + # overhangs comfortably clear alignSJDBoverhangMin. + overhang = rng.randint(read_len // 3, 2 * read_len // 3) + start = junction_offset - overhang + if start < 0 or start + read_len > len(seq): + continue + read = seq[start : start + read_len].upper() + if "N" in read: + continue + # Half the reads on the opposite strand, so both orientations + # of the annotated-junction path get exercised. + if written % 2 == 1: + read = revcomp(read) + out.write(f"@sim_{written}_{tid}_j{junction_offset}\n") + out.write(read + "\n+\n" + "I" * read_len + "\n") + written += 1 + placed_this_pass += 1 + if placed_this_pass == 0: + break + + print(f"wrote {written} junction-spanning reads to {out_path}", file=sys.stderr) + if written < n_reads: + print( + f"note: only {written} of {n_reads} requested reads were placeable", + file=sys.stderr, + ) + + +if __name__ == "__main__": + main() From 8a0a47a67fd2f36e10cbc4e888f58e6d568244ed Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 22:21:18 +0200 Subject: [PATCH 08/17] docs: drop the extend-to-end divergence note, superseded by #145 #145 handles the same case with an explicit EXTEND_TO_END_KILL sentinel checked before the score clamp, so there is no divergence left to record. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 8 ++++---- docs-old/dev/divergences.md | 22 +--------------------- 2 files changed, 5 insertions(+), 25 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d189e09..2196c27 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -27,10 +27,10 @@ Sections commonly used: Features, Bug fixes, Other changes. - Non-canonical annotated junctions are snapped onto their annotated coordinates instead of being left wherever the leftmost-flush scan put them, and the stored motif overrides the one re-derived from the genome. - - `--alignEndsType`, `--alignSoftClipAtReferenceEnds`, - `--alignInsertionFlush Right` and `--outFilterMismatchNoverReadLmax` are - implemented; `--alignEndsProtrude`, `--alignTranscriptsPerReadNmax`, - `--seedNoneLociPerWindow` and `--seedSplitMin` are accepted. + - `--alignSoftClipAtReferenceEnds`, `--alignInsertionFlush Right` and + `--outFilterMismatchNoverReadLmax` are implemented; `--alignEndsProtrude`, + `--alignTranscriptsPerReadNmax`, `--seedNoneLociPerWindow` and + `--seedSplitMin` are accepted. (`--alignEndsType` itself came from #145.) - The genomic-length penalty is applied inside the stitch recursion, where it can affect which transcripts survive, rather than only at finalisation. diff --git a/docs-old/dev/divergences.md b/docs-old/dev/divergences.md index a951cdb..6127bd6 100644 --- a/docs-old/dev/divergences.md +++ b/docs-old/dev/divergences.md @@ -21,27 +21,7 @@ test that locks the behaviour in. --- -## D-01 · `--alignEndsType` failure is signalled by rejection, not by a sentinel score - -**STAR.** When an end cannot be extended to the end of the read (the chromosome -runs out first), `extendAlign` returns a score of `-999999999` and lets that -value propagate; downstream score filters then discard the alignment. - -**rustar-aligner.** `ExtendResult::hit_boundary` says so explicitly, and -`finalize_transcript` drops the transcript. - -**Why.** `finalize_transcript` clamps the final score with `.max(0)`. A -`-999999999` sentinel would be clamped to `0` and the alignment would *survive* -as a zero-score record rather than being rejected. The sentinel only works in -STAR because STAR does not clamp. The observable outcome is the same, without -depending on arithmetic this code path does not perform. - -**Test.** `extend_to_end_reports_boundary_when_the_chromosome_runs_out` -(`src/align/stitch.rs`). - ---- - -## D-02 · Non-canonical annotated junctions carry no strand +## D-01 · Non-canonical annotated junctions carry no strand **STAR.** Tracks a per-junction strand (`sjStr`) alongside the motif, and for an annotated junction takes that strand from `sjdbStrand`. SJ.out.tab column 4 then From 27e6bbe66bbb47ca8be9812f6078572125aa4ebb Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 12:33:45 +0200 Subject: [PATCH 09/17] fix(junction): a one-base intron is a junction, and an annotated gap is never a deletion Two defects, found by measuring the annotated-junction tier against STAR rather than by reading the code. Together they took that tier from 96.957% to 99.790% exact agreement, and they compound: the first hides the second. **`genomeGenerate` dropped one-base introns.** The GTF junction extractor required `intron_end > intron_start`, which rejects an intron of length one. STAR skips a pair only when the exons touch or overlap (`GTF_transcriptGeneSJ.cpp:123`: `exonLoci[iex][exS] <= exonLoci[iex-1][exE]+1`), so `exon2.start == exon1.end + 2` is a junction it keeps. Yeast R64-1-1 has three. The index carried 361 junctions against STAR's 364, and the three missing were exactly those. `sjdbList.out.tab` is now identical to STAR's. **An annotated gap shorter than `--alignIntronMin` was scored as a deletion.** STAR resolves `sjdbInd` *before* testing `Del >= alignIntronMin` (`stitchAlignToTranscript.cpp:198-241`): on a hit it takes the annotated branch whatever the gap length, and the length test only decides for gaps the annotation does not know. Here the annotation lookup sat inside the already-decided splice branch, so a short annotated gap could never reach it. The lookup moves out; `is_splice` becomes `annotated || long enough`. The CIGAR builder re-derived the same decision from the gap length alone, so even once the score was right the operation was still `D`. It now asks the same question through a shared `gap_is_junction`, which consults the annotation for gaps below the threshold. Reads crossing yeast's one-base introns come out `74M1N76M` as STAR writes them, not `74M1D76M`. Measured on the 1000-read annotated-junction tier, against STAR 2.7.11b: CIGAR disagreements 29 -> 2 AS disagreements 29 -> 2 exact (tie-adjusted) 924/953 (96.957%) -> 951/953 (99.790%) SE and PE on the 10k yeast tiers are unchanged at 99.898% and 99.849%: those reads do not cross a one-base intron, which is why this stayed invisible until the tier existed. --- src/align/stitch.rs | 149 ++++++++++++++++++++++++++++++++++++-------- src/junction/gtf.rs | 59 +++++++++++++++++- 2 files changed, 180 insertions(+), 28 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 31ac476..5ff3596 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -1575,30 +1575,34 @@ fn stitch_align_to_transcript( gap_mm += shared_mm; d_score += shared_score; + // Look the junction up in the annotation *before* deciding what kind of + // gap this is. STAR resolves `sjdbInd` first and, on a hit, takes the + // annotated branch whatever the gap length: `Del >= alignIntronMin` + // only decides between junction and deletion for gaps the annotation + // does not know (`stitchAlignToTranscript.cpp:198-241`). An annotated + // one-base gap is a junction to STAR, not a deletion. + // + // The lookup also supplies the motif, which is what the mismatch gate + // below sees rather than the one scanned from the genome. + let junc_donor_sa = (donor_sa as i64 + jr_shift as i64) as u64; + let donor_fwd = index.sa_pos_to_forward(junc_donor_sa, cluster.is_reverse, del as usize); + let acceptor_fwd = donor_fwd + del as u64 - 1; + + // The metadata-bearing lookup. Junctions inserted by two-pass mode + // carry no motif or shift, so they are only visible through the + // annotated map below; the two lookups are therefore separate. + let sj_entry = + junction_db.and_then(|db| db.find(donor_fwd, acceptor_fwd).and_then(|i| db.entry(i))); + let is_annotated = sj_entry.is_some() + || junction_db.is_some_and(|db| { + db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 0) + || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 1) + || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 2) + }); + let is_splice = is_annotated || is_splice; + // --- Type-specific scoring and tracking --- if is_splice { - // Look the junction up in the annotation *before* gating on the - // motif, because an annotated junction supplies its own motif and - // that is what STAR's mismatch gate sees - // (`stitchAlignToTranscript.cpp`: `sjdbInd` is resolved, `jcan` is - // overwritten from `sjdbMotif`, and only then is the gate applied). - let junc_donor_sa = (donor_sa as i64 + jr_shift as i64) as u64; - let donor_fwd = - index.sa_pos_to_forward(junc_donor_sa, cluster.is_reverse, del as usize); - let acceptor_fwd = donor_fwd + del as u64 - 1; - - // The metadata-bearing lookup. Junctions inserted by two-pass mode - // carry no motif or shift, so they are only visible through the - // annotated map below; the two lookups are therefore separate. - let sj_entry = junction_db - .and_then(|db| db.find(donor_fwd, acceptor_fwd).and_then(|i| db.entry(i))); - let is_annotated = sj_entry.is_some() - || junction_db.is_some_and(|db| { - db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 0) - || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 1) - || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 2) - }); - if let Some(pj) = sj_entry { // The annotated entry's motif wins over the one scanned from // the genome (STAR: `jcan = sjdbMotif[ind]`). This is what the @@ -2079,6 +2083,28 @@ pub(crate) fn finalize_transcript( } } + // Whether a genome gap of `del` bases starting at `donor_end` (SA space) is + // a junction rather than a deletion. `--alignIntronMin` decides for gaps the + // annotation does not know; an annotated junction is a junction whatever its + // length, which is how STAR reads it (`stitchAlignToTranscript.cpp:198-241`) + // and the reason yeast's three one-base introns come out as `1N`, not `1D`. + let gap_is_junction = |donor_end: u64, del: usize| -> bool { + if del >= scorer.align_intron_min as usize && del <= scorer.align_intron_max as usize { + return true; + } + if del == 0 { + return false; + } + let donor_fwd = index.sa_pos_to_forward(donor_end, cluster.is_reverse, del); + let acceptor_fwd = donor_fwd + del as u64 - 1; + index.junction_db.find(donor_fwd, acceptor_fwd).is_some() + || (0..=2).any(|strand| { + index + .junction_db + .is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, strand) + }) + }; + // STAR finalization check: exon lengths including repeat lengths (shiftSJ) // For non-annotated junctions: exon_len >= alignSJoverhangMin + shiftSJ[side] // For annotated junctions: exon_len >= alignSJDBoverhangMin @@ -2094,7 +2120,10 @@ pub(crate) fn finalize_transcript( let genome_gap = next_exon.genome_start as i64 - exon.genome_end as i64; let read_gap = next_exon.read_start as i64 - exon.read_end as i64; let del = genome_gap - read_gap.max(0); - if del >= scorer.align_intron_min as i64 && junction_idx < wt.junction_shifts.len() { + if del > 0 + && gap_is_junction(exon.genome_end, del as usize) + && junction_idx < wt.junction_shifts.len() + { // This is a junction — check exon lengths with repeat let (shift_l, shift_r) = wt.junction_shifts[junction_idx]; let is_annotated = wt.junction_annotated[junction_idx]; @@ -2186,9 +2215,7 @@ pub(crate) fn finalize_transcript( append_match(&mut final_cigar, shared); } let del = (genome_gap - read_gap.max(0)) as usize; - if del >= scorer.align_intron_min as usize - && del <= scorer.align_intron_max as usize - { + if gap_is_junction(prev.genome_end, del) { final_cigar.push(Op::new(Kind::Skip, del)); } else { final_cigar.push(Op::new(Kind::Deletion, del)); @@ -4148,6 +4175,76 @@ mod tests { assert_eq!(unsnapped.junction_annotated, vec![true]); } + /// A gap shorter than `--alignIntronMin` is a deletion unless the + /// annotation knows it, in which case it is a junction however short. + /// + /// STAR resolves `sjdbInd` before testing `Del >= alignIntronMin` + /// (`stitchAlignToTranscript.cpp:198-241`), so the length test only decides + /// for gaps the annotation has never heard of. Yeast R64-1-1 has three + /// one-base introns; before this, reads crossing them scored as deletions + /// and came out `1D` where STAR writes `1N`. + #[test] + fn a_short_annotated_gap_is_a_junction_not_a_deletion() { + use crate::align::score::AlignmentScorer; + + let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]); + let scorer = AlignmentScorer::from_params_minimal(); + assert!( + scorer.align_intron_min > 5, + "the fixture's 5-base gap must be below the intron threshold" + ); + let cluster = sjab_cluster(); + let read = [0u8; 50]; + + // Exon A covers read[0,20) at genome[0,20); the next seed sits five + // bases further along the genome, so the gap is a 5-base deletion. + let wt = sjab_wt(20, -1); + let wa = sjab_wa(20, 20, 25, -1); + + let plain = + stitch_align_to_transcript(&wt, &wa, &read, &index, &scorer, &cluster, None, 0, "t") + .expect("unannotated stitch should succeed"); + assert_eq!(plain.n_gap, 1, "unannotated: a deletion"); + assert_eq!(plain.n_junction, 0); + + // The same gap, annotated. + let donor_end = plain.exons[0].genome_end; + let acceptor_start = plain.exons[1].genome_start; + let db = crate::junction::SpliceJunctionDb::from_prepared(vec![ + crate::junction::PreparedJunction { + chr_idx: 0, + start_pos: donor_end, + end_pos: acceptor_start - 1, + motif: 1, + shift_left: 0, + shift_right: 0, + strand: 1, + }, + ]); + let annotated = stitch_align_to_transcript( + &wt, + &wa, + &read, + &index, + &scorer, + &cluster, + Some(&db), + 0, + "t", + ) + .expect("annotated stitch should succeed"); + assert_eq!(annotated.n_junction, 1, "annotated: a junction"); + assert_eq!(annotated.n_gap, 0); + assert_eq!(annotated.junction_annotated, vec![true]); + assert!( + annotated.score > plain.score, + "the annotated form earns sjdbScore rather than a deletion penalty: \ + {} vs {}", + annotated.score, + plain.score + ); + } + #[test] fn too_large_repeat_around_annotated_junction_rejects() { use crate::align::score::AlignmentScorer; diff --git a/src/junction/gtf.rs b/src/junction/gtf.rs index 411b4f5..4411341 100644 --- a/src/junction/gtf.rs +++ b/src/junction/gtf.rs @@ -204,9 +204,16 @@ pub fn extract_junctions_configured( let intron_start_local_1b = exon1.end + 1; let intron_end_local_1b = exon2.start.saturating_sub(1); - if intron_end_local_1b <= intron_start_local_1b { + // STAR skips a pair only when the exons touch or overlap + // (`GTF_transcriptGeneSJ.cpp:123`, `exonLoci[iex][exS] <= + // exonLoci[iex-1][exE] + 1`). A one-base intron — `exon2.start == + // exon1.end + 2` — is a junction it keeps, and yeast has three of + // them. Requiring `end > start` here dropped exactly those, so the + // aligner never saw them as annotated and emitted `1D` where STAR + // emits `1N`. + if intron_end_local_1b < intron_start_local_1b { log::warn!( - "Invalid junction coordinates: {intron_start_local_1b}-{intron_end_local_1b} (possibly overlapping exons)" + "Invalid junction coordinates: {intron_start_local_1b}-{intron_end_local_1b} (overlapping exons)" ); continue; } @@ -360,6 +367,54 @@ mod tests { assert_eq!(strand, 1); } + /// A one-base intron is a junction, and touching exons are not. + /// + /// STAR skips a pair only when `exon2.start <= exon1.end + 1` + /// (`GTF_transcriptGeneSJ.cpp:123`), so `exon2.start == exon1.end + 2` + /// yields a one-base intron it keeps. Yeast R64-1-1 has three; requiring + /// `end > start` here dropped all three from the index, and reads crossing + /// them then came out with `1D` where STAR writes `1N`. + #[test] + fn a_one_base_intron_is_a_junction_and_touching_exons_are_not() { + let genome = Genome { + transform_blocks: None, + sequence: vec![0; 1000].into(), + n_genome: 1000, + n_genome_real: 1000, + n_chr_real: 1, + chr_start: vec![0, 1000], + chr_length: vec![1000], + chr_name: vec!["chr1".to_string()], + }; + let exon = |start: u64, end: u64, transcript: &str| GtfRecord { + seqname: "chr1".to_string(), + feature: "exon".to_string(), + start, + end, + strand: '+', + attributes: vec![ + ("gene_id".to_string(), "G1".to_string()), + ("transcript_id".to_string(), transcript.to_string()), + ] + .into_iter() + .collect(), + }; + + // T1: exons 100-200 and 202-300 — a one-base intron at 201. + let junctions = + extract_junctions_from_exons(vec![exon(100, 200, "T1"), exon(202, 300, "T1")], &genome) + .unwrap(); + assert_eq!(junctions.len(), 1, "the one-base intron must survive"); + let (_, start, end, _) = junctions[0]; + assert_eq!((start, end), (200, 200), "0-based, inclusive, length one"); + + // T2: exons 100-200 and 201-300 — adjacent, no intron between them. + let touching = + extract_junctions_from_exons(vec![exon(100, 200, "T2"), exon(201, 300, "T2")], &genome) + .unwrap(); + assert!(touching.is_empty(), "touching exons are not a junction"); + } + #[test] fn test_extract_junctions_multiple_transcripts() { let genome = Genome { From bcf2368c04c2ea7a8da64654e196bd6551f0070f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 12:38:19 +0200 Subject: [PATCH 10/17] fix(align): break primary ties on genomic span, as STAR does MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit STAR compares `maxScore`, then `gLength` — the alignment's genomic span — and leaves anything still tied to whichever window it reached first (`ReadAlign_stitchPieces.cpp:340`). The second key here was the junction count, which is a proxy for the span rather than the span itself: usually correlated, not always equal. Using the span directly gains 6 reads on the 10k SE tier (8791 -> 8797 exactly matching STAR) and costs none anywhere. PE and the junction tier are unchanged. What this does not fix is the residual: among the 37 cross-chromosome primary differences on the junction tier, all 37 are multimappers with identical scores, and 18 have identical CIGARs — so identical spans too. STAR's remaining discriminator is window iteration order, which this codebase deliberately does not reproduce, since windows are built per read in parallel to keep output invariant to `--runThreadN` (DIVERGENCE.md §1.1). Those are ties in both tools; only the choice among equals differs. --- src/align/read_align.rs | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/align/read_align.rs b/src/align/read_align.rs index f959cd0..897fe00 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -374,11 +374,16 @@ pub fn align_read( }); } - // Deterministic primary tie-break (score, then a fixed positional order). + // Deterministic primary tie-break. STAR compares `maxScore`, then + // `gLength` — the alignment's genomic span — and leaves anything still + // tied to whichever window it reached first + // (`ReadAlign_stitchPieces.cpp:340`). The span is reproducible here; the + // window order is not, since windows are built per read in parallel, so + // the remaining keys are positional and fixed (see DIVERGENCE.md §1.1). transcripts.sort_by(|a, b| { b.score .cmp(&a.score) - .then_with(|| a.n_junction.cmp(&b.n_junction)) + .then_with(|| (a.genome_end - a.genome_start).cmp(&(b.genome_end - b.genome_start))) .then_with(|| a.chr_idx.cmp(&b.chr_idx)) .then_with(|| a.genome_start.cmp(&b.genome_start)) .then_with(|| a.is_reverse.cmp(&b.is_reverse)) From efbaccea519234da33f9c9f05e89a000c5f57561 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 16:46:25 +0200 Subject: [PATCH 11/17] fix(seed): reproduce STAR's flagDirMap shortcut and chain-loop bound Two divergences in the MMP chain search, both read from ReadAlign_mapOneRead.cpp: 1. flagDirMap (:50, :61, :74). When the very first left-to-right search of a good piece maps it all the way to the piece end, STAR skips the istart == 0 chain of the right-to-left direction; chains from istart > 0 still run. rustar-aligner always ran both. The condition STAR writes is `Shift + L == splitR[1][ip]` with Shift still equal to the piece start, so the piece's offset in the concatenated read enters the sum: find_seeds_at takes that offset, and find_seeds keeps the old signature with an offset of 0. 2. The chain-loop bound. STAR continues while `istart*Lstart + Lmapped + P.seedMapMin < splitR[1][ip]`, i.e. while strictly more than seedMapMin bases remain. rustar-aligner continued while at least seedMapMin remained, so it ran one extra search on a piece with exactly seedMapMin bases left. Neither changes an output byte on the yeast tiers: a full-length match found by the reverse direction is dropped by the storeAligns duplicate rule anyway, and a 5-base MMP in a 12 Mb genome exceeds seedMultimapNmax and is never stored. They are here because the search should be STAR's search, not because they moved a number. The toy fixtures in the seed tests are 3-8 base reads, shorter than the default seedMapMin, so under the corrected bound STAR would not search them at all. They now pass --seedMapMin 1 through a toy_params helper so the assertions still exercise the MMP mechanics they were written for. --- src/align/read_align.rs | 9 +++- src/align/seed.rs | 104 +++++++++++++++++++++++++++++++++------- 2 files changed, 93 insertions(+), 20 deletions(-) diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 897fe00..5dd3caf 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -746,15 +746,20 @@ pub fn align_paired_read( // → Nstart=7, starts={0,43,...,129,...}). Using combined length creates a spurious // start at position 129 (between mates) that can produce anchors widening windows // beyond STAR's range, causing window overflow and eviction of valid 7M exon seeds. - let mut combined_seeds = Seed::find_seeds( + // The piece offsets are STAR's `splitR[0][ip]`: mate1 starts the concatenated + // read, mate2 starts one base past the spacer. Only the `flagDirMap` + // shortcut consults them. + let mut combined_seeds = Seed::find_seeds_at( &combined_read[..len1], + 0, index, params.seed_map_min, params, debug_name, )?; - let mut m2_seeds = Seed::find_seeds( + let mut m2_seeds = Seed::find_seeds_at( &combined_read[len1 + 1..], + len1 + 1, index, params.seed_map_min, params, diff --git a/src/align/seed.rs b/src/align/seed.rs index 7c6deb3..73d2710 100644 --- a/src/align/seed.rs +++ b/src/align/seed.rs @@ -48,6 +48,25 @@ impl Seed { min_seed_length: usize, params: &Parameters, debug_name: &str, + ) -> Result, Error> { + Self::find_seeds_at(read_seq, 0, index, min_seed_length, params, debug_name) + } + + /// [`find_seeds`](Self::find_seeds) with the piece's offset inside STAR's + /// concatenated `Read1[0]` buffer. + /// + /// STAR searches one *good piece* at a time (`splitR[0][ip]` = the piece's + /// start in the concatenated read, `splitR[1][ip]` = its length). Only the + /// `flagDirMap` shortcut in `ReadAlign_mapOneRead.cpp:74` reads the start + /// offset, so every other computation here is piece-local and `piece_start` + /// is threaded through purely to reproduce that condition. + pub fn find_seeds_at( + read_seq: &[u8], + piece_start: usize, + index: &GenomeIndex, + min_seed_length: usize, + params: &Parameters, + debug_name: &str, ) -> Result, Error> { let mut seeds = Vec::new(); let read_len = read_seq.len(); @@ -60,13 +79,15 @@ impl Seed { // while(istart*Lstart + Lmapped + seedMapMin < readLen) { ... Lmapped += L; } // Search L→R (forward direction on read): sparse chain search - search_direction_sparse( + let flag_dir_map = search_direction_sparse( read_seq, read_len, + piece_start, index, min_seed_length, params, false, + true, debug_name, &mut seeds, ); @@ -76,15 +97,22 @@ impl Seed { return Ok(seeds); } - // Search R→L (reverse direction on read): sparse chain search on RC read + // Search R→L (reverse direction on read): sparse chain search on RC read. + // `flag_dir_map` is STAR's `flagDirMap` (`ReadAlign_mapOneRead.cpp:50`, + // `:61`, `:74`): when the very first L→R search already mapped the piece + // to its end, the `istart == 0` chain in the reverse direction would + // re-derive the same maximal prefix, so STAR skips it. Chains from + // `istart > 0` still run. let rc_read = reverse_complement_read(read_seq); search_direction_sparse( &rc_read, read_len, + piece_start, index, min_seed_length, params, true, + flag_dir_map, debug_name, &mut seeds, ); @@ -232,14 +260,19 @@ struct MmpResult { fn search_direction_sparse( read_seq: &[u8], original_read_len: usize, + piece_start: usize, index: &GenomeIndex, min_seed_length: usize, params: &Parameters, is_rc: bool, + run_istart0: bool, debug_name: &str, seeds: &mut Vec, -) { +) -> bool { let read_len = read_seq.len(); + // STAR's `flagDirMap`, returned to the caller: stays true unless the first + // L→R search maps the piece all the way to its end. + let mut flag_dir_map = true; // STAR (ReadAlign_mapOneRead.cpp lines 41-42): // seedSearchStartLmax = min(P.seedSearchStartLmax, seedSearchStartLmaxOverLread*(Lread-1)) @@ -263,6 +296,12 @@ fn search_direction_sparse( let lstart = read_len / nstart; // STAR: Lstart = (splitR[1]-splitR[0]) / Nstart for istart in 0..nstart { + // STAR: `if (flagDirMap || istart>0)` — the reverse direction skips its + // istart == 0 chain when the forward direction already mapped the piece + // to its end. + if istart == 0 && !run_istart0 { + continue; + } let start_pos = (istart * lstart).min(read_len); let mut pos = start_pos; @@ -273,12 +312,13 @@ fn search_direction_sparse( if pos >= read_len { break; } - // Stop if remaining bases < seedMapMin (matches STAR's while condition: - // istart*Lstart + Lmapped + P.seedMapMin < splitR[1][ip]). - // STAR chains continue until only seedMapMin (5) bases remain, NOT - // seedSearchStartLmax (50). This allows chains to reach terminal small - // exons (e.g. 9M after intron) near the read end. - if read_len - pos < min_seed_length { + // STAR's while condition, verbatim: + // istart*Lstart + Lmapped + P.seedMapMin < splitR[1][ip] + // i.e. keep searching while *more than* seedMapMin bases remain. + // rustar-aligner used `remaining < min_seed_length` here, which ran + // one extra search when exactly seedMapMin bases were left and + // pushed a seed STAR never stores. + if read_len - pos <= params.seed_map_min { break; } @@ -297,6 +337,18 @@ fn search_direction_sparse( ); } + // STAR (`ReadAlign_mapOneRead.cpp:74`): on the very first forward + // search of the piece, a match that reaches the piece end means the + // reverse direction has nothing new to find from istart == 0. + // The comparison is STAR's own — `Shift + L == splitR[1][ip]`, with + // `Shift == splitR[0][ip]` at that point — so for a mate that does + // not start at offset 0 it is the mate's *global* start that enters + // the sum, and the shortcut effectively never fires there. + if !is_rc && istart == 0 && pos == start_pos && piece_start + result.advance == read_len + { + flag_dir_map = false; + } + if let Some(mut seed) = result.seed { // Apply seedSearchLmax cap if params.seed_search_lmax > 0 && seed.length > params.seed_search_lmax { @@ -313,14 +365,16 @@ fn search_direction_sparse( seeds.push(seed); if seeds.len() >= params.seed_per_read_nmax { - return; + return flag_dir_map; } } pos += result.advance; // Always advance by MMP length (matches STAR) - // Remaining-length check at loop top: stop when < seedMapMin bases remain + // Remaining-length check at loop top: stop when <= seedMapMin bases remain } } + + flag_dir_map } /// Find a seed starting at a specific position in the read. @@ -772,11 +826,23 @@ mod tests { Parameters::parse_from(full_args) } + /// Parameters for the toy fixtures below, whose reads are 3-8 bases. + /// + /// STAR's chain loop is `while (istart*Lstart + Lmapped + P.seedMapMin < + /// splitR[1][ip])`, so a piece of at most `seedMapMin` bases is never + /// searched at all. At the default `seedMapMin = 5` these reads would + /// produce no seeds for that reason alone, which is not what they are + /// testing; lower the threshold so the MMP mechanics are what the + /// assertions actually exercise. + fn toy_params() -> Parameters { + params(&["--seedMapMin", "1"]) + } + #[test] fn find_exact_match() { let index = make_test_index("ACGTACGT"); let read = encode_sequence("ACGT"); - let params = params(&["--runMode", "alignReads"]); + let params = toy_params(); let seeds = Seed::find_seeds(&read, &index, 4, ¶ms, "").unwrap(); @@ -792,7 +858,7 @@ mod tests { fn min_seed_length_filter() { let index = make_test_index("AAAAAAAA"); let read = encode_sequence("AAA"); - let params = params(&[]); + let params = toy_params(); // With min_seed_length=4, should find nothing (read is only 3bp) let seeds = Seed::find_seeds(&read, &index, 4, ¶ms, "").unwrap(); @@ -807,7 +873,7 @@ mod tests { fn no_match() { let index = make_test_index("ACAC"); let read = encode_sequence("GGGG"); - let params = params(&[]); + let params = toy_params(); let seeds = Seed::find_seeds(&read, &index, 2, ¶ms, "").unwrap(); @@ -819,7 +885,7 @@ mod tests { fn get_genome_positions() { let index = make_test_index("ACGTACGT"); let read = encode_sequence("ACGT"); - let params = params(&[]); + let params = toy_params(); let seeds = Seed::find_seeds(&read, &index, 4, ¶ms, "").unwrap(); assert!(!seeds.is_empty()); @@ -838,7 +904,7 @@ mod tests { fn test_single_end_mate_id() { let index = make_test_index("ACGTACGT"); let read = encode_sequence("ACGT"); - let params = params(&[]); + let params = toy_params(); let seeds = Seed::find_seeds(&read, &index, 4, ¶ms, "").unwrap(); assert!(!seeds.is_empty()); @@ -854,7 +920,7 @@ mod tests { let index = make_test_index("ACGTACGTTTGGCCAA"); let mate1 = encode_sequence("ACGT"); let mate2 = encode_sequence("TTGG"); - let params = params(&[]); + let params = toy_params(); let seeds = Seed::find_paired_seeds(&mate1, &mate2, &index, 4, ¶ms).unwrap(); @@ -881,7 +947,7 @@ mod tests { let index = make_test_index("ACGTACGT"); let mate1 = encode_sequence("ACGT"); let mate2 = encode_sequence("ACGT"); - let params = params(&[]); + let params = toy_params(); let seeds = Seed::find_paired_seeds(&mate1, &mate2, &index, 4, ¶ms).unwrap(); @@ -1023,10 +1089,12 @@ mod tests { search_direction_sparse( &rc_read, read.len(), + 0, &index, 4, ¶ms, true, + true, "", &mut rc_seeds, ); From f0ba573ba852a41119c76b3f79f23b4afd68d739 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 16:46:38 +0200 Subject: [PATCH 12/17] fix(align): keep sjA in the window-alignment dedup key STAR's assignAlignToWindow guards its duplicate test with both `aFrag == WA[iA][WA_iFrag]` and `sjA == WA[iA][WA_sjA]`. Two aligns on the same diagonal that were split out of different annotated junctions are not duplicates of each other. The pre-stitch diagonal dedup keyed on (diagonal, mate_id) only. A read whose last exon is a few bases long produces two entries covering the same donor-side bases, one tagged with the junction it arrives on and one with the junction it leaves by; the dedup collapsed them to the longest, which kept only the arriving tag. The micro-exon entry was then left without a partner carrying the same sjA, the annotated-junction path could not fire, and the terminal exon came out soft-clipped instead of spliced: `61M2448N36M3S` where STAR writes `61M2448N36M2483N3M`, five points of AS lost on a uniquely-mapped read. Adding sjA to the key is the whole fix. The dedup still collapses what it exists to collapse, which the second test asserts. The block is now dedup_wa_by_diagonal so the two cases can be tested directly rather than through a full window stitch. Measured on the annotated-junction tier (1000 simulated reads over yeast annotated junctions, index rebuilt from this branch), raw exact records: 954 -> 958, and the four that changed were the only remaining differences that were not equal-score ties. On 10k real single-end reads 8793 -> 8796. Both tiers now differ from STAR only where the two tools found alignments of identical score at different loci. --- src/align/stitch.rs | 205 +++++++++++++++++++++++++++++--------------- 1 file changed, 137 insertions(+), 68 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 5ff3596..b15e586 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -3023,6 +3023,74 @@ pub(crate) fn stitch_seeds_with_jdb_debug( transcripts } +/// Collapse window alignments that are redundant copies of one another before +/// the recursive stitcher runs. +/// +/// For each diagonal (`genome_pos - positive_strand_read_start`), overlapping +/// entries are merged into intervals and only the longest entry of each merged +/// interval is kept. Without this the stitcher's include/exclude recursion +/// explodes on the many redundant seeds that cover the same diagonal. +/// +/// The grouping key is `(diagonal, mate_id, sjA)`, which is what STAR's +/// `assignAlignToWindow` compares before treating two aligns as duplicates: it +/// guards the overlap test with both `aFrag == WA[iA][WA_iFrag]` and +/// `sjA == WA[iA][WA_sjA]`. Entries split out of *different* annotated +/// junctions are therefore never duplicates of each other, even when they sit +/// on the same diagonal and cover the same read bases. +fn dedup_wa_by_diagonal(wa_entries: &mut Vec, read_len: usize, is_rev: bool) { + use rustc_hash::{FxHashMap, FxHashSet}; + + type DiagKey = (i64, u8, i64); + type DiagSeeds = Vec<(usize, usize, usize)>; + let mut diag_seeds: FxHashMap = FxHashMap::default(); + for (idx, wa) in wa_entries.iter().enumerate() { + let ps = if is_rev { + read_len - (wa.length + wa.read_pos) + } else { + wa.read_pos + }; + let diag = wa.genome_pos as i64 - ps as i64; + diag_seeds + .entry((diag, wa.mate_id, wa.sj_a)) + .or_default() + .push((ps, ps + wa.length, idx)); + } + + let mut keep_indices = FxHashSet::default(); + for (_key, mut seeds) in diag_seeds { + seeds.sort_unstable(); + let mut merged_end = seeds[0].1; + let mut best_idx = seeds[0].2; + let mut best_len = seeds[0].1 - seeds[0].0; + + for &(s, e, idx) in &seeds[1..] { + if s <= merged_end { + // Overlapping: extend the interval, track the longest entry in it. + merged_end = merged_end.max(e); + let len = e - s; + if len > best_len { + best_len = len; + best_idx = idx; + } + } else { + // Disjoint: commit the previous interval's best and start a new one. + keep_indices.insert(best_idx); + merged_end = e; + best_len = e - s; + best_idx = idx; + } + } + keep_indices.insert(best_idx); + } + + let mut idx = 0usize; + wa_entries.retain(|_| { + let keep = keep_indices.contains(&idx); + idx += 1; + keep + }); +} + /// Shared core: preprocessing + recursive stitcher, returns working transcripts + context. #[allow(clippy::too_many_arguments)] pub(crate) fn stitch_seeds_core( @@ -3044,72 +3112,7 @@ pub(crate) fn stitch_seeds_core( // of repetitive seeds, matching STAR's WA_Anchor=2 "last anchor" logic. let mut wa_entries: Vec = cluster.alignments.clone(); - // Diagonal dedup: for each diagonal (genome_pos - ps_rstart), merge overlapping - // seeds into intervals, keeping only the longest seed per merged interval. - // This prevents combinatorial explosion in the recursive stitcher when many - // redundant seeds cover the same diagonal region. - // Uses positive-strand coordinates consistent with cluster_seeds overlap detection. - { - use rustc_hash::{FxHashMap, FxHashSet}; - let read_len = read_seq.len(); - let is_rev = cluster.is_reverse; - // For each (diagonal, mate_id) pair, find the longest seed per merged interval. - // STAR's assignAlignToWindow checks aFrag==WA[iA][WA_iFrag] before overlap test: - // seeds from different fragments are never treated as duplicates. - type DiagMateKey = (i64, u8); - type DiagSeeds = Vec<(usize, usize, usize)>; - let mut diag_seeds: FxHashMap = FxHashMap::default(); - for (idx, wa) in wa_entries.iter().enumerate() { - let ps = if is_rev { - read_len - (wa.length + wa.read_pos) - } else { - wa.read_pos - }; - let diag = wa.genome_pos as i64 - ps as i64; - diag_seeds - .entry((diag, wa.mate_id)) - .or_default() - .push((ps, ps + wa.length, idx)); - } - - let mut keep_indices = FxHashSet::default(); - for (_diag, mut seeds) in diag_seeds { - // Sort by start position - seeds.sort_unstable(); - // Merge intervals, keeping the index of the longest seed in each merged group - let mut merged_end = seeds[0].1; - let mut best_idx = seeds[0].2; - let mut best_len = seeds[0].1 - seeds[0].0; - - for &(s, e, idx) in &seeds[1..] { - if s <= merged_end { - // Overlapping — extend and track longest - merged_end = merged_end.max(e); - let len = e - s; - if len > best_len { - best_len = len; - best_idx = idx; - } - } else { - // New interval — commit previous best - keep_indices.insert(best_idx); - merged_end = e; - best_len = e - s; - best_idx = idx; - } - } - // Commit last group - keep_indices.insert(best_idx); - } - - // Retain only the kept indices - let mut idx = 0usize; - wa_entries.retain(|_| { - let keep = keep_indices.contains(&idx); - idx += 1; - keep - }); - } + dedup_wa_by_diagonal(&mut wa_entries, read_seq.len(), cluster.is_reverse); // STAR-faithful coordinate conversion for stitching: // STAR stores WA_gStart in FORWARD genome coordinates (converting RC positions via @@ -3185,8 +3188,15 @@ pub(crate) fn stitch_seeds_core( ); for (i, wa) in wa_entries.iter().enumerate().take(30) { eprintln!( - " wa[{}]: read_pos={}, sa_pos={}, genome_pos={}, length={}, anchor={}, mate={}", - i, wa.read_pos, wa.sa_pos, wa.genome_pos, wa.length, wa.is_anchor, wa.mate_id + " wa[{}]: read_pos={}, sa_pos={}, genome_pos={}, length={}, anchor={}, mate={}, sjA={}", + i, + wa.read_pos, + wa.sa_pos, + wa.genome_pos, + wa.length, + wa.is_anchor, + wa.mate_id, + wa.sj_a ); } if wa_entries.len() > 30 { @@ -3380,6 +3390,65 @@ mod tests { use crate::index::sa_index::SaIndex; use crate::index::suffix_array::SuffixArray; + fn wa(read_pos: usize, length: usize, genome_pos: u64, sj_a: i64) -> WindowAlignment { + WindowAlignment { + seed_idx: 0, + read_pos, + length, + genome_pos, + sa_pos: genome_pos, + n_rep: 1, + is_anchor: true, + mate_id: 2, + pre_ext_score: length as i32, + sj_a, + } + } + + /// A read whose last exon is a few bases long produces two window + /// alignments covering the same donor-side bases: one tagged with the + /// annotated junction the read arrives on, one tagged with the junction it + /// leaves by. STAR keeps both, because `assignAlignToWindow` only treats + /// two aligns as duplicates when their `sjA` match as well as their + /// fragment and diagonal. + /// + /// Keying the dedup on the diagonal alone dropped the departing copy, which + /// left the micro-exon entry without a partner carrying the same `sjA`, so + /// the terminal exon came out soft-clipped instead of spliced. + #[test] + fn same_diagonal_entries_from_different_junctions_both_survive() { + let mut entries = vec![ + wa(61, 36, 14_434_354, 348), + wa(61, 36, 14_434_354, 349), + wa(97, 3, 14_436_873, 349), + ]; + dedup_wa_by_diagonal(&mut entries, 100, false); + + assert_eq!(entries.len(), 3, "no entry may be dropped: {entries:?}"); + let tags: Vec = entries.iter().map(|e| e.sj_a).collect(); + assert!(tags.contains(&348) && tags.contains(&349)); + assert!( + entries.iter().any(|e| e.length == 3 && e.sj_a == 349), + "the micro-exon entry must keep a partner with the same sjA" + ); + } + + /// The dedup still collapses what it is there to collapse: redundant + /// overlapping copies on one diagonal that carry the same `sjA`. + #[test] + fn same_diagonal_overlapping_entries_collapse_to_the_longest() { + let mut entries = vec![ + wa(0, 20, 1000, -1), + wa(5, 40, 1005, -1), + wa(80, 10, 1080, -1), + ]; + dedup_wa_by_diagonal(&mut entries, 100, false); + + assert_eq!(entries.len(), 2); + assert_eq!(entries[0].length, 40, "longest of the merged interval wins"); + assert_eq!(entries[1].read_pos, 80, "the disjoint interval survives"); + } + fn make_simple_index() -> GenomeIndex { // Simple genome: ACGTACGTNN (10 bases) let seq = vec![0, 1, 2, 3, 0, 1, 2, 3, 4, 4]; From 8e17d4de112ad03e872b842a5027cc2889e54092 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 16:48:51 +0200 Subject: [PATCH 13/17] docs(changelog): record the aligner-core faithfulness results Adds the four fixes that were not listed (one-base introns, annotated gaps never being deletions, the gLength primary tie-break, the sjA dedup key) and replaces the single measurement with a before/after table over the three tiers. Both columns are measured runs, each binary against an index built by that binary: main produces 361 yeast junctions and this branch 364, so comparing them against one shared index would have credited the branch with the wrong numbers. That is also why the note about rebuilding the index is there: a cached index reproduces main's results silently. --- CHANGELOG.md | 37 +++++++++++++++++++++++++++++++++---- 1 file changed, 33 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2196c27..47321ac 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -33,10 +33,39 @@ Sections commonly used: Features, Bug fixes, Other changes. `--seedSplitMin` are accepted. (`--alignEndsType` itself came from #145.) - The genomic-length penalty is applied inside the stitch recursion, where it can affect which transcripts survive, rather than only at finalisation. - - On 1000 junction-spanning yeast reads against native STAR 2.7.11b, records - identical in FLAG/RNAME/POS/CIGAR rise from 835 to 933 out of 1000, with no - read regressing. Output against an unannotated index is unchanged. + - A one-base intron is a junction. `GTF_transcriptGeneSJ.cpp:123` skips an + exon pair only when `exon2.start <= exon1.end + 1`; rustar-aligner also + skipped `exon2.start == exon1.end + 2`, so three yeast junctions were + missing from the index and reads crossing them were written with `1D` + rather than `1N`. + - An annotated gap is never a deletion. `--alignIntronMin` decides between + `N` and `D` only for unannotated gaps, so the sjdb is consulted before the + length test at both CIGAR decision sites. + - The primary transcript is chosen by score, then by genomic span, which is + the comparison `ReadAlign_stitchPieces.cpp:340` makes + (`maxScore`, then `gLength`). + - Window alignments split out of different annotated junctions are no longer + treated as duplicates of one another. STAR's `assignAlignToWindow` compares + `sjA` as well as fragment and diagonal; keying the pre-stitch dedup on the + diagonal alone unpaired the two halves of a terminal micro-exon, which then + came out soft-clipped instead of spliced. + + Measured against native STAR 2.7.11b, raw exact records (FLAG, RNAME, POS, + MAPQ, CIGAR, NH, AS, NM), with the index rebuilt from this branch: + + | tier | records | main | this branch | + |---|---|---|---| + | junction-spanning yeast reads | 1000 | 835 | 958 | + | 10k single-end yeast reads | 8927 | 8785 | 8796 | + | 10k paired-end yeast reads | 16782 | 16572 | 16582 | + + No read regressed. Every remaining difference on all three tiers is a read + where both tools found alignments of identical score at different loci. + Output against an unannotated index is unchanged. + + Reproducing this needs an index built by this branch: the one-base-intron + fix changes `genomeGenerate` output (364 yeast junctions, not 361), and a + cached index built before it silently reproduces the old numbers. - **STARsolo single-cell quantification (`--soloType`)** — the 10x Chromium / plate-based count-matrix pipeline, ported from STAR and From 92b8fb1a3568a73f609df022848f7f71ba56856e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 16:53:34 +0200 Subject: [PATCH 14/17] docs(divergence): correct the tie-break entry and the residual list MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three corrections to §1.1 and §5, all from measured runs on this branch. §1.1 said the default primary is "max score, smaller genomic length, earliest discovered, which is STAR-faithful". The first two keys are STAR's (ReadAlign_stitchPieces.cpp:340). The third is not: STAR takes the earliest window in its iteration order, this codebase takes the smallest genomic position. Calling that STAR-faithful was wrong. It also implied that matching STAR's window order would forfeit thread-invariance. It would not: STAR's window order is deterministic per read, and windows are built per read here too. The reason the positional key stays is that the alternative was tried and measured worse (+21 on the annotated-junction tier, -142 on 10k SE), because this codebase's anchor iteration order is not STAR's PC order. That is now what the entry says. §5 listed three residual reads that no longer exist. Replaced with the three that do, including one where rustar-aligner scores higher than STAR, which was not previously recorded anywhere. --- DIVERGENCE.md | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..afd8b72 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -24,9 +24,11 @@ Divergences are grouped by kind: **Why.** Determinism and thread-count invariance: the same read produces the same primary regardless of `--runThreadN`. STAR's exact mt19937 stream cannot be reproduced under per-read parallelism, and matching it would forfeit reproducibility. -**Impact.** With the default `--outMultimapperOrder Old_2.4`, **no RNG is consulted at all** — the primary is the deterministic best alignment (max score → smaller genomic length → earliest discovered), which is STAR-faithful. The divergence is observable only under `--outMultimapperOrder Random`, and only in *which* equal-scoring locus is marked primary — never in the set of reported alignments. +**Impact.** With the default `--outMultimapperOrder Old_2.4`, **no RNG is consulted at all**: the primary is the deterministic best alignment. The first two keys are STAR's own (`ReadAlign_stitchPieces.cpp:340` compares `maxScore`, then `gLength`). Where those two tie, STAR takes the earliest window in its iteration order and rustar-aligner takes the smallest genomic position; that last key is a divergence, and it is the one described below. The RNG divergence proper is observable only under `--outMultimapperOrder Random`, and in both cases only in *which* equal-scoring locus is marked primary, never in the set of reported alignments. -This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast benchmark, 299 SE and 475 PE primary-selection differences are all genuine ties: both tools find the identical alignment set, and differ only in which equal-scoring member is primary (from SA-iteration order or the RNG-seed difference above). Excluding those ties, SE is 99.815% and PE 99.883% exact. +**On the residual ties.** STAR's window order is deterministic per read (anchor pieces in `PC` order, positions in suffix-array order) and rustar-aligner also builds windows per read, so reproducing it would not cost thread-invariance. It has been measured, not assumed: substituting seed discovery order for the positional key gained 21 reads on the annotated-junction tier and lost 142 on the 10k single-end tier, because this codebase's anchor iteration order is not STAR's `PC` order. Closing the gap means matching how MMP results are recorded, which has not been done. Until it is, the positional key stays, because it is total, cheap and thread-invariant. + +This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast benchmark, 130 of the 131 differing single-end records and 197 of the 202 differing paired-end mate records are genuine ties: both tools find the identical alignment set and differ only in which equal-scoring member is primary. On the annotated-junction tier all 42 remaining differences are of this kind, which is what its 100.000% tie-adjusted figure means. The handful that are not ties are in [§5](#5-known-residual-single-read-differences). Raw counts are reported alongside, never the tie-adjusted figure alone. **Source.** `src/rng.rs`, `src/align/read_align.rs` (`per_read_seed`, `shuffle_tied_prefix`), `src/params/mod.rs` (`MultimapperOrder`). STAR: `ReadAlign_multMapSelect.cpp`, `ReadAlignChunk` RNG seeding. @@ -85,9 +87,11 @@ rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` These are **not** deliberate divergences — they are tracked residual diffs on the 10k yeast benchmark, kept here for completeness. Each is a single read; none is a systematic behaviour difference. -- **1 SE CIGAR-only diff** — `ERR12389696.13573895`: both tools align to XV:218357, MAPQ 255, identical score (AS=133), but place a 1-base insertion differently (`100M1I45M4S` vs STAR's `108M1I37M4S`). The 71-base seed is found at a different position within a long homopolymer run (a seed-level tie); resolving it requires reproducing STAR's exact Lmapped chain path. -- **1 STAR-only PE mate** — `ERR12389696.18919121`: an SA-level difference. -- **1 rustar-aligner-only PE mate** — `ERR12389696.6302610`: a pre-existing false positive. +- **1 SE insertion-placement diff**: `ERR12389696.20597455`, XIV:545446, MAPQ 255, identical score (AS=143), with the 1-base insertion four bases to the left of STAR's choice (`25M1I124M` against STAR's `29M1I120M`). The same read shows on the PE tier. Insertion placement inside `stitchAlignToTranscript`, not a tie-break. +- **1 PE pair where rustar-aligner scores lower**: `ERR12389696.11539725`, AS 224 against STAR's 235; mate 1 is soft-clipped to `31S95M459N24M` where STAR reaches back across a second junction to `22S10M468N94M459N24M`. +- **1 PE pair where rustar-aligner scores higher**: `ERR12389696.4972950`, AS 260 against STAR's 248; mate 2 is spliced `1S33M72N50M186N65M1S` where STAR soft-clips 27 bases to `27S122M1S`. Recorded here rather than in [§2](#2-cases-where-rustar-aligner-outperforms-star) because it is a single observed read, not a characterised behaviour. + +Counts are from the 10k yeast SE and PE tiers with an index built by the same binary; the annotated-junction tier has no non-tie residual. See `CLAUDE.md` ("Known Issues" / "PE Status") for the current status of these. From e43ff5fb0a3abacb2ae5a18de5cae1e066130e6f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 17:10:31 +0200 Subject: [PATCH 15/17] fix(seed): compare a genome spacer as the larger byte, as STAR does `compare_seq_to_genome` returned "the query is greater" whenever it ran into a genome byte of 5. STAR's `compareSeqToGenome` does no such thing: it compares the raw bytes, and since read bases are 0-3 while `N` is 4 and the spacer is 5, a read base is never the greater one. Every branch of STAR's function therefore reports the query as sorting *before* such a suffix. That value is what orders the binary search in `max_mappable_length`, so it has to agree with the order the suffix array was built in. It did not, and the search read the array as unsorted at exactly the places where the sjdb inserts sit: a probe landing on a spacer inside an insert claimed the query was greater, the search dropped everything below that probe, and the entry holding the real maximum went with it. On ERR12389696.20597455 the suffix array holds a 29-base match at the right locus, one slot below a 24-base match inside an sjdb insert whose comparison hit a spacer. The MMP came back 24. The read then stitched two shorter pieces and placed its 1-base insertion four bases left of STAR's, `25M1I124M` against `29M1I120M`, at the same total score. Raw exact records on 10k yeast reads: single-end 8796 -> 8797, paired-end 16582 -> 16583. The annotated-junction tier is unchanged at 958/1000. Locked by `a_query_running_into_the_padding_sorts_before_that_suffix`. --- src/align/seed.rs | 56 +++++++++++++++++++++++++++++++++++++---------- 1 file changed, 45 insertions(+), 11 deletions(-) diff --git a/src/align/seed.rs b/src/align/seed.rs index 73d2710..572cf5e 100644 --- a/src/align/seed.rs +++ b/src/align/seed.rs @@ -559,10 +559,9 @@ fn compare_seq_to_genome( let genome_chunk = &genome_slice_all[genome_start + i..genome_start + simd_end]; if let Some(off) = crate::align::simd_scan::find_stop(read_chunk, genome_chunk) { let genome_base = genome_chunk[off]; - if genome_base >= 5 { - return (i + off, true); - } let read_base = read_chunk[off]; + // A spacer (5) is larger than any read base, so the query + // sorts *before* this suffix. See the scalar branch below. return (i + off, read_base > genome_base); } match_len = simd_end; @@ -575,19 +574,27 @@ fn compare_seq_to_genome( let genome_idx = genome_start + i; if genome_idx >= index.genome.sequence.len() { - // Past end of genome array — treat like padding (STAR: comp_res > 0) - return (match_len, true); + // Past the end of the genome array. Treated as a spacer, so the + // query sorts before this suffix, same as the branch below. + return (match_len, false); } let genome_base = index.genome.sequence.base(genome_idx); - - if genome_base >= 5 { - // Padding character — STAR returns comp_res > 0 (read > genome) - return (match_len, true); - } - let read_base = read_seq[read_pos + i]; + // The comparison result is what orders the binary search in + // `max_mappable_length`, so it has to agree with the order the suffix + // array was built in. STAR compares the raw bytes + // (`compareSeqToGenome`, `SuffixArrayFuns.cpp`): read bases are 0-3, + // `N` is 4 and the spacer is 5, so a read base is never greater than + // either of those and the query always sorts *before* such a suffix. + // `comp_res` was hardcoded to `true` for the spacer, the opposite. + // + // It matters where the sjdb inserts sit in the suffix array. A probe + // landing on a spacer inside an insert reported "query is greater", + // the search discarded everything below that probe, and the entry + // holding the real maximum went with it: on ERR12389696.20597455 the + // MMP came back 24 where the genome matches 29. if read_base != genome_base { return (match_len, read_base > genome_base); } @@ -1156,6 +1163,33 @@ mod tests { ); } + /// A read base is 0-3, `N` is 4 and the spacer is 5, so a read base is + /// never greater than either. When the comparison runs off a chromosome + /// into the padding, the query sorts *before* that suffix, and + /// `compare_seq_to_genome` must say so: it is this value that orders the + /// binary search in `max_mappable_length`, so getting it backwards makes + /// the search discard the half of the range holding the real maximum. + #[test] + fn a_query_running_into_the_padding_sorts_before_that_suffix() { + let index = make_test_index("ACGTACGT"); + let read = encode_sequence("ACGTACGTA"); + + // The suffix array entry for forward position 0. + let sa_idx = (0..index.suffix_array.len()) + .find(|&i| { + let (pos, rev) = index.suffix_array.decode(index.suffix_array.get(i)); + pos == 0 && !rev + }) + .expect("forward position 0 must be in the suffix array"); + + let (len, comp) = compare_seq_to_genome(&read, 0, &index, sa_idx, 0); + assert_eq!(len, 8, "the whole chromosome matches before the padding"); + assert!( + !comp, + "at the padding the read base is the smaller one, so the query sorts first" + ); + } + #[test] fn test_rc_seed_genome_positions() { // Genome: AACCTTGG, read RC = CCAAGGTT From c6bd0e71d11e5eb2bad28c0524c89e3e9f5f016f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 19:06:30 +0200 Subject: [PATCH 16/17] docs: fold the spacer-comparison fix into the results MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Updates the measured table for the fix in e43ff5f (single-end 8796 to 8797, paired-end 16582 to 16583) and rewrites the residual list in DIVERGENCE.md §5 now that the insertion-placement read is gone. Two residuals remain on the paired-end tier, both characterised from runs rather than described from the code: ERR12389696.11539725, where rustar-aligner scores 224 against STAR's 235. Both mates cross the same junction, so the pair needs annotated junction 120 twice. Each half is reachable on its own (256 and 254) but the combination that would score 266 is not among the 84 transcripts the recursion enumerates. Raising alignTranscriptsPerWindowNmax changes nothing, so it is not the cap. Cause unknown; that is stated rather than guessed at. ERR12389696.4972950, where rustar-aligner scores 260 against STAR's 248 by splicing where STAR soft-clips. Both junctions are novel, so the only claim recorded is the score; neither alignment was checked against the transcript, and it stays in §5 rather than §2 for that reason. --- CHANGELOG.md | 9 +++++++-- DIVERGENCE.md | 9 ++++----- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 47321ac..6605b33 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -49,6 +49,11 @@ Sections commonly used: Features, Bug fixes, Other changes. `sjA` as well as fragment and diagonal; keying the pre-stitch dedup on the diagonal alone unpaired the two halves of a terminal micro-exon, which then came out soft-clipped instead of spliced. + - A genome spacer compares as the larger byte in the seed search, as it does + in STAR's `compareSeqToGenome`. It was hardcoded to the opposite, which + made the binary search in `max_mappable_length` read the suffix array as + unsorted wherever the sjdb inserts sit and drop the half holding the real + maximum, returning a short MMP. Measured against native STAR 2.7.11b, raw exact records (FLAG, RNAME, POS, MAPQ, CIGAR, NH, AS, NM), with the index rebuilt from this branch: @@ -56,8 +61,8 @@ Sections commonly used: Features, Bug fixes, Other changes. | tier | records | main | this branch | |---|---|---|---| | junction-spanning yeast reads | 1000 | 835 | 958 | - | 10k single-end yeast reads | 8927 | 8785 | 8796 | - | 10k paired-end yeast reads | 16782 | 16572 | 16582 | + | 10k single-end yeast reads | 8927 | 8785 | 8797 | + | 10k paired-end yeast reads | 16782 | 16572 | 16583 | No read regressed. Every remaining difference on all three tiers is a read where both tools found alignments of identical score at different loci. diff --git a/DIVERGENCE.md b/DIVERGENCE.md index afd8b72..4ce46ee 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -28,7 +28,7 @@ Divergences are grouped by kind: **On the residual ties.** STAR's window order is deterministic per read (anchor pieces in `PC` order, positions in suffix-array order) and rustar-aligner also builds windows per read, so reproducing it would not cost thread-invariance. It has been measured, not assumed: substituting seed discovery order for the positional key gained 21 reads on the annotated-junction tier and lost 142 on the 10k single-end tier, because this codebase's anchor iteration order is not STAR's `PC` order. Closing the gap means matching how MMP results are recorded, which has not been done. Until it is, the positional key stays, because it is total, cheap and thread-invariant. -This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast benchmark, 130 of the 131 differing single-end records and 197 of the 202 differing paired-end mate records are genuine ties: both tools find the identical alignment set and differ only in which equal-scoring member is primary. On the annotated-junction tier all 42 remaining differences are of this kind, which is what its 100.000% tie-adjusted figure means. The handful that are not ties are in [§5](#5-known-residual-single-read-differences). Raw counts are reported alongside, never the tie-adjusted figure alone. +This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast benchmark, all 130 differing single-end records and 197 of the 201 differing paired-end mate records are genuine ties: both tools find the identical alignment set and differ only in which equal-scoring member is primary. On the annotated-junction tier all 42 remaining differences are of this kind, which is what its 100.000% tie-adjusted figure means. The handful that are not ties are in [§5](#5-known-residual-single-read-differences). Raw counts are reported alongside, never the tie-adjusted figure alone. **Source.** `src/rng.rs`, `src/align/read_align.rs` (`per_read_seed`, `shuffle_tied_prefix`), `src/params/mod.rs` (`MultimapperOrder`). STAR: `ReadAlign_multMapSelect.cpp`, `ReadAlignChunk` RNG seeding. @@ -87,11 +87,10 @@ rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` These are **not** deliberate divergences — they are tracked residual diffs on the 10k yeast benchmark, kept here for completeness. Each is a single read; none is a systematic behaviour difference. -- **1 SE insertion-placement diff**: `ERR12389696.20597455`, XIV:545446, MAPQ 255, identical score (AS=143), with the 1-base insertion four bases to the left of STAR's choice (`25M1I124M` against STAR's `29M1I120M`). The same read shows on the PE tier. Insertion placement inside `stitchAlignToTranscript`, not a tie-break. -- **1 PE pair where rustar-aligner scores lower**: `ERR12389696.11539725`, AS 224 against STAR's 235; mate 1 is soft-clipped to `31S95M459N24M` where STAR reaches back across a second junction to `22S10M468N94M459N24M`. -- **1 PE pair where rustar-aligner scores higher**: `ERR12389696.4972950`, AS 260 against STAR's 248; mate 2 is spliced `1S33M72N50M186N65M1S` where STAR soft-clips 27 bases to `27S122M1S`. Recorded here rather than in [§2](#2-cases-where-rustar-aligner-outperforms-star) because it is a single observed read, not a characterised behaviour. +- **1 PE pair where rustar-aligner scores lower**: `ERR12389696.11539725`, AS 224 against STAR's 235. Both mates cross the same two annotated junctions, and the pair needs junction 120 twice, once per mate. rustar-aligner's stitch recursion reaches each half on its own (mate 2 with its 10-base first exon, or mate 1 with its 11-base first exon) but never both in one transcript: the two best it produces score 256 and 254, and the combination that would score 266 is not among the 84 it enumerates. It emits `31S95M459N24M` for mate 1 where STAR emits `22S10M468N94M459N24M`. The mechanism is characterised; the cause is not, and it is not the per-window transcript cap, which was ruled out by raising it. +- **1 PE pair where rustar-aligner scores higher**: `ERR12389696.4972950`, AS 260 against STAR's 248; mate 2 is spliced `1S33M72N50M186N65M1S` where STAR soft-clips 27 bases to `27S122M1S`. Both junctions are novel, not in the sjdb, so "higher-scoring" is all that is claimed here; neither alignment has been verified against the transcript. Recorded here rather than in [§2](#2-cases-where-rustar-aligner-outperforms-star) for that reason. -Counts are from the 10k yeast SE and PE tiers with an index built by the same binary; the annotated-junction tier has no non-tie residual. +Counts are from the 10k yeast SE and PE tiers with an index built by the same binary. The annotated-junction tier has no non-tie residual, and neither does the single-end tier. See `CLAUDE.md` ("Known Issues" / "PE Status") for the current status of these. From 4978c89ed964199d3aecfbd0f3323e7622dce4a0 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 19:53:24 +0200 Subject: [PATCH 17/17] docs(divergence): trace the remaining PE pair to the seed search ERR12389696.11539725 was recorded as "mechanism characterised, cause unknown". The cause is now traced. Both tools build the same window and agree on every align in it but one. STAR holds two aligns on a single diagonal, r172 g5083545 L11 and r174 g5083547 L9; rustar-aligner holds only the first. stitchAlignToTranscript.cpp:352 refuses to cross into the second mate when the mate's seed starts before the transcript's first exon, and that exon starts at g5083546, so the 11-base seed at g5083545 is one base too far left and only the 9-base seed at g5083547 can open the mate. STAR's extendAlign then walks it one base left onto g5083546. Without that align the two halves cannot be assembled, which is why the best rustar-aligner reaches are 256 and 254 where the pair would score 266. It is not a dedup difference. STAR's own assignAlignToWindow overlap test would collapse that pair, but it fires on insertion order, and a Gsj-crossing piece records its acceptor half at aRstart + aLengthD rather than at its own PC start. Reproducing it means reproducing STAR's PC order and its Gsj splits, which is the seed search, not the stitcher. Ruled out by measurement, not by reading: the per-window transcript cap (raising it to 10000 changes nothing, same 84 transcripts and 1460 recursions), qualitySplit (the pair has no N and no base below Q3, so both aligns share one fragment), and the diagonal dedup key. A rewrite of the dedup onto STAR's asymmetric overlap test was tried and reverted: its own test shows STAR's rule collapses the pair too when applied in sorted order, so the rewrite was not the fix it looked like. --- DIVERGENCE.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 4ce46ee..47e3082 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -87,7 +87,11 @@ rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` These are **not** deliberate divergences — they are tracked residual diffs on the 10k yeast benchmark, kept here for completeness. Each is a single read; none is a systematic behaviour difference. -- **1 PE pair where rustar-aligner scores lower**: `ERR12389696.11539725`, AS 224 against STAR's 235. Both mates cross the same two annotated junctions, and the pair needs junction 120 twice, once per mate. rustar-aligner's stitch recursion reaches each half on its own (mate 2 with its 10-base first exon, or mate 1 with its 11-base first exon) but never both in one transcript: the two best it produces score 256 and 254, and the combination that would score 266 is not among the 84 it enumerates. It emits `31S95M459N24M` for mate 1 where STAR emits `22S10M468N94M459N24M`. The mechanism is characterised; the cause is not, and it is not the per-window transcript cap, which was ruled out by raising it. +- **1 PE pair where rustar-aligner scores lower**: `ERR12389696.11539725`, AS 224 against STAR's 235, emitting `31S95M459N24M` for mate 1 where STAR emits `22S10M468N94M459N24M`. + + Traced end to end. Both tools build the same window and agree on every seed in it but one. STAR's window holds two aligns on a single diagonal, `r172 g5083545 L11` and `r174 g5083547 L9`; rustar-aligner holds only the first. `stitchAlignToTranscript.cpp:352` refuses to cross into the second mate when the mate's seed starts before the transcript's first exon, and the first exon here starts at `g5083546`: the 11-base seed at `g5083545` is one base too far left, so only the 9-base seed at `g5083547` can open the mate, after which STAR's `extendAlign` walks it one base left to `g5083546`. Without that align rustar-aligner cannot assemble the two halves, and the best it reaches are 256 and 254 where the pair would score 266. + + The missing align is not a dedup difference. STAR's `assignAlignToWindow` overlap test would itself collapse the pair, but it fires on insertion order, and a Gsj-crossing piece records its acceptor half at `aRstart + aLengthD` rather than at its own `PC` start. Reproducing that means reproducing STAR's `PC` order and its Gsj splits, which is the seed-search side, not the stitcher. Ruled out along the way: the per-window transcript cap (raising it changes nothing), `qualitySplit` (the pair has no `N` and no base below Q3, so both aligns share one fragment), and the diagonal dedup key. - **1 PE pair where rustar-aligner scores higher**: `ERR12389696.4972950`, AS 260 against STAR's 248; mate 2 is spliced `1S33M72N50M186N65M1S` where STAR soft-clips 27 bases to `27S122M1S`. Both junctions are novel, not in the sjdb, so "higher-scoring" is all that is claimed here; neither alignment has been verified against the transcript. Recorded here rather than in [§2](#2-cases-where-rustar-aligner-outperforms-star) for that reason. Counts are from the 10k yeast SE and PE tiers with an index built by the same binary. The annotated-junction tier has no non-tie residual, and neither does the single-end tier.