Merge PR #8: SJ insertion into Genome + SA at genomeGenerate

- PreparedJunction pipeline: motif detection, shift_left/shift_right
- build_gsj(): concatenated SJ flanking-sequence buffer
- sort + cross-strand dedup for prepared junctions
- Genome::append_sjdb(): extend forward genome + rebuild RC
- sjdbInfo.txt + sjdbList.out.tab writers
- Wire sjdb insertion into GenomeIndex::build
- Extended byte-parity diff harness

Author: Psy-Fer

src/align/read_align.rs

@@ -1356,6 +1356,7 @@ mod tests {
             sa_index,
             junction_db: crate::junction::SpliceJunctionDb::empty(),
             transcriptome: None,
+            prepared_junctions: Vec::new(),
         }
     }
 

src/align/score.rs

@@ -502,64 +502,15 @@ impl AlignmentScorer {
         )
     }
 
-    /// Detect splice junction motif
-    ///
-    /// # Arguments
-    /// - `donor_pos`: Position of the donor site (first base after exon)
-    /// - `intron_len`: Length of the intron
-    /// - `genome`: Genome reference
-    ///
-    /// # Returns
-    /// The detected splice motif
+    /// Detect splice junction motif (thin wrapper over the free function
+    /// so `AlignmentScorer` callers keep working).
     pub fn detect_splice_motif(
         &self,
         donor_pos: u64,
         intron_len: u32,
         genome: &Genome,
     ) -> SpliceMotif {
-        // Read 2bp donor and 2bp acceptor from the FORWARD genome
-        // Donor: donor_pos, donor_pos+1
-        // Acceptor: donor_pos+intron_len-2, donor_pos+intron_len-1
-        // Always read forward strand — motif pattern determines the strand
-        let d1 = genome.get_base(donor_pos);
-        let d2 = genome.get_base(donor_pos + 1);
-        let a1 = genome.get_base(donor_pos + intron_len as u64 - 2);
-        let a2 = genome.get_base(donor_pos + intron_len as u64 - 1);
-
-        // Check if all bases are valid
-        // A=0, C=1, G=2, T=3
-        match (d1, d2, a1, a2) {
-            (Some(d1), Some(d2), Some(a1), Some(a2)) => {
-                // Forward-strand motifs
-                // GT-AG: (2,3,0,2)
-                if d1 == 2 && d2 == 3 && a1 == 0 && a2 == 2 {
-                    return SpliceMotif::GtAg;
-                }
-                // GC-AG: (2,1,0,2)
-                if d1 == 2 && d2 == 1 && a1 == 0 && a2 == 2 {
-                    return SpliceMotif::GcAg;
-                }
-                // AT-AC: (0,3,0,1)
-                if d1 == 0 && d2 == 3 && a1 == 0 && a2 == 1 {
-                    return SpliceMotif::AtAc;
-                }
-                // Reverse-strand motifs (reverse complement on forward genome)
-                // CT-AC: (1,3,0,1) — reverse complement of GT-AG
-                if d1 == 1 && d2 == 3 && a1 == 0 && a2 == 1 {
-                    return SpliceMotif::CtAc;
-                }
-                // CT-GC: (1,3,2,1) — reverse complement of GC-AG
-                if d1 == 1 && d2 == 3 && a1 == 2 && a2 == 1 {
-                    return SpliceMotif::CtGc;
-                }
-                // GT-AT: (2,3,0,3) — reverse complement of AT-AC
-                if d1 == 2 && d2 == 3 && a1 == 0 && a2 == 3 {
-                    return SpliceMotif::GtAt;
-                }
-                SpliceMotif::NonCanonical
-            }
-            _ => SpliceMotif::NonCanonical,
-        }
+        detect_splice_motif(donor_pos, intron_len, genome)
     }
 
     /// Score a splice junction based on motif
@@ -573,6 +524,31 @@ impl AlignmentScorer {
     }
 }
 
