feat(quant): unspliced targets in TranscriptomeSAM, splicing-status tag and GeneSplicing counts for bulk total RNA-seq

CHANGELOG.md

@@ -44,6 +44,31 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Features
 
+- **Bulk total RNA-seq: unspliced targets, splicing-status tag and
+  spliced / unspliced gene counts** (rustar-aligner extensions, opt-in, see
+  DIVERGENCE.md 1.4 and the "Bulk total RNA-seq" guide).
+  - `--quantTranscriptomeUnspliced Intron|PreMRNA` adds one `<gene_id>-I`
+    target per gene (merged introns plus `--quantTranscriptomeUnsplicedFlank`,
+    or the gene body) after the transcripts of
+    `Aligned.toTranscriptome.out.bam`; unspliced reads are projected onto
+    every compatible target, spliced and unspliced, so a read in a retained
+    intron is no longer forced onto the retained-intron isoform. Writes
+    `Aligned.toTranscriptome.targets.tsv` and, on request, the unspliced
+    sequences (`--quantTranscriptomeUnsplicedFasta Yes`).
+  - `--outSAMsplicingStatus Yes` tags genomic alignments `sp:A:S|U|A`
+    (spliced / unspliced / ambiguous).
+  - `--quantGeneSplicing Yes` writes `ReadsPerGeneSplicing.out.tab` and
+    `ReadsPerGeneSplicing.summary.tsv`, using STARsolo's spliced / unspliced
+    classification on each read or pair.
+  - These are separate rustar-aligner flags, not new values of STAR's
+    `--quantMode` or `--outSAMattributes`, so a STAR command line keeps its
+    STAR meaning.
+  - Behaviour change: `genomeGenerate` with `--quantGeneSplicing Yes` now
+    requires `--sjdbGTFfile`, as it already did for
+    `--quantMode TranscriptomeSAM`. `--quantGeneSplicing Yes`,
+    `--outSAMsplicingStatus Yes` and `--quantTranscriptomeUnspliced` all need
+    a GTF-aware index at `alignReads`.
+
 - **CLI and output parity: SAM/SJ/read-input knobs and the STAR limit
   surface** — 30 further STAR 2.7.11b parameters. (`--outSAMorder` came from #145.)
 

DIVERGENCE.md

@@ -66,6 +66,23 @@ The site that acts on those zeroed counts, however, gates on a different flag (`
 
 ---
 
+### 1.4 Bulk total-RNA options: unspliced transcriptome targets, `sp` tag, `--quantGeneSplicing`
+
+All three are rustar-aligner extensions with no STAR equivalent, and all are opt-in. Each has its own flag rather than a new value of a STAR parameter (`--quantMode`, `--outSAMattributes`), so a STAR command line never means something different here. With the defaults (`--quantTranscriptomeUnspliced None`, `--outSAMsplicingStatus No`, `--quantGeneSplicing No`) no output changes, which `tests/bulk_unspliced.rs` (`new_options_leave_existing_outputs_unchanged`) checks.
+
+**What STAR does.** `--quantMode TranscriptomeSAM` projects every alignment onto every annotated transcript whose exons contain it (`Transcriptome_quantAlign.cpp`); there are no unspliced targets, so an unspliced read inside an intron that another isoform retains is projected onto the retaining isoform only, and a read in a constitutive intron is not projected at all. STAR classifies reads as spliced / unspliced / ambiguous only for single-cell runs (STARsolo, `Transcriptome_classifyAlign.cpp` + `SoloFeature_countVelocyto.cpp`), and has no per-alignment splicing tag.
+
+**What rustar-aligner does.**
+
+- `--quantTranscriptomeUnspliced Intron|PreMRNA` appends one `<gene_id>-I` target per gene (merged introns of all isoforms plus `--quantTranscriptomeUnsplicedFlank`, or the gene body) after the annotated transcripts, and projects alignments onto them with the unchanged STAR projection, except that fragments crossing a splice junction are kept off the unspliced targets. It also writes `Aligned.toTranscriptome.targets.tsv` and, on request, the unspliced sequences.
+- `--outSAMsplicingStatus Yes` and `--quantGeneSplicing Yes` run a port of STARsolo's spliced / unspliced classification (`alignToTranscriptMinOverlap` with `minOverlapMinusOne = 6` and the 1 Mb intron cap, then the per-UMI collapse of `countVelocyto`), one read or pair standing for one UMI. Two details differ from `classifyAlign`: the containment test uses the true leftmost / rightmost aligned base of the pair, where STAR uses the first block's start and the last block's end (a `TODO` next to that line in STAR flags the case where mate 2 ends before mate 1); and blocks are sorted by position before the scan, so STAR's early exit at the last exon also holds for overlapping mates. The `sp` tag collapses over all containing transcripts regardless of gene and strand.
+
+**Why.** In ribo-depleted total RNA a large share of reads is pre-mRNA. Without unspliced targets STAR's projection gives the intronic ones to retained-intron isoforms, so isoform proportions follow the library's pre-mRNA content (see the PR for measurements, and COMBINE-lab/salmon#1229 for the same effect with Salmon decoys). Unspliced targets let the downstream EM share such reads, as splici does for single-cell data. Bulk users have no STAR option that reports how much of a library is unspliced.
+
+**Impact.** Only when the options are given. Unspliced targets add `@SQ` lines and records to `Aligned.toTranscriptome.out.bam` and change `NH` / `HI` / `MAPQ` of the reads they receive; `sp` adds one tag per record.
+
+**Source.** `src/quant/transcriptome.rs` (`with_unspliced_targets`, `unspliced_intervals`, `filter_and_project`), `src/quant/splice_status.rs`, `src/quant/mod.rs` (`QuantContext`), `src/params/sam.rs` (`SP`). STAR: `Transcriptome_quantAlign.cpp`, `Transcriptome_classifyAlign.cpp`, `SoloFeature_countVelocyto.cpp`.
+
 ## 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. They are listed here so the differential benchmark's non-exact reads are fully accounted for.

docs/astro.config.mjs

@@ -55,6 +55,7 @@ export default defineConfig({
             { label: 'Two-pass mode', slug: 'guides/two-pass' },
             { label: 'Chimeric detection', slug: 'guides/chimeric' },
             { label: 'Gene quantification', slug: 'guides/quantification' },
+            { label: 'Bulk total RNA-seq', slug: 'guides/bulk-total-rna' },
             { label: 'Migrating from STAR', slug: 'guides/migrating-from-star' },
           ],
         },

docs/src/content/docs/guides/bulk-total-rna.md

@@ -0,0 +1,165 @@
+---
+title: Bulk total RNA-seq
+description: Unspliced targets in the transcriptome BAM, a splicing-status tag and spliced / unspliced gene counts for ribo-depleted bulk libraries.
+---
+
+Ribo-depleted ("total") RNA-seq libraries contain a large share of unspliced
+pre-mRNA: intronic reads routinely make up a third to more than half of the
+fragments. rustar-aligner has three **opt-in** options for this kind of data.
+None of them exists in STAR; without them every output is the same as STAR's.
+They are separate flags, not new values of STAR's parameters, and all three
+need a GTF-aware index (`genomeGenerate` with `--quantGeneSplicing Yes`
+requires `--sjdbGTFfile`).
+
+| Option | What it adds |
+|---|---|
+| `--quantTranscriptomeUnspliced Intron` or `PreMRNA` | one unspliced target per gene, `<gene_id>-I`, in `Aligned.toTranscriptome.out.bam` |
+| `--outSAMsplicingStatus Yes` | an `sp:A` splicing-status tag on every genomic alignment |
+| `--quantGeneSplicing Yes` | spliced / unspliced / ambiguous read counts per gene |
+
+## Why pre-mRNA matters for transcript quantification
+
+`--quantMode TranscriptomeSAM` projects each genomic alignment onto every
+annotated transcript whose exons contain it. A read that lies inside an
+intron is compatible with no fully spliced isoform, but GENCODE annotates many
+`retained_intron` isoforms whose single exon covers that intron. The read is
+therefore projected onto the retained-intron isoform only, and Salmon or RSEM
+assign it there. In total RNA the number of such reads follows the pre-mRNA
+content of the library, not the abundance of the retained-intron isoform, so
+isoform proportions (and tximport's `avgTxLength` offsets) track library
+quality rather than biology. The same effect is reported for Salmon with
+genome decoys ([COMBINE-lab/salmon#1229](https://github.com/COMBINE-lab/salmon/issues/1229)).
+
+The fix used for single-cell data by
+[splici / alevin-fry](https://doi.org/10.1038/s41592-022-01408-3) is to give
+the quantifier unspliced targets next to the spliced transcripts, so that a
+read compatible with both is shared by the EM instead of being forced onto the
+retained-intron isoform. `--quantTranscriptomeUnspliced` does this for the
+transcriptome BAM.
+
+## Unspliced targets in the transcriptome BAM
+
+```bash
+rustar-aligner \
+  --genomeDir /path/to/genome_index \
+  --readFilesIn reads_1.fq.gz reads_2.fq.gz \
+  --quantMode TranscriptomeSAM \
+  --quantTranscriptomeUnspliced PreMRNA \
+  --quantTranscriptomeUnsplicedFasta Yes \
+  --outFileNamePrefix sample_
+```
+
+### Targets
+
+One target per gene, named `<gene_id>-I`, is added after the annotated
+transcripts (`@SQ` lines in gene order):
+
+- `Intron` (splici-style): the union of the introns of all the gene's
+  isoforms, merged, extended on each side by
+  `--quantTranscriptomeUnsplicedFlank` bases (default `-1` = the index's
+  `sjdbOverhang`, i.e. read length - 1 by STAR's convention), clipped to the
+  gene body, merged again and concatenated in genome order. Introns retained
+  by an isoform and introns skipped by a cassette exon are included.
+- `PreMRNA`: the whole gene body, from the first to the last annotated base.
+
+A gene is taken on the chromosome and strand of its first transcript. In
+`Intron` mode a single-exon gene has no target. Targets are built at
+alignment time from the transcript tables of the index, so any GTF-aware
+index works and the flank can change between runs.
+
+### Projection
+
+Nothing changes in how an alignment is projected: it goes to every target that
+contains all its aligned blocks, spliced and unspliced alike. So:
+
+- a read inside a retained intron is written both to the retained-intron
+  isoform and to `<gene_id>-I`; `NH`, `HI` and `MAPQ` count all targets, and
+  Salmon or RSEM decide;
+- a read inside a constitutive intron is written to `<gene_id>-I` only
+  (it used to be dropped);
+- a read or pair that **crosses a splice junction** is processed RNA and goes
+  to spliced targets only;
+- `--quantTranscriptomeSAMoutput` rules (indels, soft-clip extension,
+  single-end) apply to every target.
+
+`Intron` or `PreMRNA`? `PreMRNA` also offers a home to the exonic part of
+pre-mRNA and keeps pairs with one mate in an exon and one in an intron (in
+`Intron` mode such a pair fits no target unless an isoform retains that
+intron). `Intron` matches splici and keeps exon-only reads away from the
+unspliced targets, at the cost of counting the exonic part of pre-mRNA as
+mature.
+
+### Files
+
+- `sample_Aligned.toTranscriptome.targets.tsv`: `target_id`, `gene_id`,
+  `gene_name`, `status` (`spliced` / `unspliced`), `length`, in `@SQ` order;
+  a ready-made tx2gene table (sum by `gene_id` and `status` for spliced and
+  unspliced gene counts).
+- `sample_Aligned.toTranscriptome.unspliced.fa` (with
+  `--quantTranscriptomeUnsplicedFasta Yes`): the `<gene_id>-I` sequences, in
+  transcript orientation. It is of the order of the genome size and depends
+  only on the index and the flank, so write it once.
+
+### Salmon
+
+Salmon's alignment mode needs every `@SQ` target in its FASTA, with the same
+names. With GENCODE:
+
+```bash
+gzip -dc gencode.v50.transcripts.fa.gz | sed 's/|.*//' > transcripts.fa
+cat transcripts.fa sample_Aligned.toTranscriptome.unspliced.fa > targets.fa
+salmon quant -a sample_Aligned.toTranscriptome.out.bam -t targets.fa -l A -o salmon_out
+```
+
+Keep the `-I` targets out of isoform-level analyses, and out of tximport's
+`avgTxLength` if a spliced-only gene length is wanted.
+
+## Splicing-status tag (`sp`)
+
+`--outSAMsplicingStatus Yes` adds `sp:A:S` (spliced: compatible only with
+mature mRNA, the read need not cross a junction), `sp:A:U` (unspliced: needs
+pre-mRNA) or `sp:A:A` (ambiguous: both) to every genomic alignment record,
+from the rules below applied to every annotated transcript that contains the
+alignment, on either strand and whatever its gene. There is no tag when no
+transcript contains the alignment. Both mates of a pair carry the fragment
+status, after the `--outSAMattributes` tags. The `sp` tag name is not used by
+STAR or STARsolo.
+
+## Spliced / unspliced gene counts (`--quantGeneSplicing`)
+
+```bash
+rustar-aligner ... --quantMode GeneCounts --quantGeneSplicing Yes
+```
+
+A cheap per-library measure of pre-mRNA content, and a gene-level table. Each
+uniquely mapped read or pair is classified with the rules STARsolo uses for its
+single-cell spliced / unspliced matrices, one read standing for one molecule:
+
+1. Every annotated transcript that fully contains the alignment is tested.
+   Each aligned block is called exonic, intronic or exon/intron-spanning, with
+   a 6-base tolerance at exon boundaries. A spliced alignment that touches an
+   intron is incompatible with that transcript, and a block in an intron longer
+   than 1 Mb is not called intronic.
+2. If the compatible transcripts belong to more than one gene, the read is
+   `N_multiGene` and not counted.
+3. Otherwise: only-exonic models give **spliced**, intronic or spanning models
+   with no only-exonic model give **unspliced**, and a mix gives
+   **ambiguous** (for example a read inside an intron that another isoform
+   retains).
+
+Unmapped, too-many-loci and multimapping reads are accounted as in
+`GeneCounts`; pairs with a single mapped mate count as unmapped.
+
+- `sample_ReadsPerGeneSplicing.out.tab`: a header line, then one line per gene
+  (in `geneInfo.tab` order) with spliced, unspliced and ambiguous counts for
+  each strand convention (`unstranded_*`, `forward_*`, `reverse_*`).
+  `forward` keeps transcripts on the strand of read 1 (`htseq-count -s yes`);
+  `reverse` keeps the opposite strand (dUTP / TruSeq Stranded,
+  `-s reverse`).
+- `sample_ReadsPerGeneSplicing.summary.tsv`: `N_unmapped`, `N_multimapping`,
+  `N_noFeature`, `N_multiGene`, `N_spliced`, `N_unspliced`, `N_ambiguous` and
+  the three fractions of the assigned reads, per strand convention.
+
+The benchmark behind these options (public whole-blood total RNA, GRCh38 +
+GENCODE v50) can be reproduced with the scripts in
+`scripts/bench_bulk_unspliced/`.

docs/src/content/docs/guides/quantification.md

@@ -79,6 +79,10 @@ You can request both modes in the same run:
 
 This emits `ReadsPerGene.out.tab` *and* `Aligned.toTranscriptome.out.bam` in addition to the normal alignment output.
 
+## Total RNA (ribo-depleted) libraries
+
+For libraries with a large pre-mRNA content, `--quantTranscriptomeUnspliced` adds one unspliced `<gene_id>-I` target per gene to the transcriptome BAM (so pre-mRNA reads stop landing on retained-intron isoforms), `--outSAMsplicingStatus Yes` tags genomic alignments with their splicing status, and `--quantGeneSplicing Yes` counts spliced / unspliced / ambiguous reads per gene. All three are rustar-aligner extensions; see [Bulk total RNA-seq](/rustar-aligner/guides/bulk-total-rna/).
+
 ## Index-time vs alignment-time
 
 For best speed, supply `--sjdbGTFfile` at `--runMode genomeGenerate` time and the transcript-level data structures get persisted into the genome directory. Then at alignment time you only need `--quantMode TranscriptomeSAM` (or `GeneCounts`); rustar-aligner reuses the persisted annotations.

docs/src/content/docs/reference/cli-parameters.md

@@ -153,6 +153,11 @@ Run `rustar-aligner --help` for the full machine-generated listing.
 |-----------|---------|-------------|
 | `--quantMode` | — | `GeneCounts` and/or `TranscriptomeSAM`, space-separated. |
 | `--quantTranscriptomeSAMoutput` | `BanSingleEnd_BanIndels_ExtendSoftclip` | Variant for transcriptome BAM: `BanSingleEnd`, `BanSingleEnd_ExtendSoftclip`, or the default RSEM-compatible form. |
+| `--quantTranscriptomeUnspliced` | `None` | rustar-aligner extension. `Intron` / `PreMRNA`: add one unspliced target per gene (`<gene_id>-I`: merged introns plus flanks, or the gene body) to the transcriptome BAM. `None` is STAR's behaviour. See [Bulk total RNA-seq](/rustar-aligner/guides/bulk-total-rna/). |
+| `--quantTranscriptomeUnsplicedFlank` | `-1` | rustar-aligner extension. Flank around each merged intron for `Intron`; `-1` uses the index's `sjdbOverhang`. |
+| `--quantTranscriptomeUnsplicedFasta` | `No` | rustar-aligner extension. `Yes` writes the unspliced target sequences to `Aligned.toTranscriptome.unspliced.fa`. |
+| `--quantGeneSplicing` | `No` | rustar-aligner extension. `Yes` writes bulk spliced / unspliced / ambiguous gene counts (`ReadsPerGeneSplicing.out.tab`), see [Bulk total RNA-seq](/rustar-aligner/guides/bulk-total-rna/). |
+| `--outSAMsplicingStatus` | `No` | rustar-aligner extension. `Yes` adds an `sp:A` splicing-status tag (S / U / A) to genomic SAM/BAM records. |
 
 ## Two-pass mode
 

docs/src/content/docs/reference/output-files.md

@@ -23,7 +23,7 @@ Coordinate-sorted BAM. Written when `--outSAMtype BAM SortedByCoordinate`. The s
 
 ### `sample_Aligned.toTranscriptome.out.bam`
 
-Transcriptome-coordinate BAM. Written when `--quantMode TranscriptomeSAM` is set. Each record's reference is a transcript ID rather than a chromosome; one record is emitted per transcript that the read aligns within.
+Transcriptome-coordinate BAM. Written when `--quantMode TranscriptomeSAM` is set. Each record's reference is a transcript ID rather than a chromosome; one record is emitted per transcript that the read aligns within. With `--quantTranscriptomeUnspliced Intron|PreMRNA` (rustar-aligner extension) the references also include one `<gene_id>-I` unspliced target per gene, and `sample_Aligned.toTranscriptome.targets.tsv` (target, gene, `spliced`/`unspliced` status, length) is written next to it, plus `sample_Aligned.toTranscriptome.unspliced.fa` with `--quantTranscriptomeUnsplicedFasta Yes`. See [Bulk total RNA-seq](/rustar-aligner/guides/bulk-total-rna/).
 
 ## Log files
 
@@ -85,6 +85,10 @@ gene_id    unstranded    forward_stranded    reverse_stranded
 
 The first four rows are summary categories: `N_unmapped`, `N_multimapping`, `N_noFeature`, `N_ambiguous`. Subsequent rows are per-gene counts. Pick the column matching your library's strandedness — see the [quantification guide](/rustar-aligner/guides/quantification/).
 
+### `sample_ReadsPerGeneSplicing.out.tab` / `sample_ReadsPerGeneSplicing.summary.tsv`
+
+Written when `--quantGeneSplicing Yes` is set (rustar-aligner extension, not in STAR). The table has a header line and one line per gene with spliced, unspliced and ambiguous counts for the unstranded, forward and reverse strand conventions (nine count columns). The summary gives read accounting and the spliced / unspliced / ambiguous fractions per strand convention. See [Bulk total RNA-seq](/rustar-aligner/guides/bulk-total-rna/).
+
 ## Unmapped reads
 
 ### `sample_Unmapped.out.mate1` / `sample_Unmapped.out.mate2`