diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e6..5d85a6f2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,8 +11,55 @@ 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 +- 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/DIVERGENCE.md b/DIVERGENCE.md index f0890fbd..a693341f 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 1b2f112b..9aace6e3 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 diff --git a/examples/jrbench.rs b/examples/jrbench.rs index d987c9b5..db285e59 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 8fc7febc..382eddc7 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. @@ -493,16 +448,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 @@ -864,6 +829,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 { @@ -1133,84 +1102,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(), - )) - }); - } - - // 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() + // 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() - .filter(|(i, _)| keep[*i]) - .map(|(_, p)| p) + .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). @@ -1669,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; @@ -1736,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 899b6a4c..cd429c00 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; @@ -478,93 +444,24 @@ 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; - - // 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; - } - } - - 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, + 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 @@ -576,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 { @@ -596,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 { @@ -610,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); } } @@ -650,6 +531,202 @@ 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, &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). + 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: 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 = 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 motif = MotifWindow { + donor: 0, + d1: don[k], + d2: don[k + 1], + a1: acc[k], + a2: acc[k + 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 { + 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 }; + if bound <= max_score2 { + 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 usize); + let gd = scan.seq.base(g_dn_pos 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; + let w = match window.as_mut() { + Some(w) => { + w.slide_to(donor_sa, scan.del as u64, scan.seq); + &*w + } + 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); + 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. + if 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: 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 }; + if bound <= max_score2 { + 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( @@ -683,6 +760,24 @@ 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, + seq: SeqView<'a>, +} + /// The four bases that decide a splice motif: the intron's first two and last /// two, on the forward strand. /// @@ -777,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, @@ -792,18 +886,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, + /// 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: 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 @@ -946,6 +1056,75 @@ 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 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, + 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): @@ -1645,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 ); @@ -1689,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); @@ -1740,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 dfcc025f..cc417e72 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 09a2d799..0d4a3191 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; } @@ -137,51 +126,91 @@ fn score_region( /// staging copy is amortized, small enough to sit on the stack. const SCORE_REGION_CHUNK: usize = 256; -fn count_mismatches( +/// 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; + +/// 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. +#[allow(clippy::too_many_arguments)] +fn extend_alignment_cached( + jcache: &mut crate::align::score::JunctionScanCache, read_seq: &[u8], - cigar_ops: &[cigar::Op], - genome_start: u64, 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, -) -> 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 => {} - } + extend_to_end: bool, +) -> ExtendResult { + 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, + extend_to_end, + ); } - - n_mismatch + // 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 as usize, + 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, + extend_to_end, + ); + if jcache.ext.len() < MEMO_EXTENSIONS_MAX { + jcache + .ext + .insert(key, (r.extend_len as u32, r.max_score, r.n_mismatch)); + } + r } /// Result of extending an alignment into flanking regions @@ -204,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 @@ -212,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], @@ -225,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 { @@ -236,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. @@ -280,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, @@ -355,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; } @@ -436,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, @@ -449,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 @@ -980,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, }, ); } @@ -989,14 +1010,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; @@ -1011,12 +1038,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]) @@ -1051,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, }, ); } @@ -1090,7 +1110,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 @@ -1124,6 +1144,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 { @@ -1143,16 +1168,92 @@ impl WorkingTranscript { read_end: 0, genome_start: 0, genome_end: 0, + window_score: 0, + mapped_len: 0, + } + } +} + +/// 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, @@ -1162,8 +1263,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 @@ -1177,34 +1280,15 @@ 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 - // 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 = &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; - } - } 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: @@ -1213,33 +1297,32 @@ 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 right_ext = extend_alignment( + 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, last_exon.read_end, 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, index, - cluster.is_reverse, false, // internal stitch extension: alignEndsType applied at finalize ); if right_ext.extend_len > 0 { @@ -1288,7 +1371,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, @@ -1299,7 +1383,6 @@ fn stitch_align_to_transcript( 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 { @@ -1310,7 +1393,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 @@ -1321,7 +1404,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; @@ -1332,7 +1415,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; @@ -1346,7 +1429,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 @@ -1356,10 +1439,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; @@ -1372,7 +1454,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 @@ -1406,10 +1488,9 @@ 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); + let motif = scorer.detect_splice_motif(donor_sa, del, &index.genome); ( 0, motif, @@ -1426,8 +1507,6 @@ fn stitch_align_to_transcript( 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, ) @@ -1435,16 +1514,11 @@ fn stitch_align_to_transcript( // 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 None; + return false; } } @@ -1469,7 +1543,6 @@ fn stitch_align_to_transcript( last_exon.genome_end, junction_offset, index, - cluster.is_reverse, ); shared_mm += mm; shared_score += s; @@ -1484,7 +1557,6 @@ fn stitch_align_to_transcript( last_exon.genome_end + junction_offset as u64 + del as u64, acceptor_bases, index, - cluster.is_reverse, ); shared_mm += mm; shared_score += s; @@ -1505,7 +1577,6 @@ fn stitch_align_to_transcript( extra_genome_start, n_extra, index, - cluster.is_reverse, ); shared_mm += extra_mm; shared_score -= 2 * extra_mm as i32; @@ -1523,7 +1594,6 @@ fn stitch_align_to_transcript( extra_genome_start, n_extra, index, - cluster.is_reverse, ); shared_mm += extra_mm; shared_score -= 2 * extra_mm as i32; @@ -1540,15 +1610,13 @@ 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 // 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() @@ -1561,13 +1629,9 @@ fn stitch_align_to_transcript( // 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 @@ -1586,17 +1650,12 @@ 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 - } 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 None; // STAR -1000006 + return false; // STAR -1000006 } jr_exon = shifted; jj_l = u32::from(pj.shift_left); @@ -1663,19 +1722,13 @@ fn stitch_align_to_transcript( 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 { @@ -1715,7 +1768,7 @@ fn stitch_align_to_transcript( // 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 @@ -1771,7 +1824,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 @@ -1792,7 +1845,7 @@ fn stitch_align_to_transcript( new_wt.n_anchor += 1; } - Some(new_wt) + true } /// Stitch seeds within a cluster using recursive combinatorial stitching. @@ -1845,6 +1898,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; @@ -1944,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 { @@ -1963,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 { @@ -1984,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 { @@ -2003,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 { @@ -2199,20 +2279,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, - ); - - // Reverse CIGAR for reverse strand - if cluster.is_reverse { - final_cigar.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; // Compute total reference-consuming length from CIGAR let mut ref_len = 0u64; @@ -2229,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() { @@ -2286,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); @@ -2301,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, @@ -2562,7 +2628,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, @@ -2591,233 +2657,247 @@ 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( - 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( - 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, + 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, + 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( - 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, + 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, + 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( - 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; - } - } 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; + // 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; } - 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| wt_final_score(t, scorer)); - 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 = wt_mapped_length(existing) - overlap; - if u_new == 0 && wt_score < wt_final_score(existing, scorer) { - 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 = wt_final_score(t, scorer); - wt_score > score || (wt_score == score && wt_span < span(t)) - }) - .unwrap_or(transcripts.len()); - transcripts.insert(pos, wt); - transcripts.truncate(max_transcripts); + if let Some(mark) = mark { + wt.restore(&mark); } } return; @@ -2828,7 +2908,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, @@ -2864,10 +2945,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, @@ -2880,7 +2964,7 @@ fn stitch_recurse( ) { stitch_recurse( i_a + 1, - new_wt, + wt, wa_entries, read_seq, index, @@ -2897,6 +2981,7 @@ fn stitch_recurse( debug_name, ); } + wt.restore(&mark); } // EXCLUDE branch: skip wa_entries[i_a]. @@ -3023,12 +3108,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); @@ -3125,6 +3210,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, @@ -3141,6 +3228,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, }, )) } @@ -3261,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() }; @@ -3314,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 @@ -3460,11 +3408,13 @@ 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; stitch_recurse( 0, - WorkingTranscript::new(), + &mut WorkingTranscript::new(), &wa_entries, stitch_read, index, @@ -3480,6 +3430,8 @@ pub(crate) fn stitch_seeds_core( &mut jcache, debug_read_name, ); + jcache.clear(); + JUNCTION_SCAN_CACHE.set(jcache); if debug { eprintln!( @@ -3725,7 +3677,6 @@ mod tests { is_anchor: true, mate_id: 2, sj_a: None, - pre_ext_score: 5, }, WindowAlignment { seed_idx: 1, @@ -3737,7 +3688,6 @@ mod tests { is_anchor: true, mate_id: 2, sj_a: None, - pre_ext_score: 5, }, ]; @@ -3900,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( @@ -3943,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); @@ -3968,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); @@ -3993,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 @@ -4015,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) @@ -4031,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); @@ -4047,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); @@ -4063,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 @@ -4080,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); @@ -4097,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); @@ -4117,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, @@ -4316,6 +4265,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 @@ -4361,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 badff9d9..e42f3b84 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 7f6b5783..d6f2128b 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 be8b0403..2e55b4be 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 40a03e57..378844d2 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) /// @@ -526,8 +459,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) @@ -924,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( @@ -2295,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(); @@ -2312,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(); @@ -5370,6 +5262,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(); diff --git a/src/junction/gtf.rs b/src/junction/gtf.rs index ef40ce61..7109ec73 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/junction/mod.rs b/src/junction/mod.rs index 47557594..9df3bb4e 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, diff --git a/src/lib.rs b/src/lib.rs index 1d0c9932..7878bc46 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()); diff --git a/src/solo/count.rs b/src/solo/count.rs index 2acfb8ab..8b71d501 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 @@ -876,38 +855,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 +985,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)) @@ -1729,13 +1726,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, @@ -1763,7 +1777,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()); @@ -1906,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), @@ -2042,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 @@ -2050,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, @@ -2074,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; @@ -3003,6 +3053,54 @@ 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); + } + + /// 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 { diff --git a/src/solo/mod.rs b/src/solo/mod.rs index b91d1fbc..d12a3d24 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 diff --git a/tests/unmapped_buckets.rs b/tests/unmapped_buckets.rs index 617219c4..8ac7784c 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:?}" + ); +}