fix(quant): Aligned.toTranscriptome.out.bam record for record as STAR (primary draw, pair-level soft-clip budget, MAPQ, mate order, merged CIGAR, deletions, alignment order)

DIVERGENCE.md

@@ -143,11 +143,13 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST
 
 These differ in *how* a result is produced, not *what* is produced. They are documented so a reviewer chasing a discrepancy knows the mechanism differs by design.
 
-### 4.1 Transcriptome tables built on the fly
+### 4.1 Transcriptome tables built on the fly (GTF given at mapping time)
 
-For `--quantMode TranscriptomeSAM`, rustar-aligner builds the per-transcript exon map directly from the input GTF at run time, instead of loading STAR's persisted `transcriptInfo.tab` / `exonInfo.tab` files. The projection logic mirrors STAR's `Transcriptome_quantAlign.cpp`; the output (`Aligned.toTranscriptome.out.bam`) is intended to be equivalent.
+Like STAR, rustar-aligner loads the gene and transcript model for `--quantMode TranscriptomeSAM` / `GeneCounts` and STARsolo `Gene` / `GeneFull` / `Velocyto` / SmartSeq (and the `GX`/`GN` tags) from the annotation tables in `--genomeDir` (`transcriptInfo.tab`, `exonInfo.tab`, `geneInfo.tab`) when no `--sjdbGTFfile` is given at mapping time, and exits with STAR's "could not open ... geneInfo.tab" error when neither a GTF nor the tables exist.
 
-**Source.** `src/quant/transcriptome.rs`.
+When `--sjdbGTFfile` is given at mapping time, STAR always uses it: it re-inserts the annotated junctions on the fly, writes new tables to `_STARgenome/` and reads them from there. rustar-aligner builds the gene model for GeneCounts and solo directly from that GTF in memory (no `_STARgenome/` tables are written). For `--quantMode TranscriptomeSAM` it still reads the index tables when the index has them, and only parses the GTF if the index has none (legacy indexes); a GTF that differs from the one used at `genomeGenerate` is therefore not honoured for TranscriptomeSAM. The projection logic mirrors STAR's `Transcriptome_quantAlign.cpp`; the output (`Aligned.toTranscriptome.out.bam`) is intended to be equivalent.
+
+**Source.** `src/quant/transcriptome.rs`, `src/quant/mod.rs` (`resolve_gene_annotation`).
 
 ### 4.2 In-tree RNG generator
 

src/align/pe_overlap.rs

@@ -364,6 +364,7 @@ pub fn convert_merged_transcript_to_pe(
             n_junction,
             junction_motifs: out_junctions[i].iter().map(|(m, _)| *m).collect(),
             junction_annotated: out_junctions[i].iter().map(|(_, a)| *a).collect(),
+            star_order: 0,
         });
     }
 
@@ -478,6 +479,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let scorer = AlignmentScorer::from_params_minimal();
@@ -550,6 +552,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let scorer = AlignmentScorer::from_params_minimal();

src/align/read_align.rs

@@ -121,6 +121,11 @@ pub struct PairedAlignment {
     /// Combined pair score: sum of per-mate finalized scores (each includes genomic length penalty).
     /// Used for multi-mapper score-range ranking and mappedFilter quality check.
     pub combined_wt_score: i32,
+    /// Position of this pair in STAR's alignment order (`trMult`: window
+    /// creation order, then the order inside the window). The genomic output
+    /// sorts pairs by score and position; the transcriptome output follows
+    /// this rank, as STAR's does. Set once the pair list is final.
+    pub star_order: u32,
 }
 
 impl PairedAlignment {
@@ -164,6 +169,7 @@ impl PairedAlignment {
             n_junction: m1.n_junction + m2.n_junction,
             junction_motifs: Vec::new(),
             junction_annotated: Vec::new(),
+            star_order: 0,
         }
     }
 }
@@ -478,6 +484,10 @@ fn align_read_inner(
         });
     }
 
