fix: decode Gsj-region SA hits into split splice seeds

Pass-1 alignment was silently dropping SA hits that fall in the appended
splice-junction flanking buffer (Gsj). Those hits encode candidate
splice events for annotated junctions: when a read seed straddles the
donor/acceptor boundary of a Gsj slot, the matching genome bytes live
in two non-adjacent real-genome positions. STAR decodes these via
g1 >= P.nGenomeReal and feeds the (donor, acceptor) pair into the seed
pipeline; rustar-aligner was discarding them at position_to_chr (which
returns None for positions past the last real chromosome).

Plumb n_genome_real onto Genome (pinned at the pre-sjdb forward total)
and reload prepared junctions + sjdbOverhang from sjdbInfo.txt on the
index-load path. cluster_seeds now expands each SA hit through
decode_gsj_hit: real-genome hits pass through unchanged, single-flank
Gsj hits are skipped (the equivalent real-genome SA entry already
exists in the same range), and boundary-crossing Gsj hits split into
two virtual seeds with the donor-side read run and the acceptor-side
read run pointing at their respective real-genome positions. The
existing stitch DP then chains them via its splice branch.

Verified on the yeast SRR6357072 reproducer from the issue:
Number of splices: Total jumps from 366 to 631 (target ~720) and
GT/AG from 266 to 531 with non-canonical unchanged at 93. The
remaining gap is gated on PR #45's annotated-junction lookup landing.

Fixes #47

Co-Authored-By: Claude <noreply@anthropic.com>

Author: pinin4fjords

src/align/read_align.rs

