Refuse a junction or annotation table cut short

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
 

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

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<String, Error> {
+    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`

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<SjdbInfoTab, Error> {
-    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<SjdbInfoTab, E
         .next()
         .and_then(|s| s.parse().ok())
         .ok_or_else(|| Error::Index(format!("{}: bad header overhang", path.display())))?;
+    // Each junction takes 2*overhang+1 bytes after the chromosomes. STAR does
+    // not check (`Genome_genomeLoad.cpp:487`) and maps reads wrongly when they
+    // disagree.
+    let gsj_bytes = genome
+        .n_genome
+        .saturating_sub(genome.chr_start[genome.n_chr_real]);
+    if gsj_bytes != n_expected as u64 * (2 * sjdb_overhang as u64 + 1) {
+        return Err(Error::Index(format!(
+            "{}: {n_expected} junctions of overhang {sjdb_overhang} do not fill the {gsj_bytes} \
+             junction bytes of the Genome; regenerate the index",
+            path.display()
+        )));
+    }
 
     let mut junctions = Vec::with_capacity(n_expected);
     for (i, line) in lines.enumerate() {
@@ -1261,8 +1274,8 @@ mod tests {
         let genome = Genome {
             transform_blocks: None,
             sequence: seq.into(),
-            n_genome: 2000,
-            n_genome_real: 2000,
+            n_genome: 1000 + 2 * 201,
+            n_genome_real: 1000,
             n_chr_real: 1,
             chr_name: vec!["I".to_string()],
             chr_length: vec![1000],

src/quant/transcriptome.rs

@@ -706,6 +706,7 @@ impl TranscriptomeIndex {
         let (ex_start_rel, ex_end_rel_incl, ex_len_cum_flat) =
             read_exon_info(&dir.join("exonInfo.tab"))?;
         let (gene_ids, gene_names, gene_biotypes) = read_gene_info(&dir.join("geneInfo.tab"))?;
+        check_exon_ge_tr_info(&dir.join("exonGeTrInfo.tab"))?;
 
         let n_tr = tr_ids.len();
         let n_exons_total = ex_start_rel.len();
@@ -1073,7 +1074,7 @@ type TrInfoColumns = (
 /// 0-based-inclusive `trEnd` to rustar-aligner's 0-based-exclusive convention.
 /// Ignores the `trEmax` column (reconstructed from `tr_end` + sort order).
 fn read_transcript_info(path: &Path) -> Result<TrInfoColumns, Error> {
-    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<TrInfoColumns, Error> {
 /// 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<ExonInfoColumns, Error> {
-    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<ExonInfoColumns, Error> {
     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::<usize>().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<GeneInfoColumns, Error> {
-    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

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"));
+}