+    for (rank, t) in transcripts.iter_mut().enumerate() {
+        t.star_order = rank as u32;
+    }
+
     // Deterministic primary tie-break (score, then a fixed positional order).
     transcripts.sort_by(|a, b| {
         b.score
@@ -1034,6 +1044,7 @@ pub fn align_paired_read(
                     is_proper_pair,
                     insert_size,
                     combined_wt_score,
+                    star_order: 0,
                 });
             }
             if !converted.is_empty() {
@@ -1144,6 +1155,10 @@ pub fn align_paired_read(
         joint_pairs.retain(|pa| pa.combined_wt_score >= score_threshold);
     }
 
+    for (rank, pair) in joint_pairs.iter_mut().enumerate() {
+        pair.star_order = rank as u32;
+    }
+
     // Deterministic primary tie-break (combined score, then a fixed positional
     // order on mate1).
     joint_pairs.sort_by(|a, b| {
@@ -1353,6 +1368,7 @@ fn try_pair_transcripts(
         is_proper_pair,
         insert_size,
         combined_wt_score,
+        star_order: 0,
     })
 }
 
@@ -1650,6 +1666,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
         let pair = PairedAlignment {
             mate1_transcript: make_tr(1000, 1100, 0, 100),
@@ -1659,6 +1676,7 @@ mod tests {
             is_proper_pair: true,
             insert_size: 400,
             combined_wt_score: 200,
+            star_order: 0,
         };
         let combined = pair.combined_transcript_for_projection();
         assert_eq!(combined.exons.len(), 2);
@@ -1721,6 +1739,7 @@ mod tests {
             is_proper_pair: true,
             insert_size: 275,
             combined_wt_score: score,
+            star_order: 0,
         };
         // pairs[0] is from an earlier window: never compared.
         let mut pairs = vec![pair(66, 100), pair(75, 141), pair(66, 140), pair(66, 141)];
@@ -1755,6 +1774,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: Vec::new(),
             junction_annotated: Vec::new(),
+            star_order: 0,
         }
     }
 
@@ -1908,6 +1928,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let t2 = Transcript {
@@ -1929,6 +1950,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         // Distance = 300bp, within default limit (auto mode = unlimited)
@@ -1962,6 +1984,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let t2 = Transcript {
@@ -1983,6 +2006,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         // Distance = 400bp, exceeds limit of 100bp
@@ -2014,6 +2038,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let t2 = Transcript {
@@ -2035,6 +2060,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let tlen = calculate_insert_size(&t1, &t2);
@@ -2067,6 +2093,7 @@ mod tests {
             n_junction: 2,
             junction_motifs: vec![SpliceMotif::GtAg, SpliceMotif::CtAc], // +strand and -strand
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         // Create a transcript with consistent strand motifs (all + strand)
@@ -2089,6 +2116,7 @@ mod tests {
             n_junction: 2,
             junction_motifs: vec![SpliceMotif::GtAg, SpliceMotif::GcAg], // both + strand
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         // Note: STAR's RemoveInconsistentStrands filters transcripts where
@@ -2167,6 +2195,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let t2 = Transcript {
@@ -2188,6 +2217,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let tlen = calculate_insert_size(&t1, &t2);
@@ -2228,6 +2258,7 @@ mod tests {
             n_junction: 1,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         // Case 1: NonCanonical + unannotated → should be filtered
@@ -2313,6 +2344,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         // Test BothMapped variant
@@ -2324,6 +2356,7 @@ mod tests {
             is_proper_pair: true,
             insert_size: 200,
             combined_wt_score: 0,
+            star_order: 0,
         }));
         assert!(matches!(both, PairedAlignmentResult::BothMapped(_)));
 
@@ -2417,6 +2450,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         }
     }
 

src/align/stitch.rs

@@ -2301,6 +2301,7 @@ pub(crate) fn finalize_transcript(
         n_junction: wt.n_junction,
         junction_motifs: wt.junction_motifs.clone(),
         junction_annotated: wt.junction_annotated.clone(),
+        star_order: 0,
     })
 }
 
@@ -2387,6 +2388,27 @@ fn pe_mates_consistent(wt: &WorkingTranscript, read_seq: &[u8], scorer: &Alignme
     true
 }
 
+type NovelFilter = Option<std::rc::Rc<std::collections::HashSet<(u64, u64)>>>;
+
+thread_local! {
+    /// `--outFilterType BySJout` second stage: the novel junctions (intron
+    /// start, end; absolute 0-based) that passed the SJ filters. While set,
+    /// a transcript with an unannotated junction outside it is dropped at
+    /// finalization (STAR `stitchWindowAligns.cpp`, `outFilterBySJoutStage==2`).
+    static BYSJ_NOVEL: std::cell::RefCell<NovelFilter> = const { std::cell::RefCell::new(None) };
+}
+
+/// Run `f` with the second-stage BySJout junction filter active on this thread.
+pub(crate) fn with_bysj_novel_filter<R>(
+    novel: std::collections::HashSet<(u64, u64)>,
+    f: impl FnOnce() -> R,
+) -> R {
+    BYSJ_NOVEL.with(|c| *c.borrow_mut() = Some(std::rc::Rc::new(novel)));
+    let out = f();
+    BYSJ_NOVEL.with(|c| *c.borrow_mut() = None);
+    out
+}
+
 /// The finalization checks STAR runs before recording a transcript in its window:
 /// mate placement, exon lengths next to a junction (`alignSJoverhangMin` plus the
 /// repeat shift, `alignSJDBoverhangMin`) and the mapped length of a spliced mate
@@ -2439,6 +2461,18 @@ fn passes_finalization_filters(
             gaps.push(Gap::Other);
         }
     }
