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

[BenjaminDEMAILLE]

Problem

Ribo-depleted ("total") RNA-seq contains a lot of unspliced pre-mRNA, and --quantMode TranscriptomeSAM has nowhere to put it. The projection (Transcriptome_quantAlign.cpp) sends a genomic alignment to every transcript whose exons contain it. An unspliced read inside an intron fits no fully spliced isoform, but GENCODE annotates many retained_intron isoforms whose single exon covers that intron: the read is projected onto that isoform only and Salmon / RSEM must assign it there. The amount follows the library's pre-mRNA content, not the isoform's abundance. Bulk users also have no STAR measure of spliced vs unspliced reads (GeneFull / Velocyto exist only in STARsolo).

Motivating measurements (made by the PR author with Salmon, not with rustar-aligner): 96 whole-blood libraries, ribo-depleted total RNA, reverse-stranded, 2x51 bp. With Salmon and full GRCh38 decoys, about 58% of fragments went to decoys (intronic pre-mRNA, expected in total RNA). A pre-mRNA read compatible with a retained_intron isoform is assigned to it; library by library, the share of reads on retained_intron isoforms tracks the decoy share (Spearman rho = 0.94, n = 96), unchanged at fixed residual gDNA and unrelated to TIN. Isoform proportions, hence tximport avgTxLength offsets, did not reproduce between technical replicates (r ~ 0.05).

The public benchmark below reproduces the pattern with rustar-aligner's STAR-faithful projection: over 16 public whole-blood libraries the retained-intron share tracks the unspliced fraction with Spearman rho = 0.92.

Refs COMBINE-lab/salmon#1229. See also https://support.bioconductor.org/p/9163851/, He et al. 2022 (splici, Nat Methods, doi:10.1038/s41592-022-01408-3) and Soneson et al. 2021 (PLoS Comput Biol, doi:10.1371/journal.pcbi.1008585). I know of no published bulk benchmark, so this PR brings its own.

What this PR adds (all opt-in, all rustar-aligner extensions)

1. Unspliced targets in the transcriptome BAM: --quantTranscriptomeUnspliced None|Intron|PreMRNA

Default None = STAR. With Intron or PreMRNA, one target per gene, <gene_id>-I (splici naming), is appended after the annotated transcripts, @SQ lines in gene order:

• Intron: union of the introns of all the gene's isoforms (retained introns and cassette-exon introns included), merged, extended by --quantTranscriptomeUnsplicedFlank on each side (default -1 = the index's sjdbOverhang, i.e. read length - 1 by STAR convention), clipped to the gene body, merged again, concatenated in genome order.
• PreMRNA: the whole gene body.

Projection is the unchanged STAR logic applied to all targets: an alignment goes to every target that contains all its blocks, spliced and unspliced alike. A read in a retained intron is written to both the retained-intron isoform and <gene_id>-I (NH / HI / MAPQ count all targets) and the EM decides; a read in a constitutive intron now reaches <gene_id>-I instead of being dropped; a fragment crossing a splice junction only goes to spliced targets. --quantTranscriptomeSAMoutput rules (indels, soft-clip extension, single-end) apply to every target; minus-strand targets are read 3' to 5' like transcripts.

Also written: Aligned.toTranscriptome.targets.tsv (target_id, gene_id, gene_name, spliced|unspliced status, length, in @SQ order: a tx2gene table) and, with --quantTranscriptomeUnsplicedFasta Yes, Aligned.toTranscriptome.unspliced.fa (the -I sequences, needed next to the transcript FASTA by salmon quant -a; off by default because it is genome-sized and only depends on index + flank).

Targets are built at alignment time from the transcript tables of any GTF-aware index: no index format change, flank adjustable per run, works with indexes built before this PR. Cost measured below: none visible in wall time, +0.3 to 0.4 GiB RSS.

2. sp:A splicing-status tag in the genomic SAM/BAM: --outSAMattributes ... sp