@@ -1505,6 +1505,7 @@ mod tests {
         let genome = Genome {
             sequence,
             n_genome,
+            n_genome_real: n_genome,
             n_chr_real: 1,
             chr_name: vec!["chr1".to_string()],
             chr_length: vec![10],
@@ -1535,6 +1536,7 @@ mod tests {
             junction_db: crate::junction::SpliceJunctionDb::empty(),
             transcriptome: None,
             prepared_junctions: Vec::new(),
+            sjdb_overhang: 0,
         }
     }
 

src/align/score.rs

@@ -637,6 +637,7 @@ mod tests {
         Genome {
             sequence,
             n_genome,
+            n_genome_real: n_genome,
             n_chr_real: 1,
             chr_name: vec!["chr1".to_string()],
             chr_length: vec![seq.len() as u64],

src/align/stitch.rs

@@ -1000,6 +1000,7 @@ mod tests {
         Genome {
             sequence,
             n_genome,
+            n_genome_real: n_genome,
             n_chr_real: 2,
             chr_name: vec!["chr0".to_string(), "chr1".to_string()],
             chr_length: vec![chr_len, chr_len],
@@ -1027,6 +1028,7 @@ mod tests {
             junction_db: SpliceJunctionDb::empty(),
             transcriptome: None,
             prepared_junctions: Vec::new(),
+            sjdb_overhang: 0,
         }
     }
 

src/chimeric/detect.rs

@@ -443,6 +443,7 @@ mod tests {
         Genome {
             sequence: vec![0u8; 2048],
             n_genome: 1024,
+            n_genome_real: 1024,
             n_chr_real: 2,
             chr_name: vec!["chr9".to_string(), "chr22".to_string()],
             chr_length: vec![512, 512],

src/chimeric/output.rs

@@ -173,6 +173,7 @@ mod tests {
         Genome {
             sequence: seq,
             n_genome: 100,
+            n_genome_real: 100,
             n_chr_real: 1,
             chr_name: vec!["chr1".to_string()],
             chr_start: vec![0],

src/chimeric/score.rs

@@ -25,6 +25,12 @@ pub struct Genome {
     /// Total length of the forward (padded) genome.
     pub n_genome: u64,
 
+    /// Length of the real-genome forward strand, excluding any appended
+    /// Gsj splice-junction buffer. Equals `n_genome` before `append_sjdb`
+    /// and stays pinned at that value afterwards; matches STAR's
+    /// `mapGen.nGenomeReal`.
+    pub n_genome_real: u64,
+
     /// Number of real chromosomes (not including scaffold/contigs if excluded).
     pub n_chr_real: usize,
 
@@ -111,6 +117,7 @@ impl Genome {
         Ok(Genome {
             sequence,
             n_genome,
+            n_genome_real: n_genome,
             n_chr_real,
             chr_name,
             chr_length,

src/genome/mod.rs

@@ -11,6 +11,7 @@ use crate::index::packed_array::PackedArray;
 use crate::index::sa_index::SaIndex;
 use crate::index::suffix_array::SuffixArray;
 use crate::junction::SpliceJunctionDb;
+use crate::junction::sjdb_insert;
 use crate::params::Parameters;
 use crate::quant::transcriptome::TranscriptomeIndex;
 
@@ -99,13 +100,32 @@ impl GenomeIndex {
             );
         }
 
+        // Reload prepared junctions + sjdbOverhang from sjdbInfo.txt when
+        // present. STAR appends a Gsj flanking-sequence buffer to the
+        // genome at build time; align-time code needs the parsed junctions
+        // to decode SA hits that land inside that buffer back to real
+        // `(donor, acceptor)` genome positions.
+        let sjdb_info_path = genome_dir.join("sjdbInfo.txt");
+        let (prepared_junctions, sjdb_overhang) = if sjdb_info_path.exists() {
+            let tab = sjdb_insert::read_sjdb_info_tab(&sjdb_info_path, &genome)?;
+            log::info!(
+                "Loaded sjdbInfo.txt: {} junctions, sjdbOverhang={}",
+                tab.junctions.len(),
+                tab.sjdb_overhang,
+            );
+            (tab.junctions, tab.sjdb_overhang)
+        } else {
+            (Vec::new(), 0)
+        };
+
         Ok(GenomeIndex {
             genome,
             suffix_array,
             sa_index,
             junction_db,
             transcriptome,
-            prepared_junctions: Vec::new(),
+            prepared_junctions,
+            sjdb_overhang,
         })
     }
 }
@@ -164,7 +184,8 @@ fn load_genome(genome_dir: &Path, _params: &Parameters) -> Result<Genome, Error>
     // (real + Gsj) lives in `genomeParameters.txt` under `genomeFileSizes`.
     // Prefer that value; fall back to the chr_start boundary for indices
     // built without a GTF.
-    let n_genome = read_genome_file_size(genome_dir)?.unwrap_or(chr_start[n_chr_real]);
+    let n_genome_real = chr_start[n_chr_real];
+    let n_genome = read_genome_file_size(genome_dir)?.unwrap_or(n_genome_real);
 
     // Load Genome sequence file
     let genome_path = genome_dir.join("Genome");
@@ -192,6 +213,7 @@ fn load_genome(genome_dir: &Path, _params: &Parameters) -> Result<Genome, Error>
     Ok(Genome {
         sequence,
         n_genome,
+        n_genome_real,
         n_chr_real,
         chr_name,
         chr_length,

src/index/io.rs

@@ -27,11 +27,16 @@ pub struct GenomeIndex {
     /// `transcriptInfo.tab` + friends at `genomeGenerate`, reloaded at
     /// `alignReads` from the same files.
     pub transcriptome: Option<TranscriptomeIndex>,
-    /// Prepared splice junctions (sorted/deduped) used only on the build
-    /// path to write `sjdbInfo.txt` + `sjdbList.out.tab`. Empty on the
-    /// load path — those files have already been written and are not
-    /// needed at align time.
+    /// Prepared splice junctions in their post-dedup order (the same
+    /// order they occupy in the Gsj buffer appended to the genome).
+    /// Populated on both build and load paths when sjdb is present;
+    /// empty for indices built without a GTF. Used at align time to
+    /// decode Gsj-region SA hits back to real-genome `(donor, acceptor)`
+    /// pairs.
     pub prepared_junctions: Vec<PreparedJunction>,
+    /// `sjdbOverhang` recorded in `sjdbInfo.txt`. Zero when no sjdb
+    /// junctions are present.
+    pub sjdb_overhang: u32,
 }
 
 impl GenomeIndex {
@@ -143,13 +148,20 @@ impl GenomeIndex {
             sa_index.data.len()
         );
 
+        let sjdb_overhang = if prepared_junctions.is_empty() {
+            0
+        } else {
+            params.sjdb_overhang
+        };
+
         Ok(GenomeIndex {
             genome,
             suffix_array,
             sa_index,
             junction_db,
             transcriptome,
             prepared_junctions,
+            sjdb_overhang,
         })
     }
 

src/index/mod.rs

@@ -469,6 +469,7 @@ mod tests {
         Genome {
             sequence: vec![0, 1, 2, 3, 0, 1, 2, 3], // ACGTACGT
             n_genome: 8,
+            n_genome_real: 8,
             n_chr_real: 1,
             chr_name: vec!["chr1".to_string()],
             chr_length: vec![8],