Junctions near a genome edge, SAindex after a dropped junction, index checks at load

CHANGELOG.md

@@ -11,6 +11,12 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ## [Unreleased]
 
+### Bug fixes
+
+- A junction near the genome's start or end no longer stops the run. `sjdb flanking window ... donor underflows` was fatal, for a GTF at `genomeGenerate` and for a novel junction after pass 1 of a two-pass run, while the same flank on a later contig simply read spacers. A flank before the genome's start reads spacers, and one past its end reads what STAR reads there: the reverse strand at `genomeGenerate`, spacers at mapping. See `DIVERGENCE.md` §3.5.
+- On-the-fly insertion that drops an old junction rebuilds the SAindex. It used to keep stale values, which misdirected seeds and could panic mid-alignment.
+- A damaged index is refused at load. A short `SA` or `SAindex` used to read as zeros or panic mid-run, and garbage in `chrLength.txt` or `chrStart.txt` panicked.
+
 ### Other changes
 
 - Third-party notices completed. `LICENSES/` gains the licences of code

DIVERGENCE.md

@@ -267,6 +267,18 @@ There are none at present. The four yeast PE reads once listed here now match ST
 
 ---
 
+### 3.5 `SAindex` after a junction flank past the genome end
+
+**What STAR does.** At `genomeGenerate`, STAR keeps the reverse strand right after `G` (`Genome_genomeGenerate.cpp:176-181`) and copies junction flanks straight from `G` (`sjdbPrepare.cpp:210-211`), so an acceptor flank running past the end of the genome reads `comp(G[2n-1-p])`: the padding, then the complement of the last bases. At mapping time that strand is not in memory and the same flank reads spacers.
+
+**What rustar-aligner does.** The same (`build_gsj`, `mirror_past_end`). On a last contig of 4094 bases with exons 3000-3100 and 4000-4094, `--sjdbOverhang 100 --genomeChrBinNbits 12`, `Genome`, `SA`, `sjdbList.out.tab` and `sjdbInfo.txt` are byte-identical to STAR's whether the last base is `T` or `G`. `SAindex` is identical for `T`; for `G`, 2 of its bytes differ.
+
+**Why.** Not found: the `SAindex` difference was measured, not explained.
+
+**Impact.** Only a junction whose acceptor flank overruns the genome's end, which used to stop the run. The two bytes are SAindex entries for prefixes that cross the end of that flank.
+
+**Source.** `src/junction/sjdb_insert.rs` (`build_gsj`), locked by `build_gsj_acceptor_flank_past_the_genome_end`.
+
 ## 4. Implementation divergences (no intended output difference)
 
 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.

src/index/io.rs

@@ -149,26 +149,48 @@ impl GenomeIndex {
     }
 }
 
-/// 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> {
+/// Read `genomeFileSizes\t<n_genome> <sa_size>` from genomeParameters.txt.
+/// Either value is `None` when absent; a recorded SA size of 0 means unknown.
+fn read_genome_file_sizes(genome_dir: &Path) -> Result<(Option<u64>, 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) if e.kind() == std::io::ErrorKind::NotFound => return Ok((None, 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));
+        if let Some(rest) = line.strip_prefix("genomeFileSizes\t") {
+            let mut values = rest.split_whitespace().map(|v| v.parse::<u64>().ok());
+            let n_genome = values.next().flatten();
+            let sa_size = values.next().flatten().filter(|&v| v != 0);
+            return Ok((n_genome, sa_size));
         }
     }
-    Ok(None)
+    Ok((None, None))
+}
+
+/// One integer per line, as `chrLength.txt` and `chrStart.txt` hold.
+fn parse_u64_lines(contents: &str, path: &Path) -> Result<Vec<u64>, Error> {
+    contents
+        .lines()
+        .map(|s| {
+            s.parse().map_err(|_| {
+                Error::Index(format!(
+                    "{}: '{s}' is not a number; regenerate the index",
+                    path.display()
+                ))
+            })
+        })
+        .collect()
+}
+
+/// A short index file reads as zeros or panics mid-run, so refuse it at load.
+fn index_file_size_err(path: &Path, actual: u64, expected: u64) -> Error {
+    Error::Index(format!(
+        "{} is {actual} bytes, expected {expected}: the index is damaged or was not fully \
+         copied; copy it again or regenerate it",
+        path.display()
+    ))
 }
 
 /// The `versionGenome` an index has to carry (STAR 2.7.11b's `versionGenome`
@@ -260,20 +282,24 @@ fn load_genome(genome_dir: &Path, _params: &Parameters) -> Result<Genome, Error>
     let chr_length_path = genome_dir.join("chrLength.txt");
     let chr_length_contents =
         std::fs::read_to_string(&chr_length_path).map_err(|e| Error::io(e, &chr_length_path))?;