In no preset (default output unchanged). sp:A:S spliced (compatible only with mature mRNA), sp:A:U unspliced (needs pre-mRNA), sp:A:A ambiguous; no tag when no annotated transcript contains the alignment; both mates carry the fragment status. The name is not used by STAR / STARsolo (checked against parametersDefault); lower-case tags are reserved for local use by the SAM spec.

3. Per-gene table: --quantMode GeneSplicing (secondary, cheap)

ReadsPerGeneSplicing.out.tab (spliced / unspliced / ambiguous per gene for the unstranded / forward / reverse conventions of ReadsPerGene.out.tab) and ReadsPerGeneSplicing.summary.tsv (read accounting and fractions). It only runs when requested and reuses the classifier of the sp tag, so I kept it; its unspliced fraction is the per-library pre-mRNA measure used in the benchmark.

The classifier behind 2 and 3 is a port of STARsolo's Velocyto logic (provenance): alignToTranscriptMinOverlap from Transcriptome_classifyAlign.cpp (minOverlapMinusOne = 6, 1 Mb intron cap, spliced alignments touching an intron rejected) and the per-UMI collapse of SoloFeature_countVelocyto.cpp, one read or pair standing for one UMI.

Also in this PR

• fix(transcriptome): with the default BanSingleEnd_BanIndels_ExtendSoftclip, a CIGAR such as 1S97M2S lost its right clip when rebuilt (single-op body), producing a 98M record for a 100-base read; BAM encoding then failed and main aborts the whole run (read length-sequence length mismatch) on the public data below. STAR extends both ends (ReadAlign_quantTranscriptome.cpp); fixed, with a unit test.
• fix(params): unknown --quantMode values were silently ignored (a typo gave no counts and no error); now rejected as STAR does (Parameters.cpp).

Design history and rejected alternatives

This PR went through two earlier designs, kept in the commit history (no force-push):