+/// Detect splice junction motif from forward-strand bases at the intron
+/// boundaries. Stateless — exposed as a free function so both alignment
+/// scoring and `genomeGenerate` splice-junction insertion can share one
+/// truth table.
+///
+/// `donor_pos` is the 0-based position of the intron's first base on the
+/// forward strand; `intron_len` is the intron length in bases.
+pub fn detect_splice_motif(donor_pos: u64, intron_len: u32, genome: &Genome) -> SpliceMotif {
+    let d1 = genome.get_base(donor_pos);
+    let d2 = genome.get_base(donor_pos + 1);
+    let a1 = genome.get_base(donor_pos + intron_len as u64 - 2);
+    let a2 = genome.get_base(donor_pos + intron_len as u64 - 1);
+
+    // Base encoding: A=0, C=1, G=2, T=3.
+    match (d1, d2, a1, a2) {
+        (Some(2), Some(3), Some(0), Some(2)) => SpliceMotif::GtAg,
+        (Some(2), Some(1), Some(0), Some(2)) => SpliceMotif::GcAg,
+        (Some(0), Some(3), Some(0), Some(1)) => SpliceMotif::AtAc,
+        (Some(1), Some(3), Some(0), Some(1)) => SpliceMotif::CtAc,
+        (Some(1), Some(3), Some(2), Some(1)) => SpliceMotif::CtGc,
+        (Some(2), Some(3), Some(0), Some(3)) => SpliceMotif::GtAt,
+        _ => SpliceMotif::NonCanonical,
+    }
+}
+
 /// Splice junction motif types
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
 pub enum SpliceMotif {

src/align/stitch.rs

@@ -2976,6 +2976,7 @@ mod tests {
             sa_index,
             junction_db: crate::junction::SpliceJunctionDb::empty(),
             transcriptome: None,
+            prepared_junctions: Vec::new(),
         }
     }
 
@@ -3095,6 +3096,7 @@ mod tests {
             sa_index,
             junction_db: crate::junction::SpliceJunctionDb::empty(),
             transcriptome: None,
+            prepared_junctions: Vec::new(),
         }
     }
 

src/genome/mod.rs

@@ -118,6 +118,37 @@ impl Genome {
         })
     }
 
+    /// Append a splice-junction flanking-sequence buffer (`Gsj`) to the
+    /// forward genome and rebuild the reverse complement over the extended
+    /// forward range. Matches STAR's in-memory layout after
+    /// `sjdbBuildIndex.cpp:293` (`memcpy(G+chrStart[nChrReal], Gsj, nGsj)`):
+    /// the Gsj bytes live immediately after the forward real-genome bytes
+    /// (`chr_start[n_chr_real]` stays pinned at the pre-sjdb forward total,
+    /// matching STAR's `chrStart.txt` — verified against STAR docker output).
+    ///
+    /// `n_genome` grows to include Gsj. `chr_start` / `chr_length` /
+    /// `chr_name` / `n_chr_real` are NOT updated — Gsj lives outside the
+    /// chromosome accounting.
+    pub fn append_sjdb(&mut self, gsj: &[u8]) {
+        let old_n = self.n_genome;
+        let new_n = old_n + gsj.len() as u64;
+
+        let mut new_seq = vec![GENOME_SPACING_CHAR; (new_n * 2) as usize];
+        new_seq[..old_n as usize].copy_from_slice(&self.sequence[..old_n as usize]);
+        new_seq[old_n as usize..new_n as usize].copy_from_slice(gsj);
+
+        // Rebuild RC over the extended forward range (STAR stores Gsj_RC
+        // implicitly — ruSTAR keeps the explicit `[fwd | RC]` layout).
+        for i in 0..new_n as usize {
+            let base = new_seq[i];
+            let complement = if base < 4 { 3 - base } else { base };
+            new_seq[2 * new_n as usize - 1 - i] = complement;
+        }
+
+        self.sequence = new_seq;
+        self.n_genome = new_n;
+    }
+
     /// Access a base from the genome (forward or reverse strand).
     ///
     /// # Arguments
@@ -452,4 +483,49 @@ mod tests {
         assert_eq!(genome.position_to_chr(10), Some((1, 2)));
         assert_eq!(genome.position_to_chr(11), None); // padding
     }