-    let chr_length: Vec<u64> = chr_length_contents
-        .lines()
-        .map(|s| s.parse().unwrap())
-        .collect();
+    let chr_length = parse_u64_lines(&chr_length_contents, &chr_length_path)?;
 
     let chr_start_path = genome_dir.join("chrStart.txt");
     let chr_start_contents =
         std::fs::read_to_string(&chr_start_path).map_err(|e| Error::io(e, &chr_start_path))?;
-    let chr_start: Vec<u64> = chr_start_contents
-        .lines()
-        .map(|s| s.parse().unwrap())
-        .collect();
+    let chr_start = parse_u64_lines(&chr_start_contents, &chr_start_path)?;
 
     let n_chr_real = chr_name.len();
+    if chr_length.len() != n_chr_real || chr_start.len() != n_chr_real + 1 {
+        return Err(Error::Index(format!(
+            "{}: {n_chr_real} chromosomes but {} lengths and {} starts (expected {}); \
+             regenerate the index",
+            genome_dir.display(),
+            chr_length.len(),
+            chr_start.len(),
+            n_chr_real + 1
+        )));
+    }
 
     // `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`).
@@ -282,7 +308,9 @@ fn load_genome(genome_dir: &Path, _params: &Parameters) -> Result<Genome, Error>
     // Prefer that value; fall back to the chr_start boundary for indices
     // built without a GTF.
     let n_genome_real = chr_start[n_chr_real];
-    let n_genome = read_genome_file_size(genome_dir)?.unwrap_or(n_genome_real);
+    let n_genome = read_genome_file_sizes(genome_dir)?
+        .0
+        .unwrap_or(n_genome_real);
 
     // Memory-map the Genome sequence file (forward strand only, `n_genome`
     // bytes). The reverse-complement half is computed on access by
@@ -351,6 +379,12 @@ fn load_suffix_array(genome_dir: &Path, genome: &Genome) -> Result<SuffixArray,
     let mmap = unsafe { memmap2::Mmap::map(&file).map_err(|e| Error::io(e, &sa_path))? };
     advise_random(&mmap);
 
+    if let (_, Some(expected)) = read_genome_file_sizes(genome_dir)?
+        && mmap.len() as u64 != expected
+    {
+        return Err(index_file_size_err(&sa_path, mmap.len() as u64, expected));
+    }
+
     let gstrand_bit = SuffixArray::calculate_gstrand_bit(genome.n_genome);
     let word_length = gstrand_bit + 1;
     let gstrand_mask = (1u64 << gstrand_bit) - 1;
@@ -419,6 +453,15 @@ fn load_sa_index(genome_dir: &Path, gstrand_bit: u32) -> Result<SaIndex, Error>
     let word_length = gstrand_bit + 3;
     let num_indices = SaIndex::calculate_num_indices(nbases);
 
