From a941067393e25f6c9eedd327c7b236f9d6a61ca5 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 10:39:13 +0100 Subject: [PATCH 1/2] GeneCounts: stranded columns use each exon's strand, as STAR ReadsPerGene.out.tab columns 3 and 4 took one strand per gene (its first exon's). STAR keeps a strand per exon (exG.str, the exon's transcript strand) and a read counts for a gene in a stranded column when an exon it overlaps is on that strand; a strand-0 ('.') exon counts in both. The transcriptome tables now write strand 0 for '.' transcripts, as STAR's GTF.cpp does. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 9 ++ src/quant/mod.rs | 194 +++++++++++++++++++++++++++---------- src/quant/transcriptome.rs | 3 +- 3 files changed, 155 insertions(+), 51 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e..217ada5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -232,6 +232,15 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- **`ReadsPerGene.out.tab` stranded columns use each exon's strand.** + Columns 3 and 4 took one strand per gene, from its first exon, so a gene + ID with transcripts on both strands put every read in one column. Like + STAR (`exG.str`, `Transcriptome_geneCountsAddAlign.cpp`), a read now + counts for a gene in a stranded column when an exon it overlaps is on that + strand. A transcript with strand `.` is strand 0 in `transcriptInfo.tab` + and `exonGeTrInfo.tab` (it was 1), and its exons count in both stranded + columns, as in STAR. + - **STARsolo's stats files as STAR writes them.** `Solo.out/Barcodes.stats`, `/Features.stats` and `/UMIperCellSorted.txt` were not written, SJ and Velocyto had no `Summary.csv`, and the Gene/GeneFull diff --git a/src/quant/mod.rs b/src/quant/mod.rs index f3c39f6..31e96b4 100644 --- a/src/quant/mod.rs +++ b/src/quant/mod.rs @@ -42,9 +42,13 @@ pub struct GeneAnnotation { pub gene_names: Vec, /// Strand per gene: true = reverse/minus strand. pub gene_is_reverse: Vec, - /// Per-chromosome exon interval list, sorted by (start, end). + /// Per-chromosome exon interval list, sorted by (start, end, gene). /// Each entry: (start_0based_incl, end_0based_excl, gene_idx). pub chr_exons: Vec>, + /// Strands of each `chr_exons` entry, parallel to it: [`STRAND_PLUS`] and/or + /// [`STRAND_MINUS`], both for a strand-less (`.`) transcript. GeneCounts keys + /// strand per exon like STAR's `exG.str`, not per gene. + pub chr_exons_strand: Vec>, /// Per-chromosome **gene-body** interval list (one entry per gene: its full /// `[min exon start, max exon end)` span, covering introns), sorted by /// (start, end). Used by the STARsolo `GeneFull` feature, which counts a @@ -65,6 +69,20 @@ pub struct GeneAnnotation { chr_gene_body_seg: Vec>, } +/// Exon strand mask bits of [`GeneAnnotation::chr_exons_strand`]. +pub const STRAND_PLUS: u8 = 1; +pub const STRAND_MINUS: u8 = 2; + +/// Mask for STAR's transcript strand code (1 `+`, 2 `-`, 0 neither): a strand-0 +/// exon accepts reads of both strands (`Transcriptome_geneCountsAddAlign.cpp:34`). +fn strand_mask(star_strand: u8) -> u8 { + match star_strand { + 1 => STRAND_PLUS, + 2 => STRAND_MINUS, + _ => STRAND_PLUS | STRAND_MINUS, + } +} + impl GeneAnnotation { /// Build from GTF exon records using `gene_tag` as the gene grouping attribute. /// @@ -76,7 +94,7 @@ impl GeneAnnotation { let mut gene_is_reverse: Vec = Vec::new(); let mut gene_id_to_idx: std::collections::HashMap = std::collections::HashMap::new(); - let mut chr_exons: Vec> = vec![Vec::new(); n_chrs]; + let mut chr_exons: Vec> = vec![Vec::new(); n_chrs]; // Per-gene full span: (chr_idx, min_start, max_end). Accumulated over all // of a gene's exons to build the GeneFull gene-body intervals. let mut gene_span: Vec> = Vec::new(); @@ -120,7 +138,14 @@ impl GeneAnnotation { let start = chr_offset + exon.start.saturating_sub(1); let end = chr_offset + exon.end; - chr_exons[chr_idx].push((start, end, gene_idx)); + // shortcut: the exon's own strand, where STAR uses its transcript's + // (first exon's) strand; they differ only in a malformed GTF. + let star_strand = match exon.strand { + '+' => 1, + '-' => 2, + _ => 0, + }; + chr_exons[chr_idx].push((start, end, gene_idx, strand_mask(star_strand))); // Extend this gene's full span. (A gene's exons share one chr.) match &mut gene_span[gene_idx] { @@ -146,7 +171,7 @@ impl GeneAnnotation { let n_genes = tr.gene_ids.len(); let mut gene_is_reverse = vec![false; n_genes]; let mut gene_seen = vec![false; n_genes]; - let mut chr_exons: Vec> = vec![Vec::new(); n_chrs]; + let mut chr_exons: Vec> = vec![Vec::new(); n_chrs]; let mut gene_span: Vec> = vec![None; n_genes]; for t in 0..tr.n_transcripts() { @@ -163,7 +188,8 @@ impl GeneAnnotation { continue; } for ex in &tr.tr_exons[t] { - chr_exons[chr_idx].push((ex.genome_start, ex.genome_end, g)); + let mask = strand_mask(tr.tr_strand[t]); + chr_exons[chr_idx].push((ex.genome_start, ex.genome_end, g, mask)); match &mut gene_span[g] { Some((_, s, e)) => { *s = (*s).min(ex.genome_start); @@ -189,13 +215,28 @@ impl GeneAnnotation { gene_ids: Vec, gene_names: Vec, gene_is_reverse: Vec, - mut chr_exons: Vec>, + stranded_exons: Vec>, gene_span: &[Option<(usize, u64, u64)>], ) -> Self { - let n_chrs = chr_exons.len(); - for exons in &mut chr_exons { - exons.sort_unstable_by_key(|&(s, e, _)| (s, e)); - exons.dedup(); + let n_chrs = stranded_exons.len(); + let mut chr_exons: Vec> = Vec::with_capacity(n_chrs); + let mut chr_exons_strand: Vec> = Vec::with_capacity(n_chrs); + for mut exons in stranded_exons { + exons.sort_unstable(); + let mut iv: Vec<(u64, u64, usize)> = Vec::with_capacity(exons.len()); + let mut strands: Vec = Vec::with_capacity(exons.len()); + // One entry per (start, end, gene); its strands OR-ed, which counts + // like STAR's separate exG entries per strand. + for (s, e, g, m) in exons { + if iv.last() == Some(&(s, e, g)) { + *strands.last_mut().unwrap() |= m; + } else { + iv.push((s, e, g)); + strands.push(m); + } + } + chr_exons.push(iv); + chr_exons_strand.push(strands); } // Build the per-chromosome gene-body interval list. @@ -248,6 +289,7 @@ impl GeneAnnotation { gene_names, gene_is_reverse, chr_exons, + chr_exons_strand, chr_gene_body, gene_exons, chr_exons_seg, @@ -320,6 +362,31 @@ impl GeneAnnotation { Self::overlapping_in_into(&self.chr_exons, &self.chr_exons_seg, transcript, out); } + /// Genes whose exons overlap any of `transcripts` (the mates of one read), + /// each with the OR of the strand masks of its overlapping exons. Sorted by + /// gene, one entry per gene. Drives `ReadsPerGene.out.tab`. + pub fn overlapping_gene_strands(&self, transcripts: &[&Transcript]) -> Vec<(usize, u8)> { + let mut hits: Vec<(usize, u8)> = Vec::new(); + for t in transcripts { + let Some(strands) = self.chr_exons_strand.get(t.chr_idx) else { + continue; + }; + let chr = &self.chr_exons[t.chr_idx]; + Self::for_each_overlap(chr, &self.chr_exons_seg[t.chr_idx], t, |i| { + hits.push((chr[i].2, strands[i])); + }); + } + hits.sort_unstable(); + let mut out: Vec<(usize, u8)> = Vec::with_capacity(hits.len()); + for (g, m) in hits { + match out.last_mut() { + Some((lg, lm)) if *lg == g => *lm |= m, + _ => out.push((g, m)), + } + } + out + } + /// `overlapping_genes_full` into a caller-provided buffer. pub fn overlapping_genes_full_into(&self, transcript: &Transcript, out: &mut Vec) { Self::overlapping_in_into( @@ -353,10 +420,24 @@ impl GeneAnnotation { return; } let chr = &chr_intervals[transcript.chr_idx]; + Self::for_each_overlap(chr, &chr_segs[transcript.chr_idx], transcript, |i| { + out.push(chr[i].2); + }); + out.sort_unstable(); + out.dedup(); + } + + /// Call `f` with the index in `chr` of every interval overlapping an exon + /// of `transcript` (an interval overlapping several exons comes once per exon). + fn for_each_overlap( + chr: &[(u64, u64, usize)], + seg: &[u64], + transcript: &Transcript, + mut f: impl FnMut(usize), + ) { if chr.is_empty() { return; } - let seg = &chr_segs[transcript.chr_idx]; let sz = seg.len() / 2; // number of leaves (power of two) for exon in &transcript.exons { @@ -385,7 +466,7 @@ impl GeneAnnotation { } if nhi - nlo == 1 { // Leaf: nlo < upper (start < re) and its end > rs → overlap. - out.push(chr[nlo].2); + f(nlo); continue; } let mid = usize::midpoint(nlo, nhi); @@ -395,9 +476,6 @@ impl GeneAnnotation { sp += 1; } } - - out.sort_unstable(); - out.dedup(); } } @@ -498,15 +576,11 @@ impl GeneCounts { let pair: &PairedAlignment = both_mapped[0]; // Union of all genes overlapping either mate. - let mut all_genes = gene_ann.overlapping_genes(&pair.mate1_transcript); - let genes2 = gene_ann.overlapping_genes(&pair.mate2_transcript); - all_genes.extend_from_slice(&genes2); - all_genes.sort_unstable(); - all_genes.dedup(); + let genes = + gene_ann.overlapping_gene_strands(&[&pair.mate1_transcript, &pair.mate2_transcript]); // Strand is determined by mate1. - let read_is_rev = pair.mate1_transcript.is_reverse; - self.apply_three_columns(&all_genes, read_is_rev, gene_ann); + self.apply_three_columns(&genes, pair.mate1_transcript.is_reverse); } /// Write `ReadsPerGene.out.tab` in STAR's format. @@ -551,32 +625,28 @@ impl GeneCounts { /// Run 3 independent counting passes on a pre-computed set of overlapping genes. /// - /// `all_genes`: sorted+deduplicated gene indices overlapping the read (any strand). + /// `genes`: overlapping genes with the strand mask of their overlapping + /// exons (from [`GeneAnnotation::overlapping_gene_strands`]). /// `read_is_rev`: strand of the read (mate1 for PE). - fn apply_three_columns( - &self, - all_genes: &[usize], - read_is_rev: bool, - gene_ann: &GeneAnnotation, - ) { - // Col 1 — unstranded: use all_genes directly. - self.apply_column(all_genes, 0, &self.unstranded); - - // Col 2 — strand1: genes whose strand matches the read. - let same: Vec = all_genes - .iter() - .copied() - .filter(|&g| gene_ann.gene_is_reverse[g] == read_is_rev) - .collect(); - self.apply_column(&same, 1, &self.strand1); - - // Col 3 — strand2: genes whose strand is opposite to the read. - let opp: Vec = all_genes - .iter() - .copied() - .filter(|&g| gene_ann.gene_is_reverse[g] != read_is_rev) - .collect(); - self.apply_column(&opp, 2, &self.strand2); + /// A gene counts in a stranded column when one of its overlapping exons is + /// on that strand (STAR `Transcriptome_geneCountsAddAlign.cpp:30-45`). + fn apply_three_columns(&self, genes: &[(usize, u8)], read_is_rev: bool) { + let (same_bit, opp_bit) = if read_is_rev { + (STRAND_MINUS, STRAND_PLUS) + } else { + (STRAND_PLUS, STRAND_MINUS) + }; + let pick = |bit: u8| -> Vec { + genes + .iter() + .filter(|&&(_, m)| m & bit != 0) + .map(|&(g, _)| g) + .collect() + }; + let all: Vec = genes.iter().map(|&(g, _)| g).collect(); + self.apply_column(&all, 0, &self.unstranded); + self.apply_column(&pick(same_bit), 1, &self.strand1); + self.apply_column(&pick(opp_bit), 2, &self.strand2); } fn apply_column(&self, genes: &[usize], col: usize, per_gene: &[AtomicU64]) { @@ -600,8 +670,8 @@ impl GeneCounts { read_is_reverse: bool, gene_ann: &GeneAnnotation, ) { - let all_genes = gene_ann.overlapping_genes(transcript); - self.apply_three_columns(&all_genes, read_is_reverse, gene_ann); + let genes = gene_ann.overlapping_gene_strands(&[transcript]); + self.apply_three_columns(&genes, read_is_reverse); } } @@ -1144,6 +1214,29 @@ mod tests { assert_eq!(counts.strand2[0].load(Ordering::Relaxed), 1); } + /// STAR keys the stranded columns on each exon's strand, not on the gene's: + /// a gene with a + and a - transcript counts a forward read on either as + /// same-strand only on the + exon. A `.` exon counts in both columns. + #[test] + fn test_gene_counts_strand_per_exon() { + let genome = make_genome(); + let exons = vec![ + make_gtf_exon("chr1", 101, 200, '+', "G1"), + make_gtf_exon("chr1", 501, 600, '-', "G1"), + make_gtf_exon("chr1", 701, 800, '.', "G2"), + ]; + let ann = GeneAnnotation::from_gtf_exons(&exons, &genome); + let counts = GeneCounts::new(ann.n_genes()); + for (s, e) in [(110, 150), (510, 550), (710, 750)] { + counts.count_se_read(&[make_transcript(0, s, e, false)], 1, &ann); + } + let col = + |c: &[AtomicU64]| -> Vec { c.iter().map(|x| x.load(Ordering::Relaxed)).collect() }; + assert_eq!(col(&counts.unstranded), [2, 1]); + assert_eq!(col(&counts.strand1), [1, 1]); + assert_eq!(col(&counts.strand2), [1, 1]); + } + #[test] fn from_transcriptome_matches_from_gtf() { // The index-table gene model must count like the GTF-built one. @@ -1152,7 +1245,7 @@ mod tests { make_gtf_exon("chr1", 101, 200, '+', "G1"), make_gtf_exon("chr1", 301, 400, '+', "G1"), make_gtf_exon("chr1", 351, 450, '-', "G2"), - make_gtf_exon("chr2", 51, 150, '-', "G3"), + make_gtf_exon("chr2", 51, 150, '.', "G3"), ]; for (e, t) in exons.iter_mut().zip(["T1", "T1", "T2", "T3"]) { e.attributes @@ -1165,6 +1258,7 @@ mod tests { assert_eq!(from_tr.gene_ids, from_gtf.gene_ids); assert_eq!(from_tr.gene_is_reverse, from_gtf.gene_is_reverse); assert_eq!(from_tr.chr_exons, from_gtf.chr_exons); + assert_eq!(from_tr.chr_exons_strand, from_gtf.chr_exons_strand); assert_eq!(from_tr.chr_gene_body, from_gtf.chr_gene_body); assert_eq!(from_tr.gene_exons, from_gtf.gene_exons); } diff --git a/src/quant/transcriptome.rs b/src/quant/transcriptome.rs index 0a02059..6923fd9 100644 --- a/src/quant/transcriptome.rs +++ b/src/quant/transcriptome.rs @@ -343,8 +343,9 @@ impl TranscriptomeIndex { let total_len = ex_len_cum; let strand_u8 = match strand_char { + '+' => 1u8, '-' => 2u8, - _ => 1u8, // forward or unknown → treat as forward (STAR default) + _ => 0u8, // STAR GTF.cpp:138-144: neither strand }; let gene_id = first.attributes.get(gene_tag).cloned().unwrap_or_default(); From 823ddb9fa4408da198aa7bc527e5ffc1c460ff4e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 11:16:22 +0100 Subject: [PATCH 2/2] CHANGELOG: note the strand-0 effect on STARsolo Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 217ada5..ecdf5ec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -239,7 +239,8 @@ Sections commonly used: Features, Bug fixes, Other changes. counts for a gene in a stranded column when an exon it overlaps is on that strand. A transcript with strand `.` is strand 0 in `transcriptInfo.tab` and `exonGeTrInfo.tab` (it was 1), and its exons count in both stranded - columns, as in STAR. + columns, as in STAR. STARsolo's `tr_strand == 1` test now treats such a + transcript as reverse, also as STAR does. - **STARsolo's stats files as STAR writes them.** `Solo.out/Barcodes.stats`, `/Features.stats` and `/UMIperCellSorted.txt` were not