+
+    #[test]
+    fn append_sjdb_extends_forward_and_rebuilds_rc() {
+        // Forward "ACGT" + spacer padding (bin 3 → 8 bytes padded).
+        let mut file = NamedTempFile::new().unwrap();
+        writeln!(file, ">chr1").unwrap();
+        writeln!(file, "ACGT").unwrap();
+        let params = make_params(vec![file.path().to_path_buf()], 3);
+        let mut genome = Genome::from_fasta(&params).unwrap();
+        assert_eq!(genome.n_genome, 8);
+
+        // Gsj bytes: donor "AC" + acceptor "GT" + spacer 5 (five total).
+        let gsj: Vec<u8> = vec![0, 1, 2, 3, 5];
+        let gsj_len = gsj.len();
+        genome.append_sjdb(&gsj);
+
+        // n_genome grows by gsj_len; chr_start pinned at pre-sjdb boundary.
+        assert_eq!(genome.n_genome, 8 + gsj_len as u64);
+        assert_eq!(genome.chr_start, vec![0, 8]);
+        assert_eq!(genome.n_chr_real, 1);
+
+        // Forward is [real 0..8 | gsj 8..13].
+        assert_eq!(&genome.sequence[..4], &[0, 1, 2, 3]);
+        assert_eq!(&genome.sequence[8..13], gsj.as_slice());
+
+        // RC over the extended forward range. sequence[2n-1-i] = complement(sequence[i]).
+        let new_n = genome.n_genome as usize;
+        assert_eq!(genome.sequence[2 * new_n - 1 - 8], 3); // complement of A at fwd[8]=0
+        assert_eq!(genome.sequence[2 * new_n - 1 - 12], 5); // spacer stays 5
+        assert_eq!(genome.sequence.len(), 2 * new_n);
+    }
+
+    #[test]
+    fn append_sjdb_with_empty_gsj_is_noop() {
+        let mut file = NamedTempFile::new().unwrap();
+        writeln!(file, ">chr1").unwrap();
+        writeln!(file, "ACGT").unwrap();
+        let params = make_params(vec![file.path().to_path_buf()], 3);
+        let mut genome = Genome::from_fasta(&params).unwrap();
+        let before = genome.sequence.clone();
+        let before_n = genome.n_genome;
+        genome.append_sjdb(&[]);
+        assert_eq!(genome.n_genome, before_n);
+        assert_eq!(genome.sequence, before);
+    }
 }

src/index/io.rs

@@ -90,10 +90,33 @@ impl GenomeIndex {
             sa_index,
             junction_db,
             transcriptome,
+            prepared_junctions: Vec::new(),
         })
     }
 }
 
+/// Read `genomeFileSizes\t<n_genome> <sa_size>` from genomeParameters.txt
+/// and return the first field (total genome byte count, including Gsj if
+/// sjdb was baked in). Returns `Ok(None)` if the file or line is absent,
+/// leaving the caller to fall back to the chr_start boundary.
+fn read_genome_file_size(genome_dir: &Path) -> Result<Option<u64>, Error> {
+    let path = genome_dir.join("genomeParameters.txt");
+    let contents = match std::fs::read_to_string(&path) {
+        Ok(s) => s,
+        Err(e) if e.kind() == std::io::ErrorKind::NotFound => return Ok(None),
+        Err(e) => return Err(Error::io(e, &path)),
+    };
+    for line in contents.lines() {
+        if let Some(rest) = line.strip_prefix("genomeFileSizes\t")
+            && let Some(first) = rest.split_whitespace().next()
+            && let Ok(v) = first.parse::<u64>()
+        {
+            return Ok(Some(v));
+        }
+    }
+    Ok(None)
+}
+
 /// Load genome from disk.
 fn load_genome(genome_dir: &Path, _params: &Parameters) -> Result<Genome, Error> {
     // Read chromosome metadata
@@ -119,7 +142,14 @@ fn load_genome(genome_dir: &Path, _params: &Parameters) -> Result<Genome, Error>
         .collect();
 
     let n_chr_real = chr_name.len();
-    let n_genome = chr_start[n_chr_real]; // Last entry is total size
+
+    // `chr_start[n_chr_real]` is the forward boundary of REAL chromosomes
+    // only — it stays pinned at the pre-sjdb value in STAR (`chrStart.txt`).
+    // When sjdb has been baked into the index, the total genome size
+    // (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]);
 
     // Load Genome sequence file
     let genome_path = genome_dir.join("Genome");