+    let expected = PackedArray::data_byte_len_for(word_length, num_indices as usize);
+    if mmap.len() < expected {
+        return Err(index_file_size_err(
+            &sai_path,
+            header_len + mmap.len() as u64,
+            header_len + expected as u64,
+        ));
+    }
+
     let data = PackedArray::from_mmap(word_length, num_indices as usize, mmap);
 
     Ok(SaIndex {

src/index/mod.rs

@@ -355,8 +355,13 @@ impl GenomeIndex {
                 .collect();
             let prepared = sjdb_insert::sort_and_dedup(prepared);
 
-            let gsj =
-                sjdb_insert::build_gsj(&prepared, &genome, n_genome_real, params.sjdb_overhang)?;
+            let gsj = sjdb_insert::build_gsj(
+                &prepared,
+                &genome,
+                n_genome_real,
+                params.sjdb_overhang,
+                true,
+            )?;
             log::info!(
                 "Built Gsj buffer: {} junctions × {} bytes = {} bytes",
                 prepared.len(),

src/index/sjdb_otf.rs

@@ -239,7 +239,7 @@ pub fn insert_junctions(
         chr_length: genome.chr_length.clone(),
         chr_start: genome.chr_start.clone(),
     };
-    let gsj = sjdb_insert::build_gsj(&all, &real_genome, n_real, overhang)?;
+    let gsj = sjdb_insert::build_gsj(&all, &real_genome, n_real, overhang, false)?;
     // STAR's working `Gsj` array: forward blocks, their reverse complement,
     // then a terminating spacer (`sjdbBuildIndex.cpp:31-40`).
     let mut gsj2 = Vec::with_capacity(2 * gsj.len() + 1);
@@ -449,24 +449,34 @@ pub fn insert_junctions(
 
     let mut new_genome = real_genome;
     new_genome.append_sjdb(&gsj);
-    let mut new_sai = std::mem::replace(
-        &mut index.sa_index,
-        SaIndex {
-            nbases: 0,
-            genome_sa_index_start: Vec::new(),
-            data: PackedArray::new(1, 0),
-            word_length: 1,
-            gstrand_bit,
-            sparse_d: 1,
-        },
-    );
-    // A loaded SAindex is memory-mapped: take an owned copy to update.
-    new_sai.data = PackedArray::from_bytes(
-        new_sai.data.word_length(),
-        new_sai.data.len(),
-        new_sai.data.data().to_vec(),
-    );
-    update_sa_index(&mut new_sai, &inserts, &gsj2);
+    let new_sai = if dropped_old {
+        // `update_sa_index` only ever moves entries up. Removed suffixes move
+        // them down too, which it has no way to express, so build the table anew.
+        let mut sai = SaIndex::build(&new_genome, &new_sa, index.sa_index.nbases)?;
+        sai.sparse_d = index.sa_index.sparse_d;
+        sai
+    } else {
+        let mut new_sai = std::mem::replace(
+            &mut index.sa_index,
+            SaIndex {
+                nbases: 0,
+                genome_sa_index_start: Vec::new(),
+                data: PackedArray::new(1, 0),
+                word_length: 1,
+                gstrand_bit,
+                sparse_d: 1,
+            },
+        );
+        // A loaded SAindex is memory-mapped: take an owned copy to update.
+        new_sai.data = PackedArray::from_bytes(
+            new_sai.data.word_length(),
+            new_sai.data.len(),
+            new_sai.data.data().to_vec(),
+        );
+        update_sa_index(&mut new_sai, &inserts, &gsj2);
+
+        new_sai
+    };
 
     log::info!(
         "Inserted {} new junctions on the fly ({} total, {} SA entries)",

src/junction/sjdb_insert.rs

@@ -515,33 +515,47 @@ fn merge_cross_strand(
 /// Callers supply junctions in sorted/deduplicated order (via
 /// [`sort_and_dedup`]). `genome` is read through
 /// `Genome::sequence[..n_genome_real]`.
+///
+/// `mirror_past_end`: at `genomeGenerate` STAR keeps the reverse strand right
+/// after `G` (`Genome_genomeGenerate.cpp:176-181`), so an acceptor flank past
+/// the end reads `comp(G[2n-1-p])`. At mapping it is not there and STAR reads
+/// spacers, as for a donor flank before the start.
 pub fn build_gsj(
     junctions: &[PreparedJunction],
     genome: &Genome,
     n_genome_real: u64,
     sjdb_overhang: u32,
+    mirror_past_end: bool,
 ) -> Result<Vec<u8>, Error> {
     let overhang = sjdb_overhang as usize;
     let sjdb_length = 2 * overhang + 1;
     let forward = &genome.sequence.as_slice()[..n_genome_real as usize];
     let mut gsj = vec![GSJ_SPACING; junctions.len() * sjdb_length];
 
     for (i, pj) in junctions.iter().enumerate() {
-        let donor_start = pj
-            .original_start()
-            .checked_sub(overhang as u64)
-            .ok_or_else(|| sjdb_bounds_err(pj, overhang, "donor underflows"))?;
-        let acceptor_start = pj.original_end() + 1;
-
-        let d0 = donor_start as usize;
-        let a0 = acceptor_start as usize;
-        if d0 + overhang > forward.len() || a0 + overhang > forward.len() {
-            return Err(sjdb_bounds_err(pj, overhang, "flank overruns genome"));
+        let start = pj.original_start() as usize;
+        let a0 = pj.original_end() as usize + 1;
+        if a0 > forward.len() {
+            return Err(sjdb_bounds_err(pj, overhang, "intron overruns genome"));
         }
 
+        // A flank running past either end of the genome keeps its spacer bytes
+        // for the missing part, what STAR reads from the padding around G,
+        // except past the end at genomeGenerate (`mirror_past_end`).
         let base = i * sjdb_length;
-        gsj[base..base + overhang].copy_from_slice(&forward[d0..d0 + overhang]);
-        gsj[base + overhang..base + 2 * overhang].copy_from_slice(&forward[a0..a0 + overhang]);
+        let donor_len = overhang.min(start);
+        gsj[base + overhang - donor_len..base + overhang]
+            .copy_from_slice(&forward[start - donor_len..start]);
+        let acceptor_len = overhang.min(forward.len() - a0);
+        gsj[base + overhang..base + overhang + acceptor_len]
+            .copy_from_slice(&forward[a0..a0 + acceptor_len]);
+        if mirror_past_end {
+            let n = forward.len();
+            for k in acceptor_len..overhang.min(2 * n - a0) {
+                let b = forward[2 * n - 1 - (a0 + k)];
+                gsj[base + overhang + k] = if b < 4 { 3 - b } else { b };
+            }
+        }
         // base + 2*overhang = trailing spacer, already initialized.
     }
 
@@ -893,7 +907,7 @@ mod tests {
         f[203] = 1;
         let g = make_test_genome(f);
         let junction = pj(0, 100, 200, 1, 0, 1);
-        let gsj = build_gsj(&[junction], &g, g.n_genome, 3).unwrap();
+        let gsj = build_gsj(&[junction], &g, g.n_genome, 3, true).unwrap();
         assert_eq!(gsj, vec![0, 1, 2, 3, 0, 1, GSJ_SPACING]);
     }
 
@@ -917,7 +931,7 @@ mod tests {
         // start_pos/end_pos are SHIFTED (stored form for motif==0).
         // With shift_left=2, original_start = 102, original_end = 199.
         let junction = pj(0, 100, 197, 0, 2, 0);
-        let gsj = build_gsj(&[junction], &g, g.n_genome, 4).unwrap();
+        let gsj = build_gsj(&[junction], &g, g.n_genome, 4, true).unwrap();
         // Donor: forward[98..102] = [0,1,2,3]
         // Acceptor: forward[200..204] = [2,3,0,1]
         assert_eq!(gsj, vec![0, 1, 2, 3, 2, 3, 0, 1, GSJ_SPACING]);
@@ -941,7 +955,7 @@ mod tests {
         let g = make_test_genome(f);
         let a = pj(0, 100, 200, 1, 0, 1);
         let b = pj(0, 300, 350, 1, 0, 1);
-        let gsj = build_gsj(&[a, b], &g, g.n_genome, 2).unwrap();
+        let gsj = build_gsj(&[a, b], &g, g.n_genome, 2, true).unwrap();
         // sjdb_length = 5 (2 donor + 2 acceptor + 1 spacer) per junction.
         assert_eq!(
             gsj,
@@ -961,21 +975,35 @@ mod tests {
     }
 
     #[test]
-    fn build_gsj_errors_when_flank_underflows() {
-        // Junction too close to the left edge → donor_start underflows.
-        let g = make_test_genome(vec![0u8; 300]);
-        let bad = pj(0, 3, 200, 1, 0, 1); // original_start=3; overhang 10 underflows
-        let err = build_gsj(&[bad], &g, g.n_genome, 10).unwrap_err();
-        assert!(format!("{err}").contains("underflows"));
+    fn build_gsj_pads_a_donor_flank_before_the_genome_start() {
+        let g = make_test_genome((0..300).map(|i| (i % 4) as u8).collect());
+        let near = pj(0, 3, 200, 1, 0, 1); // original_start=3, overhang 10
+        let gsj = build_gsj(&[near], &g, g.n_genome, 10, true).unwrap();
+        // Seven spacers, then bases 0..3, then the acceptor from original_end+1.
+        assert_eq!(&gsj[..10], &[5, 5, 5, 5, 5, 5, 5, 0, 1, 2]);
+        assert_eq!(gsj[10], (201 % 4) as u8);
+    }
+
+    /// Past the end STAR reads its reverse strand at genomeGenerate, `comp(G[2n-1-p])`,
+    /// and spacers at mapping.
+    #[test]
+    fn build_gsj_acceptor_flank_past_the_genome_end() {
+        let g = make_test_genome((0..300).map(|i| (i % 4) as u8).collect());
+        let near = pj(0, 50, 294, 1, 0, 1); // acceptor starts at 295: five real bases
+        let real = [3, 0, 1, 2, 3];
+        let generate = build_gsj(std::slice::from_ref(&near), &g, g.n_genome, 10, true).unwrap();
+        assert_eq!(&generate[10..15], &real);
+        assert_eq!(&generate[15..20], &[0, 1, 2, 3, 0]); // comp of bases 299, 298, ..., 295
+        let mapping = build_gsj(&[near], &g, g.n_genome, 10, false).unwrap();
+        assert_eq!(&mapping[10..15], &real);
+        assert_eq!(&mapping[15..20], &[5, 5, 5, 5, 5]);
     }
 
     #[test]
-    fn build_gsj_errors_when_flank_overruns_genome() {
-        // Acceptor overruns the right edge.
+    fn build_gsj_errors_when_the_intron_overruns_the_genome() {
         let g = make_test_genome(vec![0u8; 300]);
-        // original_end = 295, overhang = 10 → acceptor region [296..306), exceeds n_genome=300.
-        let bad = pj(0, 100, 295, 1, 0, 1);
-        let err = build_gsj(&[bad], &g, g.n_genome, 10).unwrap_err();
+        let bad = pj(0, 50, 400, 1, 0, 1);
+        let err = build_gsj(&[bad], &g, g.n_genome, 10, true).unwrap_err();
         assert!(format!("{err}").contains("overruns"));
     }