From 460ddbae619806283d93e506887f41d96409a736 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 13:40:16 +0100 Subject: [PATCH 01/17] nM: STAR's accumulated mismatch count, not a recount STAR's trA.nMM is summed while the transcript is built and never recounted. When a junction scan moves the junction back into exon A (jR < 0), STAR counts those bases again against the acceptor without taking back the mismatches A had already counted (stitchAlignToTranscript.cpp, the commented-out nMM--). finalize_transcript now reports the stitched count plus its extensions; the two halves of a pair carry the pair's count once. Fixes nM on chr21 ERR188021.1648513/1662773 and 21 fusion pairs (DIVERGENCE.md 5). Co-Authored-By: Claude Opus 5.5 --- src/align/read_align.rs | 4 ++ src/align/stitch.rs | 103 ++++++++++++++++++---------------------- 2 files changed, 51 insertions(+), 56 deletions(-) diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 8fc7feb..b0dc5df 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -864,6 +864,10 @@ pub fn align_paired_read( ) else { continue; }; + // Both halves carry the pair's stitched mismatch count (STAR's trA.nMM, + // which also bounds each extension); keep it once, on mate1, so that + // the mates sum to STAR's nMM for the pair. + t2.n_mismatch -= m2_wt.n_mismatch; any_transcript = true; if stitch_is_reverse { diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 09a2d79..a846893 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -137,53 +137,6 @@ fn score_region( /// staging copy is amortized, small enough to sit on the stack. const SCORE_REGION_CHUNK: usize = 256; -fn count_mismatches( - read_seq: &[u8], - cigar_ops: &[cigar::Op], - genome_start: u64, - read_start: usize, - index: &GenomeIndex, - is_reverse: bool, -) -> u32 { - use cigar::op::Kind; - // Add n_genome offset for reverse-strand genome access - let genome_offset = if is_reverse { index.genome.n_genome } else { 0 }; - - let mut n_mismatch = 0u32; - let mut read_pos = read_start; - let mut genome_pos = genome_start; - - for op in cigar_ops { - match op.kind() { - Kind::Match | Kind::SequenceMatch | Kind::SequenceMismatch => { - for _i in 0..op.len() { - if read_pos < read_seq.len() { - let read_base = read_seq[read_pos]; - if let Some(genome_base) = index.genome.get_base(genome_pos + genome_offset) - && read_base != genome_base - && read_base != 4 - && genome_base != 4 - { - n_mismatch += 1; - } - } - read_pos += 1; - genome_pos += 1; - } - } - Kind::Insertion | Kind::SoftClip => { - read_pos += op.len(); - } - Kind::Deletion | Kind::Skip => { - genome_pos += op.len() as u64; - } - Kind::HardClip | Kind::Pad => {} - } - } - - n_mismatch -} - /// Result of extending an alignment into flanking regions #[derive(Debug, Clone)] struct ExtendResult { @@ -2199,15 +2152,14 @@ pub(crate) fn finalize_transcript( // Adjust score: add both left and right extensions (STAR adds both at finalization) let adjusted_score = wt.score + left_extend.max_score + right_extend.max_score; - // Count mismatches — MUST be called BEFORE CIGAR reversal - let n_mismatch = count_mismatches( - read_seq, - &final_cigar, - adjusted_genome_start, - 0, - index, - cluster.is_reverse, - ); + // STAR's nMM is accumulated while the transcript is built (stitch gaps, + // junction scans, extensions), never recounted from the final alignment. + // The two differ when a junction scan moves the junction back into exon A + // (jR < 0) over bases A had already counted as mismatches against the + // other flank: STAR counts them again against the acceptor and does not + // take the first count back (stitchAlignToTranscript.cpp, the commented-out + // `nMM--`). `wt.n_mismatch` keeps that count, as nMM does. + let n_mismatch = wt.n_mismatch + left_extend.n_mismatch + right_extend.n_mismatch; // Reverse CIGAR for reverse strand if cluster.is_reverse { @@ -4316,6 +4268,45 @@ mod tests { assert_ne!(additive_buggy, replacement_correct); } + /// `nM` is STAR's `trA.nMM` as accumulated while stitching, not a recount of + /// the final alignment: a mismatch the stitcher counted stays counted even + /// where the final exons match the genome (the junction-scan case in + /// `finalize_transcript`). + #[test] + fn finalize_keeps_the_stitched_mismatch_count() { + let index = make_simple_index(); + let scorer = AlignmentScorer::from_params_minimal(); + let cluster = SeedCluster { + alignments: vec![], + chr_idx: 0, + genome_start: 0, + genome_end: 8, + is_reverse: false, + anchor_idx: 0, + anchor_bin: 0, + }; + let read = [0u8, 1, 2, 3, 0, 1, 2, 3]; // genome[0..8], no mismatch + let mut wt = WorkingTranscript::new(); + wt.exons.push(ExonBlock { + read_start: 0, + read_end: 8, + genome_start: 0, + genome_end: 8, + mate_id: 2, + sj_a: None, + }); + wt.read_end = 8; + wt.genome_end = 8; + wt.score = 6; + wt.n_mismatch = 1; + let t = finalize_transcript( + &wt, &read, &index, &scorer, &cluster, false, false, false, 0, + ) + .expect("finalizes"); + assert_eq!(t.cigar_string(), "8M"); + assert_eq!(t.n_mismatch, 1); + } + /// STAR never stitches two seeds that sit on one diagonal (`rGap == gGap`): /// `stitchAlignToTranscript` leaves `jCan=999` there and its closing /// `jCan<0 || (jCan<7 && ...)` test returns -1000007. Adjacent seeds From 040b8477af2e3c14fd10b5e5fa09a9e002039dfb Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 13:40:16 +0100 Subject: [PATCH 02/17] Write the forward-strand mate of a pair first STAR's mate loop starts at exons[0], whose EX_iFrag is trOut.Str (ReadAlign_outputTranscriptSAM.cpp, ReadAlign_alignBAM.cpp): mate 1 when it is forward, mate 2 otherwise. rustar-aligner always wrote mate 1 first, about half of STAR's pairs; the per-read multiset comparison of the parity suite did not see it. The transcriptome path already did this. Co-Authored-By: Claude Opus 5.5 --- src/io/sam.rs | 83 +++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 81 insertions(+), 2 deletions(-) diff --git a/src/io/sam.rs b/src/io/sam.rs index 40a03e5..898bf33 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -526,8 +526,17 @@ impl SamWriter { if attrs.contains(SamAttributes::MC) { insert_mate_cigar_tags(&mut rec1, &mut rec2); } - records.push(rec1); - records.push(rec2); + // STAR writes the mate of the transcript's first exon first + // (ReadAlign_outputTranscriptSAM.cpp / ReadAlign_alignBAM.cpp: the mate loop + // starts at exons[0], whose EX_iFrag is trOut.Str): mate1 when it is on the + // forward strand, mate2 otherwise. + if paired_aln.mate1_transcript.is_reverse { + records.push(rec2); + records.push(rec1); + } else { + records.push(rec1); + records.push(rec2); + } } Ok(records) @@ -5370,6 +5379,76 @@ mod tests { assert!(!records_m2[1].flags().is_unmapped()); // mate2 is mapped } + /// STAR writes the mate on the forward strand first (the mate of the + /// transcript's first exon), so mate2 leads when mate1 is reverse. + #[test] + fn test_build_paired_records_forward_mate_first() { + use crate::align::transcript::Exon; + use cigar::op::{Kind, Op}; + + let genome = make_test_genome(); + let params = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "r1.fq", "r2.fq"]); + let tr = |start: u64, is_reverse: bool| Transcript { + output_rank: 0, + chr_idx: 0, + genome_start: start, + genome_end: start + 4, + is_reverse, + exons: vec![Exon { + genome_start: start, + genome_end: start + 4, + read_start: 0, + read_end: 4, + i_frag: 0, + }], + cigar: vec![Op::new(Kind::Match, 4)], + score: 4, + n_mismatch: 0, + n_gap: 0, + n_junction: 0, + junction_motifs: vec![], + junction_annotated: vec![], + junction_strand: vec![], + }; + let pair = |m1: Transcript, m2: Transcript| PairedAlignment { + output_rank: 0, + mate1_transcript: m1, + mate2_transcript: m2, + mate1_region: (0, 4), + mate2_region: (0, 4), + is_proper_pair: true, + insert_size: 8, + combined_wt_score: 8, + }; + let seq = vec![0, 1, 2, 3]; + let qual = vec![30, 30, 30, 30]; + let first_record = |p: PairedAlignment| { + let recs = SamWriter::build_paired_records( + "read1", + &seq, + &qual, + &seq, + &qual, + 0, + 0, + 0, + 0, + &[p], + &genome, + ¶ms, + 1, + ) + .unwrap(); + assert_eq!(recs.len(), 2); + recs[0].flags().is_first_segment() + }; + + // mate1 forward on the left, mate2 reverse: mate1 first. + assert!(first_record(pair(tr(0, false), tr(4, true)))); + // mate2 forward on the left, mate1 reverse: mate2 first. + assert!(!first_record(pair(tr(4, true), tr(0, false)))); + } + #[test] fn xs_strand_of_pair_uses_both_mates_junctions() { let spliced = spliced_gtag_transcript(); From b00182f8da22b611a5ef5abdd175c16a667e56db Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 13:40:16 +0100 Subject: [PATCH 03/17] Drop the mate-2 exon-subset dedup of finalized pairs It has no counterpart in STAR, whose window dedup (blocksOverlap in stitchWindowAligns.cpp) compares the whole pair. It collapsed two alignments with the same mate positions but a different indel or splice placement inside mate 1 (fusion fustut000009a03e, fustut000009a03f, fustut0000167ff6). Co-Authored-By: Claude Opus 5.5 --- src/align/read_align.rs | 66 ----------------------------------------- 1 file changed, 66 deletions(-) diff --git a/src/align/read_align.rs b/src/align/read_align.rs index b0dc5df..312d0f7 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -1151,72 +1151,6 @@ pub fn align_paired_read( }); } - // Post-finalization mate2-exon-subset dedup. - { - use crate::align::transcript::Exon; - let exons_subset = |b: &[Exon], a: &[Exon]| -> bool { - let total_b: u32 = b.iter().map(|e| (e.read_end - e.read_start) as u32).sum(); - if total_b == 0 { - return false; - } - let mut covered = 0u32; - for be in b { - let b_diag = be.genome_start as i64 - be.read_start as i64; - for ae in a { - let a_diag = ae.genome_start as i64 - ae.read_start as i64; - if a_diag == b_diag { - let r_start = be.read_start.max(ae.read_start); - let r_end = be.read_end.min(ae.read_end); - if r_start < r_end { - covered += (r_end - r_start) as u32; - } - } - } - } - covered == total_b - }; - let n = joint_pairs.len(); - let mut keep = vec![true; n]; - for i in 0..n { - if !keep[i] { - continue; - } - for j in 0..n { - if i == j || !keep[j] { - continue; - } - let same_pos = joint_pairs[i].mate1_transcript.chr_idx - == joint_pairs[j].mate1_transcript.chr_idx - && joint_pairs[i].mate1_transcript.genome_start - == joint_pairs[j].mate1_transcript.genome_start - && joint_pairs[i].mate1_transcript.genome_end - == joint_pairs[j].mate1_transcript.genome_end - && joint_pairs[i].mate1_transcript.is_reverse - == joint_pairs[j].mate1_transcript.is_reverse - && joint_pairs[i].mate2_transcript.genome_start - == joint_pairs[j].mate2_transcript.genome_start - && joint_pairs[i].mate2_transcript.is_reverse - == joint_pairs[j].mate2_transcript.is_reverse; - if same_pos - && joint_pairs[i].combined_wt_score < joint_pairs[j].combined_wt_score - && exons_subset( - &joint_pairs[i].mate2_transcript.exons, - &joint_pairs[j].mate2_transcript.exons, - ) - { - keep[i] = false; - break; - } - } - } - joint_pairs = joint_pairs - .into_iter() - .enumerate() - .filter(|(i, _)| keep[*i]) - .map(|(_, p)| p) - .collect(); - } - // Step 2: score-range filter (STAR's multMapSelect). if !joint_pairs.is_empty() { let best_score = joint_pairs From c94e6c3f71ccb5ca01cb6a65e934abbf2da989db Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 13:40:16 +0100 Subject: [PATCH 04/17] Drop the read when a window full of anchors receives another align STAR's assignAlignToWindow recomputes WALrec over the non-anchors of a full window before it looks at the new align; with none it sets MARKER_TOO_MANY_ANCHORS_PER_WINDOW and nW=0, so the read is 'other'. rustar-aligner skipped an incoming non-anchor instead, so the read could map (chr21 ERR188021.1789477) or count as 'too short' (yeast 1, chr21 17, fusion 73 reads). New integration test in tests/unmapped_buckets.rs; it fails on the previous code. Co-Authored-By: Claude Opus 5.5 --- src/align/stitch.rs | 18 +++++------ tests/unmapped_buckets.rs | 64 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 9 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index a846893..8846d2a 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -942,14 +942,20 @@ pub fn cluster_seeds( // Capacity check (seedPerWindowNmax) with anchor protection if window.alignments.len() >= seed_per_window_nmax { - // Find min length of non-anchor entries (STAR: recalculate WALrec) - let min_non_anchor_len = window + // Find min length of non-anchor entries (STAR: recalculate WALrec). + // With none, STAR drops the whole read before it looks at the new + // align, so a non-anchor arriving at a window full of anchors still + // triggers it: MARKER_TOO_MANY_ANCHORS_PER_WINDOW → nW=0. + let Some(min_non_anchor_len) = window .alignments .iter() .filter(|wa| !wa.is_anchor) .map(|wa| wa.length) .min() - .unwrap_or(usize::MAX); + else { + too_many_anchors = true; + break 'outer; + }; // Update persistent threshold window.wa_lrec = min_non_anchor_len; @@ -964,12 +970,6 @@ pub fn cluster_seeds( window .alignments .retain(|wa| wa.is_anchor || wa.length > min_non_anchor_len); - - if window.alignments.len() >= seed_per_window_nmax { - // STAR: MARKER_TOO_MANY_ANCHORS_PER_WINDOW → nW=0, abort entire read - too_many_anchors = true; - break 'outer; - } } // STAR's addition condition: if (aAnchor || aLength > WALrec[iW]) diff --git a/tests/unmapped_buckets.rs b/tests/unmapped_buckets.rs index 617219c..8ac7784 100644 --- a/tests/unmapped_buckets.rs +++ b/tests/unmapped_buckets.rs @@ -221,3 +221,67 @@ fn single_end_unmapped_reads_split_between_too_short_and_other() { "single-end `too short` bucket: {rows:?}" ); } + +#[test] +fn window_full_of_anchors_drops_the_read_as_other() { + // STAR's assignAlignToWindow: once a window holds seedPerWindowNmax (50) + // aligns, the next align makes it recompute the shortest non-anchor; when + // every align in the window is an anchor there is none, and STAR drops the + // whole read (MARKER_TOO_MANY_ANCHORS_PER_WINDOW, nW=0), so it is counted + // as "other". That happens whatever the incoming align is, a non-anchor + // included. Here the non-anchor is the poly-A tail each mate ends in: the + // 40 kb genome (one window) has 100 poly-A runs, so the tail maps more than + // winAnchorMultimapNmax times. STAR 2.7.11b reports both pairs as "other". + let dir = TempDir::new().unwrap(); + let root = dir.path(); + let mut genome = Vec::new(); + for k in 0..100u32 { + genome.extend_from_slice(&lcg_seq(88_888 + k, 400)); + genome.extend_from_slice(b"AAAAAAAAAAAA"); + } + let genome_dir = root.join("genome"); + build_index(&write_fasta(root, &genome), &genome_dir); + + let (p1, p2) = (root.join("r1.fq"), root.join("r2.fq")); + let mut f1 = fs::File::create(&p1).unwrap(); + let mut f2 = fs::File::create(&p2).unwrap(); + for i in [115u32, 1603] { + for (f, seed) in [(&mut f1, 700_000 + i), (&mut f2, 900_000 + i)] { + let mut seq = lcg_seq(seed, READ_LEN - 8); + seq.extend_from_slice(b"AAAAAAAA"); + writeln!(f, "@pair{i}").unwrap(); + f.write_all(&seq).unwrap(); + writeln!(f, "\n+\n{}", "I".repeat(READ_LEN)).unwrap(); + } + } + drop((f1, f2)); + + let prefix = format!("{}/anchors_", root.display()); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + p1.to_str().unwrap(), + p2.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let rows = parse_final_log(Path::new(&format!("{prefix}Log.final.out"))); + assert_eq!(get(&rows, "Number of input reads"), 2, "fixture size"); + assert_eq!( + get(&rows, "Number of reads unmapped: other"), + 2, + "a window full of anchors drops the read as `other`: {rows:?}" + ); + assert_eq!( + get(&rows, "Number of reads unmapped: too short"), + 0, + "{rows:?}" + ); +} From d8322e50b6cde77579c29a96eccfe21853478af6 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 13:40:16 +0100 Subject: [PATCH 05/17] STARsolo: STAR's EM for --soloMultiMappers, zeros written STAR's gEem starts from unique plus uniform counts, zeroes values below 0.01 at the start of each round, and stops when no gene moves by 0.01 or after 101 rounds, so it can keep a value below 0.01 (SoloFeature_collapseUMIall.cpp). The UniqueAndMult matrices also list every gene of the cell's unique or multi-gene counts, a zero included (SoloFeature_outputResults.cpp). Co-Authored-By: Claude Opus 5.5 --- src/solo/count.rs | 75 ++++++++++++++++++++++++++++++++--------------- 1 file changed, 52 insertions(+), 23 deletions(-) diff --git a/src/solo/count.rs b/src/solo/count.rs index 2acfb8a..2f8224d 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -876,38 +876,54 @@ fn distribute_multi( } } MultiMethod::Em => { - // theta_g = u_g + (multi distributed proportional to theta), iterated. - let mut theta = u.clone(); + // STAR's gEem (`SoloFeature_collapseUMIall.cpp`): start from the + // unique counts plus the uniform spread; each round zeroes the + // previous estimate's values below 0.01, then gives every molecule to + // its genes in proportion to that estimate. It stops once no gene + // moved by 0.01, or after 101 rounds, and keeps that round's estimate, + // which can still hold values below 0.01. + let genes: std::collections::BTreeSet = + molecules.iter().flatten().copied().collect(); + let mut old: HashMap = genes.iter().map(|&g| (g, 0.0)).collect(); for s in molecules { + let w = unit(s); for &g in s { - theta.entry(g).or_insert(0.0); + *old.get_mut(&g).unwrap() += w; } } - for _ in 0..100 { - let mut next = u.clone(); - for s in molecules { - for &g in s { - next.entry(g).or_insert(0.0); + for (&g, v) in &mut old { + *v += get(u, g); + } + let mut iter = 0u32; + let new = loop { + iter += 1; + let mut new: HashMap = genes.iter().map(|&g| (g, get(u, g))).collect(); + for v in old.values_mut() { + if *v < 0.01 { + *v = 0.0; } } for s in molecules { - let total: f64 = s.iter().map(|&g| get(&theta, g)).sum(); - for &g in s { - let w = if total > 0.0 { - get(&theta, g) / total - } else { - unit(s) - }; - *next.get_mut(&g).unwrap() += w; + let norm = 1.0 / s.iter().map(|g| old[g]).sum::(); + for g in s { + *new.get_mut(g).unwrap() += old[g] * norm; } } - let delta: f64 = next.iter().map(|(g, v)| (v - get(&theta, *g)).abs()).sum(); - theta = next; - if delta < 1e-6 { - break; + let max_change = genes + .iter() + .map(|g| (new[g] - old[g]).abs()) + .fold(0.0, f64::max); + if max_change < 0.01 || iter > 100 { + break new; } + old = new; + }; + // STAR keeps the multi part (gEem - gEu) and adds the unique count + // back when it writes the matrix. + for g in genes { + let ug = get(u, g); + out.insert(g, ug + (new[&g] - ug)); } - out = theta; } } out @@ -990,9 +1006,11 @@ fn build_multi_matrices( -> Result<(), Error> { let mols = cell_molecules(cb, mptr); for (k, &m) in methods.iter().enumerate() { - let counts = distribute_multi(m, u, &mols); + // Every gene of the cell's unique counts or of a multi-gene + // molecule is written, a zero included (STAR's outputResults + // walks both lists without testing the value). let mut entries: Vec<(u32, f64)> = - counts.into_iter().filter(|&(_, v)| v > 1e-9).collect(); + distribute_multi(m, u, &mols).into_iter().collect(); entries.sort_unstable_by_key(|&(g, _)| g); for (g, v) in entries { writeln!(writers[k], "{} {} {}", g + 1, cb + 1, fmt_real(v)) @@ -3003,6 +3021,17 @@ mod tests { assert!((pu0[&0] - 0.5).abs() < 1e-9 && (pu0[&1] - 0.5).abs() < 1e-9); } + /// STAR's EM zeroes values below 0.01 at the start of a round, not at the + /// end, and stops once nothing moved by 0.01: a gene with no unique count + /// halves each round from its uniform 0.5 and keeps 2^-7 when it stops. + #[test] + fn em_keeps_the_last_round_below_the_cutoff() { + let u: HashMap = [(1u32, 1.0)].into_iter().collect(); + let em = distribute_multi(MultiMethod::Em, &u, &[vec![0u32, 1u32]]); + assert_eq!(em[&0], 0.007_812_5); + assert_eq!(em[&1], 1.992_187_5); + } + #[test] fn called_cells_methods() { let mk = |cb, u| CellStat { From 824cc3f0d654267a590b6ab8e4116f418780e1c4 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 13:40:16 +0100 Subject: [PATCH 06/17] STARsolo: keep gene-ambiguous reads with a 1MM_multi barcode STAR records them with their whitelist candidates and resolves the barcode at counting time, as for a unique-gene read (SoloReadFeature_inputRecords.cpp); one it cannot resolve counts as noTooManyWLmatches. rustar-aligner dropped them, which changed GeneFull Features.stats and Summary.csv on the mouse chr19 fixture. Co-Authored-By: Claude Opus 5.5 --- src/solo/count.rs | 30 +++++++++++++++++++++++++----- src/solo/mod.rs | 39 ++++++++++++++++++++++++++------------- 2 files changed, 51 insertions(+), 18 deletions(-) diff --git a/src/solo/count.rs b/src/solo/count.rs index 2f8224d..d50c0a3 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1747,13 +1747,30 @@ pub fn write_gene_matrix( } // --soloMultiMappers: UniqueAndMult-.mtx alongside raw. - let (mut multi_counted, mut multi_dropped) = (0u64, 0u64); + let (mut multi_counted, mut multi_dropped, mut multi_posterior) = (0u64, 0u64, 0u64); if !multi_methods.is_empty() { let mut mg = recorder.multi_gene.lock().unwrap(); - let before = mg.len(); - mg.retain(|r| cb_passes_one_exact(one_exact, &exact_counts, r.cb)); + // 1MM_multi barcodes: the same posterior as unique-gene reads; a + // barcode it cannot pick is `noTooManyWLmatches`. oneExact applies + // to single corrections only. + let (pending, mut resolved): (Vec<_>, Vec<_>) = std::mem::take(&mut *mg) + .into_iter() + .partition(|r| !r.candidates.is_empty()); + let before = resolved.len(); + resolved.retain(|r| cb_passes_one_exact(one_exact, &exact_counts, r.cb)); + multi_dropped = (before - resolved.len()) as u64; + for mut r in pending { + match resolve_multi_cb(&r.candidates, &exact_counts, pseudocount) { + Some(cb) => { + r.cb = cb; + r.candidates.clear(); + resolved.push(r); + } + None => multi_posterior += 1, + } + } + *mg = resolved; multi_counted = mg.len() as u64; - multi_dropped = (before - mg.len()) as u64; build_multi_matrices( &body, &mg, @@ -1781,7 +1798,10 @@ pub fn write_gene_matrix( F::NO_MM_TO_WL_WITHOUT_EXACT, mstats.dropped_one_exact + multi_dropped, ); - st.add(F::NO_TOO_MANY_WL_MATCHES, mstats.dropped_posterior); + st.add( + F::NO_TOO_MANY_WL_MATCHES, + mstats.dropped_posterior + multi_posterior, + ); st.add(F::YES_WL_MATCH, mstats.unique_reads + multi_counted); st.add(F::YES_SUB_WL_MATCH_UNIQUE_FEATURE, mstats.unique_reads); st.add(F::YES_UMIS, mstats.cells.iter().map(|c| c.n_umis).sum()); diff --git a/src/solo/mod.rs b/src/solo/mod.rs index b91d1fb..d12a3d2 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -1208,6 +1208,10 @@ pub struct MultiGeneRecord { pub cb: u32, pub umi: u64, pub genes: Vec, + /// The whitelist candidates of a 1MM_multi barcode, resolved at counting + /// time as for a unique-gene read (`cb` is meaningless until then); empty + /// when the barcode was resolved per read. + pub candidates: Vec, } /// Thread-safe sink for the records produced during alignment. @@ -2085,19 +2089,28 @@ impl SoloContext { GeneAssignment::Gene(g) => g, GeneAssignment::Ambiguous => { // Gene-ambiguous read: record its gene set for - // --soloMultiMappers distribution (resolved CB only). - if let Some(cb) = cb_resolved { - let genes = match feature { - SoloFeature::Gene => &class.gene_multi, - SoloFeature::GeneFull => &class.gene_full_multi, - }; - if !genes.is_empty() { - fo.multi_gene = Some(MultiGeneRecord { - cb, - umi, - genes: genes.clone(), - }); - } + // --soloMultiMappers distribution. A 1MM_multi barcode + // goes along with its candidates: STAR records these + // reads too and resolves the barcode when it counts + // (`SoloReadFeature_inputRecords.cpp`). + let genes = match feature { + SoloFeature::Gene => &class.gene_multi, + SoloFeature::GeneFull => &class.gene_full_multi, + }; + let cb_and_candidates = match (cb_resolved, &cb_match) { + (Some(cb), _) => Some((cb, Vec::new())), + (None, CbMatch::Multi(cands)) => Some((0, cands.clone())), + (None, _) => None, + }; + if let Some((cb, candidates)) = cb_and_candidates + && !genes.is_empty() + { + fo.multi_gene = Some(MultiGeneRecord { + cb, + umi, + genes: genes.clone(), + candidates, + }); } // readInfo: the resolved cell and the raw UMI, whether // the read is left uncounted or goes to the multi-gene From 1fd99a320b50c9bf03e2ca5c4061c82d33399ed4 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 13:40:26 +0100 Subject: [PATCH 07/17] STARsolo: --soloUMIfiltering for the SJ feature STAR's collapseUMIall filters multi-gene UMIs on every feature. For SJ, a read crossing two junctions is one UMI on two features, a tie MultiGeneUMI_CR drops (4 UMIs on the mouse chr19 fixture). The gene path's filter-then-dedup is now shared by both. Co-Authored-By: Claude Opus 5.5 --- src/solo/count.rs | 109 +++++++++++++++++++++++++++++++++------------- 1 file changed, 79 insertions(+), 30 deletions(-) diff --git a/src/solo/count.rs b/src/solo/count.rs index d50c0a3..8b71d50 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -590,28 +590,7 @@ fn build_matrix_body( // `MultiGeneUMI_CR` decides gene ownership on *corrected* // UMIs, so it needs the correction to have happened first and // cannot go through the shared filter-then-dedup path below. - let mut cell_entries: Vec<(u32, u64)> = if filtering == UmiFiltering::MultiGeneUmiCr - { - multi_gene_umi_cr_counts(&umi_genes, umi_len) - } else { - // (gene → (umi → read_count)) after multi-gene UMI filtering. - let mut gene_umis: HashMap> = HashMap::default(); - for (&umi, genes) in &umi_genes { - for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { - *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; - } - } - - // Collapse UMIs per gene, then emit gene-ascending. - let mut entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); - for (&gene, umis) in &gene_umis { - let count = dedup_count(umis, method, umi_len); - if count > 0 { - entries.push((gene, count)); - } - } - entries - }; + let mut cell_entries = cell_feature_counts(&umi_genes, filtering, method, umi_len); cell_entries.sort_unstable_by_key(|&(g, _)| g); // STAR's readInfo: every read counted for this cell records the @@ -1944,6 +1923,7 @@ pub fn write_gene_matrix( let sj_counts = build_sj_matrix( &sj_records, &row, + filtering, method, umi_len, &sj_dir.join(&matrix_name), @@ -2080,6 +2060,41 @@ pub fn write_gene_matrix( Ok(()) } +/// One cell's `(feature, UMI count)` entries from its `umi → feature → reads` +/// map: STAR's multi-gene UMI filtering (`--soloUMIfiltering`), then UMI +/// collapsing per feature. STAR runs this on every feature, SJ included, where a +/// read crossing two junctions is one UMI on two features +/// (`SoloFeature_collapseUMIall.cpp`). +fn cell_feature_counts( + umi_genes: &HashMap>, + filtering: UmiFiltering, + method: UmiDedup, + umi_len: usize, +) -> Vec<(u32, u64)> { + // `MultiGeneUMI_CR` decides gene ownership on *corrected* UMIs, so it + // needs the correction to have happened first and cannot go through the + // shared filter-then-dedup path below. + if filtering == UmiFiltering::MultiGeneUmiCr { + return multi_gene_umi_cr_counts(umi_genes, umi_len); + } + // (gene → (umi → read_count)) after multi-gene UMI filtering. + let mut gene_umis: HashMap> = HashMap::default(); + for (&umi, genes) in umi_genes { + for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { + *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + } + } + // Collapse UMIs per gene. + let mut entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); + for (&gene, umis) in &gene_umis { + let count = dedup_count(umis, method, umi_len); + if count > 0 { + entries.push((gene, count)); + } + } + entries +} + /// Build the SJ feature matrix from (cell, UMI, junction) records, mapping each /// junction's absolute intron coords to its `SJ.out.tab` row and UMI-collapsing /// per (cell, junction). Junctions not in `row` (filtered out of SJ.out.tab) are @@ -2088,6 +2103,7 @@ pub fn write_gene_matrix( fn build_sj_matrix( records: &[crate::solo::SjCountRecord], row: &HashMap<(u64, u64), u32>, + filtering: UmiFiltering, method: UmiDedup, umi_len: usize, matrix_path: &Path, @@ -2112,22 +2128,18 @@ fn build_sj_matrix( let mut i = 0; while i < recs.len() { let cb = recs[i].cb; - // junction row → (umi → read count) for this cell. - let mut sj_umis: HashMap> = HashMap::default(); + // umi → (junction row → read count) for this cell. + let mut umi_rows: HashMap> = HashMap::default(); while i < recs.len() && recs[i].cb == cb { let r = recs[i]; if let Some(&rw) = row.get(&(r.intron_start, r.intron_end)) { - *sj_umis.entry(rw).or_default().entry(r.umi).or_insert(0) += 1; + *umi_rows.entry(r.umi).or_default().entry(rw).or_insert(0) += 1; counts.reads += 1; counts.exact += u64::from(r.exact); } i += 1; } - let mut entries: Vec<(u32, u64)> = sj_umis - .into_iter() - .map(|(rw, umis)| (rw, dedup_count(&umis, method, umi_len))) - .filter(|&(_, c)| c > 0) - .collect(); + let mut entries = cell_feature_counts(&umi_rows, filtering, method, umi_len); entries.sort_unstable_by_key(|&(rw, _)| rw); if !entries.is_empty() { counts.cells += 1; @@ -3052,6 +3064,43 @@ mod tests { assert_eq!(em[&1], 1.992_187_5); } + /// STAR filters multi-gene UMIs on every feature, SJ included: a read that + /// crosses two junctions is one UMI on two features, a tie that + /// `MultiGeneUMI_CR` drops. + #[test] + fn sj_matrix_applies_multi_gene_umi_filtering() { + let rec = |intron_start| crate::solo::SjCountRecord { + cb: 3, + umi: 42, + intron_start, + intron_end: intron_start + 100, + exact: true, + }; + let records = vec![rec(1000), rec(2000)]; + let row: HashMap<(u64, u64), u32> = + [((1000, 1100), 0), ((2000, 2100), 1)].into_iter().collect(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("matrix.mtx"); + let run = |filtering| { + build_sj_matrix( + &records, + &row, + filtering, + UmiDedup::OneMmCr, + 10, + &path, + 2, + 8, + false, + ) + .unwrap() + }; + let plain = run(UmiFiltering::None); + assert_eq!((plain.umis, plain.cells, plain.reads), (2, 1, 2)); + let cr = run(UmiFiltering::MultiGeneUmiCr); + assert_eq!((cr.umis, cr.cells, cr.reads), (0, 0, 2)); + } + #[test] fn called_cells_methods() { let mk = |cb, u| CellStat { From 0144a96ef529d24a8b3618a2b72ca33e40569f06 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 13:40:26 +0100 Subject: [PATCH 08/17] Docs: the DIVERGENCE.md 5 residuals are fixed DIVERGENCE.md 5, the README comparison tables and residuals, and CHANGELOG. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 35 +++++++++++++++++++++++++++++++++++ DIVERGENCE.md | 7 +------ README.md | 21 ++++++++------------- 3 files changed, 44 insertions(+), 19 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e..b18a76e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,41 @@ Sections commonly used: Features, Bug fixes, Other changes. ## [Unreleased] +### Bug fixes + +- The residual differences in DIVERGENCE.md §5 are fixed. On the parity + suite's yeast (10k and 100k, SE and PE, chimeric, `--outSAMunmapped + Within`), human chr21 (SE, PE) and fusion runs, alignments, `SJ.out.tab` + and `Log.final.out` now match STAR 2.7.11b, mate order included; on the + mouse chr19 STARsolo runs every `Solo.out` file does. + - `nM` is STAR's mismatch count as accumulated while stitching, not a + recount of the final alignment. The two differ when a junction scan moves + a junction back over bases already counted as mismatches against the + intron: STAR counts them again and keeps both counts. + - The two mates of a pair are written in STAR's order: the one on the + forward strand first, so mate 2 leads when mate 1 is reverse. Every + paired alignment was mate 1 first before, about half of STAR's. + - A pair whose two alignments share their mate positions but place an + indel or junction differently inside one mate is no longer collapsed to + one alignment. The extra mate-2 dedup that did it is gone; STAR's window + dedup compares the whole pair. + - A read is dropped as "other" when a window full of anchors receives one + more align, whatever that align is (STAR's + `MARKER_TOO_MANY_ANCHORS_PER_WINDOW`). A non-anchor arriving there was + skipped instead, so the read could map or count as "too short". + - STARsolo `--soloMultiMappers EM` is STAR's EM: it starts from the unique + counts plus the uniform spread, zeroes values below 0.01 at the start of + each round, and stops once no gene moves by 0.01 (or after 101 rounds), + so a gene can keep a value below 0.01. Every gene of a cell's unique or + multi-gene counts is written to `UniqueAndMult-*.mtx`, a zero included. + - A gene-ambiguous read whose barcode has several 1MM whitelist matches is + kept for `--soloMultiMappers` and its barcode resolved at counting time, + as for a unique-gene read; one it cannot resolve counts as + `noTooManyWLmatches`. It was dropped before. + - `--soloUMIfiltering` applies to the STARsolo `SJ` feature as it does to + genes: a read that crosses two junctions is one UMI on two features, a + tie `MultiGeneUMI_CR` drops. + ### Other changes - Third-party notices completed. `LICENSES/` gains the licences of code diff --git a/DIVERGENCE.md b/DIVERGENCE.md index f0890fb..a693341 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -321,12 +321,7 @@ rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` These are **not** deliberate divergences. They are the differences `test/parity_suite.sh` finds, kept here so that each one is fixed or explained before v0.2.0 (README, *Comparison with STAR 2.7.11b*). Read names are from that suite's datasets. -- **`nM` and mate order on multi-mapping pairs.** Human chr21 `ERR188021.1648513` and `ERR188021.1662773`, plus 21 pairs of the fusion set (e.g. `fustut00000271e5`, `fustut00000282b2`). The loci, CIGARs, `AS` and MAPQ match STAR's. On some alignments `nM` is 1–2 lower than STAR's, and in 11 of the 21 fusion pairs the mates are written in the opposite order. -- **Same-span alignments collapsed.** Fusion set `fustut000009a03e`, `fustut000009a03f`, `fustut0000167ff6`. STAR keeps two alignments with the same mate positions but different indel or splice placement inside one mate; rustar-aligner keeps one, so `NH` is one lower. -- **A pair STAR leaves unmapped.** Human chr21 `ERR188021.1789477`: rustar-aligner maps it to 10 loci (`AS` 136, spliced), adding 2 lines to `SJ.out.tab`. -- **Unmapped-read category.** Some unmapped reads are counted as "too short" in `Log.final.out` (and `uT:1`), where STAR has "other" (`uT:0`): yeast 100k PE 1 (`ERR12389696.13360341`), human chr21 PE 17, fusion 73. For the yeast pair, STAR builds no alignment window, while rustar-aligner builds one, aligns 31 points and rejects it as too short. -- **STARsolo `GeneFull` multi-mappers.** With `--soloMultiMappers`, STAR counts one multi-gene read whose barcode needed 1MM_multi correction, and rustar-aligner drops it (`UniqueAndMult-*.mtx`, one `Features.stats` line). -- **STARsolo `SJ` with `--soloUMIdedup 1MM_CR`.** STAR collapses 4 more SJ UMIs than rustar-aligner on the mouse chr19 fixture (SJ `matrix.mtx`, `Features.stats`, `Summary.csv`). +There are none at present. The six once listed here (`nM` and mate order, same-span alignments, a pair STAR leaves unmapped, the unmapped-read category, STARsolo `GeneFull` multi-mappers and `SJ` UMI filtering) now match STAR on the suite's yeast, chr21, fusion and STARsolo runs. --- diff --git a/README.md b/README.md index 1b2f112..9aace6e 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ A Rust reimplementation of [STAR](https://github.com/alexdobin/STAR) (Spliced Tr rustar-aligner aims to be a faithful port of STAR, matching the original behavior as closely as possible. It uses the same genome index format, accepts the same `--camelCase` command-line parameters, and produces compatible SAM/BAM output. -**Current status (v0.1, working toward v0.2.0):** single-end and paired-end alignment, splice junctions, two-pass mode, chimeric detection, `--quantMode GeneCounts` / `TranscriptomeSAM`, STARsolo (Gene / GeneFull / SJ / Velocyto, barcode correction, UMI dedup, multi-mapper resolution, cell calling), WASP filtering, mate-overlap merging, coverage tracks, adapter clipping, sparse suffix arrays and multi-threading. On the [parity suite](#comparison-with-star-2711b) the genome indices are byte-identical to STAR's, and alignments, `SJ.out.tab`, gene counts, chimeric junctions and STARsolo output match STAR 2.7.11b record for record, apart from the [residual differences](#residual-differences) listed below. Every deliberate difference is in [DIVERGENCE.md](DIVERGENCE.md); the release plan is in [ROADMAP.md](ROADMAP.md#release-plan). +**Current status (v0.1, working toward v0.2.0):** single-end and paired-end alignment, splice junctions, two-pass mode, chimeric detection, `--quantMode GeneCounts` / `TranscriptomeSAM`, STARsolo (Gene / GeneFull / SJ / Velocyto, barcode correction, UMI dedup, multi-mapper resolution, cell calling), WASP filtering, mate-overlap merging, coverage tracks, adapter clipping, sparse suffix arrays and multi-threading. On the [parity suite](#comparison-with-star-2711b) the genome indices are byte-identical to STAR's, and alignments, `SJ.out.tab`, gene counts, chimeric junctions and STARsolo output match STAR 2.7.11b record for record; [residual differences](#residual-differences) lists any that remain. Every deliberate difference is in [DIVERGENCE.md](DIVERGENCE.md); the release plan is in [ROADMAP.md](ROADMAP.md#release-plan). ## Quick Start @@ -139,17 +139,17 @@ Alignment records are compared per read: with more than one thread STAR writes r | Yeast 100k | PE, `--quantMode TranscriptomeSAM GeneCounts`, sorted BAM | 188,258 / 188,258 | 0 | identical | | Yeast 100k | PE, `--twopassMode Basic` | 188,262 / 188,262 | 0 | identical | | Yeast 100k | PE, `--outFilterType BySJout` | 187,646 / 187,646 | 0 | identical | -| Yeast 100k | PE, `--outSAMattributes All --outSAMstrandField intronMotif --outSAMunmapped Within` | 218,994 / 218,996 | 1 (unmapped pair, `uT`) | identical | +| Yeast 100k | PE, `--outSAMattributes All --outSAMstrandField intronMotif --outSAMunmapped Within` | 218,996 / 218,996 | 0 | identical | | Yeast 100k | PE, `--genomeSAsparseD 2` index | 188,436 / 188,436 | 0 | identical | | Yeast 100k | SE / PE, chimeric (`--chimSegmentMin 12`, `--chimMultimapNmax` 0 and 10) | all | 0 | identical | | Human chr21 | SE (mate 1) | 541,935 / 541,935 | 0 | identical | -| Human chr21 | PE | 902,186 / 902,190 | 3 | 3,071 / 3,073 lines shared | -| Fusion | PE, chimeric (`--chimMultimapNmax` 0 and 10) | 80,592 / 80,656 | 24 | identical | +| Human chr21 | PE | 902,190 / 902,190 | 0 | identical | +| Fusion | PE, chimeric (`--chimMultimapNmax` 0 and 10) | 80,656 / 80,656 | 0 | identical | - **Gene counts** (`ReadsPerGene.out.tab`): identical on yeast. On chr21 PE, 274 of 276 lines are identical; the other 2 come from the reads that differ. - **Transcriptome BAM** (`--quantMode TranscriptomeSAM`): every record is identical on yeast (138,084) and chr21 (28,972) except for *which* alignment of a multi-mapping read carries the primary flag. STAR draws that from a per-thread random stream, so its choice depends on how reads are split across threads. rustar-aligner reproduces STAR's `--runThreadN 1` choice at any thread count, and against STAR run single-threaded the yeast transcriptome BAM is identical record for record. See [DIVERGENCE.md §1.1](DIVERGENCE.md#11-multimapper-tie-breaking--rng). - **Chimeric junctions** (`Chimeric.out.junction`): the same lines as STAR in every configuration: yeast SE 194, PE 717 / 882 (`--chimMultimapNmax` 0 / 10); chr21 PE 301 / 3,643; fusion 31 / 675. Only the order differs, for the threading reason above. -- **`Log.final.out`**: identical on yeast SE and 10k PE. Elsewhere the mapped counts differ only by the reads in the table, and a few unmapped reads are counted as "too short" rather than "other" (yeast 100k PE 1, chr21 PE 17, fusion 73). +- **`Log.final.out`**: identical on yeast (10k and 100k, SE and PE, chimeric, `--outSAMunmapped Within`), chr21 SE and PE, and fusion. ### Single-cell (STARsolo) @@ -158,20 +158,15 @@ Mouse chr19, 10x v2, `--soloType CB_UMI_Simple --soloFeatures Gene GeneFull SJ V | Options | BAM records identical | `Solo.out` files identical | |---------|-----------------------|----------------------------| | defaults, sorted BAM with `CR UR CB UB GX GN sS sQ sM` | 19,955 / 19,955 | 36 / 36 | -| `--soloCellFilter EmptyDrops_CR --soloMultiMappers Uniform EM --soloUMIdedup 1MM_CR --soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts --soloUMIfiltering MultiGeneUMI_CR` | (no BAM) | 34 / 40 | +| `--soloCellFilter EmptyDrops_CR --soloMultiMappers Uniform EM --soloUMIdedup 1MM_CR --soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts --soloUMIfiltering MultiGeneUMI_CR` | (no BAM) | 40 / 40 | -In the default run every file matches: raw and filtered matrices, `barcodes.tsv` / `features.tsv`, `Summary.csv`, `Features.stats`, `Barcodes.stats` and `UMIperCellSorted.txt`. In the second run the EmptyDrops_CR cell calls match too, because rustar-aligner draws from the same random-number generator as STAR. The 6 files that differ are explained under [residual differences](#residual-differences). +In the default run every file matches: raw and filtered matrices, `barcodes.tsv` / `features.tsv`, `Summary.csv`, `Features.stats`, `Barcodes.stats` and `UMIperCellSorted.txt`. In the second run the EmptyDrops_CR cell calls match too, because rustar-aligner draws from the same random-number generator as STAR. ### Residual differences These aren't deliberate. Each is tracked to be fixed or explained before v0.2.0; [DIVERGENCE.md §5](DIVERGENCE.md#5-known-residual-single-read-differences) has the read names. -- **`nM` and mate order on some multi-mapping pairs** (chr21 2 pairs, fusion 21 pairs). Same loci, CIGARs, `AS` and MAPQ. On some alignments `nM` is 1–2 lower than STAR's, and in about half of these pairs the two mates are written in the opposite order. -- **Same-span alignments collapsed** (fusion, 3 pairs). STAR reports two alignments with the same mate positions but different indel or splice placement inside a mate; rustar-aligner reports one, so `NH` is one lower. -- **A pair STAR leaves unmapped** (chr21, 1 pair). rustar-aligner maps it to 10 loci (`AS` 136, spliced). This adds 2 lines to `SJ.out.tab`. -- **Unmapped-read category** ("too short" versus "other", `uT`), for the counts listed above. The yeast pair: STAR finds no alignment window (`uT:0`); rustar-aligner finds a 31-point alignment and rejects it as too short (`uT:1`). -- **STARsolo `GeneFull` multi-mappers** (`--soloMultiMappers`). STAR counts one multi-gene read whose barcode needed 1MM_multi correction, and rustar-aligner drops it. This changes `UniqueAndMult-*.mtx` and one `Features.stats` line. -- **STARsolo `SJ` with `--soloUMIdedup 1MM_CR`**. STAR collapses 4 more SJ UMIs than rustar-aligner does (SJ `matrix.mtx`, `Features.stats`, `Summary.csv`). +None at present. ### Speed and memory From 7e2212c96dbabd5d612925f4257b1d1a740b9545 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 17:45:49 +0100 Subject: [PATCH 09/17] Splice-junction database: FxHashMap instead of SipHash SpliceJunctionDb is only looked up (is_annotated), never iterated, so the hasher has no effect on output. It is consulted on every junction stitch. Co-Authored-By: Claude Opus 5.5 --- src/junction/mod.rs | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/junction/mod.rs b/src/junction/mod.rs index 4755759..9df3bb4 100644 --- a/src/junction/mod.rs +++ b/src/junction/mod.rs @@ -13,7 +13,7 @@ pub mod sjdb_insert; pub use sj_output::SpliceJunctionStats; pub(crate) use sj_output::{SjKey, encode_motif}; -use std::collections::HashMap; +use rustc_hash::FxHashMap; /// Key for junction lookup: (chr_idx, intron_start, intron_end, strand). /// @@ -40,14 +40,14 @@ pub struct JunctionInfo { #[derive(Clone)] pub struct SpliceJunctionDb { /// Map: (chr_idx, intron_start, intron_end, strand) → annotated - junctions: HashMap, + junctions: FxHashMap, } impl SpliceJunctionDb { /// Create empty database (for no-GTF mode) pub fn empty() -> Self { Self { - junctions: HashMap::new(), + junctions: FxHashMap::default(), } } @@ -57,7 +57,8 @@ impl SpliceJunctionDb { /// `TranscriptomeIndex` and the `sjdb_insert` pipeline without /// re-parsing the file. pub fn from_raw_junctions(raw: &[(usize, u64, u64, u8)]) -> Self { - let mut junctions = HashMap::with_capacity(raw.len()); + let mut junctions = + FxHashMap::with_capacity_and_hasher(raw.len(), rustc_hash::FxBuildHasher); for &(chr_idx, intron_start, intron_end, strand) in raw { let key = JunctionKey { chr_idx, From 73f7afe45f7e663927a1bd5ae197cdb0e8643991 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 17:45:49 +0100 Subject: [PATCH 10/17] Stitcher: store a recorded transcript's window score and mapped length Every leaf of the recursion compared itself with each transcript already recorded in the window, recomputing their final score (a log2) and mapped length each time. They are now computed once, when the transcript is recorded. Co-Authored-By: Claude Opus 5.5 --- src/align/stitch.rs | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 8846d2a..996c4f6 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -1077,6 +1077,11 @@ pub(crate) struct WorkingTranscript { pub(crate) read_end: usize, pub(crate) genome_start: u64, pub(crate) genome_end: u64, + /// Final score ([`wt_final_score`]) and mapped length ([`wt_mapped_length`]), + /// stored when the transcript is recorded in its window: every later + /// transcript of the window is compared against them. + pub(crate) window_score: i32, + pub(crate) mapped_len: u32, } impl WorkingTranscript { @@ -1096,6 +1101,8 @@ impl WorkingTranscript { read_end: 0, genome_start: 0, genome_end: 0, + window_score: 0, + mapped_len: 0, } } } @@ -2721,7 +2728,7 @@ fn stitch_recurse( _ => None, }; // An empty window compares against STAR's `trInit`, score 0. - let window_best = transcripts.first().map_or(0, |t| wt_final_score(t, scorer)); + let window_best = transcripts.first().map_or(0, |t| t.window_score); if !window_score_gate(score, frag, window_best, max_score_mate, scorer) { return; } @@ -2739,8 +2746,8 @@ fn stitch_recurse( let existing = &transcripts[idx]; let overlap = blocks_overlap(&wt.exons, &existing.exons); let u_new = wt_len - overlap; - let u_old = wt_mapped_length(existing) - overlap; - if u_new == 0 && wt_score < wt_final_score(existing, scorer) { + let u_old = existing.mapped_len - overlap; + if u_new == 0 && wt_score < existing.window_score { dominated = true; break; } else if u_old == 0 { @@ -2764,10 +2771,13 @@ fn stitch_recurse( let pos = transcripts .iter() .position(|t| { - let score = wt_final_score(t, scorer); + let score = t.window_score; wt_score > score || (wt_score == score && wt_span < span(t)) }) .unwrap_or(transcripts.len()); + let mut wt = wt; + wt.window_score = wt_score; + wt.mapped_len = wt_len; transcripts.insert(pos, wt); transcripts.truncate(max_transcripts); } @@ -3077,6 +3087,8 @@ pub(crate) fn split_combined_wt( read_end: m1_read_end, genome_start: m1_genome_start, genome_end: m1_genome_end, + window_score: 0, + mapped_len: 0, }, WorkingTranscript { exons: m2_exons, @@ -3093,6 +3105,8 @@ pub(crate) fn split_combined_wt( read_end: m2_read_end, genome_start: m2_genome_start, genome_end: m2_genome_end, + window_score: 0, + mapped_len: 0, }, )) } From 457c34c7848f18fa2cd2dfa1467c981cfa7d9759 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 17:45:49 +0100 Subject: [PATCH 11/17] Junction scan: stop once no later position can win score1 changes by at most +1 per position, so when score1 plus the positions left plus the best motif score cannot beat the best score so far, the remaining positions cannot be chosen. Ties count for the reverse-strand deletion case only, as in the scan itself. The chosen position, motif and repeat lengths are unchanged. Co-Authored-By: Claude Opus 5.5 --- src/align/score.rs | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/align/score.rs b/src/align/score.rs index 899b6a4..11f810c 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -488,6 +488,11 @@ impl AlignmentScorer { // runs is decided once rather than per iteration. let motif_in_range = del >= self.align_intron_min as i64 && del <= self.align_intron_max as i64; + let max_motif = self + .score_gap + .max(self.score_gap_gcag) + .max(self.score_gap_atac) + .max(self.score_gap_noncan); // The four motif bases sit at `donor`, `donor+1`, `donor+del-2` and // `donor+del-1`, and `donor` moves exactly one base per iteration, so @@ -556,6 +561,16 @@ impl AlignmentScorer { } } + // `score1` rises by at most 1 a position, so once even that cannot + // lift a later position past the best one the rest of the scan is + // settled. Ties win only in the reverse-strand deletion case. + let scan_end = r_gap as i32 + next_seed_len as i32; + let bound = score1 + (scan_end - 1 - jr1) + if motif_in_range { max_motif } else { 0 }; + let tie_wins = !motif_in_range && is_reverse; + if bound < max_score2 || (bound == max_score2 && !tie_wins) { + break; + } + jr1 += 1; // STAR: jR1 < int(rBend) - int(rAend) where rBend = rBstart + L - 1 // This equals jR1 < rGap + L, allowing the scan to extend into the B seed. From 45aa4c9a75f1411bb4bdb65a0eced12a1784e9eb Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 17:45:50 +0100 Subject: [PATCH 12/17] Stitcher: memoize flank extensions in large windows Within one window the read, strand and mismatch limits are fixed, so an extension is determined by its start, direction, length cap and the transcript's mismatches and length so far; the recursion's leaves repeat many of them. Windows with fewer than 25 seeds skip the memo, where hashing cost more than it saved (+1.7% instructions on chr21 without the threshold). Co-Authored-By: Claude Opus 5.5 --- src/align/score.rs | 9 +++++ src/align/stitch.rs | 95 ++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 98 insertions(+), 6 deletions(-) diff --git a/src/align/score.rs b/src/align/score.rs index 11f810c..3fc291f 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -813,8 +813,17 @@ struct JunctionScanKey { #[derive(Default)] pub struct JunctionScanCache { map: rustc_hash::FxHashMap, + /// Flank extensions of the same window, for `stitch::extend_alignment_cached`. + pub(crate) ext: rustc_hash::FxHashMap, + /// Whether `ext` is used for this window. + pub(crate) memo_extensions: bool, } +/// What a flank extension depends on within one window: where it starts, which +/// way and how far it may go, the transcript's mismatches and length so far, +/// and whether it must reach the read end. +pub(crate) type ExtendKey = (usize, u64, i32, usize, u32, usize, bool); + impl JunctionScanCache { pub fn new() -> Self { Self::default() diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 996c4f6..f574665 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -137,6 +137,81 @@ fn score_region( /// staging copy is amortized, small enough to sit on the stack. const SCORE_REGION_CHUNK: usize = 256; +/// Windows with fewer seeds than this have few recursion leaves, and hashing +/// the extension costs more there than it saves. +const MEMO_EXTENSIONS_MIN_SEEDS: usize = 25; + +/// [`extend_alignment`] memoized in the window's cache. Within one window the +/// read, genome strand and mismatch limits are fixed, so the [`ExtendKey`] +/// fields decide the result; the recursion's leaves repeat many of them. +#[allow(clippy::too_many_arguments)] +fn extend_alignment_cached( + jcache: &mut crate::align::score::JunctionScanCache, + read_seq: &[u8], + read_start: usize, + genome_start: u64, + direction: i32, + max_extend: usize, + n_mm_prev: u32, + len_prev: usize, + n_mm_max: u32, + p_mm_max: f64, + index: &GenomeIndex, + is_reverse: bool, + extend_to_end: bool, +) -> ExtendResult { + let key = ( + read_start, + genome_start, + direction, + max_extend, + n_mm_prev, + len_prev, + extend_to_end, + ); + if !jcache.memo_extensions { + return extend_alignment( + read_seq, + read_start, + genome_start, + direction, + max_extend, + n_mm_prev, + len_prev, + n_mm_max, + p_mm_max, + index, + is_reverse, + extend_to_end, + ); + } + if let Some(&(extend_len, max_score, n_mismatch)) = jcache.ext.get(&key) { + return ExtendResult { + extend_len, + max_score, + n_mismatch, + }; + } + let r = extend_alignment( + read_seq, + read_start, + genome_start, + direction, + max_extend, + n_mm_prev, + len_prev, + n_mm_max, + p_mm_max, + index, + is_reverse, + extend_to_end, + ); + jcache + .ext + .insert(key, (r.extend_len, r.max_score, r.n_mismatch)); + r +} + /// Result of extending an alignment into flanking regions #[derive(Debug, Clone)] struct ExtendResult { @@ -1188,7 +1263,8 @@ fn stitch_align_to_transcript( // extendToEnd=false (local mode); spacer at position 150 naturally stops the extension. // nMatch approximated as score + n_mismatch (exact for seeds-only paths, approx otherwise). let n_match_m1 = (wt.score + wt.n_mismatch as i32).max(0) as usize; - let right_ext = extend_alignment( + let right_ext = extend_alignment_cached( + jcache, read_seq, last_exon.read_end, last_exon.genome_end, @@ -1248,7 +1324,8 @@ fn stitch_align_to_transcript( } }; let n_match_for_left = n_match_m1 + right_ext.extend_len + wa.length; - let left_ext = extend_alignment( + let left_ext = extend_alignment_cached( + jcache, read_seq, wa.read_pos, wa.sa_pos, @@ -2566,7 +2643,8 @@ fn stitch_recurse( let first_ext = if do_left_first { // Left extension (5' of forward read) if wt.read_start > 0 { - extend_alignment( + extend_alignment_cached( + jcache, read_seq, wt.read_start, wt.genome_start, @@ -2586,7 +2664,8 @@ fn stitch_recurse( } else { // Right extension (5' of reverse read) if wt.read_end < read_seq.len() { - extend_alignment( + extend_alignment_cached( + jcache, read_seq, wt.read_end, wt.genome_end, @@ -2612,7 +2691,8 @@ fn stitch_recurse( // Right extension (3' of forward read) let read_end_after = wt.read_end; // right boundary unchanged by left ext if read_end_after < read_seq.len() { - extend_alignment( + extend_alignment_cached( + jcache, read_seq, read_end_after, wt.genome_end, @@ -2633,7 +2713,8 @@ fn stitch_recurse( // Left extension (3' of reverse read) let read_start_after = wt.read_start; // left boundary unchanged by right ext if read_start_after > 0 { - extend_alignment( + extend_alignment_cached( + jcache, read_seq, read_start_after, wt.genome_start, @@ -3427,6 +3508,8 @@ pub(crate) fn stitch_seeds_core( // fixed for the whole recursion below, which is what makes the six-field // key in `JunctionScanCache` a complete identifier for a scan. let mut jcache = crate::align::score::JunctionScanCache::new(); + // Memoizing extensions pays only where the recursion has many leaves. + jcache.memo_extensions = wa_entries.len() >= MEMO_EXTENSIONS_MIN_SEEDS; stitch_recurse( 0, From c44c54b5ffb8a6ef70d92fdbfb0c2454cba0f533 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 17:45:50 +0100 Subject: [PATCH 13/17] Transcript dedup: compare CIGARs without building strings The exact-duplicate filters (SE transcripts, PE pairs) formatted every CIGAR into a String and hashed it with SipHash. They now key on the positions with FxHashMap and compare the CIGAR ops of transcripts at the same position; the first occurrence is kept, as before. Co-Authored-By: Claude Opus 5.5 --- src/align/read_align.rs | 72 ++++++++++++++++++++++++++++------------- 1 file changed, 50 insertions(+), 22 deletions(-) diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 312d0f7..35e450a 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -493,16 +493,26 @@ fn align_read_inner( // to the earliest-discovered alignment, so we must not reorder here (the // old code sorted by coordinate, destroying that order). { - let mut seen = std::collections::HashSet::new(); - transcripts.retain(|t| { - seen.insert(( - t.chr_idx, - t.genome_start, - t.genome_end, - t.is_reverse, - t.cigar_string(), - )) - }); + // Keyed on the position, with the CIGARs of one position compared + // directly: no CIGAR string is built per transcript. + let mut firsts: rustc_hash::FxHashMap<(usize, u64, u64, bool), Vec> = + rustc_hash::FxHashMap::default(); + let keep: Vec = transcripts + .iter() + .enumerate() + .map(|(i, t)| { + let same = firsts + .entry((t.chr_idx, t.genome_start, t.genome_end, t.is_reverse)) + .or_default(); + let dup = same.iter().any(|&j| transcripts[j].cigar == t.cigar); + if !dup { + same.push(i); + } + !dup + }) + .collect(); + let mut keep = keep.into_iter(); + transcripts.retain(|_| keep.next().unwrap_or(true)); } // STAR's `trBest` (ReadAlign_stitchPieces.cpp), over the transcripts in window @@ -1137,18 +1147,36 @@ pub fn align_paired_read( // earliest-discovered pair, so discovery order has to survive to the // primary-selection sort below. { - let mut seen = std::collections::HashSet::new(); - joint_pairs.retain(|p| { - seen.insert(( - p.mate1_transcript.chr_idx, - p.mate1_transcript.genome_start, - p.mate1_transcript.is_reverse, - p.mate1_transcript.cigar_string(), - p.mate2_transcript.genome_start, - p.mate2_transcript.is_reverse, - p.mate2_transcript.cigar_string(), - )) - }); + // Keyed on the mate positions, with the CIGARs of one key compared + // directly: no CIGAR string is built per pair. + let mut firsts: rustc_hash::FxHashMap<(usize, u64, bool, u64, bool), Vec> = + rustc_hash::FxHashMap::default(); + let keep: Vec = joint_pairs + .iter() + .enumerate() + .map(|(i, p)| { + let (m1, m2) = (&p.mate1_transcript, &p.mate2_transcript); + let same = firsts + .entry(( + m1.chr_idx, + m1.genome_start, + m1.is_reverse, + m2.genome_start, + m2.is_reverse, + )) + .or_default(); + let dup = same.iter().any(|&j| { + let q = &joint_pairs[j]; + q.mate1_transcript.cigar == m1.cigar && q.mate2_transcript.cigar == m2.cigar + }); + if !dup { + same.push(i); + } + !dup + }) + .collect(); + let mut keep = keep.into_iter(); + joint_pairs.retain(|_| keep.next().unwrap_or(true)); } // Step 2: score-range filter (STAR's multMapSelect). From a29228ff193944d8bd3abc5bff30b7aa7d38d9b1 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 18:41:31 +0100 Subject: [PATCH 14/17] Stitcher: backtrack on one transcript instead of cloning at each branch stitch_recurse used to clone the working transcript at every include branch (and stitch_align_to_transcript cloned it again), so the exclude branch could keep the original. It now stitches in place on a single transcript and undoes the change afterwards: a stitch only pushes exons and junctions, changes the last exon and updates scalar fields, and a leaf's extensions also change the first exon, so a small WtMark restores it. A leaf takes its mark only when an extension applies, and a recorded transcript is cloned once, when it enters the window list. stitch_align_to_transcript remains, for tests, as a copy-then-stitch wrapper; ExonBlock is now Copy. Output is unchanged: alignments, SJ.out.tab, Chimeric.out.junction, Log.final.out and Solo.out are identical on the yeast, chr21, fusion and STARsolo runs, and the full fusion set still matches STAR 2.7.11b. Single-threaded instructions: fusion 30k pairs -22%, yeast 100k PE -2.7%, chr21 100k SE unchanged (+0.4%, within noise). Co-Authored-By: Claude Opus 5.5 --- src/align/stitch.rs | 633 ++++++++++++++++++++++++++------------------ 1 file changed, 377 insertions(+), 256 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index f574665..6147a1c 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -1118,7 +1118,7 @@ pub fn cluster_seeds( } /// Lightweight exon block for in-progress transcript during recursion -#[derive(Debug, Clone)] +#[derive(Debug, Clone, Copy)] pub(crate) struct ExonBlock { pub(crate) read_start: usize, // 0-based inclusive pub(crate) read_end: usize, // 0-based exclusive @@ -1182,12 +1182,86 @@ impl WorkingTranscript { } } +/// What [`WorkingTranscript::restore`] needs to undo a stitch or a leaf's +/// extensions: they only push exons and junctions (the four junction vectors +/// together), change the last exon (and, for a leaf, the first) in place, and +/// update the scalar fields. +struct WtMark { + n_exons: usize, + first: Option, + last: Option, + n_junctions: usize, + score: i32, + n_mismatch: u32, + n_gap: u32, + n_junction: u32, + n_anchor: u32, + read_start: usize, + read_end: usize, + genome_start: u64, + genome_end: u64, +} + +impl WorkingTranscript { + /// Mark before a stitch, which leaves every exon but the last alone. + fn mark(&self) -> WtMark { + WtMark { + n_exons: self.exons.len(), + first: None, + last: self.exons.last().copied(), + n_junctions: self.junction_motifs.len(), + score: self.score, + n_mismatch: self.n_mismatch, + n_gap: self.n_gap, + n_junction: self.n_junction, + n_anchor: self.n_anchor, + read_start: self.read_start, + read_end: self.read_end, + genome_start: self.genome_start, + genome_end: self.genome_end, + } + } + + /// Mark before a leaf's extensions, which change the first exon as well. + fn mark_ends(&self) -> WtMark { + WtMark { + first: self.exons.first().copied(), + ..self.mark() + } + } + + fn restore(&mut self, m: &WtMark) { + self.exons.truncate(m.n_exons); + if let Some(first) = m.first { + self.exons[0] = first; + } + if let Some(last) = m.last { + self.exons[m.n_exons - 1] = last; + } + if self.junction_motifs.len() != m.n_junctions { + self.junction_motifs.truncate(m.n_junctions); + self.junction_annotated.truncate(m.n_junctions); + self.junction_strand.truncate(m.n_junctions); + self.junction_shifts.truncate(m.n_junctions); + } + self.score = m.score; + self.n_mismatch = m.n_mismatch; + self.n_gap = m.n_gap; + self.n_junction = m.n_junction; + self.n_anchor = m.n_anchor; + self.read_start = m.read_start; + self.read_end = m.read_end; + self.genome_start = m.genome_start; + self.genome_end = m.genome_end; + } +} + /// Fill the gap between the last exon in a WorkingTranscript and the next WA entry. /// Returns None if stitching fails (mismatch limit, overhang too short, etc.). /// Matches STAR's stitchAlignToTranscript.cpp logic. #[allow(clippy::too_many_arguments)] -fn stitch_align_to_transcript( - wt: &WorkingTranscript, +fn stitch_align_in_place( + new_wt: &mut WorkingTranscript, wa: &WindowAlignment, read_seq: &[u8], index: &GenomeIndex, @@ -1197,8 +1271,10 @@ fn stitch_align_to_transcript( align_mates_gap_max: u64, jcache: &mut crate::align::score::JunctionScanCache, _debug_name: &str, -) -> Option { - let last_exon = wt.exons.last().unwrap(); +) -> bool { + // A copy: exon A is extended in place below, and its original ends are + // still read afterwards. + let last_exon = *new_wt.exons.last().unwrap(); // Mate-boundary detection: STAR canonSJ[iex] = -3 (stitchAlignToTranscript.cpp:402) // When crossing from mate1 to mate2 (or vice versa), skip junction scoring and @@ -1212,10 +1288,10 @@ fn stitch_align_to_transcript( // mates overlap (genome positions interleaved), naturally limiting WTs to 2 exons. // rustar-aligner's overlap-trimming allows continued stitching, inflating the matched-base count. // Fix: if the WT already has exons from BOTH mates, a second crossing is invalid. - let has_m0 = wt.exons.iter().any(|e| e.mate_id == 0); - let has_m1 = wt.exons.iter().any(|e| e.mate_id == 1); + let has_m0 = new_wt.exons.iter().any(|e| e.mate_id == 0); + let has_m1 = new_wt.exons.iter().any(|e| e.mate_id == 1); if has_m0 && has_m1 { - return None; + return false; } // STAR condition (stitchAlignToTranscript.cpp:352): // gBstart + trA->exons[0][EX_R] + nBasesMax >= trA->exons[0][EX_G] || EX_G < EX_R @@ -1231,15 +1307,15 @@ fn stitch_align_to_transcript( // combined_start_mate1 < len_mate2_exon - len_mate1 + (P_mate2 - P_mate1) let combined_start_mate1 = (read_seq.len() as i64) - (wa.read_pos as i64) - (wa.length as i64); - let first_exon = &wt.exons[0]; + let first_exon = &new_wt.exons[0]; let len_mate2_exon = (first_exon.read_end - first_exon.read_start) as i64; let len_mate1 = wa.length as i64; let p_diff = first_exon.genome_start as i64 - wa.sa_pos as i64; // P_mate2 - P_mate1 if combined_start_mate1 < len_mate2_exon - len_mate1 + p_diff { - return None; + return false; } } else { - let first_exon = &wt.exons[0]; + let first_exon = &new_wt.exons[0]; let ex_g = first_exon.genome_start; let ex_r = first_exon.read_start as u64; // STAR stitchAlignToTranscript.cpp line 352: @@ -1248,21 +1324,20 @@ fn stitch_align_to_transcript( // reject if EX_G >= EX_R && gBstart + EX_R < EX_G let fwd_reject = ex_g >= ex_r && wa.genome_pos + ex_r < ex_g; if fwd_reject { - return None; + return false; } } // Forward-gap check: alignMatesGapMax (disabled when 0) let genome_gap = wa.genome_pos.saturating_sub(last_exon.genome_end); if align_mates_gap_max > 0 && genome_gap > align_mates_gap_max { - return None; + return false; } - let mut new_wt = wt.clone(); // STAR stitchAlignToTranscript.cpp:374-381: right-extend mate A to fragment boundary. // extendAlign(R, G, rAend+1, gAend+1, 1, 1, DEF_readSeqLengthMax, nMatch, nMM, ...) // extendToEnd=false (local mode); spacer at position 150 naturally stops the extension. // nMatch approximated as score + n_mismatch (exact for seeds-only paths, approx otherwise). - let n_match_m1 = (wt.score + wt.n_mismatch as i32).max(0) as usize; + let n_match_m1 = (new_wt.score + new_wt.n_mismatch as i32).max(0) as usize; let right_ext = extend_alignment_cached( jcache, read_seq, @@ -1270,7 +1345,7 @@ fn stitch_align_to_transcript( last_exon.genome_end, 1, 10_000, // DEF_readSeqLengthMax (STAR); SPACER stop applies before this limit - wt.n_mismatch, + new_wt.n_mismatch, n_match_m1, scorer.n_mm_max, scorer.p_mm_max, @@ -1347,7 +1422,7 @@ fn stitch_align_to_transcript( new_wt.n_mismatch += left_ext.n_mismatch; } - return Some(new_wt); + return true; } // Overlap trimming: if new WA overlaps previous exon in read coords, shift start right @@ -1358,7 +1433,7 @@ fn stitch_align_to_transcript( if last_exon.read_end > eff_read_pos { let overlap = last_exon.read_end - eff_read_pos; if overlap >= eff_length { - return None; // Fully consumed + return false; // Fully consumed } eff_read_pos = last_exon.read_end; eff_genome_pos += overlap as u64; @@ -1369,7 +1444,7 @@ fn stitch_align_to_transcript( if last_exon.genome_end > eff_genome_pos && eff_genome_pos > last_exon.genome_start { let g_overlap = (last_exon.genome_end - eff_genome_pos) as usize; if g_overlap >= eff_length { - return None; // Fully consumed + return false; // Fully consumed } eff_read_pos += g_overlap; eff_genome_pos += g_overlap as u64; @@ -1383,7 +1458,7 @@ fn stitch_align_to_transcript( // that only starts before A's genome end (`genome_gap < 0`) is an // insertion with overlapping seeds: its score drops by the overlap below. if read_gap < 0 || eff_genome_pos + eff_length as u64 <= last_exon.genome_end { - return None; + return false; } // Seeds on one diagonal (adjacent or separated by a gap-free filler) are @@ -1393,10 +1468,9 @@ fn stitch_align_to_transcript( // end extension of whichever seed is kept, which is bounded by the // max-score and mismatch rules of `extendAlign`, not by the gap scoring. if read_gap == genome_gap { - return None; + return false; } - let mut new_wt = wt.clone(); let mut d_score: i32 = 0; let mut gap_mm: u32 = 0; @@ -1409,7 +1483,7 @@ fn stitch_align_to_transcript( // STAR: Del > alignIntronMax → reject (return -1000003) if del > scorer.align_intron_max && scorer.align_intron_max > 0 { - return None; + return false; } // No early exit on a short exon B: after the overlap trim B can be shorter than @@ -1443,7 +1517,7 @@ fn stitch_align_to_transcript( if pj.motif == 0 && (eff_length <= pj.shift_right as usize || prev_len <= pj.shift_left as usize) { - return None; // STAR -1000006: repeats around a non-canonical junction + return false; // STAR -1000006: repeats around a non-canonical junction } let donor_fwd = index.sa_pos_to_forward(donor_sa, cluster.is_reverse, del as usize); let motif = scorer.detect_splice_motif(donor_fwd, del, &index.genome); @@ -1481,7 +1555,7 @@ fn stitch_align_to_transcript( (last_exon.read_end - last_exon.read_start) as i32 + jr_shift }; if shrunk_len < 1 { - return None; + return false; } } @@ -1577,7 +1651,7 @@ fn stitch_align_to_transcript( if is_splice { // Check stitch mismatch limit if !scorer.stitch_mismatch_allowed(&motif, gap_mm) { - return None; + return false; } // Annotated strand of a database hit (STAR `sjdbStrand`); the sjdb @@ -1623,7 +1697,7 @@ fn stitch_align_to_transcript( let sl = pj.shift_left as usize; let prev_len = last_exon.read_end - last_exon.read_start; if eff_length <= sl || prev_len <= sl { - return None; // STAR -1000006 + return false; // STAR -1000006 } let shifted = if cluster.is_reverse { jr_shift - pj.shift_left as i32 @@ -1633,7 +1707,7 @@ fn stitch_align_to_transcript( // STAR: rAend+jR >= rBend, i.e. the shift runs past B. In the // reverse scan space the roles of A and B are exchanged. if shifted >= (shared + eff_length) as i32 || prev_len as i32 + shifted < 1 { - return None; // STAR -1000006 + return false; // STAR -1000006 } jr_exon = shifted; jj_l = u32::from(pj.shift_left); @@ -1808,7 +1882,7 @@ fn stitch_align_to_transcript( // `outFilterMismatchNmaxTotal`; the mismatch ratio is applied to the final // best transcript by `mappedFilter`, not while stitching. if total_mm > scorer.n_mm_max { - return None; + return false; } // Update working transcript @@ -1829,7 +1903,7 @@ fn stitch_align_to_transcript( new_wt.n_anchor += 1; } - Some(new_wt) + true } /// Stitch seeds within a cluster using recursive combinatorial stitching. @@ -1882,6 +1956,37 @@ pub fn stitch_seeds_with_jdb( /// Check if two transcripts overlap in their exon blocks. /// Returns total overlapping bases (where exons share the same read-genome diagonal). /// Used for deduplication: if new transcript is a subset of existing, drop it. +/// [`stitch_align_in_place`] on a copy of `wt`: the stitched transcript, or +/// `None` where STAR rejects the stitch. +#[cfg(test)] +#[allow(clippy::too_many_arguments)] +fn stitch_align_to_transcript( + wt: &WorkingTranscript, + wa: &WindowAlignment, + read_seq: &[u8], + index: &GenomeIndex, + scorer: &AlignmentScorer, + cluster: &SeedCluster, + junction_db: Option<&crate::junction::SpliceJunctionDb>, + align_mates_gap_max: u64, + jcache: &mut crate::align::score::JunctionScanCache, + debug_name: &str, +) -> Option { + let mut new_wt = wt.clone(); + stitch_align_in_place( + &mut new_wt, + wa, + read_seq, + index, + scorer, + cluster, + junction_db, + align_mates_gap_max, + jcache, + debug_name, + ) + .then_some(new_wt) +} fn blocks_overlap(t1_exons: &[ExonBlock], t2_exons: &[ExonBlock]) -> u32 { let mut overlap = 0u32; let mut i = 0; @@ -2598,7 +2703,7 @@ fn passes_finalization_filters( #[allow(clippy::too_many_arguments)] fn stitch_recurse( i_a: usize, - wt: WorkingTranscript, + wt: &mut WorkingTranscript, wa_entries: &[WindowAlignment], read_seq: &[u8], index: &GenomeIndex, @@ -2627,240 +2732,251 @@ fn stitch_recurse( // Base case: all entries considered if i_a >= wa_entries.len() { if !wt.exons.is_empty() { - // STAR stitchWindowAligns.cpp lines 21-81: apply left and right extensions - // inside the stitcher BEFORE the score-gate/dedup. This matches STAR's behavior - // where extensions boost the correct WT's score so the score-range filter can - // correctly eliminate spurious WTs with short first exons. - // EXTEND_ORDER=1: extend 5' of read first (left for fwd, right for rev). - let mut wt = wt; - let zero_ext = ExtendResult { - extend_len: 0, - max_score: 0, - n_mismatch: 0, - }; + // The leaf extends `wt` in place; the branches above it still need + // it as it was, so it is put back once the leaf is done (only the + // extensions change it). + let mut mark = None; + 'leaf: { + // STAR stitchWindowAligns.cpp lines 21-81: apply left and right extensions + // inside the stitcher BEFORE the score-gate/dedup. This matches STAR's behavior + // where extensions boost the correct WT's score so the score-range filter can + // correctly eliminate spurious WTs with short first exons. + // EXTEND_ORDER=1: extend 5' of read first (left for fwd, right for rev). + let zero_ext = ExtendResult { + extend_len: 0, + max_score: 0, + n_mismatch: 0, + }; - let do_left_first = !original_is_reverse; - let first_ext = if do_left_first { - // Left extension (5' of forward read) - if wt.read_start > 0 { - extend_alignment_cached( - jcache, - read_seq, - wt.read_start, - wt.genome_start, - -1, - wt.read_start, - wt.n_mismatch, - wt.read_end - wt.read_start, - scorer.n_mm_max, - scorer.p_mm_max, - index, - cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize - ) - } else { - zero_ext.clone() - } - } else { - // Right extension (5' of reverse read) - if wt.read_end < read_seq.len() { - extend_alignment_cached( - jcache, - read_seq, - wt.read_end, - wt.genome_end, - 1, - read_seq.len() - wt.read_end, - wt.n_mismatch, - wt.read_end - wt.read_start, - scorer.n_mm_max, - scorer.p_mm_max, - index, - cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize - ) + let do_left_first = !original_is_reverse; + let first_ext = if do_left_first { + // Left extension (5' of forward read) + if wt.read_start > 0 { + extend_alignment_cached( + jcache, + read_seq, + wt.read_start, + wt.genome_start, + -1, + wt.read_start, + wt.n_mismatch, + wt.read_end - wt.read_start, + scorer.n_mm_max, + scorer.p_mm_max, + index, + cluster.is_reverse, + false, // internal stitch extension: alignEndsType applied at finalize + ) + } else { + zero_ext.clone() + } } else { - zero_ext.clone() - } - }; + // Right extension (5' of reverse read) + if wt.read_end < read_seq.len() { + extend_alignment_cached( + jcache, + read_seq, + wt.read_end, + wt.genome_end, + 1, + read_seq.len() - wt.read_end, + wt.n_mismatch, + wt.read_end - wt.read_start, + scorer.n_mm_max, + scorer.p_mm_max, + index, + cluster.is_reverse, + false, // internal stitch extension: alignEndsType applied at finalize + ) + } else { + zero_ext.clone() + } + }; - let len_after_first = (wt.read_end - wt.read_start) + first_ext.extend_len; - let mm_after_first = wt.n_mismatch + first_ext.n_mismatch; - - let second_ext = if do_left_first { - // Right extension (3' of forward read) - let read_end_after = wt.read_end; // right boundary unchanged by left ext - if read_end_after < read_seq.len() { - extend_alignment_cached( - jcache, - read_seq, - read_end_after, - wt.genome_end, - 1, - read_seq.len() - read_end_after, - mm_after_first, - len_after_first, - scorer.n_mm_max, - scorer.p_mm_max, - index, - cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize - ) + let len_after_first = (wt.read_end - wt.read_start) + first_ext.extend_len; + let mm_after_first = wt.n_mismatch + first_ext.n_mismatch; + + let second_ext = if do_left_first { + // Right extension (3' of forward read) + let read_end_after = wt.read_end; // right boundary unchanged by left ext + if read_end_after < read_seq.len() { + extend_alignment_cached( + jcache, + read_seq, + read_end_after, + wt.genome_end, + 1, + read_seq.len() - read_end_after, + mm_after_first, + len_after_first, + scorer.n_mm_max, + scorer.p_mm_max, + index, + cluster.is_reverse, + false, // internal stitch extension: alignEndsType applied at finalize + ) + } else { + zero_ext + } } else { - zero_ext + // Left extension (3' of reverse read) + let read_start_after = wt.read_start; // left boundary unchanged by right ext + if read_start_after > 0 { + extend_alignment_cached( + jcache, + read_seq, + read_start_after, + wt.genome_start, + -1, + read_start_after, + mm_after_first, + len_after_first, + scorer.n_mm_max, + scorer.p_mm_max, + index, + cluster.is_reverse, + false, // internal stitch extension: alignEndsType applied at finalize + ) + } else { + zero_ext + } + }; + + // Apply the extensions to wt + if first_ext.extend_len > 0 || second_ext.extend_len > 0 { + mark = Some(wt.mark_ends()); } - } else { - // Left extension (3' of reverse read) - let read_start_after = wt.read_start; // left boundary unchanged by right ext - if read_start_after > 0 { - extend_alignment_cached( - jcache, - read_seq, - read_start_after, - wt.genome_start, - -1, - read_start_after, - mm_after_first, - len_after_first, - scorer.n_mm_max, - scorer.p_mm_max, - index, - cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize - ) + if do_left_first { + if first_ext.extend_len > 0 { + let first = wt.exons.first_mut().unwrap(); + first.read_start -= first_ext.extend_len; + first.genome_start -= first_ext.extend_len as u64; + wt.score += first_ext.max_score; + wt.n_mismatch += first_ext.n_mismatch; + wt.read_start = first.read_start; + wt.genome_start = first.genome_start; + } + if second_ext.extend_len > 0 { + let last = wt.exons.last_mut().unwrap(); + last.read_end += second_ext.extend_len; + last.genome_end += second_ext.extend_len as u64; + wt.score += second_ext.max_score; + wt.n_mismatch += second_ext.n_mismatch; + wt.read_end = last.read_end; + wt.genome_end = last.genome_end; + } } else { - zero_ext + if first_ext.extend_len > 0 { + let last = wt.exons.last_mut().unwrap(); + last.read_end += first_ext.extend_len; + last.genome_end += first_ext.extend_len as u64; + wt.score += first_ext.max_score; + wt.n_mismatch += first_ext.n_mismatch; + wt.read_end = last.read_end; + wt.genome_end = last.genome_end; + } + if second_ext.extend_len > 0 { + let first = wt.exons.first_mut().unwrap(); + first.read_start -= second_ext.extend_len; + first.genome_start -= second_ext.extend_len as u64; + wt.score += second_ext.max_score; + wt.n_mismatch += second_ext.n_mismatch; + wt.read_start = first.read_start; + wt.genome_start = first.genome_start; + } } - }; - // Apply the extensions to wt - if do_left_first { - if first_ext.extend_len > 0 { - let first = wt.exons.first_mut().unwrap(); - first.read_start -= first_ext.extend_len; - first.genome_start -= first_ext.extend_len as u64; - wt.score += first_ext.max_score; - wt.n_mismatch += first_ext.n_mismatch; - wt.read_start = first.read_start; - wt.genome_start = first.genome_start; - } - if second_ext.extend_len > 0 { - let last = wt.exons.last_mut().unwrap(); - last.read_end += second_ext.extend_len; - last.genome_end += second_ext.extend_len as u64; - wt.score += second_ext.max_score; - wt.n_mismatch += second_ext.n_mismatch; - wt.read_end = last.read_end; - wt.genome_end = last.genome_end; + // STAR's intron filters run here, after the extensions and before the + // window's dedup (`stitchWindowAligns.cpp:130-152` vs `:260-302`): a + // rejected transcript must not evict one it covers. For a pair the + // junctions of both mates are on this one transcript. + if !scorer.intron_filter.passes( + wt.junction_motifs + .iter() + .zip(wt.junction_annotated.iter()) + .zip(wt.junction_strand.iter()) + .map(|((m, a), s)| (m, a, s)), + ) { + break 'leaf; } - } else { - if first_ext.extend_len > 0 { - let last = wt.exons.last_mut().unwrap(); - last.read_end += first_ext.extend_len; - last.genome_end += first_ext.extend_len as u64; - wt.score += first_ext.max_score; - wt.n_mismatch += first_ext.n_mismatch; - wt.read_end = last.read_end; - wt.genome_end = last.genome_end; - } - if second_ext.extend_len > 0 { - let first = wt.exons.first_mut().unwrap(); - first.read_start -= second_ext.extend_len; - first.genome_start -= second_ext.extend_len as u64; - wt.score += second_ext.max_score; - wt.n_mismatch += second_ext.n_mismatch; - wt.read_start = first.read_start; - wt.genome_start = first.genome_start; - } - } - // STAR's intron filters run here, after the extensions and before the - // window's dedup (`stitchWindowAligns.cpp:130-152` vs `:260-302`): a - // rejected transcript must not evict one it covers. For a pair the - // junctions of both mates are on this one transcript. - if !scorer.intron_filter.passes( - wt.junction_motifs - .iter() - .zip(wt.junction_annotated.iter()) - .zip(wt.junction_strand.iter()) - .map(|((m, a), s)| (m, a, s)), - ) { - return; - } + // STAR rejects a transcript on its exon overhangs, spliced-mate length and + // mate placement before the window's dedup (`stitchWindowAligns.cpp:96-219`), + // so a rejected transcript never evicts the one it covers. + if !passes_finalization_filters(wt, read_seq, scorer) { + break 'leaf; + } - // STAR rejects a transcript on its exon overhangs, spliced-mate length and - // mate placement before the window's dedup (`stitchWindowAligns.cpp:96-219`), - // so a rejected transcript never evicts the one it covers. - if !passes_finalization_filters(&wt, read_seq, scorer) { - return; - } + // STAR's score gate (`stitchWindowAligns.cpp:244-250`); see + // [`window_score_gate`]. + { + let score = wt_final_score(wt, scorer); + // Single-end seeds carry mate_id 2; STAR's iFrag is 0 there. + let frag = match (wt.exons.first(), wt.exons.last()) { + (Some(f), Some(l)) if f.mate_id == l.mate_id => Some(if f.mate_id >= 2 { + 0 + } else { + f.mate_id as usize + }), + _ => None, + }; + // An empty window compares against STAR's `trInit`, score 0. + let window_best = transcripts.first().map_or(0, |t| t.window_score); + if !window_score_gate(score, frag, window_best, max_score_mate, scorer) { + break 'leaf; + } + } - // STAR's score gate (`stitchWindowAligns.cpp:244-250`); see - // [`window_score_gate`]. - { - let score = wt_final_score(&wt, scorer); - // Single-end seeds carry mate_id 2; STAR's iFrag is 0 there. - let frag = match (wt.exons.first(), wt.exons.last()) { - (Some(f), Some(l)) if f.mate_id == l.mate_id => Some(if f.mate_id >= 2 { - 0 + // Dedup like STAR (`stitchWindowAligns.cpp:260-285`). The score compared is + // the final one, with the genomic-length penalty. A transcript that + // contains an older one replaces it whatever their scores: only a + // contained transcript scoring lower than its container is dropped. + let wt_score = wt_final_score(wt, scorer); + let wt_len = wt_mapped_length(wt); + let mut dominated = false; + let mut idx = 0; + while idx < transcripts.len() { + let existing = &transcripts[idx]; + let overlap = blocks_overlap(&wt.exons, &existing.exons); + let u_new = wt_len - overlap; + let u_old = existing.mapped_len - overlap; + if u_new == 0 && wt_score < existing.window_score { + dominated = true; + break; + } else if u_old == 0 { + // STAR drops the contained transcript on the spot, even if a + // later one then dominates the new transcript. + transcripts.remove(idx); } else { - f.mate_id as usize - }), - _ => None, - }; - // An empty window compares against STAR's `trInit`, score 0. - let window_best = transcripts.first().map_or(0, |t| t.window_score); - if !window_score_gate(score, frag, window_best, max_score_mate, scorer) { - return; + idx += 1; + } } - } - // Dedup like STAR (`stitchWindowAligns.cpp:260-285`). The score compared is - // the final one, with the genomic-length penalty. A transcript that - // contains an older one replaces it whatever their scores: only a - // contained transcript scoring lower than its container is dropped. - let wt_score = wt_final_score(&wt, scorer); - let wt_len = wt_mapped_length(&wt); - let mut dominated = false; - let mut idx = 0; - while idx < transcripts.len() { - let existing = &transcripts[idx]; - let overlap = blocks_overlap(&wt.exons, &existing.exons); - let u_new = wt_len - overlap; - let u_old = existing.mapped_len - overlap; - if u_new == 0 && wt_score < existing.window_score { - dominated = true; - break; - } else if u_old == 0 { - // STAR drops the contained transcript on the spot, even if a - // later one then dominates the new transcript. - transcripts.remove(idx); - } else { - idx += 1; + if !dominated { + // STAR keeps the window's transcripts sorted by score, then by shorter + // genomic span; a full list drops its last transcript. The order + // decides ties downstream. + let span = |t: &WorkingTranscript| match (t.exons.first(), t.exons.last()) { + (Some(first), Some(last)) => last.genome_end - first.genome_start, + _ => 0, + }; + let wt_span = span(wt); + let pos = transcripts + .iter() + .position(|t| { + let score = t.window_score; + wt_score > score || (wt_score == score && wt_span < span(t)) + }) + .unwrap_or(transcripts.len()); + let mut recorded = wt.clone(); + recorded.window_score = wt_score; + recorded.mapped_len = wt_len; + transcripts.insert(pos, recorded); + transcripts.truncate(max_transcripts); } } - - if !dominated { - // STAR keeps the window's transcripts sorted by score, then by shorter - // genomic span; a full list drops its last transcript. The order - // decides ties downstream. - let span = |t: &WorkingTranscript| match (t.exons.first(), t.exons.last()) { - (Some(first), Some(last)) => last.genome_end - first.genome_start, - _ => 0, - }; - let wt_span = span(&wt); - let pos = transcripts - .iter() - .position(|t| { - let score = t.window_score; - wt_score > score || (wt_score == score && wt_span < span(t)) - }) - .unwrap_or(transcripts.len()); - let mut wt = wt; - wt.window_score = wt_score; - wt.mapped_len = wt_len; - transcripts.insert(pos, wt); - transcripts.truncate(max_transcripts); + if let Some(mark) = mark { + wt.restore(&mark); } } return; @@ -2871,7 +2987,8 @@ fn stitch_recurse( // INCLUDE branch: try stitching wa_entries[i_a] to transcript if wt.exons.is_empty() { // First seed: create initial transcript - let mut new_wt = wt.clone(); + let mark = wt.mark(); + let new_wt = &mut *wt; new_wt.exons.push(ExonBlock { read_start: wa.read_pos, read_end: wa.read_pos + wa.length, @@ -2907,10 +3024,13 @@ fn stitch_recurse( jcache, debug_name, ); + wt.restore(&mark); } else { - // Try stitching this seed onto the existing transcript - if let Some(new_wt) = stitch_align_to_transcript( - &wt, + // Try stitching this seed onto the existing transcript, in place; + // the exclude branch below gets it back as it was. + let mark = wt.mark(); + if stitch_align_in_place( + wt, wa, read_seq, index, @@ -2923,7 +3043,7 @@ fn stitch_recurse( ) { stitch_recurse( i_a + 1, - new_wt, + wt, wa_entries, read_seq, index, @@ -2940,6 +3060,7 @@ fn stitch_recurse( debug_name, ); } + wt.restore(&mark); } // EXCLUDE branch: skip wa_entries[i_a]. @@ -3066,12 +3187,12 @@ pub(crate) fn split_combined_wt( for (i, ex) in wt.exons.iter().enumerate() { if ex.mate_id == 0 { - let mut ex1 = ex.clone(); + let mut ex1 = *ex; ex1.read_start = ex1.read_start.saturating_sub(m1_offset); ex1.read_end -= m1_offset; m1_exons.push(ex1); } else if ex.mate_id == 1 { - let mut ex2 = ex.clone(); + let mut ex2 = *ex; ex2.read_start = ex2.read_start.saturating_sub(m2_offset); ex2.read_end -= m2_offset; m2_exons.push(ex2); @@ -3513,7 +3634,7 @@ pub(crate) fn stitch_seeds_core( stitch_recurse( 0, - WorkingTranscript::new(), + &mut WorkingTranscript::new(), &wa_entries, stitch_read, index, From e7853ecba7a8803258e7a4c679c2c5ab6b37acd2 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 29 Oct 2026 19:46:00 +0100 Subject: [PATCH 15/17] Junction scan: fetch the scan's bases in bulk Phase 2 of find_best_junction_position read the read base, the upstream and downstream genome bases and the four motif bases position by position, through the genome view. It now copies them into stack buffers first (SeqView::bases_into, one copy per run; the reverse strand's donor run is read backwards), computes the score1 steps in one branch-free pass, and scans the arrays: same positions, same order, same tie rules, same early stop. Positions off the read or the genome become the out-of-range sentinel, which every test there treats as STAR's skip. Scans longer than 512 positions, or whose coordinates would leave the genome, take the position-by-position path, kept as scan_phase2_scalar. buffered_junction_scan_matches_the_scalar_one checks the two on 20,000 random scans (both strands, deletions and introns, N bases, genome ends). Output is unchanged on the yeast, chr21, fusion and STARsolo runs, and the full fusion set and chr21 PE still match STAR 2.7.11b. Single-threaded instructions: chr21 100k SE -16%, yeast 100k PE -13%, fusion -1%. Co-Authored-By: Claude Opus 5.5 --- src/align/score.rs | 431 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 333 insertions(+), 98 deletions(-) diff --git a/src/align/score.rs b/src/align/score.rs index 3fc291f..7866b9a 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -478,108 +478,27 @@ impl AlignmentScorer { // jr1 is now one past where we stopped; the scan will start from jr1 // Phase 2: Scan RIGHT through all jR1 positions - score1 = 0; - let mut max_score2 = i32::MIN; - let mut best_jr: i32 = 0; - let mut best_motif = SpliceMotif::NonCanonical; - let mut best_motif_score = self.score_gap_noncan; - // `del` does not change across the scan, so whether the motif branch // runs is decided once rather than per iteration. let motif_in_range = del >= self.align_intron_min as i64 && del <= self.align_intron_max as i64; - let max_motif = self - .score_gap - .max(self.score_gap_gcag) - .max(self.score_gap_atac) - .max(self.score_gap_noncan); - - // The four motif bases sit at `donor`, `donor+1`, `donor+del-2` and - // `donor+del-1`, and `donor` moves exactly one base per iteration, so - // consecutive iterations share two of the four. Carrying the window - // instead of re-fetching it turns four genome reads per position into - // two. `window` is `None` until the first motif iteration primes it. - let mut window: Option = None; - - loop { - let ri = r_a_end_inc as i64 + jr1 as i64; - if ri >= 0 && (ri as usize) < read_seq.len() { - let read_base = read_seq[ri as usize]; - let g_up_pos = g_a_end_inc as i64 + jr1 as i64; - let g_dn_pos = g_b_start1 + jr1 as i64; - - if g_up_pos >= 0 && g_dn_pos >= 0 { - let gu = seq.base((g_up_pos as u64 + genome_offset) as usize); - let gd = seq.base((g_dn_pos as u64 + genome_offset) as usize); - - if gu < 4 && gd < 4 { - if read_base == gu && read_base != gd { - score1 += 1; - } else if read_base != gu && read_base == gd { - score1 -= 1; - } - } - } - } - - // Check splice motif at this junction position - if motif_in_range { - // Donor position in SA space: one past the last donor-exon base - let donor_sa = (g_a_end_inc as i64 + jr1 as i64 + 1) as u64; - // Convert to forward genome coordinates for motif detection - let donor_fwd = if is_reverse { - n_genome - donor_sa - del as u64 - } else { - donor_sa - }; - let w = match window.as_mut() { - Some(w) => { - w.slide_to(donor_fwd, del as u64, seq); - &*w - } - None => window.insert(MotifWindow::at(donor_fwd, del as u64, seq)), - }; - let motif = w.motif(); - let motif_score = self.score_splice_junction(motif); - let score2 = score1 + motif_score; - - if score2 > max_score2 { - max_score2 = score2; - best_jr = jr1; - best_motif = motif; - best_motif_score = motif_score; - } - } else { - // Deletion (Del < alignIntronMin) or out-of-range gap — no motif, - // pure positional score. STAR: jCan1=-1, jPen1=0, Score2=Score1. - // For reverse strand, rightmost-in-RC wins ties (= leftmost in forward). - if score1 > max_score2 || (is_reverse && score1 == max_score2) { - max_score2 = score1; - best_jr = jr1; - best_motif = SpliceMotif::NonCanonical; - best_motif_score = 0; - } - } - - // `score1` rises by at most 1 a position, so once even that cannot - // lift a later position past the best one the rest of the scan is - // settled. Ties win only in the reverse-strand deletion case. - let scan_end = r_gap as i32 + next_seed_len as i32; - let bound = score1 + (scan_end - 1 - jr1) + if motif_in_range { max_motif } else { 0 }; - let tie_wins = !motif_in_range && is_reverse; - if bound < max_score2 || (bound == max_score2 && !tie_wins) { - break; - } - - jr1 += 1; - // STAR: jR1 < int(rBend) - int(rAend) where rBend = rBstart + L - 1 - // This equals jR1 < rGap + L, allowing the scan to extend into the B seed. - // This is critical for finding canonical junctions when seeds are adjacent - // (rGap=0) but the junction needs to shift into the B seed's territory. - if jr1 >= r_gap as i32 + next_seed_len as i32 { - break; - } - } + let scan = ScanArgs { + read_seq, + r_a_end_inc: r_a_end_inc as i64, + g_a_end_inc: g_a_end_inc as i64, + g_b_start1, + jr_start: jr1, + scan_end: r_gap as i32 + next_seed_len as i32, + del, + motif_in_range, + is_reverse, + n_genome, + genome_offset, + seq, + }; + let (mut best_jr, mut best_motif, mut best_motif_score) = self + .scan_phase2_buffered(&scan) + .unwrap_or_else(|| self.scan_phase2_scalar(&scan)); // Phase 3: Repeat detection around best_jr + left-flush for non-canonical/deletion // Count matching bases between donor and acceptor sides around best_jr @@ -665,6 +584,228 @@ impl AlignmentScorer { ) } + /// The best motif score: the most a position's motif can add in the scan. + fn max_motif_score(&self) -> i32 { + self.score_gap + .max(self.score_gap_gcag) + .max(self.score_gap_atac) + .max(self.score_gap_noncan) + } + + /// Phase 2 of [`Self::find_best_junction_position`] on bases fetched in + /// bulk: the read, upstream and downstream bases (and, for an intron, the + /// donor and acceptor runs) of every position are copied into stack + /// buffers first, so the scan itself is plain array reads. Same positions, + /// same order, same ties and same early stop as [`Self::scan_phase2_scalar`]; + /// `None` hands the rare scans that do not fit the buffers to it. + fn scan_phase2_buffered(&self, scan: &ScanArgs<'_>) -> Option<(i32, SpliceMotif, i32)> { + use crate::genome::OUT_OF_RANGE; + let n = (scan.scan_end - scan.jr_start).max(1) as usize; + if n + 1 > SCAN_BUF { + return None; + } + let j0 = scan.jr_start as i64; + let up0 = scan.g_a_end_inc + j0; + let dn0 = scan.g_b_start1 + j0; + let donor_sa0 = scan.g_a_end_inc + j0 + 1; + if up0 < 0 || dn0 < 0 || donor_sa0 < 0 { + return None; + } + + let mut read = [OUT_OF_RANGE; SCAN_BUF]; + for (k, slot) in read[..n].iter_mut().enumerate() { + let ri = scan.r_a_end_inc + j0 + k as i64; + if ri >= 0 && (ri as usize) < scan.read_seq.len() { + *slot = scan.read_seq[ri as usize]; + } + } + let mut up = [OUT_OF_RANGE; SCAN_BUF]; + let mut dn = [OUT_OF_RANGE; SCAN_BUF]; + scan.seq + .bases_into(up0 as usize + scan.genome_offset as usize, &mut up[..n]); + scan.seq + .bases_into(dn0 as usize + scan.genome_offset as usize, &mut dn[..n]); + // score1 change at each position: +1 where the read matches upstream + // only, -1 where it matches downstream only, 0 otherwise or where + // either genome base is not A/C/G/T (the sentinel included). + let mut delta = [0i8; SCAN_BUF]; + for k in 0..n { + let (rb, ub, db) = (read[k], up[k], dn[k]); + let ok = ub < 4 && db < 4; + let plus = ok && rb == ub && rb != db; + let minus = ok && rb != ub && rb == db; + delta[k] = plus as i8 - minus as i8; + } + + // Donor and acceptor runs, forward strand. Forward reads move the + // donor right one base a position; reverse reads move it left, so + // their run starts at the last position's donor and is read backwards. + let mut don = [OUT_OF_RANGE; SCAN_BUF]; + let mut acc = [OUT_OF_RANGE; SCAN_BUF]; + let del = scan.del as u64; + if scan.motif_in_range { + let lo = if scan.is_reverse { + let last = donor_sa0 as u64 + (n as u64 - 1) + del; + if last > scan.n_genome { + return None; + } + scan.n_genome - last + } else { + donor_sa0 as u64 + }; + scan.seq.bases_into(lo as usize, &mut don[..=n]); + scan.seq.bases_into((lo + del - 2) as usize, &mut acc[..=n]); + } + + let max_motif = self.max_motif_score(); + let mut score1: i32 = 0; + let mut max_score2 = i32::MIN; + let mut best_jr: i32 = 0; + let mut best_motif = SpliceMotif::NonCanonical; + let mut best_motif_score = self.score_gap_noncan; + for (k, &step) in delta[..n].iter().enumerate() { + let jr1 = scan.jr_start + k as i32; + score1 += i32::from(step); + if scan.motif_in_range { + let i = if scan.is_reverse { n - 1 - k } else { k }; + let motif = MotifWindow { + donor: 0, + d1: don[i], + d2: don[i + 1], + a1: acc[i], + a2: acc[i + 1], + } + .motif(); + let motif_score = self.score_splice_junction(motif); + let score2 = score1 + motif_score; + if score2 > max_score2 { + max_score2 = score2; + best_jr = jr1; + best_motif = motif; + best_motif_score = motif_score; + } + } else if score1 > max_score2 || (scan.is_reverse && score1 == max_score2) { + max_score2 = score1; + best_jr = jr1; + best_motif = SpliceMotif::NonCanonical; + best_motif_score = 0; + } + let bound = score1 + + (scan.scan_end - 1 - jr1) + + if scan.motif_in_range { max_motif } else { 0 }; + let tie_wins = !scan.motif_in_range && scan.is_reverse; + if bound < max_score2 || (bound == max_score2 && !tie_wins) { + break; + } + } + Some((best_jr, best_motif, best_motif_score)) + } + + /// Phase 2 of [`Self::find_best_junction_position`], one position at a + /// time: the reference for [`Self::scan_phase2_buffered`]. + fn scan_phase2_scalar(&self, scan: &ScanArgs<'_>) -> (i32, SpliceMotif, i32) { + let mut jr1 = scan.jr_start; + let mut score1: i32 = 0; + let mut max_score2 = i32::MIN; + let mut best_jr: i32 = 0; + let mut best_motif = SpliceMotif::NonCanonical; + let mut best_motif_score = self.score_gap_noncan; + let max_motif = self.max_motif_score(); + + // The four motif bases sit at `donor`, `donor+1`, `donor+del-2` and + // `donor+del-1`, and `donor` moves exactly one base per iteration, so + // consecutive iterations share two of the four. Carrying the window + // instead of re-fetching it turns four genome reads per position into + // two. `window` is `None` until the first motif iteration primes it. + let mut window: Option = None; + + loop { + let ri = scan.r_a_end_inc + jr1 as i64; + if ri >= 0 && (ri as usize) < scan.read_seq.len() { + let read_base = scan.read_seq[ri as usize]; + let g_up_pos = scan.g_a_end_inc + jr1 as i64; + let g_dn_pos = scan.g_b_start1 + jr1 as i64; + + if g_up_pos >= 0 && g_dn_pos >= 0 { + let gu = scan + .seq + .base((g_up_pos as u64 + scan.genome_offset) as usize); + let gd = scan + .seq + .base((g_dn_pos as u64 + scan.genome_offset) as usize); + + if gu < 4 && gd < 4 { + if read_base == gu && read_base != gd { + score1 += 1; + } else if read_base != gu && read_base == gd { + score1 -= 1; + } + } + } + } + + // Check splice motif at this junction position + if scan.motif_in_range { + // Donor position in SA space: one past the last donor-exon base + let donor_sa = (scan.g_a_end_inc + jr1 as i64 + 1) as u64; + // Convert to forward genome coordinates for motif detection + let donor_fwd = if scan.is_reverse { + scan.n_genome - donor_sa - scan.del as u64 + } else { + donor_sa + }; + let w = match window.as_mut() { + Some(w) => { + w.slide_to(donor_fwd, scan.del as u64, scan.seq); + &*w + } + None => window.insert(MotifWindow::at(donor_fwd, scan.del as u64, scan.seq)), + }; + let motif = w.motif(); + let motif_score = self.score_splice_junction(motif); + let score2 = score1 + motif_score; + + if score2 > max_score2 { + max_score2 = score2; + best_jr = jr1; + best_motif = motif; + best_motif_score = motif_score; + } + } else { + // Deletion (Del < alignIntronMin) or out-of-range gap — no motif, + // pure positional score. STAR: jCan1=-1, jPen1=0, Score2=Score1. + // For reverse strand, rightmost-in-RC wins ties (= leftmost in forward). + if score1 > max_score2 || (scan.is_reverse && score1 == max_score2) { + max_score2 = score1; + best_jr = jr1; + best_motif = SpliceMotif::NonCanonical; + best_motif_score = 0; + } + } + + // `score1` rises by at most 1 a position, so once even that cannot + // lift a later position past the best one the rest of the scan is + // settled. Ties win only in the reverse-strand deletion case. + let scan_end = scan.scan_end; + let bound = + score1 + (scan_end - 1 - jr1) + if scan.motif_in_range { max_motif } else { 0 }; + let tie_wins = !scan.motif_in_range && scan.is_reverse; + if bound < max_score2 || (bound == max_score2 && !tie_wins) { + break; + } + + jr1 += 1; + // STAR: jR1 < int(rBend) - int(rAend) where rBend = rBstart + L - 1 + // This equals jR1 < rGap + L, allowing the scan to extend into the B seed. + // This is critical for finding canonical junctions when seeds are adjacent + // (rGap=0) but the junction needs to shift into the B seed's territory. + if jr1 >= scan.scan_end { + break; + } + } + (best_jr, best_motif, best_motif_score) + } + /// Detect splice junction motif (thin wrapper over the free function /// so `AlignmentScorer` callers keep working). pub fn detect_splice_motif( @@ -698,6 +839,27 @@ pub fn detect_splice_motif(donor_pos: u64, intron_len: u32, genome: &Genome) -> MotifWindow::at(donor_pos, intron_len as u64, genome.sequence.view()).motif() } +/// Positions a buffered junction scan holds; longer scans take the scalar path. +const SCAN_BUF: usize = 512; + +/// The inputs of a junction scan's phase 2, in STAR's inclusive coordinates. +struct ScanArgs<'a> { + read_seq: &'a [u8], + r_a_end_inc: i64, + g_a_end_inc: i64, + g_b_start1: i64, + /// First position, where phase 1 stopped. + jr_start: i32, + /// One past the last position: rGap + the B seed's length. + scan_end: i32, + del: i64, + motif_in_range: bool, + is_reverse: bool, + n_genome: u64, + genome_offset: u64, + seq: SeqView<'a>, +} + /// The four bases that decide a splice motif: the intron's first two and last /// two, on the forward strand. /// @@ -970,6 +1132,79 @@ mod tests { } } + /// The buffered phase 2 of the junction scan must pick the same position, + /// motif and motif score as the position-by-position one, on both strands, + /// for deletions and introns, near the genome's ends, with N bases in the + /// read and the genome, and with the early stop in play. + #[test] + fn buffered_junction_scan_matches_the_scalar_one() { + let mut state = 12345u64; + let mut next = move |m: u64| { + state = state + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + (state >> 33) % m + }; + // Splice motifs are rare in random sequence: seed GT..AG pairs. + let mut fwd: Vec = (0..3000).map(|_| next(4) as u8).collect(); + for _ in 0..120 { + let i = next(2900) as usize; + fwd[i] = 2; + fwd[i + 1] = 3; + } + for _ in 0..120 { + let i = next(2900) as usize; + fwd[i] = 0; + fwd[i + 1] = 2; + } + for _ in 0..30 { + fwd[next(3000) as usize] = 4; + } + let genome = make_test_genome(&fwd); + let seq = genome.sequence.view(); + let scorer = AlignmentScorer::from_params_minimal(); + let mut compared = 0; + for _ in 0..20_000 { + let is_reverse = next(2) == 1; + let read: Vec = (0..150) + .map(|_| if next(40) == 0 { 4 } else { next(4) as u8 }) + .collect(); + let del = if next(2) == 0 { + 1 + next(20) as i64 + } else { + 21 + next(1500) as i64 + }; + let r_gap = next(6) as i64; + let a = ScanArgs { + read_seq: &read, + r_a_end_inc: next(150) as i64, + g_a_end_inc: next(genome.n_genome) as i64, + g_b_start1: 0, + jr_start: -(next(40) as i32), + scan_end: (r_gap + next(80) as i64) as i32, + del, + motif_in_range: del >= scorer.align_intron_min as i64 + && del <= scorer.align_intron_max as i64, + is_reverse, + n_genome: genome.n_genome, + genome_offset: if is_reverse { genome.n_genome } else { 0 }, + seq, + }; + let a = ScanArgs { + g_b_start1: a.g_a_end_inc + del, + ..a + }; + if let Some(buffered) = scorer.scan_phase2_buffered(&a) { + assert_eq!(buffered, scorer.scan_phase2_scalar(&a)); + compared += 1; + } + } + assert!( + compared > 15_000, + "only {compared} scans took the buffered path" + ); + } + #[test] fn test_detect_gtag_motif() { // Sequence layout (0-based): From aeba6414e60b2af82ad569ac524a14eca465b4d1 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 30 Oct 2026 14:33:07 +0100 Subject: [PATCH 16/17] Stitcher: keep the window memo tables between windows The extension memo added 4 to 6 MB of private memory on read sets with many large windows, though no table exceeded 4,191 entries (about 500 KB). Creating and freeing a table per window is the likely cause, through heap fragmentation. The tables now live for the thread and are cleared between windows, and an extension entry takes 40 bytes instead of 64. As the extension table now outlives the window, it stops storing at 14,336 entries. Fusion PE, 150k pairs, 8 threads: 51 -> 39 MB of private memory (44 MB before the memo), with 0.7 % fewer instructions. Co-Authored-By: Claude Opus 5.5 --- src/align/score.rs | 19 +++++++++++++------ src/align/stitch.rs | 45 +++++++++++++++++++++++++++++++-------------- 2 files changed, 44 insertions(+), 20 deletions(-) diff --git a/src/align/score.rs b/src/align/score.rs index 7866b9a..36d3606 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -969,27 +969,34 @@ struct JunctionScanKey { /// Per-window memo table for the junction-position scan. /// -/// Create one per `stitch_seeds_core` call and pass it down the recursion; it -/// must not outlive the read, genome and strand it was filled for; the -/// per-window lifetime guarantees by construction. +/// Filled by one `stitch_seeds_core` call and passed down the recursion; it +/// must be cleared before another read, genome or strand uses it, which that +/// function does before handing the storage to the next window. #[derive(Default)] pub struct JunctionScanCache { map: rustc_hash::FxHashMap, /// Flank extensions of the same window, for `stitch::extend_alignment_cached`. - pub(crate) ext: rustc_hash::FxHashMap, + pub(crate) ext: rustc_hash::FxHashMap, /// Whether `ext` is used for this window. pub(crate) memo_extensions: bool, } /// What a flank extension depends on within one window: where it starts, which /// way and how far it may go, the transcript's mismatches and length so far, -/// and whether it must reach the read end. -pub(crate) type ExtendKey = (usize, u64, i32, usize, u32, usize, bool); +/// and whether it must reach the read end: genome start, read start, maximum +/// length, mismatches, then length with the direction and end flags packed in. +pub(crate) type ExtendKey = (u64, u32, u32, u32, u32); impl JunctionScanCache { pub fn new() -> Self { Self::default() } + + /// Forget every entry and keep the storage for the next window. + pub(crate) fn clear(&mut self) { + self.map.clear(); + self.ext.clear(); + } } /// Splice junction motif types diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 6147a1c..4fe809f 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -141,6 +141,19 @@ const SCORE_REGION_CHUNK: usize = 256; /// the extension costs more there than it saves. const MEMO_EXTENSIONS_MIN_SEEDS: usize = 25; +/// Most extensions memoized per window. The table is kept for the thread's +/// lifetime, so this bounds what it holds. 7/8 of a power of two, so the table +/// does not grow again on reaching it. Results past the limit are not stored. +const MEMO_EXTENSIONS_MAX: usize = (1 << 14) / 8 * 7; + +thread_local! { + /// The memo tables of the window being stitched, kept between windows: + /// allocating them per window fragments the heap, which costs more memory + /// than the tables themselves. + static JUNCTION_SCAN_CACHE: std::cell::Cell = + std::cell::Cell::default(); +} + /// [`extend_alignment`] memoized in the window's cache. Within one window the /// read, genome strand and mismatch limits are fixed, so the [`ExtendKey`] /// fields decide the result; the recursion's leaves repeat many of them. @@ -160,15 +173,6 @@ fn extend_alignment_cached( is_reverse: bool, extend_to_end: bool, ) -> ExtendResult { - let key = ( - read_start, - genome_start, - direction, - max_extend, - n_mm_prev, - len_prev, - extend_to_end, - ); if !jcache.memo_extensions { return extend_alignment( read_seq, @@ -185,9 +189,18 @@ fn extend_alignment_cached( extend_to_end, ); } + // Lengths are bounded by the read's, so they fit; the two flags take the low + // bits of the last field to keep the entry small. + let key = ( + genome_start, + read_start as u32, + max_extend as u32, + n_mm_prev, + (len_prev as u32) << 2 | u32::from(direction > 0) << 1 | u32::from(extend_to_end), + ); if let Some(&(extend_len, max_score, n_mismatch)) = jcache.ext.get(&key) { return ExtendResult { - extend_len, + extend_len: extend_len as usize, max_score, n_mismatch, }; @@ -206,9 +219,11 @@ fn extend_alignment_cached( is_reverse, extend_to_end, ); - jcache - .ext - .insert(key, (r.extend_len, r.max_score, r.n_mismatch)); + if jcache.ext.len() < MEMO_EXTENSIONS_MAX { + jcache + .ext + .insert(key, (r.extend_len as u32, r.max_score, r.n_mismatch)); + } r } @@ -3628,7 +3643,7 @@ pub(crate) fn stitch_seeds_core( // One memo table per window. `stitch_read`, the genome and the strand are // fixed for the whole recursion below, which is what makes the six-field // key in `JunctionScanCache` a complete identifier for a scan. - let mut jcache = crate::align::score::JunctionScanCache::new(); + let mut jcache = JUNCTION_SCAN_CACHE.take(); // Memoizing extensions pays only where the recursion has many leaves. jcache.memo_extensions = wa_entries.len() >= MEMO_EXTENSIONS_MIN_SEEDS; @@ -3650,6 +3665,8 @@ pub(crate) fn stitch_seeds_core( &mut jcache, debug_read_name, ); + jcache.clear(); + JUNCTION_SCAN_CACHE.set(jcache); if debug { eprintln!( From bcddd3164b76330abc3e7a05b499aaa1044e85e8 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 00:15:53 +0100 Subject: [PATCH 17/17] Remove dead code; stdout writers are the file writers on another output Removes, each verified unused in production: - the stitcher's scoreSeedBest block, ending in `let _ = best_pre_score`, and the `pre_ext_score` field it filled - the reverse-strand branches of the stitcher and the junction scan: `stitch_cluster.is_reverse` is always false, so `is_reverse` / `n_genome` / `genome_offset` go from those functions - BufferedBamRecords, SamWriter::write_alignment, Seed::find_paired_seeds, Seed::get_genome_positions, combined_transcript_for_projection, parse_gtf, extract_junctions_from_exons, score_gap_with_strand, word_comp_length - the transcriptome/sjdb file tail duplicated in index/mod.rs SamStdoutWriter, BamStdoutWriter and SortedBamStdoutWriter are replaced by SamWriter, BamWriter and SortedBamWriter::create_stdout. --outStd SAM gains the file writer's [SAM-MISMATCH] check. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 12 ++ examples/jrbench.rs | 2 - src/align/read_align.rs | 93 +---------- src/align/score.rs | 181 +++++----------------- src/align/seed.rs | 123 --------------- src/align/stitch.rs | 313 +++++--------------------------------- src/index/mod.rs | 136 ++++++++--------- src/index/packed_array.rs | 8 - src/io/bam.rs | 188 ++++++++--------------- src/io/sam.rs | 167 +++----------------- src/junction/gtf.rs | 41 +++-- src/lib.rs | 57 +------ 12 files changed, 265 insertions(+), 1056 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b18a76e..5d85a6f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -48,6 +48,18 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Other changes +- Dead code removed: the per-seed pre-extension scoring block in the stitcher, + which ran for every window and whose result was discarded; the reverse-strand + branches of the stitcher and junction scan, unreachable since windows are + stitched in forward coordinates; and `BufferedBamRecords`, + `SamWriter::write_alignment`, `Seed::find_paired_seeds`, + `Seed::get_genome_positions`, `combined_transcript_for_projection`, + `parse_gtf`, `extract_junctions_from_exons`, + `AlignmentScorer::score_gap_with_strand`. **API removal.** The stdout writers + are now the file writers on another output (`SamWriter`, `BamWriter`, + `SortedBamWriter::create_stdout`), so `--outStd SAM` gets the file writer's + `[SAM-MISMATCH]` check. Output is unchanged. + - Third-party notices completed. `LICENSES/` gains the licences of code ported from DegNorm (MIT) and David Elworthy's `sgt.h` via STAR (BSD 2-clause), and of libc++ (Apache 2.0 with LLVM Exceptions), whose random diff --git a/examples/jrbench.rs b/examples/jrbench.rs index d987c9b..db285e5 100644 --- a/examples/jrbench.rs +++ b/examples/jrbench.rs @@ -77,8 +77,6 @@ fn main() { 0, del as i64, &genome, - false, - N as u64, 75, iters_per_scan, ); diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 35e450a..382eddc 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -5,7 +5,7 @@ use crate::align::stitch::{ PE_SPACER_BASE, cluster_seeds, finalize_transcript, split_combined_wt, stitch_seeds_core, stitch_seeds_with_jdb_debug, }; -use crate::align::transcript::{Exon, Transcript}; +use crate::align::transcript::Transcript; use crate::error::Error; use crate::index::GenomeIndex; use crate::params::{MultimapperOrder, Parameters}; @@ -178,51 +178,6 @@ impl PairedAlignment { let (first, last) = if m1.is_reverse { (m2, m1) } else { (m1, m2) }; last.genome_end.saturating_sub(first.genome_start) } - - /// Build a STAR-style combined two-mate `Transcript` for transcriptome - /// projection. - /// - /// Matches STAR's single-`Transcript`-per-pair model: mate1 exons with - /// `i_frag = 0`, then mate2 exons rewritten to `i_frag = 1`. Only - /// meaningful for pairs on the same chromosome and strand — both are - /// invariants of a `PairedAlignment` (checked in `try_pair_transcripts`). - /// - /// The returned transcript's `cigar` is empty: transcriptome BAM - /// emission generates per-mate CIGARs from the split exon list rather - /// than consuming a combined one. - pub fn combined_transcript_for_projection(&self) -> Transcript { - let m1 = &self.mate1_transcript; - let m2 = &self.mate2_transcript; - - let mut exons: Vec = Vec::with_capacity(m1.exons.len() + m2.exons.len()); - for e in &m1.exons { - let mut ee = e.clone(); - ee.i_frag = 0; - exons.push(ee); - } - for e in &m2.exons { - let mut ee = e.clone(); - ee.i_frag = 1; - exons.push(ee); - } - - Transcript { - output_rank: 0, - chr_idx: m1.chr_idx, - genome_start: m1.genome_start.min(m2.genome_start), - genome_end: m1.genome_end.max(m2.genome_end), - is_reverse: m1.is_reverse, - exons, - cigar: Vec::new(), - score: m1.score + m2.score, - n_mismatch: m1.n_mismatch + m2.n_mismatch, - n_gap: m1.n_gap + m2.n_gap, - n_junction: m1.n_junction + m2.n_junction, - junction_motifs: Vec::new(), - junction_annotated: Vec::new(), - junction_strand: vec![], - } - } } /// Result of paired-end alignment, covering all mapping outcomes. @@ -1635,6 +1590,7 @@ pub(crate) fn pe_junctions_consistent(first: &Transcript, second: &Transcript) - mod tests { use super::*; use crate::align::score::SpliceMotif; + use crate::align::transcript::Exon; use crate::params::IntronMotifFilter; use crate::genome::Genome; @@ -1702,51 +1658,6 @@ mod tests { } } - #[test] - fn combined_transcript_for_projection_rewrites_mate2_ifrag() { - use cigar::op::{Kind, Op}; - let make_tr = |gs: u64, ge: u64, rs: usize, re: usize| Transcript { - output_rank: 0, - chr_idx: 0, - genome_start: gs, - genome_end: ge, - is_reverse: false, - exons: vec![Exon { - genome_start: gs, - genome_end: ge, - read_start: rs, - read_end: re, - i_frag: 0, - }], - cigar: vec![Op::new(Kind::Match, (ge - gs) as usize)], - score: 100, - n_mismatch: 0, - n_gap: 0, - n_junction: 0, - junction_motifs: vec![], - junction_annotated: vec![], - junction_strand: vec![], - }; - let pair = PairedAlignment { - output_rank: 0, - mate1_transcript: make_tr(1000, 1100, 0, 100), - mate2_transcript: make_tr(1300, 1400, 0, 100), - mate1_region: (0, 100), - mate2_region: (0, 100), - is_proper_pair: true, - insert_size: 400, - combined_wt_score: 200, - }; - assert_eq!(pair.star_g_length(), 400); - let combined = pair.combined_transcript_for_projection(); - assert_eq!(combined.exons.len(), 2); - assert_eq!(combined.exons[0].i_frag, 0); - assert_eq!(combined.exons[1].i_frag, 1); - assert_eq!(combined.genome_start, 1000); - assert_eq!(combined.genome_end, 1400); - assert_eq!(combined.score, 200); - } - /// Overlapping mates: STAR's `gLength` runs from the first mate's start to /// the second mate's end in transcript order, not over the pair's extent, /// and it decides which of two equal-scoring pairs is primary. The numbers diff --git a/src/align/score.rs b/src/align/score.rs index 36d3606..cd429c0 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -247,23 +247,6 @@ impl AlignmentScorer { read_gap: i64, genome_pos: u64, genome: &Genome, - ) -> (i32, GapType) { - self.score_gap_with_strand(genome_gap, read_gap, genome_pos, genome, false, 0) - } - - /// Score a gap between consecutive seeds, with strand-aware motif detection. - /// - /// For reverse-strand reads, `genome_pos` is a raw SA position in the RC genome. - /// The donor position must be converted to forward genome coordinates for motif - /// detection: `forward_donor = n_genome - rc_donor - intron_len`. - pub fn score_gap_with_strand( - &self, - genome_gap: i64, - read_gap: i64, - genome_pos: u64, - genome: &Genome, - is_reverse: bool, - n_genome: u64, ) -> (i32, GapType) { match (genome_gap, read_gap) { // Insertion: read advances but genome doesn't @@ -277,12 +260,7 @@ impl AlignmentScorer { let len = gg as u32; if len >= self.align_intron_min && len <= self.align_intron_max { // Splice junction — detect motif on forward genome - let donor = if is_reverse { - n_genome - genome_pos - len as u64 - } else { - genome_pos - }; - let motif = self.detect_splice_motif(donor, len, genome); + let motif = self.detect_splice_motif(genome_pos, len, genome); let score = self.score_splice_junction(motif); ( score, @@ -305,12 +283,7 @@ impl AlignmentScorer { 1.. => { let del_len = excess as u32; if del_len >= self.align_intron_min && del_len <= self.align_intron_max { - let rc_donor = genome_pos + rg as u64; - let donor = if is_reverse { - n_genome - rc_donor - del_len as u64 - } else { - rc_donor - }; + let donor = genome_pos + rg as u64; let motif = self.detect_splice_motif(donor, del_len, genome); let score = self.score_splice_junction(motif); ( @@ -344,8 +317,7 @@ impl AlignmentScorer { /// Memoized wrapper around [`AlignmentScorer::find_best_junction_position`]. /// /// The scan is a pure function of its arguments. Within one window's stitch - /// recursion, `read_seq`, the genome, `is_reverse` and `n_genome` are fixed, - /// so the six remaining coordinates identify a scan completely. + /// recursion, `read_seq` and the genome are fixed, so the six remaining coordinates identify a scan completely. /// `stitchWindowAligns`' include/exclude recursion reaches the same /// (exon A end, seed B) pair through many different branch paths, so without /// a memo the identical scan is repeated thousands of times per window. @@ -360,8 +332,6 @@ impl AlignmentScorer { r_gap: i64, g_gap: i64, genome: &Genome, - is_reverse: bool, - n_genome: u64, prev_exon_len: usize, next_seed_len: usize, ) -> (i32, SpliceMotif, i32, u32, u32) { @@ -383,8 +353,6 @@ impl AlignmentScorer { r_gap, g_gap, genome, - is_reverse, - n_genome, prev_exon_len, next_seed_len, ); @@ -399,7 +367,8 @@ impl AlignmentScorer { /// score how well the read matches the upstream vs downstream genome, plus /// the splice motif quality. Return the shift that maximizes the combined score. /// - /// All coordinates are in SA coordinate space (forward for fwd reads, RC genome for rev). + /// Coordinates are forward-genome: `stitch_seeds_core` reverse-complements + /// the read of a reverse-strand window before stitching. /// /// Returns: (jr_shift, best_motif, best_motif_score) /// - jr_shift: how many bases to shift the junction boundary (positive = rightward) @@ -414,8 +383,6 @@ impl AlignmentScorer { r_gap: i64, // read gap between seeds g_gap: i64, // genome gap between seeds genome: &Genome, - is_reverse: bool, - n_genome: u64, prev_exon_len: usize, next_seed_len: usize, // length of the B seed (for scan range) ) -> (i32, SpliceMotif, i32, u32, u32) { @@ -433,7 +400,6 @@ impl AlignmentScorer { // With inclusive coords: gBstart1 = g_a_end_inc + del let g_b_start1 = g_a_end_inc as i64 + del; - let genome_offset: u64 = if is_reverse { n_genome } else { 0 }; // Resolve the genome storage once: the three scans below read a base per // iteration, and `Genome::get_base` re-checks the `GenomeSeq` variant on // every one of them. @@ -457,8 +423,8 @@ impl AlignmentScorer { break; } - let g_up = seq.base((g_up_pos as u64 + genome_offset) as usize); - let g_dn = seq.base((g_dn_pos as u64 + genome_offset) as usize); + let g_up = seq.base(g_up_pos as usize); + let g_dn = seq.base(g_dn_pos as usize); if g_up >= 4 || g_dn >= 4 { break; @@ -491,9 +457,6 @@ impl AlignmentScorer { scan_end: r_gap as i32 + next_seed_len as i32, del, motif_in_range, - is_reverse, - n_genome, - genome_offset, seq, }; let (mut best_jr, mut best_motif, mut best_motif_score) = self @@ -510,8 +473,8 @@ impl AlignmentScorer { if left_pos < 0 || right_pos < 0 { break; } - let gl = seq.base((left_pos as u64 + genome_offset) as usize); - let gr = seq.base((right_pos as u64 + genome_offset) as usize); + let gl = seq.base(left_pos as usize); + let gr = seq.base(right_pos as usize); if gl < 4 && gl == gr { jj_l += 1; } else { @@ -530,8 +493,8 @@ impl AlignmentScorer { if left_pos < 0 || right_pos < 0 { break; } - let gl = seq.base((left_pos as u64 + genome_offset) as usize); - let gr = seq.base((right_pos as u64 + genome_offset) as usize); + let gl = seq.base(left_pos as usize); + let gr = seq.base(right_pos as usize); if gl < 4 && gl == gr { jj_r += 1; } else { @@ -544,33 +507,17 @@ impl AlignmentScorer { // Non-canonical or deletion: flush to be deterministic in repeat regions. // STAR operates in forward-genome space, so flush-left = leftmost in forward coords. - // Our function operates in SA space (RC genome for reverse reads), so: - // - Forward reads: flush LEFT in SA space = LEFT in forward space ✓ - // - Reverse reads: flush RIGHT in SA space = LEFT in forward space ✓ if best_motif == SpliceMotif::NonCanonical || del < self.align_intron_min as i64 { - if is_reverse { - // Reverse strand: flush RIGHT in RC space = LEFT in forward space - best_jr += jj_r; - jj_l += jj_r; - jj_r = 0; - } else { - // Forward strand: flush LEFT (STAR default) - best_jr -= jj_l; - jj_r += jj_l; - jj_l = 0; - } + best_jr -= jj_l; + jj_r += jj_l; + jj_l = 0; // STAR: `if (int(EX_L)+jR<1) return -1000005;` The caller rejects the // stitch when the flush leaves no base in the exon that shrinks // (clamping here would build a junction STAR never considers). // Re-check motif at flushed position if del >= self.align_intron_min as i64 && del <= self.align_intron_max as i64 { let donor_sa = (g_a_end_inc as i64 + best_jr as i64 + 1) as u64; - let donor_fwd = if is_reverse { - n_genome - donor_sa - del as u64 - } else { - donor_sa - }; - best_motif = self.detect_splice_motif(donor_fwd, del as u32, genome); + best_motif = self.detect_splice_motif(donor_sa, del as u32, genome); best_motif_score = self.score_splice_junction(best_motif); } } @@ -621,10 +568,8 @@ impl AlignmentScorer { } let mut up = [OUT_OF_RANGE; SCAN_BUF]; let mut dn = [OUT_OF_RANGE; SCAN_BUF]; - scan.seq - .bases_into(up0 as usize + scan.genome_offset as usize, &mut up[..n]); - scan.seq - .bases_into(dn0 as usize + scan.genome_offset as usize, &mut dn[..n]); + scan.seq.bases_into(up0 as usize, &mut up[..n]); + scan.seq.bases_into(dn0 as usize, &mut dn[..n]); // score1 change at each position: +1 where the read matches upstream // only, -1 where it matches downstream only, 0 otherwise or where // either genome base is not A/C/G/T (the sentinel included). @@ -637,22 +582,12 @@ impl AlignmentScorer { delta[k] = plus as i8 - minus as i8; } - // Donor and acceptor runs, forward strand. Forward reads move the - // donor right one base a position; reverse reads move it left, so - // their run starts at the last position's donor and is read backwards. + // Donor and acceptor runs: the donor moves right one base a position. let mut don = [OUT_OF_RANGE; SCAN_BUF]; let mut acc = [OUT_OF_RANGE; SCAN_BUF]; let del = scan.del as u64; if scan.motif_in_range { - let lo = if scan.is_reverse { - let last = donor_sa0 as u64 + (n as u64 - 1) + del; - if last > scan.n_genome { - return None; - } - scan.n_genome - last - } else { - donor_sa0 as u64 - }; + let lo = donor_sa0 as u64; scan.seq.bases_into(lo as usize, &mut don[..=n]); scan.seq.bases_into((lo + del - 2) as usize, &mut acc[..=n]); } @@ -667,13 +602,12 @@ impl AlignmentScorer { let jr1 = scan.jr_start + k as i32; score1 += i32::from(step); if scan.motif_in_range { - let i = if scan.is_reverse { n - 1 - k } else { k }; let motif = MotifWindow { donor: 0, - d1: don[i], - d2: don[i + 1], - a1: acc[i], - a2: acc[i + 1], + d1: don[k], + d2: don[k + 1], + a1: acc[k], + a2: acc[k + 1], } .motif(); let motif_score = self.score_splice_junction(motif); @@ -684,7 +618,7 @@ impl AlignmentScorer { best_motif = motif; best_motif_score = motif_score; } - } else if score1 > max_score2 || (scan.is_reverse && score1 == max_score2) { + } else if score1 > max_score2 { max_score2 = score1; best_jr = jr1; best_motif = SpliceMotif::NonCanonical; @@ -693,8 +627,7 @@ impl AlignmentScorer { let bound = score1 + (scan.scan_end - 1 - jr1) + if scan.motif_in_range { max_motif } else { 0 }; - let tie_wins = !scan.motif_in_range && scan.is_reverse; - if bound < max_score2 || (bound == max_score2 && !tie_wins) { + if bound <= max_score2 { break; } } @@ -727,12 +660,8 @@ impl AlignmentScorer { let g_dn_pos = scan.g_b_start1 + jr1 as i64; if g_up_pos >= 0 && g_dn_pos >= 0 { - let gu = scan - .seq - .base((g_up_pos as u64 + scan.genome_offset) as usize); - let gd = scan - .seq - .base((g_dn_pos as u64 + scan.genome_offset) as usize); + let gu = scan.seq.base(g_up_pos as usize); + let gd = scan.seq.base(g_dn_pos as usize); if gu < 4 && gd < 4 { if read_base == gu && read_base != gd { @@ -748,18 +677,12 @@ impl AlignmentScorer { if scan.motif_in_range { // Donor position in SA space: one past the last donor-exon base let donor_sa = (scan.g_a_end_inc + jr1 as i64 + 1) as u64; - // Convert to forward genome coordinates for motif detection - let donor_fwd = if scan.is_reverse { - scan.n_genome - donor_sa - scan.del as u64 - } else { - donor_sa - }; let w = match window.as_mut() { Some(w) => { - w.slide_to(donor_fwd, scan.del as u64, scan.seq); + w.slide_to(donor_sa, scan.del as u64, scan.seq); &*w } - None => window.insert(MotifWindow::at(donor_fwd, scan.del as u64, scan.seq)), + None => window.insert(MotifWindow::at(donor_sa, scan.del as u64, scan.seq)), }; let motif = w.motif(); let motif_score = self.score_splice_junction(motif); @@ -774,8 +697,7 @@ impl AlignmentScorer { } else { // Deletion (Del < alignIntronMin) or out-of-range gap — no motif, // pure positional score. STAR: jCan1=-1, jPen1=0, Score2=Score1. - // For reverse strand, rightmost-in-RC wins ties (= leftmost in forward). - if score1 > max_score2 || (scan.is_reverse && score1 == max_score2) { + if score1 > max_score2 { max_score2 = score1; best_jr = jr1; best_motif = SpliceMotif::NonCanonical; @@ -785,12 +707,11 @@ impl AlignmentScorer { // `score1` rises by at most 1 a position, so once even that cannot // lift a later position past the best one the rest of the scan is - // settled. Ties win only in the reverse-strand deletion case. + // settled: a tie does not move the best position. let scan_end = scan.scan_end; let bound = score1 + (scan_end - 1 - jr1) + if scan.motif_in_range { max_motif } else { 0 }; - let tie_wins = !scan.motif_in_range && scan.is_reverse; - if bound < max_score2 || (bound == max_score2 && !tie_wins) { + if bound <= max_score2 { break; } @@ -854,9 +775,6 @@ struct ScanArgs<'a> { scan_end: i32, del: i64, motif_in_range: bool, - is_reverse: bool, - n_genome: u64, - genome_offset: u64, seq: SeqView<'a>, } @@ -954,9 +872,8 @@ const fn build_motif_table() -> [SpliceMotif; 256] { /// Key for [`JunctionScanCache`]: the arguments of /// `find_best_junction_position` that vary within a single window's stitch -/// recursion. Everything else (`read_seq`, the genome, `is_reverse`, -/// `n_genome`) is loop-invariant there, so these six fields identify a scan -/// exactly. +/// recursion. Everything else (`read_seq`, the genome) is loop-invariant +/// there, so these six fields identify a scan exactly. #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] struct JunctionScanKey { r_a_end: usize, @@ -1172,7 +1089,6 @@ mod tests { let scorer = AlignmentScorer::from_params_minimal(); let mut compared = 0; for _ in 0..20_000 { - let is_reverse = next(2) == 1; let read: Vec = (0..150) .map(|_| if next(40) == 0 { 4 } else { next(4) as u8 }) .collect(); @@ -1192,9 +1108,6 @@ mod tests { del, motif_in_range: del >= scorer.align_intron_min as i64 && del <= scorer.align_intron_max as i64, - is_reverse, - n_genome: genome.n_genome, - genome_offset: if is_reverse { genome.n_genome } else { 0 }, seq, }; let a = ScanArgs { @@ -1911,15 +1824,11 @@ mod tests { let read_seq = vec![0, 1, 2, 3, 0, 0, 0, 1]; // ACGTAACG... let (jr_shift, motif, _motif_score, _, _) = scorer.find_best_junction_position( - &read_seq, - 5, // r_a_end (exclusive) + &read_seq, 5, // r_a_end (exclusive) 5, // g_a_end (exclusive) 1, // read_gap 25, // genome_gap - &genome, - false, // forward strand - genome.n_genome, - 5, // prev_exon_len + &genome, 5, // prev_exon_len 3, // next_seed_len ); @@ -1955,16 +1864,11 @@ mod tests { let read_seq = vec![0, 1, 2, 3, 0, 0, 1]; // ACGTAAC let (jr_shift, motif, _motif_score, _, _) = scorer.find_best_junction_position( - &read_seq, - 5, // r_a_end (exclusive) + &read_seq, 5, // r_a_end (exclusive) 5, // g_a_end (exclusive) 0, // read_gap = 0 (pure splice) 24, // genome_gap = intron length - &genome, - false, - genome.n_genome, - 5, - 2, // next_seed_len + &genome, 5, 2, // next_seed_len ); assert_eq!(motif, SpliceMotif::GtAg); @@ -2006,16 +1910,11 @@ mod tests { let read_seq = vec![0, 1, 2, 3, 1, 1, 1]; // the gap base let (jr_shift, motif, _, _, _) = scorer.find_best_junction_position( - &read_seq, - 5, // r_a_end + &read_seq, 5, // r_a_end 5, // g_a_end 0, // read_gap 24, // genome_gap - &genome, - false, - genome.n_genome, - 5, - 2, // next_seed_len + &genome, 5, 2, // next_seed_len ); // Non-canonical junction should be flushed left diff --git a/src/align/seed.rs b/src/align/seed.rs index dfcc025..cc417e7 100644 --- a/src/align/seed.rs +++ b/src/align/seed.rs @@ -190,65 +190,6 @@ impl Seed { Ok(seeds) } - /// Find all seeds for paired-end reads using unified seed pooling. - /// - /// This implements STAR's hybrid approach: seeds from both mates are found - /// independently, tagged with their mate origin, then pooled together for - /// unified clustering. - /// - /// # Arguments - /// * `mate1_seq` - First mate sequence (encoded) - /// * `mate2_seq` - Second mate sequence (encoded) - /// * `index` - Genome index with SA and SAindex - /// * `min_seed_length` - Minimum seed length to report - /// * `params` - Parameters including seedMultimapNmax - /// - /// # Returns - /// Vector of seeds from both mates, tagged with mate_id (0 or 1) - pub fn find_paired_seeds( - mate1_seq: &[u8], - mate2_seq: &[u8], - index: &GenomeIndex, - min_seed_length: usize, - params: &Parameters, - ) -> Result, Error> { - // Find seeds from mate1 (tag with mate_id = 0) - let lread = mate1_seq.len() + 1 + mate2_seq.len(); - let mut seeds = - Self::find_seeds_at(mate1_seq, 0, lread, index, min_seed_length, params, "")?; - for seed in &mut seeds { - seed.mate_id = 0; - } - - // Find seeds from mate2 (tag with mate_id = 1) - // IMPORTANT: read_pos is relative to mate2 start (will be adjusted during stitching) - let mut seeds2 = Self::find_seeds_at( - mate2_seq, - mate1_seq.len() + 1, - lread, - index, - min_seed_length, - params, - "", - )?; - for seed in &mut seeds2 { - seed.mate_id = 1; - } - - // Pool seeds together - seeds.extend(seeds2); - - check_seed_per_read(seeds.len(), params)?; - Ok(seeds) - } - - /// Get all genome positions for this seed. - /// - /// Expands the SA range to actual genome positions. - pub fn get_genome_positions(&self, index: &GenomeIndex) -> Vec<(u64, bool)> { - self.genome_positions(index).collect() - } - /// Iterate over genome positions for this seed without allocating. /// /// Returns an iterator that lazily decodes SA entries. @@ -1151,25 +1092,6 @@ mod tests { assert!(seeds.is_empty()); } - #[test] - fn get_genome_positions() { - let index = make_test_index("ACGTACGT"); - let read = encode_sequence("ACGT"); - let params = params(&[]); - - let seeds = Seed::find_seeds(&read, &index, 3, ¶ms, "").unwrap(); - assert!(!seeds.is_empty()); - - // Get positions for first seed - let positions = seeds[0].get_genome_positions(&index); - assert_ne!(positions, []); - - // Should have at least one valid position - for (pos, _is_reverse) in positions { - assert!(pos < index.genome.n_genome); - } - } - #[test] fn test_single_end_mate_id() { let index = make_test_index("ACGTACGT"); @@ -1185,51 +1107,6 @@ mod tests { } } - #[test] - fn test_find_paired_seeds() { - let index = make_test_index("ACGTACGTTTGGCCAA"); - let mate1 = encode_sequence("ACGT"); - let mate2 = encode_sequence("TTGG"); - let params = params(&[]); - - let seeds = Seed::find_paired_seeds(&mate1, &mate2, &index, 3, ¶ms).unwrap(); - - // Should have seeds from both mates - let mate1_seeds: Vec<_> = seeds.iter().filter(|s| s.mate_id == 0).collect(); - let mate2_seeds: Vec<_> = seeds.iter().filter(|s| s.mate_id == 1).collect(); - - assert!(!mate1_seeds.is_empty(), "Should have mate1 seeds"); - assert!(!mate2_seeds.is_empty(), "Should have mate2 seeds"); - - // Verify mate1 seeds have correct read positions - for seed in mate1_seeds { - assert!(seed.read_pos < mate1.len()); - } - - // Verify mate2 seeds have correct read positions (relative to mate2) - for seed in mate2_seeds { - assert!(seed.read_pos < mate2.len()); - } - } - - #[test] - fn test_paired_seeds_pooling() { - let index = make_test_index("ACGTACGT"); - let mate1 = encode_sequence("ACGT"); - let mate2 = encode_sequence("ACGT"); - let params = params(&[]); - - let seeds = Seed::find_paired_seeds(&mate1, &mate2, &index, 3, ¶ms).unwrap(); - - // Should have roughly double the seeds (one set from each mate) - let mate1_count = seeds.iter().filter(|s| s.mate_id == 0).count(); - let mate2_count = seeds.iter().filter(|s| s.mate_id == 1).count(); - - assert!(mate1_count > 0); - assert!(mate2_count > 0); - assert_eq!(seeds.len(), mate1_count + mate2_count); - } - #[test] fn test_reverse_complement_read() { // ACGT → RC = ACGT (palindrome) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 4fe809f..0d4a319 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -16,17 +16,15 @@ const EXTEND_TO_END_KILL: i32 = i32::MIN / 2; /// Count mismatches in an alignment by comparing read sequence to genome sequence. /// -/// The read sequence is always in forward orientation. For reverse-strand alignments, -/// the genome is accessed at `pos + n_genome` (the reverse-complement region) rather -/// than reverse-complementing the read. +/// Stitching runs in forward genome coordinates, reverse-strand reads being +/// reverse-complemented first (`stitch_seeds_core`). /// /// # Arguments /// * `read_seq` - Read sequence in forward orientation (encoded as 0=A, 1=C, 2=G, 3=T, 4=N) /// * `cigar_ops` - CIGAR operations -/// * `genome_start` - Starting position in genome (decoded SA position, WITHOUT n_genome offset) +/// * `genome_start` - Starting position in forward genome coordinates /// * `read_start` - Starting position in read /// * `index` - Genome index (contains genome sequence) -/// * `is_reverse` - Whether this is a reverse-strand alignment /// /// # Returns /// Number of mismatched bases (excluding N bases) @@ -37,9 +35,7 @@ fn count_mismatches_in_region( genome_start: u64, length: usize, index: &GenomeIndex, - is_reverse: bool, ) -> u32 { - let genome_offset = if is_reverse { index.genome.n_genome } else { 0 }; let mut n_mismatch = 0u32; for i in 0..length { @@ -48,9 +44,7 @@ fn count_mismatches_in_region( break; } let read_base = read_seq[read_pos]; - if let Some(genome_base) = index - .genome - .get_base(genome_start + i as u64 + genome_offset) + if let Some(genome_base) = index.genome.get_base(genome_start + i as u64) && read_base != genome_base && read_base != 4 && genome_base != 4 @@ -74,9 +68,7 @@ fn score_region( genome_start: u64, length: usize, index: &GenomeIndex, - is_reverse: bool, ) -> (i32, u32) { - let genome_offset = if is_reverse { index.genome.n_genome } else { 0 }; let seq = index.genome.sequence.view(); // Bound the run once instead of testing the read end per base, then walk it @@ -93,10 +85,7 @@ fn score_region( let mut done = 0usize; while done < run { let want = (run - done).min(SCORE_REGION_CHUNK); - let got = seq.bases_into( - (genome_start + (done + genome_offset as usize) as u64) as usize, - &mut gbuf[..want], - ); + let got = seq.bases_into((genome_start + done as u64) as usize, &mut gbuf[..want]); if got == 0 { break; } @@ -170,7 +159,6 @@ fn extend_alignment_cached( n_mm_max: u32, p_mm_max: f64, index: &GenomeIndex, - is_reverse: bool, extend_to_end: bool, ) -> ExtendResult { if !jcache.memo_extensions { @@ -185,7 +173,6 @@ fn extend_alignment_cached( n_mm_max, p_mm_max, index, - is_reverse, extend_to_end, ); } @@ -216,7 +203,6 @@ fn extend_alignment_cached( n_mm_max, p_mm_max, index, - is_reverse, extend_to_end, ); if jcache.ext.len() < MEMO_EXTENSIONS_MAX { @@ -247,7 +233,7 @@ struct ExtendResult { /// # Arguments /// * `read_seq` - Full read sequence (encoded) /// * `read_start` - Boundary position in read (where extension begins) -/// * `genome_start` - Corresponding genome position (WITHOUT n_genome offset) +/// * `genome_start` - Corresponding position in forward genome coordinates /// * `direction` - +1 for rightward extension, -1 for leftward /// * `max_extend` - Maximum distance to extend (to read boundary) /// * `n_mm_prev` - Mismatches already in the aligned portion @@ -255,7 +241,6 @@ struct ExtendResult { /// * `n_mm_max` - outFilterMismatchNmax (absolute max mismatches) /// * `p_mm_max` - outFilterMismatchNoverLmax (max mismatch ratio) /// * `index` - Genome index -/// * `is_reverse` - Whether this is a reverse-strand alignment #[allow(clippy::too_many_arguments)] fn extend_alignment( read_seq: &[u8], @@ -268,7 +253,6 @@ fn extend_alignment( n_mm_max: u32, p_mm_max: f64, index: &GenomeIndex, - is_reverse: bool, extend_to_end: bool, ) -> ExtendResult { if max_extend == 0 { @@ -279,7 +263,6 @@ fn extend_alignment( }; } - let genome_offset = if is_reverse { index.genome.n_genome } else { 0 }; // Resolve the genome storage once: both extension loops below read one base // per iteration, and `Genome::get_base` re-checks the `GenomeSeq` variant on // every call. @@ -323,7 +306,7 @@ fn extend_alignment( } genome_start - 1 - i as u64 }; - let genome_base = seq.base((genome_pos + genome_offset) as usize); + let genome_base = seq.base(genome_pos as usize); if genome_base == crate::genome::OUT_OF_RANGE { return ExtendResult { extend_len: 0, @@ -398,7 +381,7 @@ fn extend_alignment( }; // Get genome base (with strand offset) - let genome_base = seq.base((genome_pos + genome_offset) as usize); + let genome_base = seq.base(genome_pos as usize); if genome_base == crate::genome::OUT_OF_RANGE { break; } @@ -479,7 +462,7 @@ pub struct WindowAlignment { pub length: usize, /// Forward-strand genome position (STAR: WA_gStart) pub genome_pos: u64, - /// Raw SA position (for genome base access in DP — reverse strand uses sa_pos + n_genome) + /// Raw SA position; `stitch_seeds_core` sets it to `genome_pos` on a reverse-strand window. pub sa_pos: u64, /// SA range size of the originating seed (STAR: WA_Nrep) — for scoring pub n_rep: usize, @@ -492,10 +475,6 @@ pub struct WindowAlignment { /// only dedup against seeds with the same value, so a junction half is /// never absorbed by a genomic seed on the same diagonal. pub sj_a: Option, - /// Pre-computed extension score estimate (STAR: scoreSeedBest[iS] base case). - /// = 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, } /// A cluster of seeds mapping to the same genomic region @@ -1023,7 +1002,6 @@ pub fn cluster_seeds( is_anchor: is_anchor_seed, mate_id: seed.mate_id, sj_a: cand.sj_a, - pre_ext_score: length as i32, }, ); } @@ -1094,7 +1072,6 @@ pub fn cluster_seeds( is_anchor: is_anchor_seed, mate_id: seed.mate_id, sj_a: cand.sj_a, - pre_ext_score: length as i32, }, ); } @@ -1310,26 +1287,7 @@ fn stitch_align_in_place( } // STAR condition (stitchAlignToTranscript.cpp:352): // gBstart + trA->exons[0][EX_R] + nBasesMax >= trA->exons[0][EX_G] || EX_G < EX_R - // For forward clusters: checked in forward genome space. - // For reverse clusters: checked in STAR's encoded genome space, then converted to - // rustar-aligner's forward-position representation (Phase 16.27 convention). - if cluster.is_reverse { - // Reverse cluster: stitch_read = RC(combined) = [mate2|SPACER|RC(mate1)]. - // wt.exons[0] is the first mate2 exon; wa is the new mate1 seed. - // Recover mate1's position in the original combined read [mate1|SPACER|RC(mate2)]: - // combined_start_mate1 = combined_len - stitch_pos_mate1 - len_mate1 - // STAR's encoded-space reject condition (converted to forward-position arithmetic): - // combined_start_mate1 < len_mate2_exon - len_mate1 + (P_mate2 - P_mate1) - let combined_start_mate1 = - (read_seq.len() as i64) - (wa.read_pos as i64) - (wa.length as i64); - let first_exon = &new_wt.exons[0]; - let len_mate2_exon = (first_exon.read_end - first_exon.read_start) as i64; - let len_mate1 = wa.length as i64; - let p_diff = first_exon.genome_start as i64 - wa.sa_pos as i64; // P_mate2 - P_mate1 - if combined_start_mate1 < len_mate2_exon - len_mate1 + p_diff { - return false; - } - } else { + { let first_exon = &new_wt.exons[0]; let ex_g = first_exon.genome_start; let ex_r = first_exon.read_start as u64; @@ -1365,7 +1323,6 @@ fn stitch_align_in_place( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, false, // internal stitch extension: alignEndsType applied at finalize ); if right_ext.extend_len > 0 { @@ -1426,7 +1383,6 @@ fn stitch_align_in_place( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, false, // internal stitch extension: alignEndsType applied at finalize ); if left_ext.extend_len > 0 { @@ -1534,8 +1490,7 @@ fn stitch_align_in_place( { return false; // STAR -1000006: repeats around a non-canonical junction } - let donor_fwd = index.sa_pos_to_forward(donor_sa, cluster.is_reverse, del as usize); - let motif = scorer.detect_splice_motif(donor_fwd, del, &index.genome); + let motif = scorer.detect_splice_motif(donor_sa, del, &index.genome); ( 0, motif, @@ -1552,8 +1507,6 @@ fn stitch_align_in_place( 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, ) @@ -1561,14 +1514,9 @@ fn stitch_align_in_place( // STAR -1000005 (`stitchAlignToTranscript.cpp:168-170`): the // non-canonical / deletion flush to the left must leave at least one - // base in the exon that shrinks (A in the forward scan, B in the - // reverse-strand scan space). + // base in the exon that shrinks. if sjdb_junction.is_none() { - let shrunk_len = if cluster.is_reverse { - (eff_length + shared) as i32 - jr_shift - } else { - (last_exon.read_end - last_exon.read_start) as i32 + jr_shift - }; + let shrunk_len = (last_exon.read_end - last_exon.read_start) as i32 + jr_shift; if shrunk_len < 1 { return false; } @@ -1595,7 +1543,6 @@ fn stitch_align_in_place( last_exon.genome_end, junction_offset, index, - cluster.is_reverse, ); shared_mm += mm; shared_score += s; @@ -1610,7 +1557,6 @@ fn stitch_align_in_place( last_exon.genome_end + junction_offset as u64 + del as u64, acceptor_bases, index, - cluster.is_reverse, ); shared_mm += mm; shared_score += s; @@ -1631,7 +1577,6 @@ fn stitch_align_in_place( extra_genome_start, n_extra, index, - cluster.is_reverse, ); shared_mm += extra_mm; shared_score -= 2 * extra_mm as i32; @@ -1649,7 +1594,6 @@ fn stitch_align_in_place( extra_genome_start, n_extra, index, - cluster.is_reverse, ); shared_mm += extra_mm; shared_score -= 2 * extra_mm as i32; @@ -1672,9 +1616,7 @@ fn stitch_align_in_place( // Annotated strand of a database hit (STAR `sjdbStrand`); the sjdb // entry found below takes precedence over the GTF-derived lookup. let db_strand = junction_db.and_then(|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 donor_fwd = (donor_sa as i64 + jr_shift as i64) as u64; let acceptor_fwd = donor_fwd + del as u64 - 1; [1u8, 2, 0] .into_iter() @@ -1687,13 +1629,9 @@ fn stitch_align_in_place( // non-canonical junction to the left, where sjdbStart keeps it; an // annotated one is then shifted back by sjdbShiftLeft to its true // coordinate (the score was taken at the flushed position, which - // is equivalent on the repeat). The shift is the same in the - // forward genome; in the reverse-strand scan space it runs the - // other way. + // is equivalent on the repeat). if is_annotated && sjdb_junction.is_none() { - 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 donor_fwd = (donor_sa as i64 + jr_shift as i64) as u64; let acceptor_fwd = donor_fwd + del as u64 - 1; let prepared = &index.prepared_junctions; let hit = prepared @@ -1714,13 +1652,8 @@ fn stitch_align_in_place( if eff_length <= sl || prev_len <= sl { return false; // STAR -1000006 } - let shifted = if cluster.is_reverse { - jr_shift - pj.shift_left as i32 - } else { - jr_shift + pj.shift_left as i32 - }; - // STAR: rAend+jR >= rBend, i.e. the shift runs past B. In the - // reverse scan space the roles of A and B are exchanged. + let shifted = jr_shift + pj.shift_left as i32; + // STAR: rAend+jR >= rBend, i.e. the shift runs past B. if shifted >= (shared + eff_length) as i32 || prev_len as i32 + shifted < 1 { return false; // STAR -1000006 } @@ -1789,19 +1722,13 @@ fn stitch_align_in_place( 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; // Phase 1: Scan shared bases to find best split point // STAR: for (jR1=1; jR1<=gGap; jR1++) for jr1 in 1..=shared_usize { - let g_pos = last_exon.genome_end + (jr1 - 1) as u64 + genome_offset; + let g_pos = last_exon.genome_end + (jr1 - 1) as u64; if let Some(gb) = index.genome.get_base(g_pos) && gb < 4 { @@ -1841,7 +1768,7 @@ fn stitch_align_in_place( // Post-insertion side (skip insertion gap) last_exon.read_end + ins + ii - 1 }; - let g_pos = last_exon.genome_end + (ii - 1) as u64 + genome_offset; + let g_pos = last_exon.genome_end + (ii - 1) as u64; if let Some(gb) = index.genome.get_base(g_pos) && gb < 4 @@ -2101,7 +2028,6 @@ pub(crate) fn finalize_transcript( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, extend_left_to_end, ) } else { @@ -2120,7 +2046,6 @@ pub(crate) fn finalize_transcript( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, extend_right_to_end, ) } else { @@ -2141,7 +2066,6 @@ pub(crate) fn finalize_transcript( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, extend_right_to_end, ) } else { @@ -2160,7 +2084,6 @@ pub(crate) fn finalize_transcript( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, extend_left_to_end, ) } else { @@ -2365,11 +2288,6 @@ pub(crate) fn finalize_transcript( // `nMM--`). `wt.n_mismatch` keeps that count, as nMM does. let n_mismatch = wt.n_mismatch + left_extend.n_mismatch + right_extend.n_mismatch; - // Reverse CIGAR for reverse strand - if cluster.is_reverse { - final_cigar.reverse(); - } - // Compute total reference-consuming length from CIGAR let mut ref_len = 0u64; for op in &final_cigar { @@ -2385,20 +2303,12 @@ pub(crate) fn finalize_transcript( } } - // Guard: reverse-strand alignment must not extend past genome end - if cluster.is_reverse && adjusted_genome_start + ref_len > index.genome.n_genome { - return None; - } - - // Convert to forward genome coordinates - let forward_genome_start = - index.sa_pos_to_forward(adjusted_genome_start, cluster.is_reverse, ref_len as usize); - let forward_genome_end = forward_genome_start + ref_len; + let forward_genome_end = adjusted_genome_start + ref_len; // Build exons from CIGAR let mut exons = Vec::new(); let mut read_pos_e = 0usize; - let mut genome_pos_e = forward_genome_start; + let mut genome_pos_e = adjusted_genome_start; for op in &final_cigar { match op.kind() { @@ -2442,7 +2352,7 @@ pub(crate) fn finalize_transcript( let t_genome_start = merged_exons .first() - .map_or(forward_genome_start, |e| e.genome_start); + .map_or(adjusted_genome_start, |e| e.genome_start); let t_genome_end = merged_exons .last() .map_or(forward_genome_end, |e| e.genome_end); @@ -2457,7 +2367,7 @@ pub(crate) fn finalize_transcript( chr_idx: cluster.chr_idx, genome_start: t_genome_start, genome_end: t_genome_end, - is_reverse: cluster.is_reverse, + is_reverse: false, exons: merged_exons, cigar: final_cigar, score: final_score, @@ -2779,7 +2689,6 @@ fn stitch_recurse( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, false, // internal stitch extension: alignEndsType applied at finalize ) } else { @@ -2800,7 +2709,6 @@ fn stitch_recurse( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, false, // internal stitch extension: alignEndsType applied at finalize ) } else { @@ -2827,7 +2735,6 @@ fn stitch_recurse( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, false, // internal stitch extension: alignEndsType applied at finalize ) } else { @@ -2849,7 +2756,6 @@ fn stitch_recurse( scorer.n_mm_max, scorer.p_mm_max, index, - cluster.is_reverse, false, // internal stitch extension: alignEndsType applied at finalize ) } else { @@ -3444,14 +3350,10 @@ pub(crate) fn stitch_seeds_core( } } - // Create cluster for stitching (is_reverse=false for reverse-strand, so stitcher + // Create cluster for stitching (is_reverse=false for every window, so the stitcher // uses forward genome coords without RC genome offset) let stitch_cluster = SeedCluster { - is_reverse: if stitch_is_reverse { - false - } else { - cluster.is_reverse - }, + is_reverse: false, ..cluster.clone() }; @@ -3497,143 +3399,6 @@ pub(crate) fn stitch_seeds_core( } } - // --- STAR stitchWindowSeeds.cpp: scoreSeedBest pre-extension --- - // - // Phase 1: Pre-extend each seed left and right (base case of scoreSeedBest DP). - // Uses stitch_read (forward strand, forward genome coords after coord conversion above). - // EXTEND_ORDER=1: for forward reads, left extension first (5' of read), then right. - // for reverse reads, right extension first (5' of read in RC), then left. - // This matches stitch_recurse's `do_left_first = !original_is_reverse`. - // Mismatch budget for second ext carries over from first ext (STAR: scoreSeedBestMM[iS1]). - { - let zero_ext = ExtendResult { - extend_len: 0, - max_score: 0, - n_mismatch: 0, - }; - let do_left_first = !stitch_is_reverse; - for wa in &mut wa_entries { - let right_start = wa.read_pos + wa.length; - - let (first_ext, second_ext) = if do_left_first { - // Forward cluster: left ext first, then right - let left = if wa.read_pos > 0 { - extend_alignment( - stitch_read, - wa.read_pos, - wa.sa_pos, - -1, - wa.read_pos, - 0, - wa.length, - scorer.n_mm_max, - scorer.p_mm_max, - index, - false, - false, // pre-ext scoring: alignEndsType applied at finalize - ) - } else { - zero_ext.clone() - }; - let right_len_prev = wa.length + left.extend_len; - let right = if right_start < stitch_read.len() { - extend_alignment( - stitch_read, - right_start, - wa.sa_pos + wa.length as u64, - 1, - stitch_read.len() - right_start, - left.n_mismatch, - right_len_prev, - scorer.n_mm_max, - scorer.p_mm_max, - index, - false, - false, // pre-ext scoring: alignEndsType applied at finalize - ) - } else { - zero_ext.clone() - }; - (left, right) - } else { - // Reverse cluster: right ext first (5' of RC read), then left - let right = if right_start < stitch_read.len() { - extend_alignment( - stitch_read, - right_start, - wa.sa_pos + wa.length as u64, - 1, - stitch_read.len() - right_start, - 0, - wa.length, - scorer.n_mm_max, - scorer.p_mm_max, - index, - false, - false, // pre-ext scoring: alignEndsType applied at finalize - ) - } else { - zero_ext.clone() - }; - let left_len_prev = wa.length + right.extend_len; - let left = if wa.read_pos > 0 { - extend_alignment( - stitch_read, - wa.read_pos, - wa.sa_pos, - -1, - wa.read_pos, - right.n_mismatch, - left_len_prev, - scorer.n_mm_max, - scorer.p_mm_max, - index, - false, - false, // pre-ext scoring: alignEndsType applied at finalize - ) - } else { - zero_ext.clone() - }; - (left, right) - }; - - wa.pre_ext_score = wa.length as i32 + first_ext.max_score + second_ext.max_score; - } - } - - // Phase 2: DP chain over pre-extended scores (STAR: scoreSeedBest chain accumulation). - // dp[i] = best chain score ending at wa_entries[i]. - // For intra-read gaps: if read_gap != genome_gap this is a potential splice — 0 penalty. - // For exact-match gaps (read_gap == genome_gap): score the gap as mismatches only - // (lenient: use 0 since exact-match gaps are handled by the actual DP later). - // Both cases produce an overestimate — conservative upper bound matching STAR's approach. - let best_pre_score = { - let mut dp: Vec = wa_entries.iter().map(|wa| wa.pre_ext_score).collect(); - for i in 1..wa_entries.len() { - for j in 0..i { - let rj_end = wa_entries[j].read_pos + wa_entries[j].length; - let gj_end = wa_entries[j].sa_pos + wa_entries[j].length as u64; - // Must not overlap in read or genome - if wa_entries[i].read_pos < rj_end || wa_entries[i].sa_pos < gj_end { - continue; - } - // Gap penalty: 0 for potential splices (read_gap != genome_gap) or short gaps - // This prevents underestimating spliced reads whose per-seed pre_ext_score is - // limited by the splice junction (extension stops at the intron boundary). - let chain = dp[j] + wa_entries[i].pre_ext_score; - if chain > dp[i] { - dp[i] = chain; - } - } - } - dp.into_iter().max().unwrap_or(0) - }; - - // Note: scoreSeedBest (best_pre_score) is used by STAR for seed ordering within - // stitchWindowAligns, NOT as a hard pre-filter gate. We keep pre_ext_score on each WA - // entry for future use in seed ordering (Phase B), but do not filter here. - let _ = best_pre_score; // suppress unused warning - // Run recursive include/exclude stitcher. // Note: STAR 2.7.11b does not force the last anchor to be included (its // WlastAnchor mechanism is dead code — see stitch_recurse), so there is no @@ -3912,7 +3677,6 @@ mod tests { is_anchor: true, mate_id: 2, sj_a: None, - pre_ext_score: 5, }, WindowAlignment { seed_idx: 1, @@ -3924,7 +3688,6 @@ mod tests { is_anchor: true, mate_id: 2, sj_a: None, - pre_ext_score: 5, }, ]; @@ -4087,7 +3850,6 @@ mod tests { is_anchor: false, mate_id: 2, sj_a: None, - pre_ext_score: 15, }; let mut jcache = crate::align::score::JunctionScanCache::new(); let out = stitch_align_to_transcript( @@ -4130,7 +3892,7 @@ mod tests { 0, // no previous length 10, // n_mm_max 0.3, // p_mm_max - &index, false, false, // forward strand, local + &index, false, // forward strand, local ); assert_eq!(result.extend_len, 5); @@ -4155,7 +3917,7 @@ mod tests { 0, // no previous length 10, // n_mm_max 0.3, // p_mm_max - &index, false, false, + &index, false, ); assert_eq!(result.extend_len, 5); @@ -4180,7 +3942,7 @@ mod tests { 0, // no previous length 10, // n_mm_max 0.3, // p_mm_max - &index, false, false, + &index, false, ); // Should extend 4 bases (perfect match), then mismatches drag score down @@ -4202,7 +3964,7 @@ mod tests { 0, // genome_start 1, // rightward 5, // max_extend - 0, 0, 10, 0.3, &index, false, false, + 0, 0, 10, 0.3, &index, false, ); // Should stop at 3 bases (genome boundary) @@ -4218,7 +3980,7 @@ mod tests { // Read: A A C G (matches at 0, N skip at 1, matches at 2-3) let read_seq: Vec = vec![0, 0, 1, 2]; - let result = extend_alignment(&read_seq, 0, 0, 1, 4, 0, 0, 10, 0.3, &index, false, false); + let result = extend_alignment(&read_seq, 0, 0, 1, 4, 0, 0, 10, 0.3, &index, false); // Should extend all 4 bases: match + N(skip) + match + match = score 3 assert_eq!(result.extend_len, 4); @@ -4234,7 +3996,7 @@ mod tests { // Read: T T T T (all 3, complete mismatch) let read_seq: Vec = vec![3, 3, 3, 3]; - let result = extend_alignment(&read_seq, 0, 0, 1, 4, 0, 0, 10, 0.3, &index, false, false); + let result = extend_alignment(&read_seq, 0, 0, 1, 4, 0, 0, 10, 0.3, &index, false); // Score never goes positive, so extend_len should be 0 assert_eq!(result.extend_len, 0); @@ -4250,7 +4012,7 @@ mod tests { // M M M X M M M M M M M M M M let read_seq: Vec = vec![0, 1, 2, 0, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1]; - let result = extend_alignment(&read_seq, 0, 0, 1, 14, 0, 0, 10, 0.3, &index, false, false); + let result = extend_alignment(&read_seq, 0, 0, 1, 14, 0, 0, 10, 0.3, &index, false); // Should extend past the mismatch: 3M + 1X + 10M // Score: +3 -1 +10 = 12, best at position 14 @@ -4267,7 +4029,7 @@ mod tests { let result = extend_alignment( &read_seq, 0, 0, 1, 0, // max_extend = 0 - 0, 0, 10, 0.3, &index, false, false, + 0, 0, 10, 0.3, &index, false, ); assert_eq!(result.extend_len, 0); @@ -4284,12 +4046,12 @@ mod tests { let read_seq: Vec = vec![0, 1, 2, 1, 1]; // Local: stops at the score peak (3M), soft-clipping the 2 mismatches. - let local = extend_alignment(&read_seq, 0, 0, 1, 5, 0, 0, 10, 0.3, &index, false, false); + let local = extend_alignment(&read_seq, 0, 0, 1, 5, 0, 0, 10, 0.3, &index, false); assert_eq!(local.extend_len, 3, "local extension stops at score peak"); assert_eq!(local.max_score, 3); // EndToEnd: forces all 5 bases; score = 3 matches - 2 mismatches = 1. - let e2e = extend_alignment(&read_seq, 0, 0, 1, 5, 0, 0, 10, 0.3, &index, false, true); + let e2e = extend_alignment(&read_seq, 0, 0, 1, 5, 0, 0, 10, 0.3, &index, true); assert_eq!(e2e.extend_len, 5, "end-to-end forces full read extension"); assert_eq!(e2e.max_score, 1); assert_eq!(e2e.n_mismatch, 2); @@ -4304,7 +4066,7 @@ mod tests { // EndToEnd cannot reach the read end without crossing the chr boundary, // so it returns the kill sentinel (extend_len 0, hugely negative score). - let e2e = extend_alignment(&read_seq, 0, 0, 1, 5, 0, 0, 10, 0.3, &index, false, true); + let e2e = extend_alignment(&read_seq, 0, 0, 1, 5, 0, 0, 10, 0.3, &index, true); assert_eq!(e2e.extend_len, 0); assert!( e2e.max_score < -1_000_000, @@ -4587,7 +4349,6 @@ mod tests { is_anchor: true, mate_id: 2, sj_a: None, - pre_ext_score: 2, }; let out = stitch_align_to_transcript( &wt, diff --git a/src/index/mod.rs b/src/index/mod.rs index badff9d..e42f3b8 100644 --- a/src/index/mod.rs +++ b/src/index/mod.rs @@ -231,35 +231,9 @@ impl GenomeIndex { // path; factored into a helper. update_genome_params_sa_size(&dir.join("genomeParameters.txt"), genome.n_genome, sa_size)?; - // Transcriptome + sjdb files. Matches the tail of - // `GenomeIndex::write` byte-for-byte. - if let Some(tr) = &transcriptome { - tr.write_transcript_info(dir)?; - tr.write_exon_info(dir)?; - tr.write_gene_info(dir)?; - tr.write_exon_ge_tr_info(dir)?; - tr.write_sjdb_list_from_gtf(dir, &genome)?; - log::info!( - "Wrote transcriptome index files: {} transcripts, {} genes", - tr.n_transcripts(), - tr.gene_ids.len() - ); - } - if !prepared_junctions.is_empty() { - sjdb_insert::write_sjdb_info_tab( - &dir.join("sjdbInfo.txt"), - &prepared_junctions, - params.sjdb_overhang, - )?; - sjdb_insert::write_sjdb_list_out_tab( - &dir.join("sjdbList.out.tab"), - &prepared_junctions, - &genome, - )?; - log::info!("Wrote sjdb files: {} junctions", prepared_junctions.len()); - } - - Ok(()) + // Same tail as `GenomeIndex::write`. + write_transcript_tables(dir, transcriptome.as_ref(), &genome)?; + write_sjdb_files(dir, &prepared_junctions, &genome, params.sjdb_overhang) } /// Shared setup: load FASTA, parse GTF, build junctions, append @@ -413,8 +387,13 @@ impl GenomeIndex { /// Write index files to directory. pub fn write(&self, dir: &Path, params: &Parameters) -> Result<(), Error> { self.write_genome_files(dir, params)?; - self.write_transcript_tables(dir)?; - self.write_sjdb_files(dir, params.sjdb_overhang) + write_transcript_tables(dir, self.transcriptome.as_ref(), &self.genome)?; + write_sjdb_files( + dir, + &self.prepared_junctions, + &self.genome, + params.sjdb_overhang, + ) } /// Genome, SA, SAindex, the chromosome tables and `genomeParameters.txt`. @@ -440,46 +419,6 @@ impl GenomeIndex { ) } - /// Write transcriptome index files (STAR-compatible) when the GTF - /// was supplied. Matches STAR's `GTF_transcriptGeneSJ.cpp` outputs. - fn write_transcript_tables(&self, dir: &Path) -> Result<(), Error> { - if let Some(tr) = &self.transcriptome { - tr.write_transcript_info(dir)?; - tr.write_exon_info(dir)?; - tr.write_gene_info(dir)?; - tr.write_exon_ge_tr_info(dir)?; - tr.write_sjdb_list_from_gtf(dir, &self.genome)?; - log::info!( - "Wrote transcriptome index files: {} transcripts, {} genes", - tr.n_transcripts(), - tr.gene_ids.len() - ); - } - Ok(()) - } - - /// Write sjdbInfo.txt + sjdbList.out.tab, the files `sjdbPrepare` writes - /// whenever junctions are inserted into the genome. - fn write_sjdb_files(&self, dir: &Path, sjdb_overhang: u32) -> Result<(), Error> { - if !self.prepared_junctions.is_empty() { - sjdb_insert::write_sjdb_info_tab( - &dir.join("sjdbInfo.txt"), - &self.prepared_junctions, - sjdb_overhang, - )?; - sjdb_insert::write_sjdb_list_out_tab( - &dir.join("sjdbList.out.tab"), - &self.prepared_junctions, - &self.genome, - )?; - log::info!( - "Wrote sjdb files: {} junctions", - self.prepared_junctions.len() - ); - } - Ok(()) - } - /// STAR's run-time genome directory `_STARgenome/` /// (`Parameters.cpp:1019-1034`) after junctions were inserted at mapping: /// `sjdbPrepare` writes `sjdbInfo.txt` and `sjdbList.out.tab`, a GTF given @@ -494,9 +433,14 @@ impl GenomeIndex { let dir = params.output_path("_STARgenome"); fs::create_dir_all(&dir).map_err(|e| Error::io(e, &dir))?; if gtf_tables { - self.write_transcript_tables(&dir)?; + write_transcript_tables(&dir, self.transcriptome.as_ref(), &self.genome)?; } - self.write_sjdb_files(&dir, self.sjdb_overhang)?; + write_sjdb_files( + &dir, + &self.prepared_junctions, + &self.genome, + self.sjdb_overhang, + )?; if params.sjdb_insert_save == "All" { self.write_genome_files(&dir, &loaded_genome_params(params))?; } @@ -504,6 +448,52 @@ impl GenomeIndex { } } +/// Write transcriptome index files (STAR-compatible) when the GTF was supplied. +/// Matches STAR's `GTF_transcriptGeneSJ.cpp` outputs. +fn write_transcript_tables( + dir: &Path, + transcriptome: Option<&TranscriptomeIndex>, + genome: &Genome, +) -> Result<(), Error> { + if let Some(tr) = transcriptome { + tr.write_transcript_info(dir)?; + tr.write_exon_info(dir)?; + tr.write_gene_info(dir)?; + tr.write_exon_ge_tr_info(dir)?; + tr.write_sjdb_list_from_gtf(dir, genome)?; + log::info!( + "Wrote transcriptome index files: {} transcripts, {} genes", + tr.n_transcripts(), + tr.gene_ids.len() + ); + } + Ok(()) +} + +/// Write sjdbInfo.txt + sjdbList.out.tab, the files `sjdbPrepare` writes +/// whenever junctions are inserted into the genome. +fn write_sjdb_files( + dir: &Path, + prepared_junctions: &[PreparedJunction], + genome: &Genome, + sjdb_overhang: u32, +) -> Result<(), Error> { + if !prepared_junctions.is_empty() { + sjdb_insert::write_sjdb_info_tab( + &dir.join("sjdbInfo.txt"), + prepared_junctions, + sjdb_overhang, + )?; + sjdb_insert::write_sjdb_list_out_tab( + &dir.join("sjdbList.out.tab"), + prepared_junctions, + genome, + )?; + log::info!("Wrote sjdb files: {} junctions", prepared_junctions.len()); + } + Ok(()) +} + /// `params` as STAR holds them after loading the genome in `--genomeDir`: /// `genomeSAindexNbases`, `genomeChrBinNbits` and `genomeSAsparseD` are taken /// from its `genomeParameters.txt` (`Genome_genomeLoad.cpp:104-107`), and diff --git a/src/index/packed_array.rs b/src/index/packed_array.rs index 7f6b578..d6f2128 100644 --- a/src/index/packed_array.rs +++ b/src/index/packed_array.rs @@ -43,10 +43,6 @@ pub struct PackedArray { /// Number of bits per element word_length: u32, - /// Complement bits (64 - word_length) — kept for STAR compatibility - #[allow(dead_code)] - word_comp_length: u32, - /// Mask for extracting an element (word_length bits set) bit_rec_mask: u64, @@ -66,7 +62,6 @@ impl PackedArray { pub fn new(word_length: u32, length: usize) -> Self { assert!(word_length > 0 && word_length <= 64); - let word_comp_length = 64 - word_length; let bit_rec_mask = if word_length == 64 { u64::MAX } else { @@ -84,7 +79,6 @@ impl PackedArray { Self { word_length, - word_comp_length, bit_rec_mask, length, data, @@ -245,7 +239,6 @@ impl PackedArray { fn from_store(word_length: u32, length: usize, data: PackedBytes) -> Self { assert!(word_length > 0 && word_length <= 64); - let word_comp_length = 64 - word_length; let bit_rec_mask = if word_length == 64 { u64::MAX } else { @@ -254,7 +247,6 @@ impl PackedArray { Self { word_length, - word_comp_length, bit_rec_mask, length, data, diff --git a/src/io/bam.rs b/src/io/bam.rs index be8b040..2e55b4b 100644 --- a/src/io/bam.rs +++ b/src/io/bam.rs @@ -14,26 +14,6 @@ use std::fs::File; use std::io::{BufReader, BufWriter, Write}; use std::path::Path; -/// Buffer for BAM records built by parallel threads -#[derive(Default)] -pub struct BufferedBamRecords { - pub records: Vec, -} - -impl BufferedBamRecords { - /// Create new buffer with capacity - pub fn new() -> Self { - Self { - records: Vec::with_capacity(10000), - } - } - - /// Add a record to the buffer - pub fn push(&mut self, record: RecordBuf) { - self.records.push(record); - } -} - /// Convert STAR's `--outBAMcompression` integer (-1..9) to a noodles level. /// /// STAR mapping: -1 or 0 = uncompressed, 1-9 = deflate levels, default 1. @@ -69,16 +49,13 @@ fn bgzf_threads(params: &Parameters) -> usize { params.run_thread_n.get() } -type BgzfFile = BgzfWriter>; -type BgzfStdout = BgzfWriter>; - /// BAM file writer (streaming, unsorted) /// /// This writer streams BAM records directly to disk as they're generated, /// without buffering or sorting. The output is BGZF-compressed but unsorted. /// Users can sort the output with `samtools sort` if needed. -pub struct BamWriter { - writer: bam::io::Writer, +pub struct BamWriter> { + writer: bam::io::Writer>, header: sam::Header, } @@ -454,11 +431,12 @@ fn read_run_record( /// The header emits `SO:coordinate`. Unmapped records sort to the end. pub struct SortedBamWriter { sorter: CoordinateSorter, - output_path: std::path::PathBuf, + /// `None`: stdout. A file is only created by `finish`. + output_path: Option, } -impl BamWriter { - /// Create a BAM writer at `output_path` with the given prepared header. +impl BamWriter { + /// Create a BAM writer on `out` with the given prepared header. /// /// Uses a lenient header writer that bypasses noodles' SAM-spec /// reference-name validator. STAR accepts `(` / `)` in reference names @@ -466,25 +444,24 @@ impl BamWriter { /// validator rejects. Writing the SAM text portion of the BAM header /// manually sidesteps that validation while preserving every other byte. fn with_header( - output_path: &Path, + out: W, header: sam::Header, header_text: &[u8], compression: i32, threads: usize, ) -> Result { - let buf_writer = BufWriter::new(File::create(output_path)?); - let mut bgzf = make_bgzf_writer(buf_writer, compression, threads)?; + let mut bgzf = make_bgzf_writer(out, compression, threads)?; write_bam_header_lenient(&mut bgzf, &header, header_text)?; let writer = bam::io::Writer::from(bgzf); Ok(Self { writer, header }) } - /// Create a new BAM writer with header from genome index. - pub fn create(output_path: &Path, genome: &Genome, params: &Parameters) -> Result { + /// Create a BAM writer on `out`, with header from genome index. + pub fn with_output(out: W, genome: &Genome, params: &Parameters) -> Result { let header = crate::io::sam::build_sam_header(genome, params)?; let text = crate::io::sam::render_header_text(&header, params, HeaderKind::Main)?; Self::with_header( - output_path, + out, header, &text, params.out_bam_compression, @@ -492,6 +469,31 @@ impl BamWriter { ) } + /// Write batch of buffered records (for parallel processing) + /// + /// # Arguments + /// * `batch` - Slice of records to write + pub fn write_batch(&mut self, batch: &[RecordBuf]) -> Result<(), Error> { + for record in batch { + self.writer.write_alignment_record(&self.header, record)?; + } + Ok(()) + } + + /// Flush and close BAM file + pub fn finish(&mut self) -> Result<(), Error> { + self.writer.get_mut().finish()?; + log::info!("BAM file written successfully"); + Ok(()) + } +} + +impl BamWriter { + /// Create a new BAM file writer with header from genome index. + pub fn create(output_path: &Path, genome: &Genome, params: &Parameters) -> Result { + Self::with_output(BufWriter::new(File::create(output_path)?), genome, params) + } + /// Create a BAM writer whose @SQ header lists the transcripts (not /// chromosomes) from `tr_idx`. Used for /// `Aligned.toTranscriptome.out.bam`. @@ -508,31 +510,13 @@ impl BamWriter { let header = crate::io::sam::build_sam_header_from_refs(refs, params)?; let text = crate::io::sam::render_header_text(&header, params, HeaderKind::Transcriptome)?; Self::with_header( - output_path, + BufWriter::new(File::create(output_path)?), header, &text, params.out_bam_compression, bgzf_threads(params), ) } - - /// Write batch of buffered records (for parallel processing) - /// - /// # Arguments - /// * `batch` - Slice of records to write - pub fn write_batch(&mut self, batch: &[RecordBuf]) -> Result<(), Error> { - for record in batch { - self.writer.write_alignment_record(&self.header, record)?; - } - Ok(()) - } - - /// Flush and close BAM file - pub fn finish(&mut self) -> Result<(), Error> { - self.writer.get_mut().finish()?; - log::info!("BAM file written successfully"); - Ok(()) - } } impl SortedBamWriter { @@ -542,12 +526,28 @@ impl SortedBamWriter { output_path: &Path, genome: &crate::genome::Genome, params: &Parameters, + ) -> Result { + Self::with_output(Some(output_path.to_path_buf()), genome, params) + } + + /// [`create`](Self::create) writing to stdout on `finish`. + pub fn create_stdout( + genome: &crate::genome::Genome, + params: &Parameters, + ) -> Result { + Self::with_output(None, genome, params) + } + + fn with_output( + output_path: Option, + genome: &crate::genome::Genome, + params: &Parameters, ) -> Result { let header = crate::io::sam::build_sam_header(genome, params)?; let text = crate::io::sam::render_header_text(&header, params, HeaderKind::MainSorted)?; Ok(Self { sorter: CoordinateSorter::new(header, text, params), - output_path: output_path.to_path_buf(), + output_path, }) } @@ -573,13 +573,14 @@ impl SortedBamWriter { &mut self, finalize: Option<&mut dyn FnMut(&mut RecordBuf)>, ) -> Result<(), Error> { - let file = File::create(&self.output_path) - .map_err(|source| Error::io(source, &self.output_path))?; - self.sorter.write_sorted( - BufWriter::new(file), - &self.output_path.display().to_string(), - finalize, - ) + let Some(path) = &self.output_path else { + return self + .sorter + .write_sorted(BufWriter::new(std::io::stdout()), "stdout", finalize); + }; + let file = File::create(path).map_err(|source| Error::io(source, path))?; + self.sorter + .write_sorted(BufWriter::new(file), &path.display().to_string(), finalize) } } @@ -642,71 +643,6 @@ fn write_bam_header_lenient( Ok(()) } -/// Streaming unsorted BAM writer that writes to stdout. -pub struct BamStdoutWriter { - writer: bam::io::Writer, - header: sam::Header, -} - -impl BamStdoutWriter { - pub fn create(genome: &crate::genome::Genome, params: &Parameters) -> Result { - let header = crate::io::sam::build_sam_header(genome, params)?; - let header_text = crate::io::sam::render_header_text(&header, params, HeaderKind::Main)?; - let mut bgzf = make_bgzf_writer( - BufWriter::new(std::io::stdout()), - params.out_bam_compression, - bgzf_threads(params), - )?; - write_bam_header_lenient(&mut bgzf, &header, &header_text)?; - let writer = bam::io::Writer::from(bgzf); - Ok(Self { writer, header }) - } - - pub fn write_batch(&mut self, batch: &[RecordBuf]) -> Result<(), Error> { - for record in batch { - self.writer.write_alignment_record(&self.header, record)?; - } - Ok(()) - } - - pub fn finish(&mut self) -> Result<(), Error> { - self.writer.get_mut().finish()?; - Ok(()) - } -} - -/// Coordinate-sorted BAM writer that writes to stdout on `finish()`, with the -/// same bounded-memory spill/merge behavior as [`SortedBamWriter`]. -pub struct SortedBamStdoutWriter { - sorter: CoordinateSorter, -} - -impl SortedBamStdoutWriter { - pub fn create(genome: &crate::genome::Genome, params: &Parameters) -> Result { - let header = crate::io::sam::build_sam_header(genome, params)?; - let text = crate::io::sam::render_header_text(&header, params, HeaderKind::MainSorted)?; - Ok(Self { - sorter: CoordinateSorter::new(header, text, params), - }) - } - - pub fn write_batch(&mut self, batch: &[RecordBuf]) -> Result<(), Error> { - self.sorter.push_batch(batch) - } - - pub fn finish(&mut self) -> Result<(), Error> { - self.sorter - .write_sorted(BufWriter::new(std::io::stdout()), "stdout", None) - } - - /// [`finish`](Self::finish), passing each record through `finalize` as it - /// is written (see [`SortedBamWriter::finish_with`]). - pub fn finish_with(&mut self, finalize: &mut dyn FnMut(&mut RecordBuf)) -> Result<(), Error> { - self.sorter - .write_sorted(BufWriter::new(std::io::stdout()), "stdout", Some(finalize)) - } -} - #[cfg(test)] mod tests { use super::*; diff --git a/src/io/sam.rs b/src/io/sam.rs index 898bf33..378844d 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -24,7 +24,7 @@ use noodles::sam::header::record::value::{ use std::collections::HashSet; use std::fmt::Write as FmtWrite; use std::fs::File; -use std::io::{BufWriter, Write as _}; +use std::io::BufWriter; use std::num::NonZeroUsize; use std::path::Path; @@ -133,24 +133,21 @@ pub(crate) fn insert_mate_cigar_tags(rec1: &mut RecordBuf, rec2: &mut RecordBuf) } /// SAM file writer -pub struct SamWriter { - writer: sam::io::Writer>, +pub struct SamWriter> { + writer: sam::io::Writer, header: sam::Header, } -impl SamWriter { - /// Create a new SAM writer with header from genome index +impl SamWriter { + /// Create a SAM writer on `out`, with header from genome index /// /// # Arguments - /// * `output_path` - Path to output SAM file + /// * `out` - Output (a buffered file or stdout) /// * `genome` - Genome index with chromosome information /// * `params` - Parameters (for @PG header) - pub fn create(output_path: &Path, genome: &Genome, params: &Parameters) -> Result { - let file = File::create(output_path)?; - let buf_writer = BufWriter::new(file); - + pub fn with_output(out: W, genome: &Genome, params: &Parameters) -> Result { let header = build_sam_header(genome, params)?; - let mut writer = sam::io::Writer::new(buf_writer); + let mut writer = sam::io::Writer::new(out); writer .get_mut() @@ -159,77 +156,6 @@ impl SamWriter { Ok(Self { writer, header }) } - /// Write alignment record(s) for a read - /// - /// # Arguments - /// * `read_name` - Read identifier - /// * `read_seq` - Read sequence (encoded) - /// * `read_qual` - Quality scores - /// * `transcripts` - Alignment transcripts (1 or more for multi-mappers) - /// * `genome` - Genome index - /// * `params` - Parameters - #[allow(clippy::too_many_arguments)] - pub fn write_alignment( - &mut self, - read_name: &str, - read_seq: &[u8], - read_qual: &[u8], - transcripts: &[Transcript], - genome: &Genome, - params: &Parameters, - n_for_mapq: usize, - ) -> Result<(), Error> { - if transcripts.is_empty() { - return Ok(()); - } - - let n_alignments = transcripts.len(); - let max_output = if params.out_sam_mult_nmax < 0 { - n_alignments - } else { - (params.out_sam_mult_nmax as usize).min(n_alignments) - }; - let effective_n = n_alignments.max(n_for_mapq); - let mapq = calculate_mapq(effective_n, params.out_sam_mapq_unique); - let attrs = params.out_sam_attributes; - let rg_id_owned = params.primary_rg_id()?; - let rg_id = rg_id_owned.as_deref(); - - let best_score = transcripts - .iter() - .map(|t| t.score) - .max() - .unwrap_or(i32::MIN); - // STAR writes the alignments in its output order (`output_rank`); the primary - // (`transcripts[0]`) is not necessarily the first. - for (hit_index, (is_primary, transcript)) in output_order(transcripts, |t| t.output_rank) - .take(max_output) - .enumerate() - { - let mut record = transcript_to_record( - transcript, - read_name, - read_seq, - read_qual, - genome, - mapq, - max_output, // NH = number of reported alignments - hit_index + 1, // 1-based - params.out_sam_attr_ih_start, - attrs, - )?; - set_secondary(&mut record, !is_primary); - maybe_insert_rg_tag(&mut record, rg_id); - order_record_tags(&mut record, &attrs); - apply_sam_flag_or_and(&mut record, params); - apply_primary_flag(&mut record, transcript.score, best_score, params); - - self.writer.write_alignment_record(&self.header, &record)?; - } - - Ok(()) - } - /// Write batch of buffered records (for parallel processing) /// /// # Arguments @@ -262,6 +188,13 @@ impl SamWriter { } Ok(()) } +} + +impl SamWriter { + /// Create a new SAM file writer. + pub fn create(output_path: &Path, genome: &Genome, params: &Parameters) -> Result { + Self::with_output(BufWriter::new(File::create(output_path)?), genome, params) + } /// Build unmapped record (without writing) /// @@ -933,30 +866,6 @@ impl SamWriter { } } -/// SAM writer that streams to stdout. -pub struct SamStdoutWriter { - writer: sam::io::Writer>, - header: sam::Header, -} - -impl SamStdoutWriter { - pub fn create(genome: &Genome, params: &Parameters) -> Result { - let header = build_sam_header(genome, params)?; - let mut writer = sam::io::Writer::new(BufWriter::new(std::io::stdout())); - writer - .get_mut() - .write_all(&render_header_text(&header, params, HeaderKind::Main)?)?; - Ok(Self { writer, header }) - } - - pub fn write_batch(&mut self, batch: &[RecordBuf]) -> Result<(), Error> { - for record in batch { - self.writer.write_alignment_record(&self.header, record)?; - } - Ok(()) - } -} - /// Build paired SAM header from genome pub fn build_sam_header(genome: &Genome, params: &Parameters) -> Result { build_sam_header_from_refs( @@ -2304,8 +2213,6 @@ mod tests { #[test] fn test_sam_output_includes_rg_header_and_tag() { - use crate::align::transcript::Exon; - use cigar::op::{Kind, Op}; use std::io::Read; let genome = make_test_genome(); @@ -2321,40 +2228,16 @@ mod tests { let tmpfile = NamedTempFile::new().unwrap(); let mut writer = SamWriter::create(tmpfile.path(), &genome, ¶ms).unwrap(); - let transcript = Transcript { - output_rank: 0, - chr_idx: 0, - genome_start: 0, - genome_end: 4, - is_reverse: false, - exons: vec![Exon { - genome_start: 0, - genome_end: 4, - read_start: 0, - read_end: 4, - i_frag: 0, - }], - cigar: vec![Op::new(Kind::Match, 4)], - score: 100, - n_mismatch: 0, - n_gap: 0, - n_junction: 0, - junction_motifs: vec![], - junction_annotated: vec![], - junction_strand: vec![], - }; - - writer - .write_alignment( - "read1", - &[0, 1, 2, 3], - &[30, 30, 30, 30], - &[transcript], - &genome, - ¶ms, - 1, - ) - .unwrap(); + let record = SamWriter::build_unmapped_record( + "read1", + &[0, 1, 2, 3], + &[30, 30, 30, 30], + ¶ms, + UnmappedReason::Other, + BestTr::default(), + ) + .unwrap(); + writer.write_batch(&[record]).unwrap(); drop(writer); let mut contents = String::new(); diff --git a/src/junction/gtf.rs b/src/junction/gtf.rs index ef40ce6..7109ec7 100644 --- a/src/junction/gtf.rs +++ b/src/junction/gtf.rs @@ -97,12 +97,6 @@ pub fn parse_gtf_configured( Ok(exons) } -/// Parse GTF file and extract exon features (default: feature `"exon"`, no chr prefix). -#[allow(dead_code)] -pub fn parse_gtf(path: &Path) -> Result, Error> { - parse_gtf_configured(path, "exon", "") -} - /// Parse a single GTF line fn parse_gtf_line(line: &str) -> Result { let fields: Vec<&str> = line.split('\t').collect(); @@ -264,15 +258,6 @@ pub fn extract_junctions_configured( Ok(junctions) } -/// Extract junctions grouping by `"transcript_id"` (default, backward-compatible wrapper). -#[allow(dead_code)] -pub fn extract_junctions_from_exons( - exons: Vec, - genome: &Genome, -) -> Result, Error> { - extract_junctions_configured(exons, genome, "transcript_id") -} - #[cfg(test)] mod tests { use super::*; @@ -307,8 +292,16 @@ mod tests { std::fs::write(&gz, bytes).unwrap(); let key = |r: &GtfRecord| (r.seqname.clone(), r.start, r.end, r.strand); - let want: Vec<_> = parse_gtf(&plain).unwrap().iter().map(key).collect(); - let got: Vec<_> = parse_gtf(&gz).unwrap().iter().map(key).collect(); + let want: Vec<_> = parse_gtf_configured(&plain, "exon", "") + .unwrap() + .iter() + .map(key) + .collect(); + let got: Vec<_> = parse_gtf_configured(&gz, "exon", "") + .unwrap() + .iter() + .map(key) + .collect(); assert_eq!(want.len(), 2); assert_eq!(got, want); assert_eq!(read_annotation_to_string(&gz).unwrap(), gtf); @@ -361,7 +354,7 @@ mod tests { ) .unwrap(); - let exons = parse_gtf(file.path()).unwrap(); + let exons = parse_gtf_configured(file.path(), "exon", "").unwrap(); assert_eq!(exons.len(), 1); } @@ -375,7 +368,7 @@ mod tests { ) .unwrap(); - let exons = parse_gtf(file.path()).unwrap(); + let exons = parse_gtf_configured(file.path(), "exon", "").unwrap(); assert_eq!(exons.len(), 1); // Only exon, not gene assert_eq!(exons[0].feature, "exon"); } @@ -424,7 +417,7 @@ mod tests { }, ]; - let junctions = extract_junctions_from_exons(exons, &genome).unwrap(); + let junctions = extract_junctions_configured(exons, &genome, "transcript_id").unwrap(); assert_eq!(junctions.len(), 1); let (chr_idx, start, end, strand) = junctions[0]; @@ -505,7 +498,7 @@ mod tests { }, ]; - let junctions = extract_junctions_from_exons(exons, &genome).unwrap(); + let junctions = extract_junctions_configured(exons, &genome, "transcript_id").unwrap(); assert_eq!(junctions.len(), 1); assert_eq!(junctions[0], (0, 200, 298, 1)); @@ -539,7 +532,7 @@ mod tests { .collect(), }]; - let junctions = extract_junctions_from_exons(exons, &genome).unwrap(); + let junctions = extract_junctions_configured(exons, &genome, "transcript_id").unwrap(); assert_eq!(junctions.len(), 0); // No junctions } @@ -587,7 +580,7 @@ mod tests { }, ]; - let junctions = extract_junctions_from_exons(exons, &genome).unwrap(); + let junctions = extract_junctions_configured(exons, &genome, "transcript_id").unwrap(); assert_eq!(junctions.len(), 0); // Skipped due to unknown chr } @@ -636,7 +629,7 @@ mod tests { }, ]; - let junctions = extract_junctions_from_exons(exons, &genome).unwrap(); + let junctions = extract_junctions_configured(exons, &genome, "transcript_id").unwrap(); assert_eq!(junctions.len(), 1); let (_chr_idx, start, end, _strand) = junctions[0]; diff --git a/src/lib.rs b/src/lib.rs index 1d0c993..7878bc4 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -291,7 +291,7 @@ impl AlignmentWriter for NullWriter { } } -impl AlignmentWriter for crate::io::sam::SamWriter { +impl AlignmentWriter for crate::io::sam::SamWriter { fn write_batch( &mut self, batch: &[noodles::sam::alignment::record_buf::RecordBuf], @@ -300,7 +300,7 @@ impl AlignmentWriter for crate::io::sam::SamWriter { } } -impl AlignmentWriter for crate::io::bam::BamWriter { +impl AlignmentWriter for crate::io::bam::BamWriter { fn write_batch( &mut self, batch: &[noodles::sam::alignment::record_buf::RecordBuf], @@ -333,48 +333,6 @@ impl AlignmentWriter for crate::io::bam::SortedBamWriter { } } -impl AlignmentWriter for crate::io::sam::SamStdoutWriter { - fn write_batch( - &mut self, - batch: &[noodles::sam::alignment::record_buf::RecordBuf], - ) -> Result<(), error::Error> { - self.write_batch(batch) - } -} - -impl AlignmentWriter for crate::io::bam::BamStdoutWriter { - fn write_batch( - &mut self, - batch: &[noodles::sam::alignment::record_buf::RecordBuf], - ) -> Result<(), error::Error> { - self.write_batch(batch) - } - fn finish(&mut self) -> Result<(), error::Error> { - self.finish() - } -} - -impl AlignmentWriter for crate::io::bam::SortedBamStdoutWriter { - fn write_batch( - &mut self, - batch: &[noodles::sam::alignment::record_buf::RecordBuf], - ) -> Result<(), error::Error> { - self.write_batch(batch) - } - fn finish(&mut self) -> Result<(), error::Error> { - self.finish() - } - fn finish_with( - &mut self, - finalize: &mut dyn FnMut(&mut noodles::sam::alignment::record_buf::RecordBuf), - ) -> Result<(), error::Error> { - self.finish_with(finalize) - } - fn holds_records(&self) -> bool { - true - } -} - fn align_reads(params: &Parameters) -> anyhow::Result<()> { use crate::index::GenomeIndex; @@ -941,7 +899,8 @@ fn run_single_pass( OutSamFormat::Sam => { if params.out_std == OutStd::Sam { info!("Writing SAM to stdout (--outStd SAM)"); - Box::new(crate::io::sam::SamStdoutWriter::create( + Box::new(SamWriter::with_output( + std::io::BufWriter::new(std::io::stdout()), &index.genome, params, )?) @@ -959,7 +918,8 @@ fn run_single_pass( if out_type.bam_unsorted { outputs.push(if params.out_std == OutStd::BamUnsorted { info!("Writing unsorted BAM to stdout (--outStd BAM_Unsorted)"); - Box::new(crate::io::bam::BamStdoutWriter::create( + Box::new(BamWriter::with_output( + std::io::BufWriter::new(std::io::stdout()), &index.genome, params, )?) @@ -977,10 +937,7 @@ fn run_single_pass( info!( "Writing coordinate-sorted BAM to stdout (--outStd BAM_SortedByCoordinate)" ); - Box::new(crate::io::bam::SortedBamStdoutWriter::create( - &index.genome, - params, - )?) + Box::new(SortedBamWriter::create_stdout(&index.genome, params)?) } else { let output_path = params.output_path("Aligned.sortedByCoord.out.bam"); info!("Writing BAM to {}", output_path.display());