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/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/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 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/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 1c25ecc..b6d3b2d 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,258 @@ 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, + // 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; + } + 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, + false, // as above: a chaining-time extension is always local + ); + 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 { + // 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; + } + + ( + vec![wt], + stitch_cluster, + stitch_is_reverse, + stitch_read_owned, + ) +} + #[cfg(test)] mod tests { use super::*; @@ -3532,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, @@ -3576,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, @@ -3698,4 +4006,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")] 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()