Merge PR #7: --quantMode TranscriptomeSAM output

- Transcriptome index built from GTF at genomeGenerate, persisted to disk
- Genome alignment projection onto transcripts (align_to_transcripts)
- Aligned.toTranscriptome.out.bam pipeline
- fix(params): don't require --sjdbGTFfile at alignReads when index has it
- fix(io/bam): lenient BAM header validator for STAR tRNA names
- fix(genome): match STAR's genomeParameters.txt byte-for-byte
- fix(quant): STAR-exact geneInfo/exonGeTrInfo/sjdbList tab formats
- Exon gains i_frag: u8 field (mirrors STAR's EX_iFrag)
- try_pair_transcripts simplified (single caller always supplies score)

Author: Psy-Fer

src/align/read_align.rs

@@ -2,10 +2,10 @@
 use crate::align::score::{AlignmentScorer, SpliceMotif};
 use crate::align::seed::Seed;
 use crate::align::stitch::{
-    cluster_seeds, finalize_transcript, split_combined_wt, stitch_seeds_core,
-    stitch_seeds_with_jdb_debug, PE_SPACER_BASE,
+    PE_SPACER_BASE, cluster_seeds, finalize_transcript, split_combined_wt, stitch_seeds_core,
+    stitch_seeds_with_jdb_debug,
 };
-use crate::align::transcript::Transcript;
+use crate::align::transcript::{Exon, Transcript};
 use crate::error::Error;
 use crate::index::GenomeIndex;
 use crate::params::{IntronMotifFilter, IntronStrandFilter, Parameters};
@@ -20,7 +20,7 @@ use std::hash::{DefaultHasher, Hash, Hasher};
 /// via rayon, so we instead fold the read name into the seed — this keeps tie
 /// breaks reproducible regardless of thread count while still honoring the
 /// user's `--runRNGseed` value.
-fn per_read_seed(run_rng_seed: u64, read_name: &str) -> u64 {
+pub(crate) fn per_read_seed(run_rng_seed: u64, read_name: &str) -> u64 {
     let mut hasher = DefaultHasher::new();
     read_name.hash(&mut hasher);
     run_rng_seed.wrapping_mul(hasher.finish().wrapping_add(1))
@@ -74,6 +74,52 @@ pub struct PairedAlignment {
     pub combined_n_match: u32,
 }
 
+impl PairedAlignment {
+    /// 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<Exon> = 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 {
+            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(),
+            read_seq: Vec::new(),
+        }
+    }
+}
+
 /// Result of paired-end alignment, covering all mapping outcomes.
 #[derive(Debug, Clone)]
 pub enum PairedAlignmentResult {
@@ -634,8 +680,13 @@ pub fn align_paired_read(
     // PE_SPACER_BASE=11 stops MMP search at the boundary — seeds never span mates.
     // For 301bp: Nstart=7 (301/50+1). Position 129 gives ≤21bp seeds (spacer at 150)
     // → less likely unique for repetitive (rDNA) sequences → fewer N² cross-copy pairings.
-    let mut combined_seeds =
-        Seed::find_seeds(&combined_read, index, params.seed_map_min, params, debug_name)?;
+    let mut combined_seeds = Seed::find_seeds(
+        &combined_read,
+        index,
+        params.seed_map_min,
+        params,
+        debug_name,
+    )?;
     // After find_seeds, ALL seeds (search_rc=false and search_rc=true) have read_pos
     // in combined_read coordinates. The RC conversion in search_direction_sparse:
     //   seed.read_pos = combined_len - p - L (where p is position in RC(combined_read))
@@ -670,13 +721,7 @@ pub fn align_paired_read(
         )?;
 
         for wt in &wts {
-            match split_combined_wt(
-                wt,
-                len1,
-                len2,
-                stitch_is_reverse,
-                scorer.align_intron_min,
-            ) {
+            match split_combined_wt(wt, len1, len2, stitch_is_reverse, scorer.align_intron_min) {
                 Some((m1_wt, m2_wt)) => {
                     let (m1_read_slice, m1_orig_rev, m2_read_slice, m2_orig_rev) =
                         if stitch_is_reverse {
@@ -706,8 +751,8 @@ pub fn align_paired_read(
                         &scorer,
                         &stitch_cluster,
                         m1_orig_rev,
-                        m1_orig_rev,   // no_left_ext = inner for reverse (orig_is_rev=true)
-                        !m1_orig_rev,  // no_right_ext = inner for forward (orig_is_rev=false)
+                        m1_orig_rev, // no_left_ext = inner for reverse (orig_is_rev=true)
+                        !m1_orig_rev, // no_right_ext = inner for forward (orig_is_rev=false)
                     ) else {
                         continue;
                     };
@@ -718,8 +763,8 @@ pub fn align_paired_read(
                         &scorer,
                         &stitch_cluster,
                         m2_orig_rev,
-                        m2_orig_rev,   // no_left_ext = inner for reverse (orig_is_rev=true)
-                        !m2_orig_rev,  // no_right_ext = inner for forward (orig_is_rev=false)
+                        m2_orig_rev, // no_left_ext = inner for reverse (orig_is_rev=true)
+                        !m2_orig_rev, // no_right_ext = inner for forward (orig_is_rev=false)
                     ) else {
                         continue;
                     };
@@ -736,9 +781,7 @@ pub fn align_paired_read(
 
                     let combined_span =
                         t1.genome_end.max(t2.genome_end) - t1.genome_start.min(t2.genome_start);
-                    let combined_score_override = Some(
-                        wt.score + scorer.genomic_length_penalty(combined_span),
-                    );
+                    let combined_wt_score = wt.score + scorer.genomic_length_penalty(combined_span);
 
                     if let Some(pair) = try_pair_transcripts(
                         &t1,
@@ -747,8 +790,7 @@ pub fn align_paired_read(
                         len2,
                         params,
                         combined_score_threshold,
-                        &scorer,
-                        combined_score_override,
+                        combined_wt_score,
                     ) {
                         joint_pairs.push(pair);
                     }
@@ -803,7 +845,6 @@ pub fn align_paired_read(
         }
     }
 
-
     // --- Decision tree: dedup, score-filter, quality-filter, then half-mapped fallback ---
 
     // Step 1: position dedup — remove exact (chr, mate1_pos, mate2_pos, strand, CIGAR) duplicates.
@@ -961,19 +1002,17 @@ pub fn align_paired_read(
     }
 
     // Half-mapped fallback: report the best-scoring single-mate transcript.
-    // STAR applies quality filter to the COMBINED read (Lread-1 = len1+len2), so use the
-    // same combined_score_threshold here. This matches STAR's behavior where a single-mate
-    // alignment that falls below the combined threshold is not output at all.
-    let thresh1 = combined_score_threshold.max(params.out_filter_score_min);
-    let thresh2 = combined_score_threshold.max(params.out_filter_score_min);
+    // STAR applies the quality filter to the COMBINED read (Lread-1 = len1+len2), so we
+    // use the same threshold for each mate here.
+    let single_mate_threshold = combined_score_threshold.max(params.out_filter_score_min);
 
     let best_m1 = single_mate1_transcripts
         .into_iter()
-        .filter(|t| t.score >= thresh1)
+        .filter(|t| t.score >= single_mate_threshold)
         .max_by_key(|t| t.score);
     let best_m2 = single_mate2_transcripts
         .into_iter()
-        .filter(|t| t.score >= thresh2)
+        .filter(|t| t.score >= single_mate_threshold)
         .max_by_key(|t| t.score);
 
     match (best_m1, best_m2) {
@@ -1031,8 +1070,7 @@ fn try_pair_transcripts(
     len2: usize,
     params: &Parameters,
     combined_score_threshold: i32,
-    scorer: &AlignmentScorer,
-    combined_wt_score_override: Option<i32>,
+    combined_wt_score: i32,
 ) -> Option<PairedAlignment> {
     // Must be same chromosome
     if t1.chr_idx != t2.chr_idx {
@@ -1067,20 +1105,6 @@ fn try_pair_transcripts(
         return None;
     }
 
-    // Combined score: use override if provided (combined-read path supplies the original
-    // stitch_recurse score + combined-span penalty, which correctly includes mismatch
-    // penalties). Fallback formula (per-mate path) undoes per-mate span penalties and
-    // applies a single combined-span penalty.
-    let combined_span = right.genome_end - left.genome_start;
-    let combined_wt_score = combined_wt_score_override.unwrap_or_else(|| {
-        let t1_span = t1.genome_end - t1.genome_start;
-        let t2_span = t2.genome_end - t2.genome_start;
-        let p1 = scorer.genomic_length_penalty(t1_span);
-        let p2 = scorer.genomic_length_penalty(t2_span);
-        let combined_p = scorer.genomic_length_penalty(combined_span);
-        t1.score + t2.score - p1 - p2 + combined_p
-    });
-
     // SCORE-GATE: reject pairs where score is below the absolute floor
     if combined_wt_score + params.out_filter_multimap_score_range < combined_score_threshold {
         return None;
@@ -1331,9 +1355,54 @@ mod tests {
             suffix_array,
             sa_index,
             junction_db: crate::junction::SpliceJunctionDb::empty(),
+            transcriptome: None,
         }
     }
 
+    #[test]
+    fn combined_transcript_for_projection_rewrites_mate2_ifrag() {
+        use crate::align::transcript::CigarOp;
+
+        let make_tr = |gs: u64, ge: u64, rs: usize, re: usize| Transcript {
+            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![CigarOp::Match((ge - gs) as u32)],
+            score: 100,
+            n_mismatch: 0,
+            n_gap: 0,
+            n_junction: 0,
+            junction_motifs: vec![],
+            junction_annotated: vec![],
+            read_seq: vec![],
+        };
+        let pair = PairedAlignment {
+            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,
+            combined_n_match: 200,
+        };
+        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);
+    }
+
     #[test]
     fn test_align_read_no_seeds() {
         let index = make_test_index();
@@ -1424,6 +1493,7 @@ mod tests {
                 genome_end: 1100,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1445,6 +1515,7 @@ mod tests {
                 genome_end: 1300,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1477,6 +1548,7 @@ mod tests {
                 genome_end: 1100,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1498,6 +1570,7 @@ mod tests {
                 genome_end: 1400,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1528,6 +1601,7 @@ mod tests {
                 genome_end: 1100,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1549,6 +1623,7 @@ mod tests {
                 genome_end: 1300,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1580,6 +1655,7 @@ mod tests {
                 genome_end: 1300,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1602,6 +1678,7 @@ mod tests {
                 genome_end: 1300,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1679,6 +1756,7 @@ mod tests {
                 genome_end: 1300,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1700,6 +1778,7 @@ mod tests {
                 genome_end: 1100,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1739,6 +1818,7 @@ mod tests {
                 genome_end: 1200,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,
@@ -1823,6 +1903,7 @@ mod tests {
                 genome_end: 1100,
                 read_start: 0,
                 read_end: 100,
+                i_frag: 0,
             }],
             cigar: vec![CigarOp::Match(100)],
             score: 100,