[BenjaminDEMAILLE]

Part of #354 (findings 6 and 7). Independent of the other open PRs apart from CHANGELOG.md and DIVERGENCE.md.

Problem

extract_junctions_configured (src/junction/gtf.rs) turns GTF exon lines into junctions. Both the genome build (src/index/mod.rs) and mapping-time insertion (src/index/io.rs) call it. It failed in four ways, in both directions.

An annotation that matches nothing exits 0. A transcript on an unknown chromosome was skipped with a log::warn, and zero junctions was treated as "no annotation". An Ensembl GTF (X) on a UCSC FASTA (chrX) gave exit 0, no sjdbList.out.tab, and at mapping, with --quantMode GeneCounts, every read in N_noFeature. --sjdbGTFfile /dev/null and a GTF with no exon lines did the same.

GFF3 could not be read. parse_attributes split only on a space, so Parent=T was dropped and --sjdbGTFtagExonParentTranscript Parent stopped with Exon missing Parent attribute.

One exon with no parent stopped the build. Exon missing transcript_id attribute was fatal, while the transcriptome builder in the same run skips such an exon (src/quant/transcriptome.rs:245-248), so the two paths disagreed.

A transcript ID on two chromosomes invented junctions. Exons were grouped by ID alone, sorted by chromosome-local start, and all placed on the first exon's chromosome.

What STAR does

Run with Homebrew STAR 2.7.11b, with --sjdbOverhang 50 on two 20 kb contigs c1 and c2:

GTF	STAR
Only X exons	exit 104, Fatal INPUT FILE error, no valid exon lines in the GTF file, with a solution line about chromosome naming
GFF3, Parent=T	exit 0, c1 1201 1999
One exon without transcript_id	exit 0, that exon ignored
T on c1 (1000-1200, 2000-2200) and c2 (1500-1700, 2500-2700)	exit 0, c1 1201 1999, c2 1701 2499, and c1 2201 263643

I did not read STAR's source for this PR; these are observed outputs. The file reference in the DIVERGENCE.md entry comes from existing comments in gtf.rs.

Change

All in src/junction/gtf.rs:

• parse_attributes splits on the first space or =, and the value is its first word, as STAR reads it (GTF.cpp turns =, quotes and tabs into spaces): Parent = rna-1 gives rna-1, "T=1" and T=2 both give T, gene_name "A B" gives A.
• extract_junctions_configured counts every exon on a known chromosome, with or without a parent, as STAR's exonN does. With none it returns Error::Gtf, using STAR's wording. This also covers an empty file and a file with no exon lines. A file whose exons on the genome all lack the parent tag is accepted, as in STAR.
• An exon on the genome without the parent tag gets a log::warn and STAR's invented transcript ID, tr_<chr>_<start>_<end>_<n> (name_orphan_exons, called before both the junction extraction and the transcriptome builder, at genome generation and at mapping). It is a transcript of its own, so it adds no junction but appears in transcriptInfo.tab as in STAR. An exon with no gene ID belongs to MissingGeneID, as in STAR.
• Exons are grouped by (seqname, transcript_id).

Deliberate divergence: DIVERGENCE.md §1.9

STAR's third junction in the last row above, c1 2201 263643, runs off the end of c1. It comes from joining the two chromosomes' exons in one chain. This PR does not reproduce it, and the entry says why. This needs maintainer sign-off, as CONTRIBUTING asks of any divergence. If parity wins, the alternative is to sort by absolute position as STAR does.

Not changed

TranscriptomeIndex::from_gtf_exons_configured still groups by ID alone and drops a transcript with inconsistent chromosomes (src/quant/transcriptome.rs:289), so such a transcript has its junctions but is not in the transcript tables. I could not compare that with STAR here: Homebrew's macOS binary exits with "could not open input file /geneInfo.tab" on --quantMode runs, as #345 describes.

Measured

Same fixtures and arguments, before and after this PR:

GTF	main	This branch
Only X exons	exit 0, no junctions	exit 1, STAR's message
/dev/null	exit 0, no junctions	exit 1, STAR's message
GFF3, Parent=T	exit 1	exit 0, c1 1201 1999, as STAR
One exon without transcript_id	exit 1	exit 0, c1 1201 1999, as STAR
T on two chromosomes	exit 0, false c1 1201-1499, 1701-1999, 2201-2499	exit 0, c1 1201 1999, c2 1701 2499

A GTF with no ID shared across chromosomes and all exons parented produces the same junctions: the grouping only splits an ID that spans chromosomes. I did not run test/parity_suite.sh.

A GTF that used to run with no matching exon now stops at startup. That is the point, but it is a behaviour change for anyone relying on it.

Tests

In src/junction/gtf.rs: same_transcript_id_on_two_chromosomes_is_two_chains, no_exon_on_a_known_chromosome_is_fatal, exon_without_parent_is_its_own_transcript_and_gff3_attributes_parse, only_parentless_exons_on_the_genome_is_not_fatal, attribute_value_is_its_first_word_as_in_star. test_extract_junctions_unknown_chromosome gains one exon on a known chromosome, since a GTF with none is now refused.

cargo test --release, cargo clippy --all-targets and cargo fmt --check are clean on macOS.

Second commit, after review

Review against STAR 2.7.11b found three cases where the first commit was not what STAR does: a file whose only exons on the genome had no transcript ID was refused (STAR exits 0), Parent = rna-X split one transcript in two, and IDs containing = were kept whole. b1c16e8 reads them as STAR does. On a test GTF (IDs "T=1", T=2, exons with no transcript_id, an exon on an unknown chromosome) and a test GFF3 (Parent = rna-1, an exon with only ID=, no gene IDs), against a Linux build of STAR 2.7.11b:

File	sjdbList.out.tab, transcriptInfo.tab, exonInfo.tab, geneInfo.tab, exonGeTrInfo.tab	ReadsPerGene.out.tab, 60 reads
GTF	byte-identical	identical
GFF3	byte-identical	column 2 identical; columns 3 and 4 differ

The GFF3 difference is older than this PR: MissingGeneID holds exons on both strands, and rustar-aligner gives a gene one strand for stranded counting where STAR uses each exon's. Any gene ID spanning both strands is affected the same way.

🤖 Generated with Claude Code