From 0cdd00e004421c33c386b2e245e06f8a9dc9ab3d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 10:45:35 +0100 Subject: [PATCH] Refuse a junction or annotation table cut short sjdbInfo.txt, sjdbList.out.tab, transcriptInfo.tab, exonInfo.tab, geneInfo.tab and exonGeTrInfo.tab must end with a newline and hold the rows their header (or sjdbInfo.txt) says. sjdbInfo.txt's junctions must fill the Genome after the chromosomes. sjdbList.out.tab is checked only when junctions are inserted, the one time STAR reads it. STAR reads most of these silently and maps or counts wrongly, or hangs: DIVERGENCE.md 1.9. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 1 + DIVERGENCE.md | 18 ++++++++ src/index/io.rs | 52 ++++++++++++++++++++--- src/junction/sjdb_insert.rs | 19 +++++++-- src/quant/transcriptome.rs | 28 ++++++++++-- tests/index_load.rs | 85 +++++++++++++++++++++++++++++++++++++ 6 files changed, 192 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 06f9821..a3e342b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ Sections commonly used: Features, Bug fixes, Other changes. - 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. +- A junction or annotation table cut short is refused too. A last line cut mid-value in `sjdbInfo.txt`, `exonInfo.tab` or `geneInfo.tab` used to load, a short `sjdbList.out.tab` or `exonGeTrInfo.tab` went unnoticed, and `sjdbInfo.txt` was not checked against the Genome. STAR reads most of these silently and then maps or counts wrongly, or hangs. See `DIVERGENCE.md` §1.9. ### Other changes diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 1c7ff7e..d546032 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -171,6 +171,24 @@ A rustar-aligner extension with no STAR equivalent, off by default (`No`). With **Source.** `src/solo/count.rs` (`write_cellranger_summary`), `src/solo/gene.rs` (`Region`). +### 1.9 A damaged junction or annotation table is refused + +**Needs maintainer sign-off.** + +**What STAR does.** It reads the index's text tables with `>>` and does not check that they are complete, except `transcriptInfo.tab` (`Transcriptome.cpp:51-56`, exit 105). On an index whose table lost its last lines or was cut mid-line: + +- `sjdbInfo.txt` (`Genome_genomeLoad.cpp:487-519`): missing rows read as zeros; with 10 rows missing STAR segfaults. A header whose count or overhang does not match the Genome is used as is: with the overhang off by one, 210 of 400 reads map uniquely instead of 400. +- `sjdbList.out.tab`, read only to insert junctions at mapping or in two-pass mode (`sjdbInsertJunctions.cpp:15-18`): with one row missing STAR runs at 100% CPU and never finishes. +- `exonInfo.tab`, `geneInfo.tab`, `exonGeTrInfo.tab` (`Transcriptome.cpp:18-29, 62-97`): exit 0 with other numbers. With 10 rows missing, `exonInfo.tab` puts 360 reads in `Aligned.toTranscriptome.out.bam` instead of 400, `exonGeTrInfo.tab` turns 40 counted reads into `N_noFeature`, and `geneInfo.tab` changes `ReadsPerGene.out.tab`. + +**What rustar-aligner does.** It refuses such an index with a message naming the file. Each table has to end with a newline and hold as many rows as its header says. `sjdbInfo.txt`'s junctions, each `2*sjdbOverhang+1` bytes, have to fill the Genome after the chromosomes. With junction insertion, `sjdbList.out.tab` has to hold as many junctions as `sjdbInfo.txt`. `exonGeTrInfo.tab`, which rustar-aligner does not otherwise read, is checked when present. + +**Why.** STAR's results on such an index are wrong or never come, and nothing says why. An index cut short by a failed copy is a common mistake. + +**Impact.** Only damaged indices. An intact one, from STAR or from rustar-aligner, loads as before. A table that was edited by hand and lost its final newline is refused too. + +**Source.** `src/index/io.rs` (`read_index_table`, `check_sjdb_list`), `src/junction/sjdb_insert.rs` (`read_sjdb_info_tab`), `src/quant/transcriptome.rs` (`check_exon_ge_tr_info`), locked by `table_with_a_cut_last_line_is_refused`, `short_sjdb_list_is_refused_when_junctions_are_inserted` and `sjdb_info_that_disagrees_with_the_genome_is_refused` in `tests/index_load.rs`. + --- ## 2. Cases where rustar-aligner outperforms STAR diff --git a/src/index/io.rs b/src/index/io.rs index b4fb7a4..8f7c581 100644 --- a/src/index/io.rs +++ b/src/index/io.rs @@ -57,6 +57,9 @@ impl GenomeIndex { tab.junctions.len(), tab.sjdb_overhang, ); + if sjdb_insert(params) { + check_sjdb_list(genome_dir, tab.junctions.len())?; + } (tab.junctions, tab.sjdb_overhang) } else { (Vec::new(), 0) @@ -236,11 +239,9 @@ fn check_genome_version(genome_dir: &Path, params: &Parameters) -> Result<(), Er } Some(_) => log::info!("Genome version is compatible with current STAR"), } - // `sjdbInsert.yes` (`Parameters.cpp:1002-1010`). - let sjdb_insert = params.sjdb_gtf_file.is_some() - || !params.sjdb_file_chr_start_end.is_empty() - || params.twopass_mode == crate::params::TwopassMode::Basic; - if sjdb_insert && genome_dir.join("sjdbInfo.txt").exists() && value("sjdbInsertSave").is_none() + if sjdb_insert(params) + && genome_dir.join("sjdbInfo.txt").exists() + && value("sjdbInsertSave").is_none() { return Err(Error::Index( "EXITING because of FATAL ERROR: old Genome is INCOMPATIBLE with on the fly junction \ @@ -252,6 +253,47 @@ fn check_genome_version(genome_dir: &Path, params: &Parameters) -> Result<(), Er Ok(()) } +/// `sjdbInsert.yes` (`Parameters.cpp:1002-1010`). +fn sjdb_insert(params: &Parameters) -> bool { + params.sjdb_gtf_file.is_some() + || !params.sjdb_file_chr_start_end.is_empty() + || params.twopass_mode == crate::params::TwopassMode::Basic +} + +/// Read a text table of the index. Each one ends with a newline, so a missing +/// one means its last line was cut, which can change a value without any +/// parse error. +pub(crate) fn read_index_table(path: &Path) -> Result { + let contents = std::fs::read_to_string(path).map_err(|e| Error::io(e, path))?; + if !contents.ends_with('\n') { + return Err(Error::Index(format!( + "{} is empty or its last line is cut: the index is damaged or was not fully \ + copied; copy it again or regenerate it", + path.display() + ))); + } + Ok(contents) +} + +/// Junction insertion reloads the index's junctions from `sjdbList.out.tab` +/// (`sjdbInsertJunctions.cpp:15-18`). STAR hangs on a short list, so refuse +/// one that disagrees with `sjdbInfo.txt`. +fn check_sjdb_list(genome_dir: &Path, n_junctions: usize) -> Result<(), Error> { + let path = genome_dir.join("sjdbList.out.tab"); + let n = read_index_table(&path)? + .lines() + .filter(|l| !l.is_empty()) + .count(); + if n != n_junctions { + return Err(Error::Index(format!( + "{} has {n} junctions but sjdbInfo.txt has {n_junctions}: the index is damaged \ + or was not fully copied; copy it again or regenerate it", + path.display() + ))); + } + Ok(()) +} + /// The suffix-array stride the index was built with (`genomeSAsparseD` in /// genomeParameters.txt), 1 when absent. The seed search needs it: with a /// sparse SA only every `D`-th position is a suffix, so STAR searches `D` diff --git a/src/junction/sjdb_insert.rs b/src/junction/sjdb_insert.rs index 217bf2b..709fb62 100644 --- a/src/junction/sjdb_insert.rs +++ b/src/junction/sjdb_insert.rs @@ -37,7 +37,7 @@ pub struct SjdbInfoTab { /// `genome.chr_start[..n_chr_real]`. Positions falling in chromosome /// padding return an `Error::Index`. pub fn read_sjdb_info_tab(path: &Path, genome: &Genome) -> Result { - let contents = std::fs::read_to_string(path).map_err(|e| Error::io(e, path))?; + let contents = crate::index::io::read_index_table(path)?; let mut lines = contents.lines(); let header = lines @@ -52,6 +52,19 @@ pub fn read_sjdb_info_tab(path: &Path, genome: &Genome) -> Result Result { - let body = std::fs::read_to_string(path).map_err(|e| Error::io(e, path))?; + let body = crate::index::io::read_index_table(path)?; let mut lines = body.lines(); let header = lines @@ -1146,7 +1147,7 @@ fn read_transcript_info(path: &Path) -> Result { /// Read `exonInfo.tab` into flat arrays. Start/end are transcript-relative, /// 0-based; exEnd is INCLUSIVE (STAR's convention). fn read_exon_info(path: &Path) -> Result { - let body = std::fs::read_to_string(path).map_err(|e| Error::io(e, path))?; + let body = crate::index::io::read_index_table(path)?; let mut lines = body.lines(); let header = lines @@ -1199,9 +1200,30 @@ fn read_exon_info(path: &Path) -> Result { Ok((starts, ends, len_cums)) } +/// rustar-aligner rebuilds the exon-gene table from the others, but STAR reads +/// `exonGeTrInfo.tab` for gene counts (`Transcriptome.cpp:78-122`) and counts +/// nothing for the exons a short one lacks. Refuse it, as STAR would give other +/// numbers on the same index. +fn check_exon_ge_tr_info(path: &Path) -> Result<(), Error> { + if !path.exists() { + return Ok(()); + } + let body = crate::index::io::read_index_table(path)?; + let mut lines = body.lines(); + let header = lines.next().unwrap_or(""); + let rows = lines.count(); + if header.trim().parse::().ok() != Some(rows) { + return Err(Error::Index(format!( + "{}: header says '{header}' exons but found {rows}", + path.display() + ))); + } + Ok(()) +} + /// Read `geneInfo.tab`. Each record: geneID \t geneName \t geneBiotype. fn read_gene_info(path: &Path) -> Result { - let body = std::fs::read_to_string(path).map_err(|e| Error::io(e, path))?; + let body = crate::index::io::read_index_table(path)?; let mut lines = body.lines(); let header = lines diff --git a/tests/index_load.rs b/tests/index_load.rs index 590f395..2111187 100644 --- a/tests/index_load.rs +++ b/tests/index_load.rs @@ -184,3 +184,88 @@ fn insertion_that_drops_an_old_junction_keeps_a_valid_saindex() { ); } } + +/// An index built with a GTF of one three-exon transcript: two junctions and +/// the transcript tables. +fn annotated_index() -> TempDir { + let dir = TempDir::new().unwrap(); + let seq = lcg_seq(11, 5000); + write_genome(dir.path(), &seq); + write_reads(dir.path(), &seq); + let gtf: String = [(1001, 1200), (2001, 2200), (3001, 3200)] + .iter() + .map(|(s, e)| { + format!("c1\tt\texon\t{s}\t{e}\t.\t+\t.\tgene_id \"G\"; transcript_id \"T\";\n") + }) + .collect(); + fs::write(dir.path().join("a.gtf"), gtf).unwrap(); + let gtf = dir.path().join("a.gtf"); + generate( + dir.path(), + &[ + "--sjdbGTFfile", + gtf.to_str().unwrap(), + "--sjdbOverhang", + "50", + ], + ) + .success(); + dir +} + +/// A table cut mid-line still parses, with a wrong last value; STAR reads it +/// silently, except `transcriptInfo.tab`, so every one is refused. +#[test] +fn table_with_a_cut_last_line_is_refused() { + let dir = annotated_index(); + let new = dir.path().join("new.tab"); + fs::write(&new, "c1\t4001\t4100\t+\n").unwrap(); + for name in [ + "sjdbInfo.txt", + "sjdbList.out.tab", + "transcriptInfo.tab", + "exonInfo.tab", + "geneInfo.tab", + "exonGeTrInfo.tab", + ] { + let path = dir.path().join("idx").join(name); + let intact = fs::read(&path).unwrap(); + truncate(&path, 2); + // sjdbList.out.tab is read only when junctions are inserted. + align( + dir.path(), + &["--sjdbFileChrStartEnd", new.to_str().unwrap()], + ) + .failure() + .stderr(contains(format!("{name} is empty or its last line is cut"))); + fs::write(&path, intact).unwrap(); + } + align(dir.path(), &[]).success(); +} + +/// STAR reloads the index's junctions from sjdbList.out.tab to insert new ones +/// and hangs when one is missing; without insertion it never reads the file. +#[test] +fn short_sjdb_list_is_refused_when_junctions_are_inserted() { + let dir = annotated_index(); + let path = dir.path().join("idx/sjdbList.out.tab"); + let list = fs::read_to_string(&path).unwrap(); + fs::write(&path, list.lines().next().unwrap().to_string() + "\n").unwrap(); + align(dir.path(), &[]).success(); + align(dir.path(), &["--twopassMode", "Basic"]) + .failure() + .stderr(contains("has 1 junctions but sjdbInfo.txt has 2")); +} + +/// The junction count must fill the junction part of the Genome; STAR maps +/// reads wrongly when it does not. +#[test] +fn sjdb_info_that_disagrees_with_the_genome_is_refused() { + let dir = annotated_index(); + let path = dir.path().join("idx/sjdbInfo.txt"); + let info = fs::read_to_string(&path).unwrap(); + fs::write(&path, info.replacen("\t50\n", "\t49\n", 1)).unwrap(); + align(dir.path(), &[]) + .failure() + .stderr(contains("do not fill")); +}