[BenjaminDEMAILLE]

Part of #354 (GeneCounts stranded columns: strand per exon, not per gene).

Problem

--quantMode GeneCounts gave each gene one strand, taken from its first exon (gene_is_reverse in src/quant/mod.rs). Columns 3 and 4 of ReadsPerGene.out.tab filtered genes on that strand. So when one gene ID has exons on both strands, every read went to the column of the first exon's strand. This happens with a GTF that reuses a gene_id on a + and a - transcript, and with GFF3 exons that all fall into STAR's MissingGeneID (#356).

A second, smaller case: a transcript with strand . was stored as + (strand 1). STAR stores it as strand 0, and strand-0 exons count in both stranded columns.

What STAR does (2.7.11b source)

• GTF.cpp:138-144: a transcript's strand is 1 for +, 2 for -, 0 for anything else, from its first exon line.
• GTF_transcriptGeneSJ.cpp:40: each exon in exonGeTrInfo.tab carries its transcript's strand.
• Transcriptome.cpp:79-91: exG.str is loaded per exon from exonGeTrInfo.tab. There is no per-gene strand for counting.
• Transcriptome_geneCountsAddAlign.cpp:30-45: for each exon overlapping a block of the read, str1 = exG.str - 1. Column "same strand" skips the exon when a.Str != str1 && str1 < 2; column "reverse" skips it when a.Str == str1 && str1 < 2. So the strand test is per exon, and a strand-0 exon (str1 wraps to a large value) is never skipped.
• A read overlapping a + exon and a - exon of the same gene therefore counts for that gene in both stranded columns. It is not ambiguous, because ambiguity is only about different gene IDs (line 42).

Change

• GeneAnnotation gains chr_exons_strand, parallel to chr_exons: a strand mask per exon (+, -, or both for strand 0). Exons with the same start, end and gene are merged and their masks OR-ed, which counts the same as STAR's separate entries.
• From a GTF, the mask comes from the exon's strand. From the index tables, it comes from tr_strand (what STAR writes into exonGeTrInfo.tab).
• New overlapping_gene_strands returns each overlapping gene with the OR of the masks of the exons the read overlaps (both mates for PE). GeneCounts::apply_three_columns keeps a gene in a stranded column when that mask has the bit for the column's strand. The overlap walk is shared with overlapping_genes (for_each_overlap), so the segment tree is unchanged.
• TranscriptomeIndex stores strand 0 for a . transcript (was 1). This makes transcriptInfo.tab and exonGeTrInfo.tab byte-identical to STAR for such transcripts. Other code only tests tr_strand == 2 or == 1; the == 1 test in solo (src/solo/gene.rs) is STAR's trStr==1 ? aG.Str : 1-aG.Str (Transcriptome_classifyAlign.cpp:208), so a . transcript now behaves there as in STAR too.
• gene_is_reverse stays: STARsolo still uses it.
• chr_exons is now sorted by (start, end, gene) rather than (start, end). The old dedup could miss duplicates that were not adjacent.

Measured against STAR 2.7.11b

Synthetic fixture: 20 kb random chromosome, 105 single-end 50 bp reads, Linux STAR in docker, each tool with its own index. G1 has transcript T1 on + and T2 on - (30 forward reads on T1, 30 reverse reads on T2). G2 is a + gene with 20 reverse reads. G3 is a . transcript with 10 forward and 15 reverse reads.

ReadsPerGene.out.tab	STAR	rustar main	rustar this PR
N_noFeature	0 20 60	0 65 40	0 20 60
G1 (+ and - transcripts)	60 60 0	60 30 30	60 60 0
G2 (+)	20 0 20	20 0 20	20 0 20
G3 (.)	25 25 25	25 10 15	25 25 25

Both rustar gene-model sources give the same table: the index tables, and --sjdbGTFfile at mapping. transcriptInfo.tab and exonGeTrInfo.tab of this fixture are now byte-identical to STAR's (main wrote strand 1 for T4).

Yeast parity config (ERR12389696, 100k PE, --quantMode TranscriptomeSAM GeneCounts --outSAMtype BAM SortedByCoordinate --runThreadN 8): ReadsPerGene.out.tab is byte-identical to STAR's, on main and with this PR. That annotation has no gene on both strands, so this is a no-regression check only.

Tests

• test_gene_counts_strand_per_exon: a gene with a + and a - exon, and a . gene, each hit by one forward read. Fails on main (strand1 is [2, 1] instead of [1, 1]), passes here.
• from_transcriptome_matches_from_gtf now has a . transcript and also compares chr_exons_strand. It fails with main's transcriptome.rs (strand 1 for .), passes here.
• cargo test --release, cargo clippy --all-targets -- -D warnings and cargo fmt --check are clean on macOS.

Not checked

• From a GTF given at mapping time, the mask uses each exon line's own strand. STAR uses the strand of the transcript's first exon line. They differ only for a transcript whose exons disagree on strand; TranscriptomeIndex already skips such transcripts with a warning, but the GeneCounts model does not. Marked shortcut: in the code.
• The GTF parser: refuse an annotation that matches nothing, accept GFF3 #356 GFF3 MissingGeneID case was not rerun here (that branch is not merged). It goes through the same code path.
• STARsolo's own strand handling (gene_is_reverse) was not changed or compared.
• sjdbList.fromGTF.out.tab now writes . for junctions of a . transcript, as STAR's strandChar (GTF_transcriptGeneSJ.cpp:142) does; my fixture has no spliced . transcript, so this was not compared.

🤖 Generated with Claude Code