From 893f3aa5bdbb80578e263339daee93475bcfb87e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Mon, 2 Nov 2026 23:48:51 +0100 Subject: [PATCH 1/2] Junctions near a genome edge, SAindex after a dropped junction, index checks A flank past either end of the genome keeps its spacer bytes instead of stopping the run. An on-the-fly insertion that drops an old junction rebuilds the SAindex. The loader refuses short SA and SAindex files and bad chromosome tables instead of reading zeros or panicking. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 6 ++ DIVERGENCE.md | 12 +++ src/index/io.rs | 85 ++++++++++++---- src/index/sjdb_otf.rs | 46 +++++---- src/junction/sjdb_insert.rs | 49 ++++++---- tests/index_load.rs | 186 ++++++++++++++++++++++++++++++++++++ 6 files changed, 324 insertions(+), 60 deletions(-) create mode 100644 tests/index_load.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e..80b7729 100644 --- a/CHANGELOG.md +++ b/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. 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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index f0890fb..cf446d0 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -267,6 +267,18 @@ There are none at present. The four yeast PE reads once listed here now match ST --- +### 3.5 `Genome` file: one byte of a flank past the genome end + +**What STAR does.** A junction's acceptor flank that runs past the end of the genome reads STAR's padding after `G`. With a last chromosome ending at the bin boundary minus 2, `--sjdbOverhang 100` and an acceptor 97 bases from the end, the flank's last three bases lie past the end, and STAR's `Genome` holds `5 5 0` for them. + +**What rustar-aligner does.** It writes spacers for the whole missing part: `5 5 5`. + +**Why.** The `0` is not explained by STAR's output, and its source was not read for this. The padding is spacers everywhere else, and a donor flank past the genome's start matches STAR byte for byte. + +**Impact.** Only a junction whose acceptor flank overruns the genome's end by at least three bases, which used to stop the run. That byte also moves `SA` and `SAindex` entries for the suffixes crossing it. A read would have to span the end of that flank to notice. + +**Source.** `src/junction/sjdb_insert.rs` (`build_gsj`), locked by `build_gsj_pads_an_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. diff --git a/src/index/io.rs b/src/index/io.rs index a274bb3..b4fb7a4 100644 --- a/src/index/io.rs +++ b/src/index/io.rs @@ -149,26 +149,48 @@ impl GenomeIndex { } } -/// Read `genomeFileSizes\t ` 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, Error> { +/// Read `genomeFileSizes\t ` 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, Option), 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::() - { - return Ok(Some(v)); + if let Some(rest) = line.strip_prefix("genomeFileSizes\t") { + let mut values = rest.split_whitespace().map(|v| v.parse::().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, 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 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 = 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 = 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 // 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 Result 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 { diff --git a/src/index/sjdb_otf.rs b/src/index/sjdb_otf.rs index 082bff5..484d47c 100644 --- a/src/index/sjdb_otf.rs +++ b/src/index/sjdb_otf.rs @@ -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)", diff --git a/src/junction/sjdb_insert.rs b/src/junction/sjdb_insert.rs index 5781b34..30516d7 100644 --- a/src/junction/sjdb_insert.rs +++ b/src/junction/sjdb_insert.rs @@ -527,21 +527,21 @@ pub fn build_gsj( 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. 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]); // base + 2*overhang = trailing spacer, already initialized. } @@ -961,20 +961,27 @@ mod tests { } #[test] - fn build_gsj_errors_when_flank_underflows() { - // Junction too close to the left edge → donor_start 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).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); + } + + #[test] + fn build_gsj_pads_an_acceptor_flank_past_the_genome_end() { 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")); + let near = pj(0, 50, 294, 1, 0, 1); // acceptor starts at 295: five real bases + let gsj = build_gsj(&[near], &g, g.n_genome, 10).unwrap(); + assert_eq!(&gsj[10..20], &[0, 0, 0, 0, 0, 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 bad = pj(0, 50, 400, 1, 0, 1); let err = build_gsj(&[bad], &g, g.n_genome, 10).unwrap_err(); assert!(format!("{err}").contains("overruns")); } diff --git a/tests/index_load.rs b/tests/index_load.rs new file mode 100644 index 0000000..590f395 --- /dev/null +++ b/tests/index_load.rs @@ -0,0 +1,186 @@ +//! A damaged index is refused at load, and junctions near a genome edge or +//! dropped by an on-the-fly insertion leave a working index. + +use assert_cmd::cargo::cargo_bin_cmd; +use predicates::str::contains; +use std::fs; +use std::path::Path; +use tempfile::TempDir; + +fn lcg_seq(seed: u32, length: usize) -> Vec { + let mut state = seed; + (0..length) + .map(|_| { + state = state.wrapping_mul(1_103_515_245).wrapping_add(12345); + b"ACGT"[((state >> 16) & 3) as usize] + }) + .collect() +} + +fn write_genome(dir: &Path, seq: &[u8]) { + let mut fa = b">c1\n".to_vec(); + fa.extend_from_slice(seq); + fa.push(b'\n'); + fs::write(dir.join("g.fa"), fa).unwrap(); +} + +fn generate(dir: &Path, extra: &[&str]) -> assert_cmd::assert::Assert { + cargo_bin_cmd!("rustar-aligner") + .args(["--runMode", "genomeGenerate", "--genomeSAindexNbases", "8"]) + .arg("--genomeDir") + .arg(dir.join("idx")) + .arg("--genomeFastaFiles") + .arg(dir.join("g.fa")) + .arg("--outFileNamePrefix") + .arg(dir.join("gg_")) + .args(extra) + .assert() +} + +fn write_reads(dir: &Path, seq: &[u8]) { + let mut out = String::new(); + for i in 0..50 { + let s = std::str::from_utf8(&seq[100 + i * 40..200 + i * 40]).unwrap(); + out += &format!("@r{i}\n{s}\n+\n{}\n", "I".repeat(100)); + } + fs::write(dir.join("r.fq"), out).unwrap(); +} + +fn align(dir: &Path, extra: &[&str]) -> assert_cmd::assert::Assert { + cargo_bin_cmd!("rustar-aligner") + .arg("--genomeDir") + .arg(dir.join("idx")) + .arg("--readFilesIn") + .arg(dir.join("r.fq")) + .args(["--outSAMtype", "None"]) + .arg("--outFileNamePrefix") + .arg(dir.join("a_")) + .args(extra) + .assert() +} + +fn built_index() -> TempDir { + let dir = TempDir::new().unwrap(); + let seq = lcg_seq(3, 5000); + write_genome(dir.path(), &seq); + write_reads(dir.path(), &seq); + generate(dir.path(), &[]).success(); + align(dir.path(), &[]).success(); + dir +} + +fn truncate(path: &Path, drop_bytes: u64) { + let f = fs::OpenOptions::new().write(true).open(path).unwrap(); + f.set_len(f.metadata().unwrap().len() - drop_bytes).unwrap(); +} + +#[test] +fn short_sa_is_refused() { + let dir = built_index(); + truncate(&dir.path().join("idx/SA"), 40); + align(dir.path(), &[]).failure().stderr(contains("SA is")); +} + +#[test] +fn short_saindex_is_refused() { + let dir = built_index(); + truncate(&dir.path().join("idx/SAindex"), 40); + align(dir.path(), &[]) + .failure() + .stderr(contains("SAindex is")); +} + +#[test] +fn garbage_chr_length_is_refused_without_a_panic() { + let dir = built_index(); + fs::write(dir.path().join("idx/chrLength.txt"), "not-a-number\n").unwrap(); + align(dir.path(), &[]) + .failure() + .stderr(contains("not a number")); +} + +#[test] +fn short_chr_start_is_refused_without_a_panic() { + let dir = built_index(); + fs::write(dir.path().join("idx/chrStart.txt"), "0\n").unwrap(); + align(dir.path(), &[]).failure().stderr(contains("starts")); +} + +/// An exon at the very start of the genome, overhang 100: the donor flank runs +/// off the start. STAR reads its padding there and so builds the index. +#[test] +fn junction_at_the_genome_start_builds_and_maps() { + let dir = TempDir::new().unwrap(); + let seq = lcg_seq(5, 5000); + write_genome(dir.path(), &seq); + write_reads(dir.path(), &seq); + fs::write(dir.path().join("sj.tab"), "c1\t51\t499\t+\n").unwrap(); + let sj = dir.path().join("sj.tab"); + generate( + dir.path(), + &[ + "--sjdbFileChrStartEnd", + sj.to_str().unwrap(), + "--sjdbOverhang", + "100", + ], + ) + .success(); + align(dir.path(), &[]).success(); +} + +/// Two listed junctions that are one repeat collapse when the index is reloaded, +/// so an on-the-fly insertion drops one old block from the SA. The SAindex has to +/// follow: it must equal the one a fresh build of the surviving junctions writes. +#[test] +fn insertion_that_drops_an_old_junction_keeps_a_valid_saindex() { + let dir = TempDir::new().unwrap(); + let mut seq = lcg_seq(9, 40_000); + // A GT..AG intron at 5001-6000, whose flanking base repeats, so that + // 5000-5999 is the same junction shifted by one. + seq[5000..5002].copy_from_slice(b"GT"); + seq[5998..6000].copy_from_slice(b"AG"); + seq[4999] = b'G'; + write_genome(dir.path(), &seq); + write_reads(dir.path(), &seq); + + let old = dir.path().join("old.tab"); + fs::write(&old, "c1\t5001\t6000\t.\nc1\t5000\t5999\t+\n").unwrap(); + generate( + dir.path(), + &["--sjdbFileChrStartEnd", old.to_str().unwrap()], + ) + .success(); + + let new = dir.path().join("new.tab"); + fs::write(&new, "c1\t20001\t21000\t+\n").unwrap(); + align( + dir.path(), + &[ + "--sjdbFileChrStartEnd", + new.to_str().unwrap(), + "--sjdbInsertSave", + "All", + ], + ) + .success(); + + let fresh = TempDir::new().unwrap(); + write_genome(fresh.path(), &seq); + let merged = fresh.path().join("merged.tab"); + fs::write(&merged, "c1\t5001\t6000\t+\nc1\t20001\t21000\t+\n").unwrap(); + generate( + fresh.path(), + &["--sjdbFileChrStartEnd", merged.to_str().unwrap()], + ) + .success(); + + let inserted = dir.path().join("a__STARgenome"); + for name in ["SA", "SAindex"] { + assert!( + fs::read(inserted.join(name)).unwrap() + == fs::read(fresh.path().join("idx").join(name)).unwrap(), + "{name} differs from a fresh build" + ); + } +} From c18748fa1ce71f176aa5d2b19a4b366ef49bf5f1 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 09:13:45 +0100 Subject: [PATCH 2/2] Junction flank past the genome end: read STAR's reverse strand At genomeGenerate STAR keeps the reverse strand right after G (Genome_genomeGenerate.cpp:176-181) and copies flanks straight from G (sjdbPrepare.cpp:210-211), so an acceptor flank past the end reads comp(G[2n-1-p]), not spacers. That is the "5 5 0" byte DIVERGENCE.md 3.5 could not explain. build_gsj now does the same at genomeGenerate and keeps spacers at mapping, where STAR has no reverse strand in memory. Last contig of 4094 bases, exons 3000-3100 and 4000-4094, --sjdbOverhang 100: Genome, SA, sjdbList.out.tab and sjdbInfo.txt are byte-identical to STAR's with the last base T or G. SAindex is identical for T and differs in 2 bytes for G; 3.5 now describes only that. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 2 +- DIVERGENCE.md | 12 +++++------ src/index/mod.rs | 9 ++++++-- src/index/sjdb_otf.rs | 2 +- src/junction/sjdb_insert.rs | 41 ++++++++++++++++++++++++++++--------- 5 files changed, 46 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 80b7729..06f9821 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,7 +13,7 @@ Sections commonly used: Features, Bug fixes, Other changes. ### 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. See `DIVERGENCE.md` §3.5. +- 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. diff --git a/DIVERGENCE.md b/DIVERGENCE.md index cf446d0..1c7ff7e 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -267,17 +267,17 @@ There are none at present. The four yeast PE reads once listed here now match ST --- -### 3.5 `Genome` file: one byte of a flank past the genome end +### 3.5 `SAindex` after a junction flank past the genome end -**What STAR does.** A junction's acceptor flank that runs past the end of the genome reads STAR's padding after `G`. With a last chromosome ending at the bin boundary minus 2, `--sjdbOverhang 100` and an acceptor 97 bases from the end, the flank's last three bases lie past the end, and STAR's `Genome` holds `5 5 0` for them. +**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.** It writes spacers for the whole missing part: `5 5 5`. +**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.** The `0` is not explained by STAR's output, and its source was not read for this. The padding is spacers everywhere else, and a donor flank past the genome's start matches STAR byte for byte. +**Why.** Not found: the `SAindex` difference was measured, not explained. -**Impact.** Only a junction whose acceptor flank overruns the genome's end by at least three bases, which used to stop the run. That byte also moves `SA` and `SAindex` entries for the suffixes crossing it. A read would have to span the end of that flank to notice. +**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_pads_an_acceptor_flank_past_the_genome_end`. +**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) diff --git a/src/index/mod.rs b/src/index/mod.rs index badff9d..82bda34 100644 --- a/src/index/mod.rs +++ b/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(), diff --git a/src/index/sjdb_otf.rs b/src/index/sjdb_otf.rs index 484d47c..32a6c40 100644 --- a/src/index/sjdb_otf.rs +++ b/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); diff --git a/src/junction/sjdb_insert.rs b/src/junction/sjdb_insert.rs index 30516d7..217bf2b 100644 --- a/src/junction/sjdb_insert.rs +++ b/src/junction/sjdb_insert.rs @@ -515,11 +515,17 @@ 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, Error> { let overhang = sjdb_overhang as usize; let sjdb_length = 2 * overhang + 1; @@ -534,7 +540,8 @@ pub fn build_gsj( } // 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. + // 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; let donor_len = overhang.min(start); gsj[base + overhang - donor_len..base + overhang] @@ -542,6 +549,13 @@ pub fn build_gsj( 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, @@ -964,25 +978,32 @@ mod tests { 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).unwrap(); + 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_pads_an_acceptor_flank_past_the_genome_end() { - let g = make_test_genome(vec![0u8; 300]); + 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 gsj = build_gsj(&[near], &g, g.n_genome, 10).unwrap(); - assert_eq!(&gsj[10..20], &[0, 0, 0, 0, 0, 5, 5, 5, 5, 5]); + 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_the_intron_overruns_the_genome() { let g = make_test_genome(vec![0u8; 300]); let bad = pj(0, 50, 400, 1, 0, 1); - let err = build_gsj(&[bad], &g, g.n_genome, 10).unwrap_err(); + let err = build_gsj(&[bad], &g, g.n_genome, 10, true).unwrap_err(); assert!(format!("{err}").contains("overruns")); }