• A Velocyto-named bulk mode (--quantMode GeneVelocyto): renamed, Velocyto being a single-cell notion (refactor(quant): rename GeneVelocyto to GeneSplicing).
• Dropping reads instead of adding targets (--quantTranscriptomePreMRNA BanRetainedIntron: do not project an unspliced read overlapping an intron that another isoform retains): removed. It also deleted the only evidence for genuinely expressed retained-intron isoforms (simulation: true share 0.333 estimated as 0.000) and pushed the exonic pre-mRNA reads onto the spliced isoforms (+45% on them). Giving the EM an unspliced target is the principled fix and is what splici does.
• Separate spliced / unspliced counting only (Velocyto-style): kept as the secondary GeneSplicing table, but it does not fix isoform quantification, which needs both kinds of targets in one reference.
• Other alternatives considered: one -I record per disjoint intron interval (splici's -I1, -I2...; loses pairs whose mates fall in two introns, so the intervals are concatenated into one target per gene); using GENCODE's transcript_type (annotation-specific; the rules here use structure only); tagging projections instead of adding targets (Salmon / RSEM ignore tags); building targets at genomeGenerate (would change the index format and freeze the flank).

Compatibility / divergence from STAR

• Without the new options, outputs are unchanged:
  • tests/bulk_unspliced.rs::new_options_leave_existing_outputs_unchanged: GeneSplicing + --quantTranscriptomeUnspliced None leave Aligned.out.sam (minus @PG/@CO), ReadsPerGene.out.tab, SJ.out.tab, Log.final.out (minus times) and the transcriptome BAM records identical; Intron leaves the genomic outputs identical and only appends @SQ lines after the annotated ones.
  • splicing_status_tag_in_genomic_sam: records with sp minus the tag equal the default records.
  • Against a main build on public data (scripts/bench_bulk_unspliced/compare_main.sh, 2 libraries x 2 M pairs, GRCh38 + GENCODE v50): Aligned.out.bam records and header (minus @PG/@CO), transcriptome BAM records, SJ.out.tab and Log.final.out identical, for the branch with and without GeneSplicing --quantTranscriptomeUnspliced None (4.4 M genomic and 33 to 37 M transcriptome records). This check uses --quantTranscriptomeSAMoutput BanSingleEnd because main aborts on this data with the default mode (fix above); with the default mode the branch differs from main exactly where main fails.
• All three options are non-STAR and documented in DIVERGENCE.md 1.4, including two small departures of the classifier port from classifyAlign (containment test on the true pair span; blocks sorted before the scan). Needs maintainer sign-off per CONTRIBUTING.
• Shared-file changes: QuantContext holds optional GeneCounts / GeneSplicing / sp parts (4 call sites in lib.rs); filter_and_project takes a fragment-level junction flag; TranscriptomeIndex gains an optional unspliced marker; SamAttributes gains SP (bit 15). Default hot path: unchanged work (new code behind Options).

Tests

• src/quant/transcriptome.rs: unspliced_mode_from_str, unspliced_intervals_intron_and_premrna (union of isoform introns incl. retained ones, flanks clipped and merged, gene body), single_exon_gene_has_no_intron_target, unspliced_targets_are_appended_in_gene_order (annotated order preserved), unspliced_target_projection (retained intron -> RI isoform + -I, spliced fragment -> spliced only, constitutive exon untouched, intron/exon boundary needs the flank, cassette exon -> inclusion isoform + -I), unspliced_target_on_minus_strand_is_flipped, targets_table_and_unspliced_fasta, softclip_rebuild_folds_both_clips_into_a_single_op.
• src/quant/splice_status.rs (19 tests): per-block calls, STAR's collapse rules, full classification on a synthetic GTF (retained intron, antisense gene inside an intron, boundary exons, protruding read, single-exon chrM gene, GENCODE _PAR_Y copies, pooled pair blocks, indel merging), table / summary writers, read_status_ignores_gene_boundaries, sp_tag_on_records.
• src/params/mod.rs: quant_gene_splicing_is_opt_in, sp_attribute_is_opt_in, unknown_quant_mode_is_rejected.
• tests/bulk_unspliced.rs (synthetic genome, 3 genes on both strands, each with a spliced isoform A and a retained-intron isoform R; reverse-stranded library of 300 mature A, 0 / 150 / 0 mature R and 400 pre-mRNA reads per gene):

  • simulated_total_rna_mixture: every pre-mRNA read inside the retained intron is written to both R and <gene>-I. Minimal effective-length EM on the transcriptome BAM:

    gene (truth RI share)	STAR projection A / R / -I, RI share	Intron	PreMRNA
    G1 (0.000)	259 / 306 / 0, 0.541	389 / 79 / 232, 0.169	269 / 47 / 384, 0.148
    G2 (0.333)	295 / 412 / 0, 0.582	426 / 172 / 252, 0.288	303 / 132 / 415, 0.304
    G3 (0.000)	275 / 294 / 0, 0.517	400 / 74 / 226, 0.156	284 / 43 / 373, 0.132

    Summed RI-share error 1.308 / 0.370 / 0.310; summed relative error on A 0.236 / 1.049 / 0.164; PreMRNA recovers the 400 simulated pre-mRNA reads per gene on -I within 15%. Intron cannot place the exonic part of pre-mRNA anywhere but the spliced isoforms (splici has the same limitation), which the test states.

  • splicing_status_tag_in_genomic_sam: pre-mRNA reads in exon 1 / the constitutive intron / the retained intron get S / U / A.

cargo fmt --check, cargo clippy --all-targets (0 warnings) and cargo test pass locally and in CI.

Public benchmark

Scripts and raw results: scripts/bench_bulk_unspliced/ (results_PRJEB57727.txt). GRCh38 primary assembly + GENCODE v50 comprehensive annotation (full genome index, no subset); ENA PRJEB57727 (whole blood in PAXgene, rRNA + globin depleted, reverse-stranded, 2x100 bp), 8 runs, first 4 M pairs of each split into two 2 M-pair pseudo-replicates (16 libraries). Each library aligned three times (STAR projection, Intron, PreMRNA), each transcriptome BAM quantified with salmon quant -a (salmon 2.8.0). macOS aarch64, 16 threads.

	STAR projection	Intron	PreMRNA
fragments in transcriptome BAM per input pair	0.617	0.917	0.924
Salmon mapping rate	0.986	0.983	0.985
Salmon fragments on -I targets	0	0.350	0.470
retained_intron share of spliced-target fragments (mean)	0.0537	0.0440	0.0512
slope of RI share vs unspliced fraction	0.151	0.107	0.111
Spearman rho, RI share vs unspliced fraction (n = 16; rep1 only n = 8)	0.92 (0.93)	0.92 (0.93)	0.92 (0.93)
pseudo-replicate r of within-gene isoform proportions (median of 8)	0.817	0.819	0.828
pseudo-replicate r of log avgTxLength (median of 8)	0.975	0.974	0.972
wall time per 2 M pairs (mean, noisy shared machine)	49 s	48 s	55 s
peak RSS	29.2 GiB	29.5 GiB	29.6 GiB

GeneSplicing (reverse column, picked from the data): unspliced fraction 0.29 to 0.58 across the 16 libraries (mean 0.39). genomeGenerate (full GRCh38 + GENCODE v50, 644 k transcripts): 926 s, peak RSS 20.7 GiB.

What this says, honestly:

• The bias the PR targets is real on public data: with STAR's projection the retained-intron share tracks the unspliced fraction almost perfectly across libraries.
• Unspliced targets reduce the retained-intron share (-18% with Intron) and its dependence on pre-mRNA content (slope -26% to -29%), and raise the fragments reaching the transcriptome BAM from 0.62 to 0.92 per pair, but they do not remove the rank correlation (rho stays 0.92). With Salmon's length-normalised EM, a short retained-intron isoform still explains local intronic coverage better than a long gene-body target; part of the correlation may also be biological (libraries richer in nascent RNA also carry more genuine retention). This is the same limit splici has.
• Pseudo-replicates share library prep, so they only measure sampling noise; isoform proportions improve slightly with PreMRNA (0.817 -> 0.828). The technical-replicate effect in the motivation (r ~ 0.05) cannot be tested on this dataset; the PR author will run the private benchmark and comment the numbers here.
• Limits of this benchmark: 2 M pairs per library; one tissue; the index was built with --sjdbOverhang 49 (ENA metadata suggested 2x50 bp, the reads are 2x100 bp), which is valid but not STAR's recommended value and sets the Intron flank to 49; wall times are noisy (other jobs on the machine).

Limits

• Intron vs PreMRNA: Intron counts the exonic part of pre-mRNA as mature and drops pairs with one mate in an exon and one in an intron (unless an isoform retains that intron); PreMRNA avoids both but makes every exonic read also compatible with -I, leaving the split to the EM.
• Overlapping genes: each gene gets its own -I; an intronic read of gene A inside gene B's body goes to both -I targets (EM decides). GeneSplicing drops multi-gene reads (N_multiGene, as STAR's Velocyto does); with GENCODE v50 that is 12% to 25% (mean 18%) of uniquely mapped pairs on the reverse column in the benchmark, which I attribute to overlapping lncRNA / readthrough annotation (not measured). The sp tag ignores gene boundaries and strand.
• Genes on several chromosomes (shared gene_id): the -I target uses the chromosome and strand of the gene's first transcript.
• Genuine intron retention cannot be separated from pre-mRNA read by read; the targets only let the EM try.
• Ambiguity rules of the classifier are STAR's, including the 6-base tolerance and "span in one model + only-exonic in another = spliced".
• Cost: target construction is a single pass over the transcript tables at startup; the -I FASTA is genome-sized, hence opt-in.

This PR was prepared with Claude Code and reviewed by the user.

🤖 Generated with Claude Code