From 16b7df3355cbc8e7b39dd42f25059cf6e39bbbed Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 24 Jul 2026 09:31:13 +0200 Subject: [PATCH 1/5] feat: STARlong binary + long-read window-coverage filter Add rustar-aligner-long binary (STAR's STARlong: identical CLI, forces the chaining-DP long-read stitcher via Parameters::long_read, matching how native STAR selects it by binary rather than a flag). Wire params.long_read through stitch_seeds_with_jdb_debug into a new window_read_coverage (ReadAlign_stitchPieces.cpp's WC accumulation, faithful to upstream's rLast=0 under-count for a seed starting at read position 0) and gate the winReadCoverageRelativeMin/winReadCoverageBasesMin cluster filter on long_read, since standard STAR does not filter clusters by seed coverage. Refactor stitch_seeds_core's preamble into prepare_stitch_input, shared with the (upcoming) STARlong chaining-DP stitcher stitch_seeds_long_read. Co-Authored-By: Claude Sonnet 5 --- Cargo.toml | 6 + src/align/read_align.rs | 103 +++++++- src/align/stitch.rs | 437 +++++++++++++++++++++++++++++++-- src/bin/rustar_aligner_long.rs | 19 ++ src/params/mod.rs | 10 + 5 files changed, 553 insertions(+), 22 deletions(-) create mode 100644 src/bin/rustar_aligner_long.rs diff --git a/Cargo.toml b/Cargo.toml index 8a4638f..23e4b06 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -25,6 +25,12 @@ exclude = [ name = "rustar-aligner" path = "src/main.rs" +[[bin]] +# STAR's `STARlong` (`-DCOMPILE_FOR_LONG_READS`): identical CLI, but forces the +# chaining-DP long-read stitcher (see params::Parameters::long_read). +name = "rustar-aligner-long" +path = "src/bin/rustar_aligner_long.rs" + [dependencies] clap = { version = "4", features = ["derive"] } anyhow = "1" diff --git a/src/align/read_align.rs b/src/align/read_align.rs index b4d2f20..91ab2ec 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -2,8 +2,8 @@ use crate::align::score::{AlignmentScorer, SpliceMotif}; use crate::align::seed::Seed; use crate::align::stitch::{ - PE_SPACER_BASE, cluster_seeds, finalize_transcript, split_combined_wt, stitch_seeds_core, - stitch_seeds_with_jdb_debug, + PE_SPACER_BASE, SeedCluster, cluster_seeds, finalize_transcript, split_combined_wt, + stitch_seeds_core, stitch_seeds_with_jdb_debug, }; use crate::align::transcript::{Exon, Transcript}; use crate::error::Error; @@ -156,6 +156,30 @@ pub enum PairedAlignmentResult { /// - n_for_mapq: effective alignment count for MAPQ calculation (max of transcript count /// and valid cluster count, to avoid undercounting from coordinate dedup on tandem repeats) /// - unmapped_reason: `Some(reason)` if no alignments produced, `None` if mapped +/// +/// STARlong-only per-window read-space seed coverage (`ReadAlign_stitchPieces.cpp`'s +/// `WC` accumulation): the union, over the window's seeds sorted by read position, of newly +/// covered read bases. Faithful to upstream's `rLast=0` initialization, which under-counts a +/// seed starting at read position 0 by one base. +fn window_read_coverage(cluster: &SeedCluster) -> u64 { + let mut wa: Vec<&crate::align::stitch::WindowAlignment> = cluster.alignments.iter().collect(); + wa.sort_by_key(|a| a.read_pos); + let mut cov = 0u64; + let mut r_last = 0u64; + for a in &wa { + let (r1, l1) = (a.read_pos as u64, a.length as u64); + if r1 + l1 > r_last + 1 { + cov += if r1 > r_last { + l1 + } else { + r1 + l1 - (r_last + 1) + }; + r_last = r1 + l1 - 1; + } + } + cov +} + pub fn align_read( read_seq: &[u8], read_name: &str, @@ -272,9 +296,23 @@ pub fn align_read( let mut clusters = clusters; clusters.truncate(params.align_windows_per_read_nmax); - // NOTE: STAR's winReadCoverageRelativeMin filter is long-reads-only - // (#ifdef COMPILE_FOR_LONG_READS in stitchPieces.cpp). Standard STAR - // does NOT filter clusters by seed coverage. Removed to match STAR. + // STARlong-only (`-DCOMPILE_FOR_LONG_READS`): drop windows whose read-space seed + // coverage falls below both a relative (of the best window) and an absolute floor. + // Inert for the short-read `rustar-aligner` binary (winReadCoverageRelativeMin's + // default 0.5 only applies once `params.long_read` is set). + if params.long_read { + let coverage: Vec = clusters.iter().map(window_read_coverage).collect(); + let cov_max = coverage.iter().copied().max().unwrap_or(0); + let cov_min = cov_max as f64 * params.win_read_coverage_relative_min; + clusters = clusters + .into_iter() + .zip(coverage) + .filter(|(_, cov)| { + (*cov as f64) >= cov_min && *cov >= params.win_read_coverage_bases_min + }) + .map(|(c, _)| c) + .collect(); + } // Step 2b: Detect chimeric alignments from multi-cluster seeds (Tier 2) let mut chimeric_alignments = Vec::new(); @@ -307,6 +345,7 @@ pub fn align_read( &scorer, junction_db, params.align_transcripts_per_window_nmax, + params.long_read, debug_name, ); if debug_read { @@ -2300,4 +2339,58 @@ mod tests { shuffle_tied_prefix(&mut items, |t| t.0, 42); assert_eq!(items, before); } + + fn wa(read_pos: usize, length: usize) -> crate::align::stitch::WindowAlignment { + crate::align::stitch::WindowAlignment { + seed_idx: 0, + read_pos, + length, + genome_pos: read_pos as u64, + sa_pos: read_pos as u64, + n_rep: 1, + is_anchor: true, + mate_id: 2, + pre_ext_score: 0, + } + } + + fn cluster_with(alignments: Vec) -> SeedCluster { + SeedCluster { + alignments, + chr_idx: 0, + genome_start: 0, + genome_end: 0, + is_reverse: false, + anchor_idx: 0, + anchor_bin: 0, + } + } + + #[test] + fn window_read_coverage_single_seed() { + // One seed covering read[0..10): the `rLast=0` init under-counts a seed starting + // at read position 0 by one base (faithful to STAR's own quirk). + let cluster = cluster_with(vec![wa(0, 10)]); + assert_eq!(window_read_coverage(&cluster), 9); + } + + #[test] + fn window_read_coverage_seed_not_at_zero() { + // A seed starting after position 0 is not subject to the rLast=0 undercount. + let cluster = cluster_with(vec![wa(5, 10)]); + assert_eq!(window_read_coverage(&cluster), 10); + } + + #[test] + fn window_read_coverage_union_of_overlapping_seeds() { + // Two overlapping seeds: coverage is the union of read bases, not the sum. + let cluster = cluster_with(vec![wa(5, 10), wa(10, 10)]); + assert_eq!(window_read_coverage(&cluster), 15); + } + + #[test] + fn window_read_coverage_disjoint_seeds() { + let cluster = cluster_with(vec![wa(0, 5), wa(20, 5)]); + assert_eq!(window_read_coverage(&cluster), 9); + } } diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 1c25ecc..54ab616 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -1125,6 +1125,7 @@ impl WorkingTranscript { #[allow(clippy::too_many_arguments)] fn stitch_align_to_transcript( wt: &WorkingTranscript, + last_exon: &ExonBlock, wa: &WindowAlignment, read_seq: &[u8], index: &GenomeIndex, @@ -1134,8 +1135,6 @@ fn stitch_align_to_transcript( align_mates_gap_max: u64, _debug_name: &str, ) -> Option { - let last_exon = 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 // check alignMatesGapMax instead. @@ -1701,6 +1700,7 @@ pub fn stitch_seeds_with_jdb( scorer, junction_db, max_transcripts_per_window, + false, "", ) } @@ -2437,6 +2437,7 @@ fn stitch_recurse( // Try stitching this seed onto the existing transcript if let Some(new_wt) = stitch_align_to_transcript( &wt, + wt.exons.last().unwrap(), wa, read_seq, index, @@ -2658,6 +2659,7 @@ pub(crate) fn split_combined_wt( /// Uses STAR's recursive combinatorial stitcher (stitchWindowAligns) instead of /// forward DP. For each seed, explores include/exclude branches, allowing the /// algorithm to skip spurious short seeds that would create false splices. +#[allow(clippy::too_many_arguments)] pub(crate) fn stitch_seeds_with_jdb_debug( cluster: &SeedCluster, read_seq: &[u8], @@ -2665,18 +2667,31 @@ pub(crate) fn stitch_seeds_with_jdb_debug( scorer: &AlignmentScorer, junction_db: Option<&crate::junction::SpliceJunctionDb>, max_transcripts_per_window: usize, + long_read: bool, debug_read_name: &str, ) -> Vec { - let (working_transcripts, stitch_cluster, stitch_is_reverse, stitch_read) = stitch_seeds_core( - cluster, - read_seq, - index, - scorer, - junction_db, - max_transcripts_per_window, - 0, - debug_read_name, - ); + let (working_transcripts, stitch_cluster, stitch_is_reverse, stitch_read) = if long_read { + stitch_seeds_long_read( + cluster, + read_seq, + index, + scorer, + junction_db, + 0, + debug_read_name, + ) + } else { + stitch_seeds_core( + cluster, + read_seq, + index, + scorer, + junction_db, + max_transcripts_per_window, + 0, + debug_read_name, + ) + }; // Finalize working transcripts → Transcript (filtering by overhang+repeat check) let mut transcripts: Vec = Vec::with_capacity(working_transcripts.len()); @@ -2786,16 +2801,20 @@ pub(crate) fn stitch_seeds_with_jdb_debug( /// Shared core: preprocessing + recursive stitcher, returns working transcripts + context. #[allow(clippy::too_many_arguments)] -pub(crate) fn stitch_seeds_core( +/// Shared preamble for both the recursive short-read stitcher ([`stitch_seeds_core`]) and the +/// STARlong chaining DP ([`stitch_seeds_long_read`]): diagonal-dedup the cluster's seeds, convert +/// reverse-strand coordinates into the stitcher's positive-strand/forward-genome frame, pre-extend +/// each seed's base score (STAR's `scoreSeedBest` base case), sort by read position, and cap the +/// seed count. Returns the prepared WA entries (pre-extended; empty if the cluster had none after +/// dedup/capping), the stitch-frame cluster, whether the original cluster was reverse-strand, and +/// the (possibly RC'd) stitch-frame read. +fn prepare_stitch_input( cluster: &SeedCluster, read_seq: &[u8], index: &GenomeIndex, scorer: &AlignmentScorer, - junction_db: Option<&crate::junction::SpliceJunctionDb>, - max_transcripts_per_window: usize, - align_mates_gap_max: u64, debug_read_name: &str, -) -> (Vec, SeedCluster, bool, Vec) { +) -> (Vec, SeedCluster, bool, Vec) { let debug = !debug_read_name.is_empty(); // Include ALL seeds (anchor and non-anchor) in the stitcher. @@ -3059,6 +3078,41 @@ pub(crate) fn stitch_seeds_core( } } + ( + wa_entries, + stitch_cluster, + stitch_is_reverse, + stitch_read_owned, + ) +} + +/// Inner implementation of the recursive short-read stitcher, sharing [`prepare_stitch_input`] +/// with [`stitch_seeds_long_read`]. +#[allow(clippy::too_many_arguments)] +pub(crate) fn stitch_seeds_core( + cluster: &SeedCluster, + read_seq: &[u8], + index: &GenomeIndex, + scorer: &AlignmentScorer, + junction_db: Option<&crate::junction::SpliceJunctionDb>, + max_transcripts_per_window: usize, + align_mates_gap_max: u64, + debug_read_name: &str, +) -> (Vec, SeedCluster, bool, Vec) { + let debug = !debug_read_name.is_empty(); + let (wa_entries, stitch_cluster, stitch_is_reverse, stitch_read_owned) = + prepare_stitch_input(cluster, read_seq, index, scorer, debug_read_name); + let stitch_read: &[u8] = &stitch_read_owned; + + if wa_entries.is_empty() { + return ( + Vec::new(), + stitch_cluster, + stitch_is_reverse, + stitch_read_owned, + ); + } + // 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. @@ -3133,6 +3187,250 @@ pub(crate) fn stitch_seeds_core( ) } +/// STARlong's `stitchWindowSeeds` (`-DCOMPILE_FOR_LONG_READS`): an O(n^2) chaining DP over a +/// window's seeds that assembles the single best-scoring seed chain into ONE working transcript, +/// replacing [`stitch_seeds_core`]'s recursive include/exclude enumeration (exponential in seed +/// count, intractable for long-read seed counts). Reuses [`stitch_align_to_transcript`] for +/// per-edge scoring; end extension and CIGAR assembly are left to the shared +/// [`finalize_transcript`] pass, exactly as for the short-read path — this port does not +/// duplicate that logic with STARlong's own `L_PREV`-widened end extension. +#[allow(clippy::too_many_arguments)] +pub(crate) fn stitch_seeds_long_read( + cluster: &SeedCluster, + read_seq: &[u8], + index: &GenomeIndex, + scorer: &AlignmentScorer, + junction_db: Option<&crate::junction::SpliceJunctionDb>, + align_mates_gap_max: u64, + debug_read_name: &str, +) -> (Vec, SeedCluster, bool, Vec) { + let (wa_entries, stitch_cluster, stitch_is_reverse, stitch_read_owned) = + prepare_stitch_input(cluster, read_seq, index, scorer, debug_read_name); + let stitch_read: &[u8] = &stitch_read_owned; + + let n = wa_entries.len(); + if n == 0 { + return ( + Vec::new(), + stitch_cluster, + stitch_is_reverse, + stitch_read_owned, + ); + } + + // `extendAlign`'s ratio-based mismatch budget collapses to a flat `n_mm_max` when + // `len_prev` is this large (STAR `stitchWindowSeeds.cpp`'s `L_PREV` constant). + const L_PREV: usize = 100_000; + + let mut score_best = vec![0i32; n]; + let mut mm_best = vec![0u32; n]; + // Predecessor seed index, -1 = none. A chain start points to itself. + let mut pred: Vec = vec![-1; n]; + + let exon_block = |wa: &WindowAlignment| ExonBlock { + read_start: wa.read_pos, + read_end: wa.read_pos + wa.length, + genome_start: wa.sa_pos, + genome_end: wa.sa_pos + wa.length as u64, + mate_id: wa.mate_id, + }; + + // ---- Phase A: forward DP (best chain score ending at each seed) ---- + for i1 in 0..n { + for i2 in 0..=i1 { + if i2 < i1 { + let a = &wa_entries[i2]; + let b = &wa_entries[i1]; + let a_exon = exon_block(a); + let wt_a = WorkingTranscript { + exons: vec![a_exon.clone()], + score: a.length as i32, + n_mismatch: mm_best[i2], + n_gap: 0, + n_junction: 0, + junction_motifs: Vec::new(), + junction_annotated: Vec::new(), + junction_shifts: Vec::new(), + n_anchor: u32::from(a.is_anchor), + read_start: a_exon.read_start, + read_end: a_exon.read_end, + genome_start: a_exon.genome_start, + genome_end: a_exon.genome_end, + }; + let Some(new_wt) = stitch_align_to_transcript( + &wt_a, + &a_exon, + b, + stitch_read, + index, + scorer, + &stitch_cluster, + junction_db, + align_mates_gap_max, + debug_read_name, + ) else { + continue; + }; + // Marginal score of the A->B edge (STAR's `score2`): the DP accumulates + // these, it does not re-add A's own base score at every step. + let d_score = new_wt.score - wt_a.score; + let is_annotated = new_wt.junction_annotated.last().copied().unwrap_or(false); + let overhang = if is_annotated { + scorer.align_sjdb_overhang_min + } else { + scorer.align_sj_overhang_min + } as usize; + let exon_a_len = new_wt.exons[0].read_end - new_wt.exons[0].read_start; + if exon_a_len >= overhang + && d_score > 0 + && score_best[i2] + d_score > score_best[i1] + { + score_best[i1] = score_best[i2] + d_score; + mm_best[i1] = new_wt.n_mismatch; + pred[i1] = i2 as i64; + } + } else { + // Chain start: seed length plus a left extension to the read start. + let s = &wa_entries[i1]; + let mut score2 = s.length as i32; + let mut ext_len = 0usize; + if s.read_pos > 0 { + let e = extend_alignment( + stitch_read, + s.read_pos, + s.sa_pos, + -1, + s.read_pos, + 0, + L_PREV, + scorer.n_mm_max, + scorer.p_mm_max, + index, + false, + ); + score2 += e.max_score; + ext_len = e.extend_len; + } + let exon_long_enough = + (s.length + ext_len) >= scorer.align_sj_overhang_min as usize; + if exon_long_enough && score2 > score_best[i1] { + score_best[i1] = score2; + pred[i1] = i1 as i64; + // mm_best[i1] intentionally left at 0 here: STARlong's own chain-start + // case never records the left extension's mismatches into + // `scoreSeedBestMM` either (faithful to the upstream quirk). + } + } + } + } + + // ---- Phase B: right-extend each chain, track the global best ---- + let mut global_best = 0i32; + let mut best_ind = 0usize; + let mut found = false; + for i1 in 0..n { + let s = &wa_entries[i1]; + let t_r2 = s.read_pos + s.length; + let mut ext_len = 0usize; + if t_r2 < stitch_read.len() { + let e = extend_alignment( + stitch_read, + t_r2, + s.sa_pos + s.length as u64, + 1, + stitch_read.len() - t_r2, + mm_best[i1], + L_PREV, + scorer.n_mm_max, + scorer.p_mm_max, + index, + false, + ); + score_best[i1] += e.max_score; + ext_len = e.extend_len; + } + let exon_long_enough = (s.length + ext_len) >= scorer.align_sj_overhang_min as usize; + if exon_long_enough && (!found || score_best[i1] > global_best) { + global_best = score_best[i1]; + best_ind = i1; + found = true; + } + } + if !found { + return ( + Vec::new(), + stitch_cluster, + stitch_is_reverse, + stitch_read_owned, + ); + } + + // ---- Phase C: traceback the best chain (rightmost -> leftmost) ---- + let mut chain: Vec = Vec::new(); + let mut cur = best_ind; + loop { + chain.push(cur); + let p = pred[cur]; + if p >= 0 && (cur as i64) > p { + cur = p as usize; + } else { + break; + } + } + + // ---- Phase D: assemble the transcript by replaying the chain, left to right ---- + let first = &wa_entries[chain[chain.len() - 1]]; + let first_exon = exon_block(first); + let mut wt = WorkingTranscript { + exons: vec![first_exon.clone()], + score: first.length as i32, + n_mismatch: 0, + n_gap: 0, + n_junction: 0, + junction_motifs: Vec::new(), + junction_annotated: Vec::new(), + junction_shifts: Vec::new(), + n_anchor: u32::from(first.is_anchor), + read_start: first_exon.read_start, + read_end: first_exon.read_end, + genome_start: first_exon.genome_start, + genome_end: first_exon.genome_end, + }; + for isc in (1..chain.len()).rev() { + let last_exon = wt.exons.last().unwrap().clone(); + let b = &wa_entries[chain[isc - 1]]; + let Some(new_wt) = stitch_align_to_transcript( + &wt, + &last_exon, + b, + stitch_read, + index, + scorer, + &stitch_cluster, + junction_db, + align_mates_gap_max, + debug_read_name, + ) else { + // The traceback replays exactly the edges Phase A already validated; a + // failure here means this window yields no long-read alignment. + return ( + Vec::new(), + stitch_cluster, + stitch_is_reverse, + stitch_read_owned, + ); + }; + wt = new_wt; + } + + ( + vec![wt], + stitch_cluster, + stitch_is_reverse, + stitch_read_owned, + ) +} + #[cfg(test)] mod tests { use super::*; @@ -3698,4 +3996,109 @@ mod tests { let replacement_correct = baseline + scorer.sjdb_score; assert_ne!(additive_buggy, replacement_correct); } + + #[test] + fn long_read_single_seed_chain() { + // STARlong chaining DP with exactly one seed spanning the whole read: Phase A/B/C/D + // should hand back that one seed's own exon, untouched (end extension is deferred to + // the shared `finalize_transcript` pass, not duplicated here). + let genome_seq: Vec = (0..40u32).map(|i| ((i * 7 + 3) % 4) as u8).collect(); + let index = make_index_with_seq(&genome_seq); + let read_seq = genome_seq[10..20].to_vec(); + + let cluster = SeedCluster { + alignments: vec![WindowAlignment { + seed_idx: 0, + read_pos: 0, + length: 10, + genome_pos: 10, + sa_pos: 10, + n_rep: 1, + is_anchor: true, + mate_id: 2, + pre_ext_score: 0, + }], + chr_idx: 0, + genome_start: 10, + genome_end: 20, + is_reverse: false, + anchor_idx: 0, + anchor_bin: 0, + }; + let scorer = AlignmentScorer::from_params_minimal(); + + let (transcripts, _cluster, is_reverse, _read) = + stitch_seeds_long_read(&cluster, &read_seq, &index, &scorer, None, 0, ""); + + assert_eq!(transcripts.len(), 1); + assert!(!is_reverse); + let wt = &transcripts[0]; + assert_eq!(wt.exons.len(), 1); + assert_eq!(wt.exons[0].read_start, 0); + assert_eq!(wt.exons[0].read_end, 10); + assert_eq!(wt.exons[0].genome_start, 10); + assert_eq!(wt.exons[0].genome_end, 20); + assert_eq!(wt.score, 10); + } + + #[test] + fn long_read_chains_two_seeds_across_deletion() { + // Two seeds bridged by a short (3bp) deletion — below alignIntronMin, so + // stitchAlignToTranscript's deletion path is taken, not the splice path. The + // chaining DP should still assemble both seeds into a single 2-exon transcript. + let genome_seq: Vec = (0..40u32).map(|i| ((i * 7 + 3) % 4) as u8).collect(); + let index = make_index_with_seq(&genome_seq); + let mut read_seq = genome_seq[0..10].to_vec(); + read_seq.extend_from_slice(&genome_seq[13..23]); + + let cluster = SeedCluster { + alignments: vec![ + WindowAlignment { + seed_idx: 0, + read_pos: 0, + length: 10, + genome_pos: 0, + sa_pos: 0, + n_rep: 1, + is_anchor: true, + mate_id: 2, + pre_ext_score: 0, + }, + WindowAlignment { + seed_idx: 1, + read_pos: 10, + length: 10, + genome_pos: 13, + sa_pos: 13, + n_rep: 1, + is_anchor: true, + mate_id: 2, + pre_ext_score: 0, + }, + ], + chr_idx: 0, + genome_start: 0, + genome_end: 23, + is_reverse: false, + anchor_idx: 0, + anchor_bin: 0, + }; + let scorer = AlignmentScorer::from_params_minimal(); + + let (transcripts, _cluster, _is_reverse, _read) = + stitch_seeds_long_read(&cluster, &read_seq, &index, &scorer, None, 0, ""); + + assert_eq!(transcripts.len(), 1); + let wt = &transcripts[0]; + assert_eq!( + wt.exons.len(), + 2, + "both seeds must chain into one transcript" + ); + assert_eq!(wt.n_gap, 1, "the 3bp gap is a deletion, not a junction"); + assert_eq!(wt.n_junction, 0); + assert_eq!(wt.exons[0].read_start, 0); + assert_eq!(wt.exons[1].read_end, 20); + assert_eq!(wt.exons[1].genome_end, 23); + } } diff --git a/src/bin/rustar_aligner_long.rs b/src/bin/rustar_aligner_long.rs new file mode 100644 index 0000000..304b5e1 --- /dev/null +++ b/src/bin/rustar_aligner_long.rs @@ -0,0 +1,19 @@ +use rustar_aligner::cpu; +use rustar_aligner::params::Parameters; + +/// Global allocator override — see `main.rs` / `Cargo.toml`'s `mimalloc` comment. +#[global_allocator] +static GLOBAL: mimalloc::MiMalloc = mimalloc::MiMalloc; + +/// STAR's `STARlong` binary: identical CLI to `rustar-aligner`, but forces the +/// chaining-DP long-read stitcher (STAR's `-DCOMPILE_FOR_LONG_READS`), matching +/// how native STAR selects it by which binary you run, not a CLI flag. +fn main() -> anyhow::Result<()> { + env_logger::Builder::from_env(env_logger::Env::default().default_filter_or("info")).init(); + + cpu::check_cpu_compat()?; + + let mut params = Parameters::parse(); + params.long_read = true; + rustar_aligner::run(¶ms) +} diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..c37fd22 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -481,6 +481,12 @@ pub struct Parameters { #[arg(long = "tempDir")] pub temp_dir: Option, + /// STAR's `STARlong` mode (chaining-DP stitcher for long reads): not a CLI + /// flag, matching native STAR (a separate binary, `-DCOMPILE_FOR_LONG_READS`); + /// set by the `rustar-aligner-long` binary entry point after parsing. + #[arg(skip)] + pub long_read: bool, + // ── Genome ────────────────────────────────────────────────────────── /// Path to genome index directory #[arg(long = "genomeDir", default_value = "./GenomeDir")] @@ -915,6 +921,10 @@ pub struct Parameters { #[arg(long = "alignTranscriptsPerWindowNmax", default_value_t = 100)] pub align_transcripts_per_window_nmax: usize, + /// STARlong only: min absolute covered bases per window + #[arg(long = "winReadCoverageBasesMin", default_value_t = 0)] + pub win_read_coverage_bases_min: u64, + // ── Splice junction database ──────────────────────────────────────── /// GTF file for splice junction annotations #[arg(long = "sjdbGTFfile")] From e76435143666bbbdf7b1dacddfb85fd3062a7626 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 00:35:34 +0200 Subject: [PATCH 2/5] fix(starlong): adopt the extend_to_end parameter added by #145 `extend_alignment` gained an `extend_to_end` argument when --alignEndsType landed. The two long-read chaining call sites pass `false`: they score a seed's own extension potential during chaining, and --alignEndsType applies at the read ends, not there. Co-Authored-By: Claude Opus 5 (1M context) --- src/align/stitch.rs | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 54ab616..cd3f779 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -3307,6 +3307,10 @@ pub(crate) fn stitch_seeds_long_read( scorer.p_mm_max, index, false, + // Long-read chaining scores a seed's own extension + // potential; --alignEndsType applies at the read ends, + // not here. + false, ); score2 += e.max_score; ext_len = e.extend_len; @@ -3345,6 +3349,7 @@ pub(crate) fn stitch_seeds_long_read( scorer.p_mm_max, index, false, + false, // as above: a chaining-time extension is always local ); score_best[i1] += e.max_score; ext_len = e.extend_len; From e2fa5a2f7960b8e4ad16c3f241f3b1a0526e83e7 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 20:51:20 +0200 Subject: [PATCH 3/5] fix(starlong): do not discard a window when a traceback edge fails The chaining DP scores an edge in Phase A with one function and the traceback re-stitches it with stitch_align_to_transcript, so the two can disagree. On a disagreement this returned an empty transcript vector, discarding the whole window. STAR does not. stitchWindowSeeds (ReadAlign_stitchWindowSeeds.cpp:130-138) calls stitchAlignToTranscript in its traceback and adds whatever it returns unconditionally; the `return` on failure is commented out and marked "this should not happen". Keep the chain assembled so far and let the score filters decide, which is what STAR's arithmetic amounts to. No effect on this branch alone: 500 simulated long reads against STARlong 2.7.11b give the same 419 primaries and 360 exact records before and after. It matters once the annotated-junction work in #146 is in, because that adds rejection paths the traceback can hit. Measured on a local merge of the two branches: window discarded 407 primaries, 388 exact, 25 STAR-only chain kept 422 primaries, 388 exact, 10 STAR-only STAR aligns 430. Adds test/simulate_long_reads.py, which draws whole spliced transcripts as reads from a seeded RNG, so the fixture is byte-reproducible; the bundled 150 bp tiers never reach this path. --- src/align/stitch.rs | 19 +++--- test/simulate_long_reads.py | 120 ++++++++++++++++++++++++++++++++++++ 2 files changed, 131 insertions(+), 8 deletions(-) create mode 100644 test/simulate_long_reads.py diff --git a/src/align/stitch.rs b/src/align/stitch.rs index cd3f779..fb87917 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -3416,14 +3416,17 @@ pub(crate) fn stitch_seeds_long_read( align_mates_gap_max, debug_read_name, ) else { - // The traceback replays exactly the edges Phase A already validated; a - // failure here means this window yields no long-read alignment. - return ( - Vec::new(), - stitch_cluster, - stitch_is_reverse, - stitch_read_owned, - ); + // STAR does not abandon the chain here. `stitchWindowSeeds` + // (`ReadAlign_stitchWindowSeeds.cpp:130-138`) calls + // `stitchAlignToTranscript` in the traceback and adds whatever it + // returns unconditionally; the `return` on failure is commented out + // and marked "this should not happen". Phase A scores an edge with + // one function and the traceback re-stitches it with another, so + // the two can disagree, and treating that as "this window has no + // alignment" throws away chains that are otherwise the best in + // their window. Keep what has been assembled and let the score + // filters decide. + break; }; wt = new_wt; } diff --git a/test/simulate_long_reads.py b/test/simulate_long_reads.py new file mode 100644 index 0000000..a6a6a2d --- /dev/null +++ b/test/simulate_long_reads.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python3 +"""Draw long reads from spliced transcripts, for the STARlong differential. + +The bundled tiers are 150 bp Illumina reads, which never reach the long-read +code path: `--winReadCoverageRelativeMin` and `--winReadCoverageBasesMin` only +filter windows when the long-read stitcher runs, and a 150 bp read produces one +window with near-total coverage. Reads have to be long enough, and span enough +junctions, for the coverage filter to be able to discard anything. + +Each read is the full spliced sequence of one transcript, so it crosses every +junction that transcript has. Transcripts shorter than `--min-length` are +skipped; longer ones are emitted whole. Errors are substitutions only, at a +fixed rate, drawn from a seeded RNG so the FASTQ is byte-reproducible. + + simulate_long_reads.py [n] [min_len] +""" + +import random +import sys +from collections import defaultdict + +COMPLEMENT = str.maketrans("ACGTNacgtn", "TGCANtgcan") + + +def read_fasta(path): + seqs, name, chunks = {}, None, [] + with open(path) as fh: + for line in fh: + if line.startswith(">"): + if name is not None: + seqs[name] = "".join(chunks) + name = line[1:].split()[0] + chunks = [] + else: + chunks.append(line.strip()) + if name is not None: + seqs[name] = "".join(chunks) + return seqs + + +def read_exons(path): + """transcript_id -> (chrom, strand, [(start, end), ...]) in 1-based inclusive GTF coords.""" + exons = defaultdict(list) + meta = {} + with open(path) as fh: + for line in fh: + if line.startswith("#"): + continue + f = line.rstrip("\n").split("\t") + if len(f) < 9 or f[2] != "exon": + continue + attrs = f[8] + key = 'transcript_id "' + i = attrs.find(key) + if i < 0: + continue + j = attrs.find('"', i + len(key)) + tid = attrs[i + len(key) : j] + exons[tid].append((int(f[3]), int(f[4]))) + meta[tid] = (f[0], f[6]) + return exons, meta + + +def main(): + if len(sys.argv) < 4: + sys.exit(__doc__) + genome_path, gtf_path, out_path = sys.argv[1:4] + n_reads = int(sys.argv[4]) if len(sys.argv) > 4 else 500 + min_len = int(sys.argv[5]) if len(sys.argv) > 5 else 1000 + error_rate = 0.02 + + rng = random.Random(100) + genome = read_fasta(genome_path) + exons, meta = read_exons(gtf_path) + + # Multi-exon transcripts only: a single-exon read crosses no junction and + # exercises none of the long-read window logic. + candidates = [] + for tid, blocks in exons.items(): + if len(blocks) < 2: + continue + chrom, strand = meta[tid] + if chrom not in genome: + continue + blocks = sorted(blocks) + length = sum(e - s + 1 for s, e in blocks) + if length >= min_len: + candidates.append((tid, chrom, strand, blocks, length)) + + candidates.sort() + if not candidates: + sys.exit(f"no multi-exon transcript of at least {min_len} bp found") + rng.shuffle(candidates) + + written = 0 + with open(out_path, "w") as out: + while written < n_reads: + tid, chrom, strand, blocks, _ = candidates[written % len(candidates)] + chrom_seq = genome[chrom] + spliced = "".join(chrom_seq[s - 1 : e] for s, e in blocks).upper() + if strand == "-": + spliced = spliced.translate(COMPLEMENT)[::-1] + if "N" in spliced: + # An N run would be clipped rather than aligned and tells us + # nothing about the coverage filter. + written += 1 + continue + bases = list(spliced) + for i in range(len(bases)): + if rng.random() < error_rate: + bases[i] = rng.choice([b for b in "ACGT" if b != bases[i]]) + seq = "".join(bases) + out.write(f"@sim_{written}_{tid}_{len(blocks)}ex\n{seq}\n+\n{'I' * len(seq)}\n") + written += 1 + + print(f"wrote {written} reads to {out_path}", file=sys.stderr) + + +if __name__ == "__main__": + main() From f6f03d8b8283c45b9fba380164d964cc0e48024a Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 21:15:23 +0200 Subject: [PATCH 4/5] fix(starlong): STARlong does not apply the per-motif stitch mismatch cap stitchAlignToTranscript.cpp:308-315 gates the junction differently in the two binaries: #ifdef COMPILE_FOR_LONG_READS if ( (trA->nMM + nMM) <= outFilterMismatchNmaxTotal ) #else if ( (trA->nMM + nMM) <= outFilterMismatchNmaxTotal && ( jCan<0 || (jCan<7 && nMM <= P.alignSJstitchMismatchNmax[(jCan+1)/2]) ) ) #endif The alignSJstitchMismatchNmax term is compiled out of STARlong. rustar applied it in both, so a long read that reaches a junction with a mismatch or two beside it lost the edge. The cost is not one junction: the chaining DP then splits into a chain on each side of the junction, and neither half clears --outFilterScoreMinOverLread, so the read comes out unmapped. That is what it was doing. On 500 simulated spliced long reads against STARlong 2.7.11b, 13 reads were unmapped that STAR maps; 11 of them were multi-exon mitochondrial transcripts, and one is an annotated 35-base intron (Mito 24871-24905) rejected as non-canonical with a single mismatch nearby. before 419 primaries, 360 exact, 405 same-locus, 13 STAR-only after 432 primaries, 372 exact, 423 same-locus, 1 STAR-only STAR maps 430. The three rustar-only reads are 6-exon mitochondrial transcripts it now aligns and STAR does not. --- src/align/score.rs | 78 +++++++++++++++++++++++++++++++++++++++++++++ src/align/stitch.rs | 2 ++ 2 files changed, 80 insertions(+) diff --git a/src/align/score.rs b/src/align/score.rs index 7a34f55..d742cee 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -28,6 +28,12 @@ pub struct AlignmentScorer { /// Max mismatches for stitching SJs: [non-canonical, GT-AG, GC-AG, AT-AC] /// -1 means unlimited pub align_sj_stitch_mismatch_nmax: [i32; 4], + /// Whether this run is STARlong. STAR compiles the per-motif stitch + /// mismatch cap out of the long-read binary entirely + /// (`stitchAlignToTranscript.cpp:308-315`: the + /// `alignSJstitchMismatchNmax` term sits in the `#else` of + /// `COMPILE_FOR_LONG_READS`), leaving only the total mismatch budget. + pub long_read: bool, /// Max absolute number of mismatches for alignment extension (outFilterMismatchNmax) pub n_mm_max: u32, /// Max ratio of mismatches to total alignment length (outFilterMismatchNoverLmax) @@ -69,6 +75,7 @@ impl AlignmentScorer { align_intron_min: 21, sjdb_score: 2, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], + long_read: false, n_mm_max: 10, p_mm_max: 0.3, align_sj_overhang_min: 5, @@ -102,6 +109,7 @@ impl AlignmentScorer { params.align_sj_stitch_mismatch_nmax[2], params.align_sj_stitch_mismatch_nmax[3], ], + long_read: params.long_read, n_mm_max: params.out_filter_mismatch_nmax, p_mm_max: params.out_filter_mismatch_nover_lmax, align_sj_overhang_min: params.align_sj_overhang_min, @@ -143,6 +151,23 @@ impl AlignmentScorer { /// Check if the number of mismatches is allowed for this junction motif type pub fn stitch_mismatch_allowed(&self, motif: &SpliceMotif, n_mismatch: u32) -> bool { + // STARlong does not apply this cap. In STAR's source the condition + // reads + // + // #ifdef COMPILE_FOR_LONG_READS + // if ( (trA->nMM + nMM) <= outFilterMismatchNmaxTotal ) + // #else + // if ( (trA->nMM + nMM) <= outFilterMismatchNmaxTotal + // && ( jCan<0 || (jCan<7 && nMM <= P.alignSJstitchMismatchNmax[(jCan+1)/2]) ) ) + // #endif + // + // so the per-motif term is compiled out of STARlong and only the total + // mismatch budget survives. Applying it anyway kills junctions that a + // long read reaches with a mismatch or two nearby, which on spliced + // long reads means losing the read outright rather than one junction. + if self.long_read { + return true; + } let idx = match motif { SpliceMotif::NonCanonical => 0, SpliceMotif::GtAg | SpliceMotif::CtAc => 1, @@ -624,6 +649,47 @@ mod tests { use super::*; + /// STAR compiles the per-motif stitch mismatch cap out of STARlong: + /// `stitchAlignToTranscript.cpp:308-315` keeps only + /// `(trA->nMM + nMM) <= outFilterMismatchNmaxTotal` under + /// `COMPILE_FOR_LONG_READS`, and puts the `alignSJstitchMismatchNmax` + /// term in the `#else`. + /// + /// Applying it in the long-read path costs whole reads, not single + /// junctions: a spliced long read that hits one mismatch beside a + /// junction loses the edge, the chaining DP then splits into two chains + /// on either side, and neither half clears + /// `--outFilterScoreMinOverLread`. + #[test] + fn the_stitch_mismatch_cap_does_not_apply_to_long_reads() { + let mut scorer = AlignmentScorer::from_params_minimal(); + assert_eq!(scorer.align_sj_stitch_mismatch_nmax[0], 0); + + // Short reads: a non-canonical junction tolerates no mismatch, and + // GT/AG tolerates any number. + assert!(scorer.stitch_mismatch_allowed(&SpliceMotif::NonCanonical, 0)); + assert!(!scorer.stitch_mismatch_allowed(&SpliceMotif::NonCanonical, 1)); + assert!(scorer.stitch_mismatch_allowed(&SpliceMotif::GtAg, 50)); + + // STARlong: the cap is gone for every motif. + scorer.long_read = true; + assert!(scorer.stitch_mismatch_allowed(&SpliceMotif::NonCanonical, 1)); + assert!(scorer.stitch_mismatch_allowed(&SpliceMotif::GcAg, 9)); + assert!(scorer.stitch_mismatch_allowed(&SpliceMotif::AtAc, 1000)); + } + + /// `--soloType`-style plumbing check: the flag reaches the scorer from + /// the parameters rather than being set by hand at one call site. + #[test] + fn the_scorer_takes_long_read_from_the_parameters() { + let short = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "r.fq"]); + assert!(!AlignmentScorer::from_params(&short).long_read); + + let mut long = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "r.fq"]); + long.long_read = true; + assert!(AlignmentScorer::from_params(&long).long_read); + } + fn make_test_genome(seq: &[u8]) -> Genome { // Create simple genome with one chromosome let n_genome = ((seq.len() as u64 + 1) / 64 + 1) * 64; // Pad to 64-byte boundary @@ -681,6 +747,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -726,6 +793,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -770,6 +838,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -812,6 +881,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -847,6 +917,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -880,6 +951,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -921,6 +993,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -960,6 +1033,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -1000,6 +1074,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -1108,6 +1183,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -1188,6 +1264,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -1240,6 +1317,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, diff --git a/src/align/stitch.rs b/src/align/stitch.rs index fb87917..b6d3b2d 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -3838,6 +3838,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, @@ -3882,6 +3883,7 @@ mod tests { score_ins_base: -2, align_intron_min: 21, sjdb_score: 2, + long_read: false, align_sj_stitch_mismatch_nmax: [0, -1, 0, 0], n_mm_max: 10, p_mm_max: 0.3, From bbda1319129c4beeacac027edeefb712873ce2d1 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 21:36:26 +0200 Subject: [PATCH 5/5] docs: record the STARlong faithfulness fixes and the reads they recover MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds the CHANGELOG entry the branch was missing and DIVERGENCE.md §2.2 for the three long reads rustar-aligner maps and STARlong does not. Those three are recorded as a divergence rather than fixed. STAR reports them as "unmapped: other", which is a failure to find a window rather than a decision to reject, and rustar-aligner places them on the mitochondrial transcript they were simulated from with the same exon structure STAR produces for other reads of that same transcript. Making them match would mean discarding correct alignments. Rebased onto current main, which restores DIVERGENCE.md: the branch predated it, so the diff read as deleting the file. --- CHANGELOG.md | 30 ++++++++++++++++++++++++++++++ DIVERGENCE.md | 19 +++++++++++++++++++ 2 files changed, 49 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..69f251d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,36 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- **`rustar-aligner-long`, mirroring STAR's STARlong**, with the long-read + window-coverage filter (`--winReadCoverageRelativeMin`, and the new + `--winReadCoverageBasesMin`). STAR selects the long-read path by binary + rather than by flag, and so does this. + + Two places where the long-read path diverged from STARlong and now does not: + + - **The per-motif stitch mismatch cap does not apply.** + `stitchAlignToTranscript.cpp:308-315` keeps only the total mismatch budget + under `COMPILE_FOR_LONG_READS` and puts `alignSJstitchMismatchNmax` in the + `#else`, so it is compiled out of STARlong. Applying it cost whole reads, + not single junctions: the chaining DP split into a chain either side of the + rejected junction and neither half cleared + `--outFilterScoreMinOverLread`. + - **A traceback edge that fails to re-stitch no longer discards its window.** + STAR's `stitchWindowSeeds` adds whatever `stitchAlignToTranscript` returns + without checking, its failure `return` commented out and marked "this + should not happen". + + On 500 simulated spliced long reads against STARlong 2.7.11b (STAR maps 430 + of 500): + + | | primaries | exact records | same locus | unmapped by rustar | + |---|---|---|---|---| + | before | 419 | 360 | 405 | 13 | + | after | 432 | 372 | 423 | 1 | + + Three of the extra primaries are reads STARlong leaves unmapped and + rustar-aligner places correctly; see `DIVERGENCE.md` §2.2. + - **STARsolo single-cell quantification (`--soloType`)** — the 10x Chromium / plate-based count-matrix pipeline, ported from STAR and verified against real STARsolo (#90). diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..6ab3bfe 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -47,6 +47,25 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST **Source.** See `CLAUDE.md` (PE status) and `STAR-RS-COMPARISON.md`. +### 2.2 Three long reads STARlong does not map at all + +On 500 simulated spliced long reads (`test/simulate_long_reads.py`, mean 2.1 kb, every read multi-exon), `rustar-aligner-long` maps three that STARlong 2.7.11b leaves unmapped: `sim_371_Q0105_mRNA_6ex`, `sim_448_Q0105_mRNA_6ex` and `sim_476_Q0115_mRNA_3ex`. + +**What STAR does.** It reports all three as `unmapped: other`, not as filtered: `Log.final.out` shows 0 too-short, 0 too-many-mismatches, 3 other. That marker means no window produced a usable transcript, so this is a failure to find the alignment rather than a decision to reject it. + +**What rustar-aligner does.** It places them on the mitochondrial transcript they were drawn from, with the same exon structure STAR itself produces for other reads of that transcript: + +``` +STAR sim_140_Q0105 Mito:36540 415M768N14M1404N77M1623N250M1417N46M738N356M +rustar sim_371_Q0105 Mito:36540 415M768N28M1404N63M1623N250M1417N50M738N352M +``` + +**Why.** These reads are simulated from `Q0105` and `Q0115`, so the locus and the junction set are known, and both match. The alignments are correct. + +**Impact.** Three reads on this tier count against raw exactness while being improvements. Not corrected: matching STAR here would mean discarding correct alignments to reproduce a failure. + +**Source.** `test/simulate_long_reads.py` for the fixture; STARlong `Log.final.out` for the unmapped classification. + --- ## 3. Output-file metadata divergences