From dc2bd3a68f2caa8b19aef71c36d8be6a9bd3570c Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Mon, 2 Nov 2026 23:34:59 +0100 Subject: [PATCH 1/2] GTF parser: refuse an annotation that matches nothing, accept GFF3 No exon on a known chromosome is fatal, as in STAR. Attributes split on a space or `=`. An exon without a parent transcript is skipped. Exons are grouped by chromosome and transcript ID, so an ID on two chromosomes no longer produces false junctions. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 6 +++ DIVERGENCE.md | 12 +++++ src/junction/gtf.rs | 110 ++++++++++++++++++++++++++++++++++++++------ 3 files changed, 115 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e..530e781 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ## [Unreleased] +### Bug fixes + +- A GTF with no exon on a known chromosome is now fatal, as in STAR, with STAR's hint about chromosome naming. An Ensembl GTF on a UCSC FASTA used to exit 0 with an unannotated index and every read in `N_noFeature`. +- GTF attributes accept GFF3's `key=value`, so `--sjdbGTFtagExonParentTranscript Parent` works on GFF3 as in STAR's manual. An exon with no parent transcript is skipped with a warning rather than stopping the build. +- A transcript ID used on two chromosomes no longer produces false junctions and loses the real ones. See `DIVERGENCE.md` §1.9. + ### Other changes - Third-party notices completed. `LICENSES/` gains the licences of code diff --git a/DIVERGENCE.md b/DIVERGENCE.md index f0890fb..6056f6f 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -173,6 +173,18 @@ A rustar-aligner extension with no STAR equivalent, off by default (`No`). With --- +### 1.9 The same transcript ID on two chromosomes + +**What STAR does.** `GTF_transcriptGeneSJ.cpp` sorts a transcript's exons by absolute position and joins neighbours, so an ID used on two chromosomes also yields one junction joining them. With `T` on `c1` (1000-1200, 2000-2200) and `c2` (1500-1700, 2500-2700), `sjdbList.out.tab` holds the real `c1 1201-1999` and `c2 1701-2499`, plus `c1 2201 263643`, which runs off the end of `c1`. + +**What rustar-aligner does.** It treats each chromosome's exons as their own chain: the two real junctions, without the third. + +**Why.** The third junction is an artefact of the sort, not annotation, and a junction crossing chromosomes cannot be aligned to. + +**Impact.** Only GTFs that reuse a transcript ID across chromosomes (the PAR genes of refGene-style files), and only by one junction per such ID. Coordinates elsewhere in the sjdb are unaffected. + +**Source.** `src/junction/gtf.rs` (`extract_junctions_configured`), locked by `same_transcript_id_on_two_chromosomes_is_two_chains`. + ## 2. Cases where rustar-aligner outperforms STAR These are not chosen divergences and not bugs: rustar-aligner reports a **higher-scoring, correct** alignment that STAR misses. diff --git a/src/junction/gtf.rs b/src/junction/gtf.rs index ef40ce6..b1469e7 100644 --- a/src/junction/gtf.rs +++ b/src/junction/gtf.rs @@ -152,14 +152,13 @@ fn parse_attributes(attr_str: &str) -> HashMap { continue; } - // Split on first space to separate key and value - let parts: Vec<&str> = pair.splitn(2, ' ').collect(); - if parts.len() != 2 { + // Split on the first space, or `=` for GFF3 (`Parent=T1`). + let Some((key, value)) = pair.split_once([' ', '=']) else { continue; // Skip malformed attributes - } + }; - let key = parts[0].trim().to_string(); - let value = parts[1] + let key = key.trim().to_string(); + let value = value .trim() .trim_matches('"') // Remove quotes .to_string(); @@ -184,16 +183,39 @@ pub fn extract_junctions_configured( genome: &Genome, transcript_tag: &str, ) -> Result, Error> { - let mut transcripts: HashMap> = HashMap::new(); + // Keyed by chromosome too: the same ID on two chromosomes (PAR genes in + // refGene-style GTFs) is two exon chains, and one list would mix their + // chromosome-local coordinates into false introns. + let mut transcripts: HashMap<(String, String), Vec> = HashMap::new(); + let mut n_valid = 0usize; for exon in exons { - let transcript_id = exon - .attributes - .get(transcript_tag) - .ok_or_else(|| Error::Gtf(format!("Exon missing {transcript_tag} attribute")))? - .clone(); + // As STAR: an exon without a parent transcript is skipped, not fatal. + let Some(transcript_id) = exon.attributes.get(transcript_tag).cloned() else { + log::warn!( + "Skipping exon at {}:{} without a {transcript_tag} attribute", + exon.seqname, + exon.start + ); + continue; + }; + if genome.chr_name.contains(&exon.seqname) { + n_valid += 1; + } + transcripts + .entry((exon.seqname.clone(), transcript_id)) + .or_default() + .push(exon); + } - transcripts.entry(transcript_id).or_default().push(exon); + // STAR's fatal error (exit 104): an annotation that matches nothing is a + // mistake, nearly always `X` vs `chrX`, not an unannotated genome. + if n_valid == 0 { + return Err(Error::Gtf( + "no valid exon lines in the GTF file. Check its formatting. One likely cause \ + is the difference in chromosome naming between GTF and FASTA file." + .to_string(), + )); } // Extract junctions from each transcript @@ -585,6 +607,17 @@ mod tests { .into_iter() .collect(), }, + // One exon on a known chromosome, or the whole GTF would be refused. + GtfRecord { + seqname: "chr1".to_string(), + feature: "exon".to_string(), + start: 100, + end: 200, + strand: '+', + attributes: vec![("transcript_id".to_string(), "T2".to_string())] + .into_iter() + .collect(), + }, ]; let junctions = extract_junctions_from_exons(exons, &genome).unwrap(); @@ -723,4 +756,55 @@ mod tests { assert_eq!(junctions.len(), 1); assert_eq!(junctions[0], (0, 200, 298, 1)); } + + fn two_chr_genome() -> Genome { + Genome { + transform_blocks: None, + sequence: vec![0; 8000].into(), + n_genome: 8000, + n_genome_real: 8000, + n_chr_real: 2, + chr_start: vec![0, 4000, 8000], + chr_length: vec![4000, 4000], + chr_name: vec!["c1".to_string(), "c2".to_string()], + } + } + + fn exon(line: &str) -> GtfRecord { + parse_gtf_line(line).unwrap() + } + + #[test] + fn same_transcript_id_on_two_chromosomes_is_two_chains() { + let a = |c: &str, s, e| format!("{c}\tt\texon\t{s}\t{e}\t.\t+\t.\ttranscript_id \"T\";"); + let exons = vec![ + exon(&a("c1", 1000, 1200)), + exon(&a("c1", 2000, 2200)), + exon(&a("c2", 1500, 1700)), + exon(&a("c2", 2500, 2700)), + ]; + let got = extract_junctions_configured(exons, &two_chr_genome(), "transcript_id").unwrap(); + assert_eq!(got, vec![(0, 1200, 1998, 1), (1, 5700, 6498, 1)]); + } + + #[test] + fn no_exon_on_a_known_chromosome_is_fatal() { + let exons = vec![exon( + "X\tt\texon\t1000\t1200\t.\t+\t.\ttranscript_id \"T\";", + )]; + let err = extract_junctions_configured(exons, &two_chr_genome(), "transcript_id"); + assert!(err.unwrap_err().to_string().contains("chromosome naming")); + assert!(extract_junctions_configured(vec![], &two_chr_genome(), "transcript_id").is_err()); + } + + #[test] + fn exon_without_parent_is_skipped_and_gff3_attributes_parse() { + let exons = vec![ + exon("c1\tt\texon\t1000\t1200\t.\t+\t.\tParent=T;Name=x"), + exon("c1\tt\texon\t2000\t2200\t.\t+\t.\tParent=T"), + exon("c1\tt\texon\t3000\t3200\t.\t+\t.\tID=orphan"), + ]; + let got = extract_junctions_configured(exons, &two_chr_genome(), "Parent").unwrap(); + assert_eq!(got, vec![(0, 1200, 1998, 1)]); + } } From b1c16e87fdb37efb7945f5f1add8fdbc24e5f5d1 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 08:58:59 +0100 Subject: [PATCH 2/2] GTF: read attributes and parentless exons as STAR does A value is its first word, as STAR turns `=`, quotes and tabs into spaces: `Parent = rna-1` was read as `= rna-1` and split a transcript in two, and `"T=1"` and `T=2` are both transcript `T` in STAR. An exon on the genome with no parent transcript now gets STAR's invented ID and counts towards "no valid exon lines", so a GTF whose only exons on the genome lack a transcript ID builds as in STAR instead of exiting 1. An exon with no gene ID belongs to MissingGeneID. On a test GTF and GFF3 covering these cases, sjdbList.out.tab, transcriptInfo.tab, exonInfo.tab, geneInfo.tab and exonGeTrInfo.tab are byte-identical to STAR 2.7.11b's. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 2 +- src/index/io.rs | 7 +++- src/index/mod.rs | 7 +++- src/junction/gtf.rs | 77 +++++++++++++++++++++++++++++++------- src/quant/mod.rs | 9 +++-- src/quant/transcriptome.rs | 15 +++++++- 6 files changed, 95 insertions(+), 22 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 530e781..8a754dc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,7 +14,7 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes - A GTF with no exon on a known chromosome is now fatal, as in STAR, with STAR's hint about chromosome naming. An Ensembl GTF on a UCSC FASTA used to exit 0 with an unannotated index and every read in `N_noFeature`. -- GTF attributes accept GFF3's `key=value`, so `--sjdbGTFtagExonParentTranscript Parent` works on GFF3 as in STAR's manual. An exon with no parent transcript is skipped with a warning rather than stopping the build. +- GTF attributes are read as STAR reads them: GFF3's `key=value` works, so `--sjdbGTFtagExonParentTranscript Parent` works on GFF3 as in STAR's manual, and a value is its first word (`Parent = rna-1` gives `rna-1`, `"T=1"` gives `T`). An exon with no parent transcript becomes a transcript of its own named as STAR names it (`tr____`) rather than stopping the build, and one with no gene ID belongs to `MissingGeneID`, so `transcriptInfo.tab`, `geneInfo.tab` and `ReadsPerGene.out.tab` match STAR's. - A transcript ID used on two chromosomes no longer produces false junctions and loses the real ones. See `DIVERGENCE.md` §1.9. ### Other changes diff --git a/src/index/io.rs b/src/index/io.rs index a274bb3..851319b 100644 --- a/src/index/io.rs +++ b/src/index/io.rs @@ -126,11 +126,16 @@ impl GenomeIndex { prepared.extend(index.prepare_raw_junctions(&raw, 10)); } if let Some(ref gtf_path) = params.sjdb_gtf_file { - let exons = crate::junction::gtf::parse_gtf_configured( + let mut exons = crate::junction::gtf::parse_gtf_configured( gtf_path, ¶ms.sjdb_gtf_feature_exon, ¶ms.sjdb_gtf_chr_prefix, )?; + crate::junction::gtf::name_orphan_exons( + &mut exons, + &index.genome, + ¶ms.sjdb_gtf_tag_exon_parent_transcript, + ); let raw = crate::junction::gtf::extract_junctions_configured( exons, &index.genome, diff --git a/src/index/mod.rs b/src/index/mod.rs index badff9d..6fff261 100644 --- a/src/index/mod.rs +++ b/src/index/mod.rs @@ -287,11 +287,16 @@ impl GenomeIndex { let mut transcriptome = None; if let Some(ref gtf_path) = params.sjdb_gtf_file { - let exons = crate::junction::gtf::parse_gtf_configured( + let mut exons = crate::junction::gtf::parse_gtf_configured( gtf_path, ¶ms.sjdb_gtf_feature_exon, ¶ms.sjdb_gtf_chr_prefix, )?; + crate::junction::gtf::name_orphan_exons( + &mut exons, + &genome, + ¶ms.sjdb_gtf_tag_exon_parent_transcript, + ); log::debug!("Parsed {} exon features from GTF", exons.len()); let tr = TranscriptomeIndex::from_gtf_exons_configured( diff --git a/src/junction/gtf.rs b/src/junction/gtf.rs index b1469e7..d791209 100644 --- a/src/junction/gtf.rs +++ b/src/junction/gtf.rs @@ -157,18 +157,43 @@ fn parse_attributes(attr_str: &str) -> HashMap { continue; // Skip malformed attributes }; - let key = key.trim().to_string(); - let value = value - .trim() - .trim_matches('"') // Remove quotes - .to_string(); + // As STAR (`GTF.cpp`), which turns `=`, quotes and tabs into spaces: the + // value is its first word, so `Parent = T1` gives `T1` and `"T=1"` gives `T`. + let Some(value) = value.split([' ', '=', '"', '\t']).find(|w| !w.is_empty()) else { + continue; + }; - attributes.insert(key, value); + attributes.insert(key.trim().to_string(), value.to_string()); } attributes } +/// Give each exon on the genome that has no `transcript_tag` the transcript ID +/// STAR invents for it (`GTF.cpp`): `tr____`, `n` counting +/// the exons on the genome before it. Such an exon is a transcript of its own. +pub fn name_orphan_exons(exons: &mut [GtfRecord], genome: &Genome, transcript_tag: &str) { + let mut n = 0usize; + for exon in exons { + let Some(chr) = genome.chr_name.iter().position(|c| *c == exon.seqname) else { + continue; + }; + if exon.end > genome.chr_length[chr] { + continue; + } + if !exon.attributes.contains_key(transcript_tag) { + log::warn!( + "Exon at {}:{} has no {transcript_tag} attribute; it is a transcript of its own", + exon.seqname, + exon.start + ); + let id = format!("tr_{}_{}_{}_{n}", exon.seqname, exon.start, exon.end); + exon.attributes.insert(transcript_tag.to_string(), id); + } + n += 1; + } +} + /// Extract junctions from exon records, grouping by `transcript_tag`. /// /// `transcript_tag` is the GTF attribute key for the parent transcript @@ -190,7 +215,11 @@ pub fn extract_junctions_configured( let mut n_valid = 0usize; for exon in exons { - // As STAR: an exon without a parent transcript is skipped, not fatal. + // STAR counts every exon on the genome, parent or not. + if genome.chr_name.contains(&exon.seqname) { + n_valid += 1; + } + // Callers name orphans first (`name_orphan_exons`); one left is skipped. let Some(transcript_id) = exon.attributes.get(transcript_tag).cloned() else { log::warn!( "Skipping exon at {}:{} without a {transcript_tag} attribute", @@ -199,9 +228,6 @@ pub fn extract_junctions_configured( ); continue; }; - if genome.chr_name.contains(&exon.seqname) { - n_valid += 1; - } transcripts .entry((exon.seqname.clone(), transcript_id)) .or_default() @@ -798,13 +824,38 @@ mod tests { } #[test] - fn exon_without_parent_is_skipped_and_gff3_attributes_parse() { - let exons = vec![ + fn exon_without_parent_is_its_own_transcript_and_gff3_attributes_parse() { + let mut exons = vec![ + exon("X\tt\texon\t10\t20\t.\t+\t.\tID=off_genome"), exon("c1\tt\texon\t1000\t1200\t.\t+\t.\tParent=T;Name=x"), - exon("c1\tt\texon\t2000\t2200\t.\t+\t.\tParent=T"), + exon("c1\tt\texon\t2000\t2200\t.\t+\t.\tParent = T"), exon("c1\tt\texon\t3000\t3200\t.\t+\t.\tID=orphan"), ]; + name_orphan_exons(&mut exons, &two_chr_genome(), "Parent"); + assert!(!exons[0].attributes.contains_key("Parent")); + assert_eq!(exons[3].attributes["Parent"], "tr_c1_3000_3200_2"); let got = extract_junctions_configured(exons, &two_chr_genome(), "Parent").unwrap(); assert_eq!(got, vec![(0, 1200, 1998, 1)]); } + + /// STAR's exit 0 when only parentless exons are on the genome. + #[test] + fn only_parentless_exons_on_the_genome_is_not_fatal() { + let exons = vec![ + exon("X\tt\texon\t1000\t1200\t.\t+\t.\ttranscript_id \"T\";"), + exon("c1\tt\texon\t1000\t1200\t.\t+\t.\tgene_id \"G\";"), + ]; + let got = extract_junctions_configured(exons, &two_chr_genome(), "transcript_id").unwrap(); + assert_eq!(got, vec![]); + } + + /// STAR keeps the first word: `"T=1"` and `T=2` are both transcript `T`. + #[test] + fn attribute_value_is_its_first_word_as_in_star() { + let attrs = parse_attributes(r#"transcript_id "T=1"; gene_name "A B"; gene_id G"#); + assert_eq!(attrs["transcript_id"], "T"); + assert_eq!(attrs["gene_name"], "A"); + assert_eq!(attrs["gene_id"], "G"); + assert_eq!(parse_attributes("transcript_id=T=2")["transcript_id"], "T"); + } } diff --git a/src/quant/mod.rs b/src/quant/mod.rs index f3c39f6..839f538 100644 --- a/src/quant/mod.rs +++ b/src/quant/mod.rs @@ -82,10 +82,11 @@ impl GeneAnnotation { let mut gene_span: Vec> = Vec::new(); for exon in exons { - let gene_id = match exon.attributes.get(gene_tag) { - Some(id) => id.clone(), - None => continue, - }; + let gene_id = exon + .attributes + .get(gene_tag) + .cloned() + .unwrap_or_else(|| transcriptome::MISSING_GENE_ID.to_string()); let gene_idx = if let Some(&idx) = gene_id_to_idx.get(&gene_id) { idx diff --git a/src/quant/transcriptome.rs b/src/quant/transcriptome.rs index 0a02059..1a7e9ec 100644 --- a/src/quant/transcriptome.rs +++ b/src/quant/transcriptome.rs @@ -27,6 +27,8 @@ const GTF_ATTR_GENE_ID: &str = "gene_id"; /// STAR's fallback for GTF records missing `gene_biotype` /// (`source/GTF.cpp`). const MISSING_GENE_TYPE: &str = "MissingGeneType"; +/// STAR's gene ID for an exon without one (`GTF.cpp`). +pub(crate) const MISSING_GENE_ID: &str = "MissingGeneID"; /// STAR's `find_attr_list` (`GTF.cpp`): the last of `tags`, in list order, that has an entry in /// `attrs`. `--sjdbGTFtagExonParentGeneName`/`GeneType` accept several candidate keys; when more @@ -347,7 +349,11 @@ impl TranscriptomeIndex { _ => 1u8, // forward or unknown → treat as forward (STAR default) }; - let gene_id = first.attributes.get(gene_tag).cloned().unwrap_or_default(); + let gene_id = first + .attributes + .get(gene_tag) + .cloned() + .unwrap_or_else(|| MISSING_GENE_ID.to_string()); // STAR-faithful fallbacks: when the GTF record omits gene_name, // STAR's GTF.cpp writes geneAttr[ig][0] = gene_id (not empty). // When gene_biotype is omitted, it writes `MISSING_GENE_TYPE`. @@ -438,11 +444,16 @@ impl TranscriptomeIndex { gtf_path: &Path, genome: &Genome, ) -> Result { - let exons = crate::junction::gtf::parse_gtf_configured( + let mut exons = crate::junction::gtf::parse_gtf_configured( gtf_path, ¶ms.sjdb_gtf_feature_exon, ¶ms.sjdb_gtf_chr_prefix, )?; + crate::junction::gtf::name_orphan_exons( + &mut exons, + genome, + ¶ms.sjdb_gtf_tag_exon_parent_transcript, + ); Ok(Self::from_gtf_exons_configured( &exons, genome,