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