+    // BySJout 2nd stage: unannotated junctions have to be in the filtered set.
+    let novel_ok = BYSJ_NOVEL.with(|c| {
+        c.borrow().as_ref().is_none_or(|set| {
+            wt.exons.windows(2).zip(&gaps).all(|(w, gap)| {
+                !matches!(gap, Gap::Junction(false, _))
+                    || set.contains(&(w[0].genome_end, w[1].genome_start - 1))
+            })
+        })
+    });
+    if !novel_ok {
+        return false;
+    }
     // As in the finalization of `stitchWindowAligns.cpp:109-123`: an annotated
     // junction's short exon only counts when the exon sits at a mate end or next to
     // an unannotated junction.

src/align/transcript.rs

@@ -29,6 +29,10 @@ pub struct Transcript {
     pub junction_motifs: Vec<crate::align::score::SpliceMotif>,
     /// Whether each junction is annotated in the GTF (for jM +20 offset)
     pub junction_annotated: Vec<bool>,
+    /// Rank in STAR's alignment order (`trMult`: window creation order); the
+    /// genomic output sorts by score and position, the transcriptome output
+    /// follows this rank. Set when the alignment list is final.
+    pub star_order: u32,
 }
 
 /// An exon segment in a transcript.
@@ -180,6 +184,7 @@ mod tests {
             n_junction: 1,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         assert_eq!(transcript.cigar_string(), "50M100N50M");
@@ -227,6 +232,7 @@ mod tests {
             n_junction: 1,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         // Query: 45 + 3 + 2 + 50 = 100
@@ -275,6 +281,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let [left, right] = transcript.count_soft_clips();
@@ -298,6 +305,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let [left, right] = transcript.count_soft_clips();
@@ -321,6 +329,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![],
             junction_annotated: vec![],
+            star_order: 0,
         };
 
         let [left, right] = transcript.count_soft_clips();

src/chimeric/detect.rs

@@ -1394,6 +1394,7 @@ mod tests {
             n_junction: 0,
             junction_motifs: vec![SpliceMotif::GtAg],
             junction_annotated: vec![false],
+            star_order: 0,
         }
     }