diff --git a/CHANGELOG.md b/CHANGELOG.md index d66f6b9a..7d5ba48b 100644 --- a/CHANGELOG.md +++ b/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 `-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.) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 39f0e17e..690825bd 100644 --- a/DIVERGENCE.md +++ b/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 `-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. diff --git a/docs/astro.config.mjs b/docs/astro.config.mjs index 548198d1..b2479d75 100644 --- a/docs/astro.config.mjs +++ b/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' }, ], }, diff --git a/docs/src/content/docs/guides/bulk-total-rna.md b/docs/src/content/docs/guides/bulk-total-rna.md new file mode 100644 index 00000000..b3b45cf7 --- /dev/null +++ b/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, `-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 `-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 `-I`; `NH`, `HI` and `MAPQ` count all targets, and + Salmon or RSEM decide; +- a read inside a constitutive intron is written to `-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 `-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/`. diff --git a/docs/src/content/docs/guides/quantification.md b/docs/src/content/docs/guides/quantification.md index 1c23856e..a80dc2ff 100644 --- a/docs/src/content/docs/guides/quantification.md +++ b/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 `-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. diff --git a/docs/src/content/docs/reference/cli-parameters.md b/docs/src/content/docs/reference/cli-parameters.md index bbdd0f48..2f2daea3 100644 --- a/docs/src/content/docs/reference/cli-parameters.md +++ b/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 (`-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 diff --git a/docs/src/content/docs/reference/output-files.md b/docs/src/content/docs/reference/output-files.md index 0d62a67b..cf742412 100644 --- a/docs/src/content/docs/reference/output-files.md +++ b/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 `-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` diff --git a/scripts/bench_bulk_unspliced/README.md b/scripts/bench_bulk_unspliced/README.md new file mode 100644 index 00000000..36a2666f --- /dev/null +++ b/scripts/bench_bulk_unspliced/README.md @@ -0,0 +1,34 @@ +# Bulk total-RNA benchmark + +Public-data benchmark for `--quantTranscriptomeUnspliced` and +`--quantGeneSplicing Yes` (see the "Bulk total RNA-seq" guide). + +- Reference: GRCh38 primary assembly + GENCODE v50 comprehensive annotation + (reference chromosomes) + GENCODE v50 transcript FASTA. +- Reads: ENA project PRJEB57727 (whole blood in PAXgene, rRNA + globin + depleted total RNA, reverse-stranded, 2x100 bp), runs ERR10501066, + ERR10501099, ERR10501182, ERR10501185, ERR10501210, ERR10501288, + ERR10501299, ERR10501321. From each run the first 4 M read pairs are split + into two pseudo-replicates of 2 M pairs (`rep1`, `rep2`). Pseudo-replicates + only share library prep and sequencing, so they measure sampling noise, not + technical-replicate variation. +- The index is built with `--sjdbOverhang 49` (the ENA metadata suggested + 2x50 bp; the reads are 2x100 bp). This is valid but not STAR's recommended + read length - 1, and it sets the default `Intron` flank to 49 bases. +- Tools used for the committed results: salmon 2.8.0, samtools 1.24, + Python 3 with numpy + scipy; macOS aarch64, 16 threads, 128 GB RAM. + +```bash +export DATA=/path/to/bench +bash scripts/bench_bulk_unspliced/fetch_data.sh # ~6 GB download +cargo build --release +BIN=target/release/rustar-aligner THREADS=16 bash scripts/bench_bulk_unspliced/run.sh +python3 scripts/bench_bulk_unspliced/analyze.py "$DATA" # results_PRJEB57727.txt +BIN_MAIN=/path/to/main/rustar-aligner LIB=ERR10501185_rep1 \ + bash scripts/bench_bulk_unspliced/compare_main.sh # outputs identical to main +``` + +`run.sh` aligns every library three times (STAR projection, `Intron`, +`PreMRNA`) and runs `salmon quant -a` on each transcriptome BAM. +`results_PRJEB57727.txt` is the output of `analyze.py` for the run behind the +PR. diff --git a/scripts/bench_bulk_unspliced/analyze.py b/scripts/bench_bulk_unspliced/analyze.py new file mode 100644 index 00000000..0d2fc7fb --- /dev/null +++ b/scripts/bench_bulk_unspliced/analyze.py @@ -0,0 +1,226 @@ +#!/usr/bin/env python3 +"""Summarise the bulk total-RNA benchmark produced by run.sh. + + python3 scripts/bench_bulk_unspliced/analyze.py "$DATA" + +Per library (run x pseudo-replicate) and TranscriptomeSAM mode (star = STAR +projection, intron / premrna = --quantTranscriptomeUnspliced Intron / PreMRNA): + - share of Salmon-assigned spliced-target fragments on GENCODE + retained_intron isoforms (RI share); + - share of Salmon-assigned fragments on the -I targets; + - unspliced / ambiguous fractions from ReadsPerGeneSplicing.summary.tsv + (star runs; strand column picked from the data: most assigned reads); + - fragments in Aligned.toTranscriptome.out.bam per input pair, and the + Salmon mapping rate; + - wall time and peak RSS of the rustar-aligner run. +Across libraries: Spearman rho and least-squares slope of RI share against +unspliced fraction. +Between the two pseudo-replicates of each run: Pearson r of within-gene +isoform proportions (annotated isoforms only) and of log tximport-style +average transcript lengths, for genes with >= MIN_READS spliced fragments. +Needs numpy and scipy only. +""" +import json +import re +import sys +from pathlib import Path + +import numpy as np +from scipy.stats import pearsonr, spearmanr + +MIN_READS = 100 +MODES = ("star", "intron", "premrna") + +data = Path(sys.argv[1]) +out = data / "out" + +# transcript_id -> (gene_id, transcript_type) +tx = {} +with open(data / "ref" / "gencode.v50.annotation.gtf") as f: + for line in f: + if line.startswith("#"): + continue + c = line.split("\t", 9) + if c[2] != "transcript": + continue + a = c[8] + tid = re.search(r'transcript_id "([^"]+)"', a).group(1) + gid = re.search(r'gene_id "([^"]+)"', a).group(1) + tt = re.search(r'transcript_type "([^"]+)"', a).group(1) + tx[tid] = (gid, tt) + + +def quant(path): + names, efflen, reads = [], [], [] + with open(path) as f: + next(f) + for line in f: + n, _l, el, _tpm, nr = line.rstrip("\n").split("\t") + names.append(n) + efflen.append(float(el)) + reads.append(float(nr)) + return names, np.array(efflen), np.array(reads) + + +def time_log(path): + txt = path.read_text() + m = re.search(r"([\d.]+) real", txt) + wall = float(m.group(1)) if m else float("nan") + m = re.search(r"(\d+)\s+maximum resident set size", txt) # macOS, bytes + if m: + rss = int(m.group(1)) / 2**30 + else: + m = re.search(r"Maximum resident set size \(kbytes\): (\d+)", txt) # GNU + rss = int(m.group(1)) / 2**20 if m else float("nan") + m = re.search(r"Elapsed \(wall clock\) time.*: ([\d:.]+)", txt) + if m: + parts = [float(p) for p in m.group(1).split(":")] + wall = sum(p * 60**i for i, p in enumerate(reversed(parts))) + return wall, rss + + +def log_final(path): + d = {} + for line in path.read_text().splitlines(): + if "|" in line: + k, v = line.split("|", 1) + d[k.strip()] = v.strip() + return d + + +def splicing(path): + rows = {} + for line in path.read_text().splitlines()[1:]: + k, *v = line.split("\t") + rows[k] = v + assigned = [sum(int(rows[k][i]) for k in ("N_spliced", "N_unspliced", "N_ambiguous")) for i in range(3)] + col = 1 if assigned[1] > assigned[2] else 2 + return { + "strand": ["unstranded", "forward", "reverse"][col], + "unspliced": float(rows["fraction_unspliced"][col]), + "ambiguous": float(rows["fraction_ambiguous"][col]), + } + + +rows = [] +q = {} +for libdir in sorted(p for p in out.iterdir() if p.is_dir()): + lib = libdir.name + for mode in MODES: + d = libdir / mode + sf = d / "salmon" / "quant.sf" + if not sf.exists(): + continue + names, efflen, reads = quant(sf) + unspl = np.array([n.endswith("-I") and n not in tx for n in names]) + ri = np.array([tx.get(n, ("", ""))[1] == "retained_intron" for n in names]) + q[(lib, mode)] = (names, efflen, reads, unspl) + meta = json.loads((d / "salmon" / "aux_info" / "meta_info.json").read_text()) + lf = log_final(d / "Log.final.out") + n_in = int(lf["Number of input reads"]) + wall, rss = time_log(d / "time.log") + spliced_total = reads[~unspl].sum() + r = { + "lib": lib, + "mode": mode, + "input_pairs": n_in, + "unique_pct": lf["Uniquely mapped reads %"], + "trsam_frag_per_pair": meta["num_processed"] / n_in, + "salmon_rate": meta["percent_mapped"] / 100, + "ri_share": reads[ri].sum() / spliced_total, + "unspliced_target_share": reads[unspl].sum() / reads.sum(), + "wall_s": wall, + "rss_gib": rss, + } + vs = d / "ReadsPerGeneSplicing.summary.tsv" + if vs.exists(): + r.update(splicing(vs)) + rows.append(r) + +print("| library | mode | pairs | unique % | trSAM frag/pair | Salmon mapped | RI share | -I share | unspliced (GeneSplicing) | wall s | RSS GiB |") +print("|---|---|---|---|---|---|---|---|---|---|---|") +for r in rows: + print( + f"| {r['lib']} | {r['mode']} | {r['input_pairs']} | {r['unique_pct']} | " + f"{r['trsam_frag_per_pair']:.3f} | {r['salmon_rate']:.3f} | {r['ri_share']:.4f} | " + f"{r['unspliced_target_share']:.4f} | {r.get('unspliced', float('nan')):.4f} | " + f"{r['wall_s']:.0f} | {r['rss_gib']:.1f} |" + ) + +by = {m: {r["lib"]: r for r in rows if r["mode"] == m} for m in MODES} +libs = sorted(set.intersection(*(set(v) for v in by.values()))) +unspliced = {l: by["star"][l]["unspliced"] for l in libs} +print() +print("GeneSplicing strand column:", sorted({by["star"][l].get("strand") for l in libs})) +for subset, label in ((libs, "all libraries"), ([l for l in libs if l.endswith("rep1")], "rep1 only")): + if len(subset) < 3: + continue + u = [unspliced[l] for l in subset] + parts = [] + for m in MODES: + s = [by[m][l]["ri_share"] for l in subset] + slope = np.polyfit(u, s, 1)[0] + parts.append( + f"{m}: RI share mean {np.mean(s):.4f}, rho vs unspliced {spearmanr(s, u)[0]:.2f}, " + f"slope {slope:.3f}" + ) + print(f"{label} (n={len(subset)}), unspliced fraction mean {np.mean(u):.4f}: " + "; ".join(parts)) +for m in MODES: + print( + f"{m}: mean trSAM frag/pair {np.mean([by[m][l]['trsam_frag_per_pair'] for l in libs]):.3f}, " + f"Salmon mapped {np.mean([by[m][l]['salmon_rate'] for l in libs]):.3f}, " + f"-I share {np.mean([by[m][l]['unspliced_target_share'] for l in libs]):.4f}, " + f"wall {np.mean([by[m][l]['wall_s'] for l in libs]):.0f} s, " + f"peak RSS {np.max([by[m][l]['rss_gib'] for l in libs]):.1f} GiB" + ) + + +def rep_stats(mode): + iso_r, len_r, n_genes = [], [], [] + runs = sorted({l.rsplit("_rep", 1)[0] for l in libs}) + for run in runs: + a, b = q.get((run + "_rep1", mode)), q.get((run + "_rep2", mode)) + if not a or not b: + continue + names, efflen, ra, unspl = a + names_b, efflen_b, rb, _ = b + assert names == names_b + genes = {} + for i, n in enumerate(names): + if not unspl[i]: + genes.setdefault(tx.get(n, (n, "?"))[0], []).append(i) + pa, pb, la, lb = [], [], [], [] + for idx in genes.values(): + if len(idx) < 2: + continue + idx = np.array(idx) + ga, gb = ra[idx].sum(), rb[idx].sum() + if ga < MIN_READS or gb < MIN_READS: + continue + pa.extend(ra[idx] / ga) + pb.extend(rb[idx] / gb) + # tximport avgTxLength: abundance-weighted effective length + ta, tb = ra[idx] / efflen[idx], rb[idx] / efflen_b[idx] + la.append(np.log((ta / ta.sum() * efflen[idx]).sum())) + lb.append(np.log((tb / tb.sum() * efflen_b[idx]).sum())) + iso_r.append(pearsonr(pa, pb)[0]) + len_r.append(pearsonr(la, lb)[0]) + n_genes.append(len(la)) + return iso_r, len_r, n_genes + + +print() +for m in MODES: + iso_r, len_r, n_genes = rep_stats(m) + if iso_r: + print( + f"pseudo-replicates, {m}: isoform-proportion r median {np.median(iso_r):.3f} " + f"(range {min(iso_r):.3f}-{max(iso_r):.3f}); log avgTxLength r median {np.median(len_r):.3f} " + f"(range {min(len_r):.3f}-{max(len_r):.3f}); multi-isoform genes with >= {MIN_READS} " + f"fragments: median {int(np.median(n_genes))}" + ) + +idx_time = data / "index" / "time.log" +if idx_time.exists(): + w, m = time_log(idx_time) + print(f"\ngenomeGenerate: wall {w:.0f} s, peak RSS {m:.1f} GiB") diff --git a/scripts/bench_bulk_unspliced/compare_main.sh b/scripts/bench_bulk_unspliced/compare_main.sh new file mode 100644 index 00000000..755741e8 --- /dev/null +++ b/scripts/bench_bulk_unspliced/compare_main.sh @@ -0,0 +1,52 @@ +#!/usr/bin/env bash +# Check that this branch, run WITHOUT the new options, reproduces a main +# build's outputs on real data (and that adding --quantGeneSplicing / +# --quantTranscriptomeUnspliced None changes none of them). +# +# DATA=... BIN_MAIN=/path/to/main/rustar-aligner BIN=target/release/rustar-aligner \ +# LIB=ERR10501185_rep1 bash scripts/bench_bulk_unspliced/compare_main.sh +# +# Compares Aligned.out.bam records and header (minus @PG/@CO, which carry the +# command line), transcriptome BAM records, SJ.out.tab and Log.final.out +# (minus wall-clock lines). Needs samtools. The transcriptome BAM uses +# --quantTranscriptomeSAMoutput BanSingleEnd: with the default soft-clip +# extension, main aborts on reads such as 1S97M2S (fixed on this branch). +set -euo pipefail + +DATA=${DATA:?}; BIN_MAIN=${BIN_MAIN:?}; BIN=${BIN:-target/release/rustar-aligner} +LIB=${LIB:-ERR10501185_rep1}; THREADS=${THREADS:-16} +OUT=$DATA/compare_main/$LIB +r1=$DATA/reads/${LIB}_1.fq.gz; r2=$DATA/reads/${LIB}_2.fq.gz + +run() { # bin outdir extra... + local bin=$1 o=$2; shift 2 + mkdir -p "$o" + "$bin" --runThreadN "$THREADS" --genomeDir "$DATA/index" --readFilesIn "$r1" "$r2" \ + --outSAMtype BAM Unsorted --quantTranscriptomeSAMoutput BanSingleEnd \ + --outFileNamePrefix "$o/" "$@" > "$o/stdout.log" 2>&1 +} +run "$BIN_MAIN" "$OUT/main" --quantMode TranscriptomeSAM +run "$BIN" "$OUT/branch" --quantMode TranscriptomeSAM +run "$BIN" "$OUT/branch_new" --quantMode TranscriptomeSAM --quantGeneSplicing Yes \ + --quantTranscriptomeUnspliced None + +digest() { # dir + local d=$1 + { + samtools view -H "$d/Aligned.out.bam" | grep -v -e '^@PG' -e '^@CO' + samtools view "$d/Aligned.out.bam" + } | shasum | cut -c1-16 + samtools view "$d/Aligned.toTranscriptome.out.bam" | shasum | cut -c1-16 + shasum < "$d/SJ.out.tab" | cut -c1-16 + grep -v -e Started -e Finished -e 'Mapping speed' "$d/Log.final.out" | shasum | cut -c1-16 +} +status=0 +for v in branch branch_new; do + if diff <(digest "$OUT/main") <(digest "$OUT/$v") > /dev/null; then + echo "$LIB: main == $v (Aligned.out.bam, toTranscriptome, SJ, Log.final.out)" + else + echo "$LIB: main != $v"; status=1 + fi +done +echo "records: $(samtools view -c "$OUT/main/Aligned.out.bam") genome, $(samtools view -c "$OUT/main/Aligned.toTranscriptome.out.bam") transcriptome" +exit $status diff --git a/scripts/bench_bulk_unspliced/fetch_data.sh b/scripts/bench_bulk_unspliced/fetch_data.sh new file mode 100644 index 00000000..135f1041 --- /dev/null +++ b/scripts/bench_bulk_unspliced/fetch_data.sh @@ -0,0 +1,45 @@ +#!/usr/bin/env bash +# Fetch the public data for the bulk total-RNA benchmark +# (--quantTranscriptomeUnspliced, --quantGeneSplicing Yes). +# +# DATA=/path/to/bench bash scripts/bench_bulk_unspliced/fetch_data.sh +# +# Reference: GRCh38 primary assembly + GENCODE v50 comprehensive annotation +# (reference chromosomes) + matching transcript FASTA. +# Reads: 8 libraries of ENA project PRJEB57727 (whole blood in PAXgene, +# rRNA + globin depleted total RNA, 2x100 bp). Only the first +# 2 x NPAIRS read pairs of each run are streamed, then split into two +# pseudo-replicates of NPAIRS pairs each (first half / second half). +set -euo pipefail + +DATA=${DATA:?set DATA to a writable directory} +NPAIRS=${NPAIRS:-2000000} +RUNS=${RUNS:-"ERR10501185 ERR10501288 ERR10501210 ERR10501321 ERR10501066 ERR10501182 ERR10501299 ERR10501099"} +GENCODE=https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_50 + +mkdir -p "$DATA/ref" "$DATA/reads" +cd "$DATA/ref" +[ -n "${SKIP_REF:-}" ] || for f in GRCh38.primary_assembly.genome.fa.gz gencode.v50.annotation.gtf.gz gencode.v50.transcripts.fa.gz; do + [ -s "$f" ] || curl -sfSLO "$GENCODE/$f" +done +if [ -z "${SKIP_REF:-}" ]; then + [ -s GRCh38.primary_assembly.genome.fa ] || gunzip -k GRCh38.primary_assembly.genome.fa.gz + [ -s gencode.v50.annotation.gtf ] || gunzip -k gencode.v50.annotation.gtf.gz +fi + +cd "$DATA/reads" +lines=$((NPAIRS * 4)) +for run in $RUNS; do + [ -s "${run}_rep2_2.fq.gz" ] && continue + urls=$(curl -sf "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=$run&result=read_run&fields=fastq_ftp" | tail -1 | cut -f2) + mate=1 + for u in ${urls//;/ }; do + # head closes the pipe early; ignore curl's resulting write error. + (curl -sfL "https://$u" || true) | gzip -dc 2>/dev/null | head -n $((2 * lines)) > "${run}_${mate}.fq" || true + head -n $lines "${run}_${mate}.fq" | gzip -1 > "${run}_rep1_${mate}.fq.gz" + tail -n +$((lines + 1)) "${run}_${mate}.fq" | gzip -1 > "${run}_rep2_${mate}.fq.gz" + rm "${run}_${mate}.fq" + mate=$((mate + 1)) + done +done +ls -la "$DATA/reads" diff --git a/scripts/bench_bulk_unspliced/results_PRJEB57727.txt b/scripts/bench_bulk_unspliced/results_PRJEB57727.txt new file mode 100644 index 00000000..30849766 --- /dev/null +++ b/scripts/bench_bulk_unspliced/results_PRJEB57727.txt @@ -0,0 +1,63 @@ +| library | mode | pairs | unique % | trSAM frag/pair | Salmon mapped | RI share | -I share | unspliced (GeneSplicing) | wall s | RSS GiB | +|---|---|---|---|---|---|---|---|---|---|---| +| ERR10501066_rep1 | star | 2000000 | 91.12% | 0.663 | 0.988 | 0.0534 | 0.0000 | 0.3334 | 35 | 28.9 | +| ERR10501066_rep1 | intron | 2000000 | 91.12% | 0.920 | 0.984 | 0.0453 | 0.2972 | nan | 42 | 29.3 | +| ERR10501066_rep1 | premrna | 2000000 | 91.12% | 0.927 | 0.985 | 0.0522 | 0.4136 | nan | 47 | 29.5 | +| ERR10501066_rep2 | star | 2000000 | 91.12% | 0.661 | 0.988 | 0.0545 | 0.0000 | 0.3337 | 46 | 28.8 | +| ERR10501066_rep2 | intron | 2000000 | 91.12% | 0.918 | 0.984 | 0.0464 | 0.2986 | nan | 54 | 29.5 | +| ERR10501066_rep2 | premrna | 2000000 | 91.12% | 0.925 | 0.986 | 0.0534 | 0.4141 | nan | 59 | 29.4 | +| ERR10501099_rep1 | star | 2000000 | 89.97% | 0.647 | 0.988 | 0.0456 | 0.0000 | 0.3570 | 65 | 29.0 | +| ERR10501099_rep1 | intron | 2000000 | 89.97% | 0.917 | 0.985 | 0.0381 | 0.3131 | nan | 46 | 29.5 | +| ERR10501099_rep1 | premrna | 2000000 | 89.97% | 0.924 | 0.986 | 0.0442 | 0.4345 | nan | 38 | 29.5 | +| ERR10501099_rep2 | star | 2000000 | 90.06% | 0.645 | 0.988 | 0.0449 | 0.0000 | 0.3574 | 36 | 29.1 | +| ERR10501099_rep2 | intron | 2000000 | 90.06% | 0.916 | 0.985 | 0.0377 | 0.3146 | nan | 44 | 29.4 | +| ERR10501099_rep2 | premrna | 2000000 | 90.06% | 0.923 | 0.986 | 0.0439 | 0.4347 | nan | 60 | 29.5 | +| ERR10501182_rep1 | star | 2000000 | 91.80% | 0.602 | 0.984 | 0.0595 | 0.0000 | 0.3990 | 52 | 29.0 | +| ERR10501182_rep1 | intron | 2000000 | 91.80% | 0.915 | 0.980 | 0.0492 | 0.3664 | nan | 61 | 29.5 | +| ERR10501182_rep1 | premrna | 2000000 | 91.80% | 0.921 | 0.982 | 0.0561 | 0.4704 | nan | 54 | 29.6 | +| ERR10501182_rep2 | star | 2000000 | 91.82% | 0.600 | 0.984 | 0.0602 | 0.0000 | 0.4002 | 50 | 29.1 | +| ERR10501182_rep2 | intron | 2000000 | 91.82% | 0.913 | 0.981 | 0.0494 | 0.3680 | nan | 47 | 29.4 | +| ERR10501182_rep2 | premrna | 2000000 | 91.82% | 0.920 | 0.982 | 0.0565 | 0.4716 | nan | 52 | 29.5 | +| ERR10501185_rep1 | star | 2000000 | 88.81% | 0.648 | 0.987 | 0.0480 | 0.0000 | 0.3871 | 57 | 28.7 | +| ERR10501185_rep1 | intron | 2000000 | 88.81% | 0.922 | 0.983 | 0.0401 | 0.3172 | nan | 42 | 29.2 | +| ERR10501185_rep1 | premrna | 2000000 | 88.81% | 0.929 | 0.984 | 0.0493 | 0.4665 | nan | 64 | 29.3 | +| ERR10501185_rep2 | star | 2000000 | 88.91% | 0.645 | 0.987 | 0.0496 | 0.0000 | 0.3878 | 35 | 28.8 | +| ERR10501185_rep2 | intron | 2000000 | 88.91% | 0.919 | 0.983 | 0.0416 | 0.3188 | nan | 39 | 29.2 | +| ERR10501185_rep2 | premrna | 2000000 | 88.91% | 0.926 | 0.984 | 0.0507 | 0.4674 | nan | 62 | 29.4 | +| ERR10501210_rep1 | star | 2000000 | 87.67% | 0.734 | 0.990 | 0.0346 | 0.0000 | 0.2939 | 39 | 28.7 | +| ERR10501210_rep1 | intron | 2000000 | 87.67% | 0.931 | 0.987 | 0.0294 | 0.2256 | nan | 41 | 29.1 | +| ERR10501210_rep1 | premrna | 2000000 | 87.67% | 0.936 | 0.988 | 0.0367 | 0.3960 | nan | 58 | 29.3 | +| ERR10501210_rep2 | star | 2000000 | 87.75% | 0.732 | 0.990 | 0.0342 | 0.0000 | 0.2936 | 61 | 28.7 | +| ERR10501210_rep2 | intron | 2000000 | 87.75% | 0.929 | 0.987 | 0.0290 | 0.2264 | nan | 52 | 29.3 | +| ERR10501210_rep2 | premrna | 2000000 | 87.75% | 0.934 | 0.988 | 0.0359 | 0.3958 | nan | 77 | 29.3 | +| ERR10501288_rep1 | star | 2000000 | 93.08% | 0.519 | 0.981 | 0.0713 | 0.0000 | 0.4772 | 61 | 28.9 | +| ERR10501288_rep1 | intron | 2000000 | 93.08% | 0.910 | 0.981 | 0.0565 | 0.4598 | nan | 42 | 29.3 | +| ERR10501288_rep1 | premrna | 2000000 | 93.08% | 0.917 | 0.983 | 0.0625 | 0.5397 | nan | 62 | 29.5 | +| ERR10501288_rep2 | star | 2000000 | 92.48% | 0.513 | 0.981 | 0.0717 | 0.0000 | 0.4773 | 63 | 29.2 | +| ERR10501288_rep2 | intron | 2000000 | 92.48% | 0.900 | 0.981 | 0.0573 | 0.4607 | nan | 78 | 29.5 | +| ERR10501288_rep2 | premrna | 2000000 | 92.48% | 0.908 | 0.983 | 0.0637 | 0.5402 | nan | 54 | 29.5 | +| ERR10501299_rep1 | star | 2000000 | 89.87% | 0.694 | 0.990 | 0.0386 | 0.0000 | 0.3196 | 54 | 28.8 | +| ERR10501299_rep1 | intron | 2000000 | 89.87% | 0.926 | 0.987 | 0.0325 | 0.2670 | nan | 50 | 29.4 | +| ERR10501299_rep1 | premrna | 2000000 | 89.87% | 0.933 | 0.988 | 0.0385 | 0.4088 | nan | 80 | 29.3 | +| ERR10501299_rep2 | star | 2000000 | 89.98% | 0.692 | 0.990 | 0.0385 | 0.0000 | 0.3199 | 55 | 28.9 | +| ERR10501299_rep2 | intron | 2000000 | 89.98% | 0.925 | 0.986 | 0.0325 | 0.2684 | nan | 42 | 29.2 | +| ERR10501299_rep2 | premrna | 2000000 | 89.98% | 0.932 | 0.988 | 0.0387 | 0.4091 | nan | 42 | 29.4 | +| ERR10501321_rep1 | star | 2000000 | 93.09% | 0.440 | 0.977 | 0.0775 | 0.0000 | 0.5750 | 38 | 28.8 | +| ERR10501321_rep1 | intron | 2000000 | 93.09% | 0.911 | 0.980 | 0.0605 | 0.5516 | nan | 45 | 29.3 | +| ERR10501321_rep1 | premrna | 2000000 | 93.09% | 0.918 | 0.982 | 0.0693 | 0.6298 | nan | 39 | 29.1 | +| ERR10501321_rep2 | star | 2000000 | 92.47% | 0.434 | 0.977 | 0.0763 | 0.0000 | 0.5752 | 33 | 28.8 | +| ERR10501321_rep2 | intron | 2000000 | 92.47% | 0.901 | 0.980 | 0.0593 | 0.5525 | nan | 40 | 29.2 | +| ERR10501321_rep2 | premrna | 2000000 | 92.47% | 0.909 | 0.982 | 0.0675 | 0.6303 | nan | 35 | 29.2 | + +GeneSplicing strand column: ['reverse'] +all libraries (n=16), unspliced fraction mean 0.3930: star: RI share mean 0.0537, rho vs unspliced 0.92, slope 0.151; intron: RI share mean 0.0440, rho vs unspliced 0.92, slope 0.107; premrna: RI share mean 0.0512, rho vs unspliced 0.92, slope 0.111 +rep1 only (n=8), unspliced fraction mean 0.3928: star: RI share mean 0.0536, rho vs unspliced 0.93, slope 0.152; intron: RI share mean 0.0439, rho vs unspliced 0.93, slope 0.108; premrna: RI share mean 0.0511, rho vs unspliced 0.93, slope 0.113 +star: mean trSAM frag/pair 0.617, Salmon mapped 0.986, -I share 0.0000, wall 49 s, peak RSS 29.2 GiB +intron: mean trSAM frag/pair 0.917, Salmon mapped 0.983, -I share 0.3504, wall 48 s, peak RSS 29.5 GiB +premrna: mean trSAM frag/pair 0.924, Salmon mapped 0.985, -I share 0.4702, wall 55 s, peak RSS 29.6 GiB + +pseudo-replicates, star: isoform-proportion r median 0.817 (range 0.786-0.827); log avgTxLength r median 0.975 (range 0.967-0.978); multi-isoform genes with >= 100 fragments: median 1616 +pseudo-replicates, intron: isoform-proportion r median 0.819 (range 0.790-0.833); log avgTxLength r median 0.974 (range 0.963-0.976); multi-isoform genes with >= 100 fragments: median 1525 +pseudo-replicates, premrna: isoform-proportion r median 0.828 (range 0.803-0.840); log avgTxLength r median 0.972 (range 0.963-0.976); multi-isoform genes with >= 100 fragments: median 1462 + +genomeGenerate: wall 926 s, peak RSS 20.7 GiB diff --git a/scripts/bench_bulk_unspliced/run.sh b/scripts/bench_bulk_unspliced/run.sh new file mode 100644 index 00000000..92b3d4cc --- /dev/null +++ b/scripts/bench_bulk_unspliced/run.sh @@ -0,0 +1,72 @@ +#!/usr/bin/env bash +# Bulk total-RNA benchmark for --quantTranscriptomeUnspliced and +# --quantGeneSplicing Yes. Run fetch_data.sh first. +# +# DATA=/path/to/bench BIN=target/release/rustar-aligner THREADS=16 \ +# bash scripts/bench_bulk_unspliced/run.sh +# python3 scripts/bench_bulk_unspliced/analyze.py "$DATA" +# +# For every pseudo-replicate library it runs rustar-aligner three times: +# star/ : --quantMode TranscriptomeSAM --quantGeneSplicing Yes (STAR projection) +# intron/ : TranscriptomeSAM + --quantTranscriptomeUnspliced Intron +# premrna/ : TranscriptomeSAM + --quantTranscriptomeUnspliced PreMRNA +# then quantifies each transcriptome BAM with `salmon quant -a`, against the +# GENCODE transcripts (+ the -I sequences written by rustar-aligner). +# Wall time and peak RSS come from /usr/bin/time (-l on macOS, -v on Linux). +set -euo pipefail + +DATA=${DATA:?set DATA} +BIN=$(cd "$(dirname "${BIN:-target/release/rustar-aligner}")" && pwd)/$(basename "${BIN:-target/release/rustar-aligner}") +THREADS=${THREADS:-16} +MODES=${MODES:-"star intron premrna"} +REF=$DATA/ref +IDX=$DATA/index +OUT=$DATA/out +if [ "$(uname)" = Darwin ]; then TIMEFLAG=-l; else TIMEFLAG=-v; fi + +if [ ! -s "$IDX/SA" ]; then + mkdir -p "$IDX" + /usr/bin/time $TIMEFLAG "$BIN" --runMode genomeGenerate --runThreadN "$THREADS" \ + --genomeDir "$IDX" --genomeFastaFiles "$REF/GRCh38.primary_assembly.genome.fa" \ + --sjdbGTFfile "$REF/gencode.v50.annotation.gtf" --sjdbOverhang 49 \ + --outFileNamePrefix "$IDX/" > "$IDX/stdout.log" 2> "$IDX/time.log" +fi + +# Salmon's alignment mode needs FASTA names equal to the BAM @SQ names +# (GENCODE transcript_id); drop the |-separated suffix of GENCODE headers. +TXFA=$REF/gencode.v50.transcripts.ids.fa +[ -s "$TXFA" ] || gzip -dc "$REF/gencode.v50.transcripts.fa.gz" | sed 's/|.*//' > "$TXFA" +mkdir -p "$DATA/targets" + +for r1 in "$DATA"/reads/*_rep?_1.fq.gz; do + lib=$(basename "$r1" _1.fq.gz) + r2=${r1%_1.fq.gz}_2.fq.gz + for mode in $MODES; do + o=$OUT/$lib/$mode + [ -s "$o/salmon/quant.sf" ] && continue + mkdir -p "$o" + fa=$DATA/targets/$mode.fa + case $mode in + star) extra=(--quantMode TranscriptomeSAM --quantGeneSplicing Yes) ;; + intron) extra=(--quantMode TranscriptomeSAM --quantTranscriptomeUnspliced Intron) ;; + premrna) extra=(--quantMode TranscriptomeSAM --quantTranscriptomeUnspliced PreMRNA) ;; + esac + # The unspliced sequences depend only on the index and the flank: + # write them on the first run of each mode. + if [ "$mode" != star ] && [ ! -s "$fa" ]; then + extra+=(--quantTranscriptomeUnsplicedFasta Yes) + fi + /usr/bin/time $TIMEFLAG "$BIN" --runThreadN "$THREADS" --genomeDir "$IDX" \ + --readFilesIn "$r1" "$r2" \ + --outSAMtype None "${extra[@]}" \ + --outFileNamePrefix "$o/" > "$o/stdout.log" 2> "$o/time.log" + if [ "$mode" = star ]; then + [ -s "$fa" ] || ln -sf "$TXFA" "$fa" + elif [ ! -s "$fa" ]; then + cat "$TXFA" "$o/Aligned.toTranscriptome.unspliced.fa" > "$fa" + rm "$o/Aligned.toTranscriptome.unspliced.fa" + fi + salmon quant -a "$o/Aligned.toTranscriptome.out.bam" -t "$fa" -l A -p "$THREADS" \ + -o "$o/salmon" > "$o/salmon.log" 2>&1 + done +done diff --git a/src/lib.rs b/src/lib.rs index 1b0024f3..4fc0ed0b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -322,48 +322,106 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { let mut params = params.clone(); params.redefine_window_params(index.genome.n_genome); - // Build gene-count context if --quantMode GeneCounts was requested. - // GTF requirement is already validated in params.validate(). - let quant_ctx: Option> = - if params.quant_gene_counts() { - let gtf_path = params.sjdb_gtf_file.as_ref().unwrap(); - info!( - "quantMode GeneCounts: building gene annotation from {}", - gtf_path.display() - ); - let ctx = crate::quant::QuantContext::build( - gtf_path, - &index.genome, - ¶ms.sjdb_gtf_feature_exon, - ¶ms.sjdb_gtf_chr_prefix, - ¶ms.sjdb_gtf_tag_exon_parent_gene, - )?; - Some(std::sync::Arc::new(ctx)) - } else { - None - }; - // Use the transcriptome index loaded alongside the genome (populated // from transcriptInfo.tab / exonInfo.tab / geneInfo.tab at load time // — see GenomeIndex::load). Only wire it through to the pipeline when - // `--quantMode TranscriptomeSAM` is requested. - let tr_idx: Option> = - if params.quant_transcriptome_sam() { + // `--quantMode TranscriptomeSAM`, `--quantGeneSplicing Yes` or + // `--outSAMsplicingStatus Yes` is requested. + let tr_idx_all: Option> = + if params.quant_transcriptome_sam() + || params.quant_gene_splicing() + || params.out_sam_splicing_status() + { let tr = index.transcriptome.as_ref().ok_or_else(|| { anyhow::anyhow!( - "--quantMode TranscriptomeSAM requires a GTF-aware index; \ + "--quantMode TranscriptomeSAM, --quantGeneSplicing and --outSAMsplicingStatus \ + require a GTF-aware index; \ re-run genomeGenerate with --sjdbGTFfile or pass --sjdbGTFfile \ at alignReads so transcriptInfo.tab can be (re)built" ) })?; info!( - "quantMode TranscriptomeSAM: using {} transcripts from genome index", + "TranscriptomeSAM / GeneSplicing / sp tag: using {} transcripts from genome index", tr.n_transcripts() ); Some(std::sync::Arc::new(tr.clone())) } else { None }; + let tr_idx = tr_idx_all + .as_ref() + .filter(|_| params.quant_transcriptome_sam()) + .map(std::sync::Arc::clone); + + // --quantTranscriptomeUnspliced: append one `-I` target per gene + // to the TranscriptomeSAM index (the GeneSplicing classifier keeps using + // the annotated transcripts only), and describe every target. + let tr_idx = match tr_idx { + Some(tr) + if params.quant_transcriptome_unspliced + != crate::quant::transcriptome::QuantTranscriptomeUnspliced::None => + { + let flank = u64::try_from(params.quant_transcriptome_unspliced_flank) + .unwrap_or(u64::from(index.sjdb_overhang)); + let ext = tr.with_unspliced_targets(params.quant_transcriptome_unspliced, flank); + info!( + "quantTranscriptomeUnspliced {:?}: {} unspliced targets (flank {flank})", + params.quant_transcriptome_unspliced, + ext.n_transcripts() - tr.n_transcripts() + ); + let path = params.output_path("Aligned.toTranscriptome.targets.tsv"); + ext.write_targets_tsv(&path)?; + info!("Wrote {}", path.display()); + if params.quant_transcriptome_unspliced_fasta == "Yes" { + let path = params.output_path("Aligned.toTranscriptome.unspliced.fa"); + ext.write_unspliced_fasta(&path, &index.genome)?; + info!("Wrote {}", path.display()); + } + Some(std::sync::Arc::new(ext)) + } + other => other, + }; + + // Build the per-read quantification context if --quantMode GeneCounts + // and/or --quantGeneSplicing / --outSAMsplicingStatus was requested. GeneCounts' GTF requirement is + // already validated in params.validate(). + let quant_ctx: Option> = if params + .quant_gene_counts() + || params.quant_gene_splicing() + || params.out_sam_splicing_status() + { + let gene = if params.quant_gene_counts() { + let gtf_path = params.sjdb_gtf_file.as_ref().unwrap(); + info!( + "quantMode GeneCounts: building gene annotation from {}", + gtf_path.display() + ); + Some(crate::quant::GeneQuant::build( + gtf_path, + &index.genome, + ¶ms.sjdb_gtf_feature_exon, + ¶ms.sjdb_gtf_chr_prefix, + ¶ms.sjdb_gtf_tag_exon_parent_gene, + )?) + } else { + None + }; + let splicing = tr_idx_all + .as_ref() + .filter(|_| params.quant_gene_splicing()) + .map(|tr| crate::quant::SplicingQuant::new(std::sync::Arc::clone(tr))); + let splice_tag = tr_idx_all + .as_ref() + .filter(|_| params.out_sam_splicing_status()) + .map(std::sync::Arc::clone); + Some(std::sync::Arc::new(crate::quant::QuantContext { + gene, + splicing, + splice_tag, + })) + } else { + None + }; // SmartSeq has no barcodes/UMIs — a dedicated manifest-driven path. if params.solo_type == params::SoloType::SmartSeq { @@ -448,11 +506,12 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { crate::io::log::write_log_progress_out(&log_progress_path, &stats, time_start, time_finish)?; info!("Wrote {}", log_progress_path.display()); - // Write ReadsPerGene.out.tab if quantMode GeneCounts was requested. + // Write ReadsPerGene.out.tab / ReadsPerGeneSplicing.* for the requested + // --quantMode values. if let Some(ref ctx) = quant_ctx { - let quant_path = params.output_path("ReadsPerGene.out.tab"); - ctx.counts.write_output(&quant_path, &ctx.gene_ann)?; - info!("Wrote {}", quant_path.display()); + for path in ctx.write_outputs(|name| params.output_path(name))? { + info!("Wrote {}", path.display()); + } } info!("Alignment complete!"); @@ -1212,7 +1271,14 @@ fn build_transcriptome_records_se( for aln in transcripts { let bases: &[u8] = if aln.is_reverse { &rc } else { read_seq }; projected_all.extend(filter_and_project( - aln, bases, genome, tr_idx, lread, mode, params, + aln, + bases, + genome, + tr_idx, + lread, + mode, + params, + aln.n_junction > 0, )); } @@ -1281,8 +1347,11 @@ where let m2 = &pair.mate2_transcript; let m1_bases: &[u8] = if m1.is_reverse { &m1_rc } else { m1_seq }; let m2_bases: &[u8] = if m2.is_reverse { &m2_rc } else { m2_seq }; - let proj_m1 = filter_and_project(m1, m1_bases, genome, tr_idx, lread1, mode, params); - let proj_m2 = filter_and_project(m2, m2_bases, genome, tr_idx, lread2, mode, params); + let spliced = m1.n_junction + m2.n_junction > 0; + let proj_m1 = + filter_and_project(m1, m1_bases, genome, tr_idx, lread1, mode, params, spliced); + let proj_m2 = + filter_and_project(m2, m2_bases, genome, tr_idx, lread2, mode, params, spliced); let mut by_tr1: HashMap> = HashMap::new(); for p in &proj_m1 { @@ -1872,7 +1941,7 @@ fn align_reads_single_end( stats.record_alignment(0, max_multimaps); stats.record_unmapped_reason(crate::stats::UnmappedReason::Other); if let Some(ref q) = quant { - q.counts.count_se_read(&[], 0, &q.gene_ann); + q.count_se_read(&[], 0); } if output_unmapped { // Unmapped reads keep the full original read (STAR: clipped @@ -1938,8 +2007,7 @@ fn align_reads_single_end( // Gene-level quantification (lock-free atomic counts) if let Some(ref q) = quant { - q.counts - .count_se_read(&transcripts, n_for_mapq, &q.gene_ann); + q.count_se_read(&transcripts, n_for_mapq); } // Record junction statistics (per-read dedup, fix A) @@ -2013,6 +2081,15 @@ fn align_reads_single_end( params.out_sam_attributes, )?; } + // --outSAMattributes sp: splicing status tag. + if let Some(tx) = quant.as_ref().and_then(|q| q.splice_tag.as_ref()) + { + crate::quant::splice_status::tag_records_se( + &mut records, + &transcripts, + tx, + ); + } for record in records { buffer.push(record); } @@ -3241,7 +3318,7 @@ fn align_reads_paired_end( stats.record_alignment(0, max_multimaps); stats.record_unmapped_reason(crate::stats::UnmappedReason::Other); if let Some(ref q) = quant { - q.counts.count_pe_read(&[], true, false, &q.gene_ann); + q.count_pe_read(&[], true, false); } if output_unmapped { // Full original mates for unmapped pairs (STAR convention). @@ -3342,12 +3419,7 @@ fn align_reads_paired_end( // Dereference Box to get &PairedAlignment slice. let bm_deref: Vec<&crate::align::read_align::PairedAlignment> = both_mapped.iter().map(AsRef::as_ref).collect(); - q.counts.count_pe_read( - &bm_deref, - results.is_empty(), - has_half_mapped, - &q.gene_ann, - ); + q.count_pe_read(&bm_deref, results.is_empty(), has_half_mapped); } // Record junction statistics @@ -3511,6 +3583,14 @@ fn align_reads_paired_end( params.out_sam_attributes, )?; } + // --outSAMattributes sp: splicing status tag. + if let Some(tx) = quant.as_ref().and_then(|q| q.splice_tag.as_ref()) { + crate::quant::splice_status::tag_records_pe( + &mut records, + &paired_alns, + tx, + ); + } for record in records { buffer.push(record); } diff --git a/src/params/mod.rs b/src/params/mod.rs index a3539e08..2b896187 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1133,6 +1133,51 @@ pub struct Parameters { )] pub quant_transcriptome_sam_output: crate::quant::transcriptome::QuantTranscriptomeSAMoutput, + /// rustar-aligner extension (not in STAR), for `--quantMode + /// TranscriptomeSAM` on total RNA-seq: add one unspliced target per gene, + /// named `-I`, after the annotated transcripts. + /// * `None` (default): STAR behaviour, annotated transcripts only + /// * `Intron`: the gene's merged annotated introns plus flanks + /// * `PreMRNA`: the whole gene body + /// + /// Unspliced reads and pairs are projected onto every compatible target, + /// spliced and unspliced; fragments crossing a junction only onto + /// spliced ones. + #[arg(long = "quantTranscriptomeUnspliced", default_value = "None")] + pub quant_transcriptome_unspliced: crate::quant::transcriptome::QuantTranscriptomeUnspliced, + + /// Flank (bases) added on each side of every merged intron for + /// `--quantTranscriptomeUnspliced Intron`; -1 (default) uses the index's + /// sjdbOverhang, i.e. read length - 1 by STAR's convention. + #[arg( + long = "quantTranscriptomeUnsplicedFlank", + default_value_t = -1, + allow_hyphen_values = true + )] + pub quant_transcriptome_unspliced_flank: i64, + + /// rustar-aligner extension: `Yes` also writes the unspliced target + /// sequences to `Aligned.toTranscriptome.unspliced.fa` (Salmon's alignment + /// mode needs them next to the transcript FASTA). `No` by default, as + /// the file is of the order of the genome size. + #[arg(long = "quantTranscriptomeUnsplicedFasta", default_value = "No")] + pub quant_transcriptome_unspliced_fasta: String, + + /// rustar-aligner extension (not in STAR): `Yes` writes per-gene spliced / + /// unspliced / ambiguous counts for bulk data + /// (`ReadsPerGeneSplicing.out.tab` and `.summary.tsv`), using STARsolo's + /// spliced / unspliced rules. Needs a GTF-aware index. + #[arg(long = "quantGeneSplicing", default_value = "No")] + pub quant_gene_splicing: String, + + /// rustar-aligner extension (not in STAR): `Yes` adds an `sp:A` tag to the + /// genomic SAM/BAM records with the alignment's splicing status against + /// the annotation (S spliced, U unspliced, A ambiguous). Needs a + /// GTF-aware index. Kept out of `--outSAMattributes` so that STAR's + /// parameter keeps STAR's meaning. + #[arg(long = "outSAMsplicingStatus", default_value = "No")] + pub out_sam_splicing_status: String, + // ── Two-pass ──────────────────────────────────────────────────────── /// Two-pass mode: None or Basic #[arg(long = "twopassMode", default_value = "None")] @@ -1779,12 +1824,45 @@ impl Parameters { // for alignReads, GenomeIndex::load checks for the on-disk files // and surfaces a clear error if neither source is available. if params.run_mode() == RunMode::GenomeGenerate - && params.quant_transcriptome_sam() + && (params.quant_transcriptome_sam() || params.quant_gene_splicing()) && params.sjdb_gtf_file.is_none() { return Err(command.error( ErrorKind::MissingRequiredArgument, - "--quantMode TranscriptomeSAM requires --sjdbGTFfile at genomeGenerate", + "--quantMode TranscriptomeSAM and --quantGeneSplicing Yes require \ + --sjdbGTFfile at genomeGenerate", + )); + } + + for (flag, value) in [ + ("--quantGeneSplicing", ¶ms.quant_gene_splicing), + ("--outSAMsplicingStatus", ¶ms.out_sam_splicing_status), + ] { + if !matches!(value.as_str(), "Yes" | "No") { + return Err(command.error( + ErrorKind::InvalidValue, + format!("{flag} must be Yes or No, got '{value}'"), + )); + } + } + + // --quantTranscriptomeUnspliced* only make sense with TranscriptomeSAM. + if !matches!( + params.quant_transcriptome_unspliced_fasta.as_str(), + "Yes" | "No" + ) { + return Err(command.error( + ErrorKind::InvalidValue, + "--quantTranscriptomeUnsplicedFasta must be Yes or No", + )); + } + if params.quant_transcriptome_unspliced + != crate::quant::transcriptome::QuantTranscriptomeUnspliced::None + && !params.quant_transcriptome_sam() + { + return Err(command.error( + ErrorKind::MissingRequiredArgument, + "--quantTranscriptomeUnspliced requires --quantMode TranscriptomeSAM", )); } @@ -2153,6 +2231,18 @@ impl Parameters { self.quant_mode.iter().any(|m| m == "TranscriptomeSAM") } + /// Returns true if `--quantGeneSplicing Yes` (rustar-aligner extension: + /// bulk spliced / unspliced / ambiguous gene counts) was requested. + pub fn quant_gene_splicing(&self) -> bool { + self.quant_gene_splicing == "Yes" + } + + /// Returns true if `--outSAMsplicingStatus Yes` (rustar-aligner extension: + /// `sp:A` splicing-status tag on genomic records) was requested. + pub fn out_sam_splicing_status(&self) -> bool { + self.out_sam_splicing_status == "Yes" + } + /// True when a single-cell run is requested (`--soloType` != None). pub fn solo_enabled(&self) -> bool { self.solo_type != SoloType::None @@ -2841,6 +2931,45 @@ mod tests { assert!(try_parse(&["--readFilesIn", "r.fq", "--quantMode", "TranscriptomeSAM"]).is_ok()); } + #[test] + fn splicing_status_tag_is_opt_in() { + let p = try_parse(&["--readFilesIn", "r.fq", "--outSAMattributes", "All"]).unwrap(); + assert!(!p.out_sam_splicing_status()); + let p = try_parse(&["--readFilesIn", "r.fq", "--outSAMsplicingStatus", "Yes"]).unwrap(); + assert!(p.out_sam_splicing_status()); + assert!(p.out_sam_attributes.contains(SamAttributes::NH)); + // Not a value of STAR's --outSAMattributes. + assert!( + try_parse(&[ + "--readFilesIn", + "r.fq", + "--outSAMattributes", + "Standard", + "sp" + ]) + .is_err() + ); + assert!(try_parse(&["--readFilesIn", "r.fq", "--outSAMsplicingStatus", "yes"]).is_err()); + } + + #[test] + fn quant_gene_splicing_is_opt_in() { + let p = try_parse(&["--readFilesIn", "r.fq", "--quantMode", "TranscriptomeSAM"]).unwrap(); + assert!(!p.quant_gene_splicing()); + let p = try_parse(&[ + "--readFilesIn", + "r.fq", + "--quantMode", + "TranscriptomeSAM", + "--quantGeneSplicing", + "Yes", + ]) + .unwrap(); + assert!(p.quant_gene_splicing()); + assert!(p.quant_transcriptome_sam()); + assert!(try_parse(&["--readFilesIn", "r.fq", "--quantGeneSplicing", "1"]).is_err()); + } + #[test] fn quant_transcriptome_sam_enabled() { let p = try_parse(&["--readFilesIn", "r.fq", "--quantMode", "TranscriptomeSAM"]).unwrap(); diff --git a/src/params/sam.rs b/src/params/sam.rs index 0f363f61..de17a909 100644 --- a/src/params/sam.rs +++ b/src/params/sam.rs @@ -117,7 +117,7 @@ impl clap::Args for SamAttributes { .default_values(["Standard"]) .help( "SAM optional tags: Standard, All, None, or any combination of \ - NH HI AS NM nM MD jM jI XS RG vW vA vG.", + NH HI AS NM nM MD jM jI XS RG vW vA vG, and sp (rustar-aligner: splicing status).", ), ) } diff --git a/src/quant/mod.rs b/src/quant/mod.rs index 10451612..2705e3a1 100644 --- a/src/quant/mod.rs +++ b/src/quant/mod.rs @@ -1,3 +1,4 @@ +pub mod splice_status; /// Gene-level read quantification (`--quantMode GeneCounts`). /// /// Implements HTSeq-style "union" counting: a read counts toward a gene if any @@ -541,13 +542,31 @@ impl GeneCounts { // QuantContext — top-level bundle (passed as Arc to alignment loops) // --------------------------------------------------------------------------- -/// Bundles GeneAnnotation + GeneCounts for cheap Arc sharing across threads. +/// Bundles the per-read quantifications requested by `--quantMode` +/// (`GeneCounts`, `--quantGeneSplicing`, `--outSAMsplicingStatus`) for cheap Arc sharing across threads. pub struct QuantContext { + /// `--quantMode GeneCounts`: exon-union gene counts. + pub gene: Option, + /// `--quantGeneSplicing Yes`: spliced / unspliced / ambiguous per gene. + pub splicing: Option, + /// `--outSAMattributes sp`: annotated transcripts for the per-alignment + /// splicing-status tag. + pub splice_tag: Option>, +} + +/// GeneAnnotation + GeneCounts (`ReadsPerGene.out.tab`). +pub struct GeneQuant { pub gene_ann: GeneAnnotation, pub counts: GeneCounts, } -impl QuantContext { +/// Transcript models + counters for `ReadsPerGeneSplicing.out.tab`. +pub struct SplicingQuant { + pub transcriptome: std::sync::Arc, + pub counts: splice_status::SplicingCounts, +} + +impl GeneQuant { /// Build from a GTF file. Call once before alignment. pub fn build( gtf_path: &Path, @@ -561,7 +580,70 @@ impl QuantContext { let n = gene_ann.n_genes(); log::info!("quantMode GeneCounts: {n} genes loaded from GTF"); let counts = GeneCounts::new(n); - Ok(QuantContext { gene_ann, counts }) + Ok(GeneQuant { gene_ann, counts }) + } +} + +impl SplicingQuant { + /// Counters over the genes of an already-loaded transcriptome index. + pub fn new(transcriptome: std::sync::Arc) -> Self { + let counts = splice_status::SplicingCounts::new(transcriptome.gene_ids.len()); + SplicingQuant { + transcriptome, + counts, + } + } +} + +impl QuantContext { + /// Count a single-end read in every enabled quantification. + pub fn count_se_read(&self, transcripts: &[Transcript], n_for_mapq: usize) { + if let Some(g) = &self.gene { + g.counts.count_se_read(transcripts, n_for_mapq, &g.gene_ann); + } + if let Some(v) = &self.splicing { + v.counts.count_se_read(transcripts, &v.transcriptome); + } + } + + /// Count a read pair in every enabled quantification. + pub fn count_pe_read( + &self, + both_mapped: &[&PairedAlignment], + unmapped: bool, + half_mapped: bool, + ) { + if let Some(g) = &self.gene { + g.counts + .count_pe_read(both_mapped, unmapped, half_mapped, &g.gene_ann); + } + if let Some(v) = &self.splicing { + v.counts + .count_pe_read(both_mapped, unmapped, &v.transcriptome); + } + } + + /// Write the output files of every enabled quantification; returns the + /// paths written. + pub fn write_outputs( + &self, + output_path: impl Fn(&str) -> std::path::PathBuf, + ) -> Result, Error> { + let mut written = Vec::new(); + if let Some(g) = &self.gene { + let path = output_path("ReadsPerGene.out.tab"); + g.counts.write_output(&path, &g.gene_ann)?; + written.push(path); + } + if let Some(v) = &self.splicing { + let path = output_path("ReadsPerGeneSplicing.out.tab"); + v.counts.write_table(&path, &v.transcriptome)?; + written.push(path); + let path = output_path("ReadsPerGeneSplicing.summary.tsv"); + v.counts.write_summary(&path)?; + written.push(path); + } + Ok(written) } } diff --git a/src/quant/splice_status.rs b/src/quant/splice_status.rs new file mode 100644 index 00000000..c0b5cc2a --- /dev/null +++ b/src/quant/splice_status.rs @@ -0,0 +1,1179 @@ +//! Bulk spliced / unspliced / ambiguous gene quantification +//! (`--quantMode GeneSplicing`, a rustar-aligner extension). +//! +//! Provenance: the classification is a port of STARsolo's Velocyto feature; +//! only the STAR source file names below still carry that name. +//! +//! STARsolo's `--soloFeatures Velocyto` classifies every read against every +//! annotated transcript that contains it, then collapses those per-transcript +//! calls into one of three categories per gene. STAR only exposes this for +//! single-cell runs. This module ports the same two steps so they can run on +//! bulk data, where each read (or read pair) plays the role of one UMI: +//! +//! 1. Per transcript: [`align_to_transcript_min_overlap`], a port of STAR's +//! `alignToTranscriptMinOverlap` (`Transcriptome_classifyAlign.cpp`), with +//! the hard-coded `minOverlapMinusOne = 6` (velocyto.py's `MIN_FLANK = 5`) +//! and the 1 Mb intron cap. The transcript loop mirrors +//! `Transcriptome::classifyAlign`: only transcripts that fully contain the +//! alignment are tested, and the strand filter follows `--soloStrand` +//! semantics (`Forward` = read 1 on the transcript strand). +//! 2. Per gene: [`collapse_gene_category`], a port of the per-UMI collapse in +//! `SoloFeature_countVelocyto.cpp`: all transcripts must belong to one gene +//! (else the read is multi-gene and dropped, as STAR does); only-exonic +//! models give `spliced`, only-intronic / spanning models give +//! `unspliced`, anything mixed gives `ambiguous`. +//! +//! As in STAR, "spliced" means "compatible only with mature mRNA" (the read +//! need not cross a junction), "unspliced" means "requires pre-mRNA", and +//! "ambiguous" means "compatible with both", which is where reads falling in +//! a retained intron that is exonic in another isoform end up. +//! +//! The per-gene table has the three categories for each of the three strand +//! conventions of `ReadsPerGene.out.tab` (unstranded, forward, reverse), so a +//! single run serves any library type. +use std::io::Write as _; +use std::path::Path; +use std::sync::atomic::{AtomicU64, Ordering}; + +use noodles::sam::alignment::record::cigar::op::Kind; + +use crate::align::read_align::PairedAlignment; +use crate::align::transcript::Transcript; +use crate::error::Error; +use crate::quant::transcriptome::{TrExon, TranscriptomeIndex}; + +/// STAR's hard-coded `minOverlapMinusOne` for Velocyto (`classifyAlign`: +/// "6 is the hard code minOverlapMinusOne, to agree with velocyto's +/// MIN_FLANK=5"). +pub const MIN_OVERLAP_MINUS_ONE: u64 = 6; + +/// STAR's intron size cap in `alignToTranscriptMinOverlap`: a read that is +/// intronic in an intron longer than this is not called intronic for that +/// transcript ("prevents large introns from swallowing small genes"). +pub const MAX_INTRONIC_INTRON: u64 = 1_000_000; + +/// STAR's `AlignVsTranscript` enum (`AlignVsTranscript.h`). The discriminants +/// are the bit positions used in the per-transcript type bitset. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[repr(u8)] +pub enum AlignVsTranscript { + /// Purely intronic. + Intron = 0, + /// Exonic and intronic blocks, none crossing an exon/intron boundary. + ExonIntron = 1, + /// At least one block crosses an exon/intron boundary. + ExonIntronSpan = 2, + /// Purely exonic. + Concordant = 3, +} + +impl AlignVsTranscript { + /// STAR's `reAnn1` bitset for one transcript: the status bit, plus the + /// `Intron` and `Concordant` bits for a span ("span is also considered + /// intronic ... also considered exonic"). + pub fn type_bits(self) -> u8 { + let mut bits = 1u8 << (self as u8); + if self == Self::ExonIntronSpan { + bits |= 1 << (Self::Intron as u8); + bits |= 1 << (Self::Concordant as u8); + } + bits + } +} + +/// Splicing status of one read for one gene. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum SpliceStatus { + /// Compatible only with exonic (mature) models. + Spliced = 0, + /// Requires an intronic (pre-mRNA) model. + Unspliced = 1, + /// Compatible with both. + Ambiguous = 2, +} + +/// Outcome of classifying one read under one strand convention. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum ReadSplicing { + /// No transcript on the selected strand contains the read (or every + /// containing transcript rejected it). + NoFeature, + /// Containing transcripts belong to more than one gene. + MultiGene, + /// One gene, one category. + Gene(u32, SpliceStatus), +} + +/// Strand conventions, in output column order (same as `ReadsPerGene.out.tab`). +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum StrandMode { + /// Strand ignored. + Unstranded = 0, + /// Read 1 on the transcript strand (STARsolo `--soloStrand Forward`). + Forward = 1, + /// Read 1 opposite to the transcript strand (`--soloStrand Reverse`, + /// dUTP / TruSeq Stranded libraries). + Reverse = 2, +} + +impl StrandMode { + /// All three modes in output order. + pub const ALL: [StrandMode; 3] = [Self::Unstranded, Self::Forward, Self::Reverse]; + + fn name(self) -> &'static str { + match self { + Self::Unstranded => "unstranded", + Self::Forward => "forward", + Self::Reverse => "reverse", + } + } + + /// STAR's filter in `classifyAlign`: + /// `(trStr==1 ? aG.Str : 1-aG.Str) != pSolo.strand` rejects the transcript. + fn keeps(self, tr_is_reverse: bool, read_is_reverse: bool) -> bool { + match self { + Self::Unstranded => true, + Self::Forward => tr_is_reverse == read_is_reverse, + Self::Reverse => tr_is_reverse != read_is_reverse, + } + } +} + +/// A genomic alignment reduced to what the classifier reads: its +/// aligned blocks with indels merged (STAR expands a block across +/// `canonSJ` -1/-2), whether it has a splice junction (`sjYes`), and the +/// strand of read 1 (`aG.Str`). +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct AlignBlocks { + /// Chromosome index. + pub chr_idx: usize, + /// Blocks as absolute 0-based inclusive `(start, end)`, sorted by start. + pub blocks: Vec<(u64, u64)>, + /// Any mate crosses a splice junction. + pub has_junction: bool, + /// Strand of read 1. + pub is_reverse: bool, +} + +impl AlignBlocks { + /// Blocks of a single alignment (SE read, or one mate). + pub fn from_transcript(t: &Transcript) -> Self { + let mut blocks = Vec::new(); + push_blocks(t, &mut blocks); + blocks.sort_unstable(); + AlignBlocks { + chr_idx: t.chr_idx, + blocks, + has_junction: t.n_junction > 0, + is_reverse: t.is_reverse, + } + } + + /// Blocks of a read pair, both mates pooled like STAR's combined + /// paired-end `Transcript`; the strand is mate 1's. + pub fn from_pair(pair: &PairedAlignment) -> Self { + let mut blocks = Vec::new(); + push_blocks(&pair.mate1_transcript, &mut blocks); + push_blocks(&pair.mate2_transcript, &mut blocks); + blocks.sort_unstable(); + AlignBlocks { + chr_idx: pair.mate1_transcript.chr_idx, + blocks, + has_junction: pair.mate1_transcript.n_junction + pair.mate2_transcript.n_junction > 0, + is_reverse: pair.mate1_transcript.is_reverse, + } + } + + /// Leftmost aligned base (0-based). + fn start(&self) -> u64 { + self.blocks.first().map_or(0, |b| b.0) + } + + /// Rightmost aligned base (0-based, inclusive). + fn end_incl(&self) -> u64 { + self.blocks.iter().map(|b| b.1).max().unwrap_or(0) + } +} + +/// Append the indel-merged blocks of `t`: a new block starts only at a splice +/// junction (`N`). Falls back to the exon list when no CIGAR is attached. +fn push_blocks(t: &Transcript, out: &mut Vec<(u64, u64)>) { + let Some(first) = t.exons.first() else { + return; + }; + if t.cigar.is_empty() { + // No CIGAR (synthetic transcripts): exons are the blocks; merge + // across insertions (read gap, no genome gap). + let mut cur = (first.genome_start, first.genome_end); + for e in &t.exons[1..] { + if e.genome_start == cur.1 { + cur.1 = e.genome_end; + } else { + out.push((cur.0, cur.1 - 1)); + cur = (e.genome_start, e.genome_end); + } + } + out.push((cur.0, cur.1 - 1)); + return; + } + let mut pos = first.genome_start; + let mut block_start = pos; + for op in &t.cigar { + match op.kind() { + Kind::Match | Kind::SequenceMatch | Kind::SequenceMismatch | Kind::Deletion => { + pos += op.len() as u64; + } + Kind::Skip => { + if pos > block_start { + out.push((block_start, pos - 1)); + } + pos += op.len() as u64; + block_start = pos; + } + Kind::Insertion | Kind::SoftClip | Kind::HardClip | Kind::Pad => {} + } + } + if pos > block_start { + out.push((block_start, pos - 1)); + } +} + +/// Port of STAR's `alignToTranscriptMinOverlap` +/// (`Transcriptome_classifyAlign.cpp`). `exons` are the transcript's exons in +/// ascending order; the alignment must be contained in the transcript span +/// (the caller checks, as in STAR). Returns `None` for STAR's `-1`: a spliced +/// alignment touching an intron, or an intronic call inside an intron longer +/// than [`MAX_INTRONIC_INTRON`]. +pub fn align_to_transcript_min_overlap( + blocks: &[(u64, u64)], + has_junction: bool, + exons: &[TrExon], + min_overlap_minus_one: u64, +) -> Option { + let m = min_overlap_minus_one; + let n_ex = exons.len(); + let mut intronic = false; + let mut exonic = false; + let mut span = false; + + for &(bs, be) in blocks { + // STAR `binarySearch1` over the interleaved exon start/end array, + // halved: the last exon whose start is <= the block start. + let ex1 = exons + .partition_point(|e| e.genome_start <= bs) + .checked_sub(1)?; + if ex1 == n_ex - 1 { + // Reached the last exon: with the alignment inside the transcript + // (and blocks sorted), everything from here on is exonic. + exonic = true; + break; + } + if be - bs < m { + continue; // block too short to call + } + let e_e = exons[ex1].genome_end - 1; // exon1 end (inclusive) + let en_s = exons[ex1 + 1].genome_start; // exon2 start + let en_e = exons[ex1 + 1].genome_end - 1; // exon2 end (inclusive) + + if bs + m <= e_e { + // start is certainly in exon1 + if be <= e_e + m { + exonic = true; + } else { + span = true; + } + } else if bs + m < en_s { + // start is in intron1 + if be >= en_s + m { + span = true; + } else if be > e_e + m { + if en_s - e_e > MAX_INTRONIC_INTRON { + return None; + } + intronic = true; + } + } else if be > en_e + m { + // start too close to exon2 start; end certainly in intron2 + span = true; + } else if be >= en_s + m { + exonic = true; + } + + if has_junction && (intronic || span) { + return None; // a spliced alignment cannot overlap an intron + } + } + + Some(if span { + AlignVsTranscript::ExonIntronSpan + } else if !intronic { + AlignVsTranscript::Concordant + } else if exonic { + AlignVsTranscript::ExonIntron + } else { + AlignVsTranscript::Intron + }) +} + +/// Per-transcript type bits for every transcript that contains the +/// alignment, regardless of strand: `(transcript index, type bits)`. Mirrors +/// the transcript walk of STAR's `classifyAlign` (the strand filter is applied +/// later, per output column, by [`classify_read`]). +pub fn transcript_types(align: &AlignBlocks, idx: &TranscriptomeIndex, out: &mut Vec<(usize, u8)>) { + out.clear(); + if align.blocks.is_empty() || idx.n_transcripts() == 0 { + return; + } + let a_start = align.start(); + let a_end_excl = align.end_incl() + 1; + let upper = idx.tr_starts_sorted.partition_point(|&s| s <= a_start); + if upper == 0 { + return; + } + let mut i = upper - 1; + loop { + if idx.tr_end_max_sorted[i] < a_end_excl { + break; + } + let tr = idx.tr_order[i]; + if idx.tr_chr_idx[tr] == align.chr_idx + && idx.tr_start[tr] <= a_start + && idx.tr_end[tr] >= a_end_excl + && let Some(status) = align_to_transcript_min_overlap( + &align.blocks, + align.has_junction, + &idx.tr_exons[tr], + MIN_OVERLAP_MINUS_ONE, + ) + { + out.push((tr, status.type_bits())); + } + if i == 0 { + break; + } + i -= 1; + } +} + +/// Port of the per-UMI collapse in STAR's `SoloFeature_countVelocyto.cpp` +/// for one read: `types` are the `(transcript, type bits)` kept under one +/// strand convention. +pub fn collapse_gene_category( + types: impl IntoIterator, + idx: &TranscriptomeIndex, +) -> ReadSplicing { + let mut gene: Option = None; + let mut models = Models::default(); + for (tr, ty) in types { + let g = idx.tr_gene_idx[tr]; + match gene { + None => gene = Some(g), + Some(g0) if g0 != g => return ReadSplicing::MultiGene, + Some(_) => {} + } + models.add(ty); + } + match gene { + None => ReadSplicing::NoFeature, + Some(g) => ReadSplicing::Gene(g, models.status()), + } +} + +/// The model flags of STAR's `countVelocyto` collapse. +#[allow(clippy::struct_excessive_bools)] // mirrors STAR's four flags +struct Models { + exon: bool, + intron: bool, + span: bool, + mixed: bool, +} + +impl Default for Models { + fn default() -> Self { + // `spanModel` starts true and is and-ed over the transcripts. + Models { + exon: false, + intron: false, + span: true, + mixed: false, + } + } +} + +impl Models { + fn add(&mut self, ty: u8) { + const INTRON: u8 = 1 << AlignVsTranscript::Intron as u8; + const EXON_INTRON: u8 = 1 << AlignVsTranscript::ExonIntron as u8; + const SPAN: u8 = 1 << AlignVsTranscript::ExonIntronSpan as u8; + const CONCORDANT: u8 = 1 << AlignVsTranscript::Concordant as u8; + let has = |bit: u8| ty & bit != 0; + self.mixed |= ((has(INTRON) && has(CONCORDANT)) || has(EXON_INTRON)) && !has(SPAN); + self.span &= has(SPAN); + self.exon |= has(CONCORDANT) && !has(INTRON) && !has(EXON_INTRON); + self.intron |= has(INTRON) && !has(EXON_INTRON) && !has(CONCORDANT); + } + + fn status(&self) -> SpliceStatus { + if self.exon && !self.intron && !self.mixed { + SpliceStatus::Spliced + } else if self.span || ((self.intron || self.mixed) && !self.exon) { + SpliceStatus::Unspliced + } else { + SpliceStatus::Ambiguous + } + } +} + +/// Read-level splicing status for the `sp` SAM tag: the same collapse as +/// [`collapse_gene_category`], but over every annotated transcript that +/// contains the alignment on either strand, whatever its gene (a tag +/// describes the read, not a gene assignment). `None` when no transcript +/// contains it. +pub fn read_status(align: &AlignBlocks, idx: &TranscriptomeIndex) -> Option { + let mut types = Vec::new(); + transcript_types(align, idx, &mut types); + if types.is_empty() { + return None; + } + let mut models = Models::default(); + for (_, ty) in types { + models.add(ty); + } + Some(models.status()) +} + +impl SpliceStatus { + /// Value of the `sp:A` SAM tag. + pub fn tag_char(self) -> u8 { + match self { + SpliceStatus::Spliced => b'S', + SpliceStatus::Unspliced => b'U', + SpliceStatus::Ambiguous => b'A', + } + } +} + +/// Add `sp:A:{S,U,A}` to single-end records (one record per alignment, in +/// `transcripts` order, as built by `SamWriter::build_alignment_records`). +pub fn tag_records_se( + records: &mut [noodles::sam::alignment::RecordBuf], + transcripts: &[Transcript], + idx: &TranscriptomeIndex, +) { + for (rec, t) in records.iter_mut().zip(transcripts) { + if let Some(st) = read_status(&AlignBlocks::from_transcript(t), idx) { + insert_sp(rec, st); + } + } +} + +/// Add `sp:A:{S,U,A}` to paired-end records (mate 1 and mate 2 records per +/// pair, in `pairs` order, as built by `SamWriter::build_paired_records`); +/// both mates carry the fragment's status. +pub fn tag_records_pe( + records: &mut [noodles::sam::alignment::RecordBuf], + pairs: &[PairedAlignment], + idx: &TranscriptomeIndex, +) { + for (recs, pair) in records.chunks_mut(2).zip(pairs) { + if let Some(st) = read_status(&AlignBlocks::from_pair(pair), idx) { + for rec in recs { + insert_sp(rec, st); + } + } + } +} + +fn insert_sp(rec: &mut noodles::sam::alignment::RecordBuf, st: SpliceStatus) { + use noodles::sam::alignment::record::data::field::Tag; + use noodles::sam::alignment::record_buf::data::field::Value; + rec.data_mut() + .insert(Tag::new(b's', b'p'), Value::Character(st.tag_char())); +} + +/// Classify one uniquely mapped read (or pair) under the three strand +/// conventions, in [`StrandMode::ALL`] order. +pub fn classify_read(align: &AlignBlocks, idx: &TranscriptomeIndex) -> [ReadSplicing; 3] { + let mut types = Vec::new(); + transcript_types(align, idx, &mut types); + // STAR's per-UMI intersection works on transcript-sorted lists; for a + // single read the order does not change the result, but sort anyway so + // the multi-gene check sees the same first gene as STAR. + types.sort_unstable_by_key(|&(tr, _)| tr); + StrandMode::ALL.map(|mode| { + collapse_gene_category( + types + .iter() + .copied() + .filter(|&(tr, _)| mode.keeps(idx.tr_strand[tr] == 2, align.is_reverse)), + idx, + ) + }) +} + +// --------------------------------------------------------------------------- +// Counters + output +// --------------------------------------------------------------------------- + +/// Thread-safe counters for `ReadsPerGeneSplicing.out.tab`. +pub struct SplicingCounts { + /// Per gene: `[strand mode][category]`, flattened as `mode * 3 + cat`. + per_gene: Vec<[AtomicU64; 9]>, + /// Reads that did not map (including too-many-loci), as in GeneCounts. + pub n_unmapped: AtomicU64, + /// Multi-mapping reads (not classified, as STAR's Velocyto requires a + /// single alignment). + pub n_multimapping: AtomicU64, + /// Per strand mode: unique reads contained in no transcript. + pub n_no_feature: [AtomicU64; 3], + /// Per strand mode: unique reads whose transcripts span several genes. + pub n_multi_gene: [AtomicU64; 3], +} + +impl SplicingCounts { + /// Zeroed counters for `n_genes` genes. + pub fn new(n_genes: usize) -> Self { + SplicingCounts { + per_gene: (0..n_genes) + .map(|_| std::array::from_fn(|_| AtomicU64::new(0))) + .collect(), + n_unmapped: AtomicU64::new(0), + n_multimapping: AtomicU64::new(0), + n_no_feature: std::array::from_fn(|_| AtomicU64::new(0)), + n_multi_gene: std::array::from_fn(|_| AtomicU64::new(0)), + } + } + + fn record(&self, classes: &[ReadSplicing; 3]) { + for (mode, class) in classes.iter().enumerate() { + match *class { + ReadSplicing::NoFeature => { + self.n_no_feature[mode].fetch_add(1, Ordering::Relaxed); + } + ReadSplicing::MultiGene => { + self.n_multi_gene[mode].fetch_add(1, Ordering::Relaxed); + } + ReadSplicing::Gene(g, cat) => { + self.per_gene[g as usize][mode * 3 + cat as usize] + .fetch_add(1, Ordering::Relaxed); + } + } + } + } + + /// Count a single-end read, with GeneCounts' unmapped / multimapping + /// rules (`GeneCounts::count_se_read`). + pub fn count_se_read(&self, transcripts: &[Transcript], idx: &TranscriptomeIndex) { + match transcripts.len() { + 0 => { + self.n_unmapped.fetch_add(1, Ordering::Relaxed); + } + 1 => self.record(&classify_read( + &AlignBlocks::from_transcript(&transcripts[0]), + idx, + )), + _ => { + self.n_multimapping.fetch_add(1, Ordering::Relaxed); + } + } + } + + /// Count a read pair, with GeneCounts' rules (`GeneCounts::count_pe_read`: + /// half-mapped pairs count as unmapped). + pub fn count_pe_read( + &self, + both_mapped: &[&PairedAlignment], + unmapped: bool, + idx: &TranscriptomeIndex, + ) { + if unmapped || both_mapped.is_empty() { + self.n_unmapped.fetch_add(1, Ordering::Relaxed); + } else if both_mapped.len() > 1 { + self.n_multimapping.fetch_add(1, Ordering::Relaxed); + } else { + self.record(&classify_read(&AlignBlocks::from_pair(both_mapped[0]), idx)); + } + } + + fn get(&self, g: usize, mode: usize, cat: usize) -> u64 { + self.per_gene[g][mode * 3 + cat].load(Ordering::Relaxed) + } + + /// Totals per `[mode][category]`. + fn totals(&self) -> [[u64; 3]; 3] { + let mut t = [[0u64; 3]; 3]; + for g in 0..self.per_gene.len() { + for (mode, row) in t.iter_mut().enumerate() { + for (cat, v) in row.iter_mut().enumerate() { + *v += self.get(g, mode, cat); + } + } + } + t + } + + /// Write `ReadsPerGeneSplicing.out.tab`: a header line, then one line per + /// gene (`geneInfo.tab` order) with spliced / unspliced / ambiguous + /// counts for each strand convention. + pub fn write_table(&self, path: &Path, idx: &TranscriptomeIndex) -> Result<(), Error> { + let mut out = String::new(); + out.push_str("gene_id"); + for mode in StrandMode::ALL { + for cat in ["spliced", "unspliced", "ambiguous"] { + out.push('\t'); + out.push_str(mode.name()); + out.push('_'); + out.push_str(cat); + } + } + out.push('\n'); + for (g, gene_id) in idx.gene_ids.iter().enumerate() { + out.push_str(gene_id); + for mode in 0..3 { + for cat in 0..3 { + out.push('\t'); + out.push_str(&self.get(g, mode, cat).to_string()); + } + } + out.push('\n'); + } + std::fs::write(path, out).map_err(|e| Error::io(e, path)) + } + + /// Write `ReadsPerGeneSplicing.summary.tsv`: read accounting and the + /// spliced / unspliced / ambiguous shares for each strand convention. + pub fn write_summary(&self, path: &Path) -> Result<(), Error> { + let mut f = std::fs::File::create(path).map_err(|e| Error::io(e, path))?; + let io = |e| Error::io(e, path); + let nu = self.n_unmapped.load(Ordering::Relaxed); + let nm = self.n_multimapping.load(Ordering::Relaxed); + let t = self.totals(); + let nf: [u64; 3] = std::array::from_fn(|m| self.n_no_feature[m].load(Ordering::Relaxed)); + let ng: [u64; 3] = std::array::from_fn(|m| self.n_multi_gene[m].load(Ordering::Relaxed)); + + writeln!(f, "metric\tunstranded\tforward\treverse").map_err(io)?; + let row = |f: &mut std::fs::File, name: &str, v: [u64; 3]| { + writeln!(f, "{name}\t{}\t{}\t{}", v[0], v[1], v[2]) + }; + row(&mut f, "N_unmapped", [nu; 3]).map_err(io)?; + row(&mut f, "N_multimapping", [nm; 3]).map_err(io)?; + row(&mut f, "N_noFeature", nf).map_err(io)?; + row(&mut f, "N_multiGene", ng).map_err(io)?; + for (cat, name) in ["N_spliced", "N_unspliced", "N_ambiguous"] + .iter() + .enumerate() + { + row(&mut f, name, std::array::from_fn(|m| t[m][cat])).map_err(io)?; + } + for (cat, name) in [ + "fraction_spliced", + "fraction_unspliced", + "fraction_ambiguous", + ] + .iter() + .enumerate() + { + let frac: [String; 3] = std::array::from_fn(|m| { + let assigned: u64 = t[m].iter().sum(); + if assigned == 0 { + "NA".to_string() + } else { + format!("{:.4}", t[m][cat] as f64 / assigned as f64) + } + }); + writeln!(f, "{name}\t{}\t{}\t{}", frac[0], frac[1], frac[2]).map_err(io)?; + } + Ok(()) + } +} + +// --------------------------------------------------------------------------- +// Tests +// --------------------------------------------------------------------------- + +#[cfg(test)] +mod tests { + use super::*; + use crate::align::transcript::Exon; + use crate::genome::Genome; + use crate::junction::gtf::GtfRecord; + use noodles::sam::alignment::record::cigar::Op; + use std::collections::HashMap; + + const V_I: u8 = 1 << AlignVsTranscript::Intron as u8; + const V_EI: u8 = 1 << AlignVsTranscript::ExonIntron as u8; + const V_C: u8 = 1 << AlignVsTranscript::Concordant as u8; + + fn exons(v: &[(u64, u64)]) -> Vec { + // (start, end_exclusive) + let mut cum = 0u32; + v.iter() + .map(|&(s, e)| { + let ex = TrExon { + genome_start: s, + genome_end: e, + ex_len_cum: cum, + }; + cum += (e - s) as u32; + ex + }) + .collect() + } + + /// Three-exon model: [100,200) [300,400) [500,600). + fn model() -> Vec { + exons(&[(100, 200), (300, 400), (500, 600)]) + } + + fn call(blocks: &[(u64, u64)], sj: bool) -> Option { + align_to_transcript_min_overlap(blocks, sj, &model(), MIN_OVERLAP_MINUS_ONE) + } + + #[test] + fn min_overlap_exonic_intronic_span() { + use AlignVsTranscript::*; + assert_eq!(call(&[(120, 169)], false), Some(Concordant)); + assert_eq!(call(&[(220, 269)], false), Some(Intron)); + // Crosses the exon1 / intron1 boundary by more than 6 bases. + assert_eq!(call(&[(170, 219)], false), Some(ExonIntronSpan)); + // Crosses intron1 / exon2. + assert_eq!(call(&[(270, 319)], false), Some(ExonIntronSpan)); + // Exonic block plus an intronic block (pair), no span. + assert_eq!(call(&[(120, 169), (220, 269)], false), Some(ExonIntron)); + // Overhang of <= 6 bases into the intron is tolerated (MIN_FLANK). + assert_eq!(call(&[(156, 205)], false), Some(Concordant)); + // Reaching the last exon short-circuits to exonic. + assert_eq!(call(&[(520, 569)], false), Some(Concordant)); + } + + #[test] + fn min_overlap_spliced_alignment_touching_intron_is_rejected() { + // Spliced read whose second block runs into intron2. + assert_eq!(call(&[(170, 199), (300, 449)], true), None); + // Spliced read with both blocks exonic stays concordant. + assert_eq!( + call(&[(150, 199), (300, 349)], true), + Some(AlignVsTranscript::Concordant) + ); + } + + #[test] + fn min_overlap_huge_intron_is_not_intronic() { + let ex = exons(&[(0, 100), (2_000_000, 2_000_100)]); + assert_eq!( + align_to_transcript_min_overlap(&[(500_000, 500_049)], false, &ex, 6), + None + ); + } + + #[test] + fn span_type_bits_include_intron_and_concordant() { + let b = AlignVsTranscript::ExonIntronSpan.type_bits(); + assert_eq!(b, 0b1101); + assert_eq!(AlignVsTranscript::Intron.type_bits(), V_I); + } + + // --- Gene-level collapse and full classification on a GTF --------------- + + fn genome() -> Genome { + // chr1 (10 kb), chrM (1 kb), chrX (2 kb), chrY (2 kb) + let lens = [10_000u64, 1_000, 2_000, 2_000]; + let mut starts = vec![0u64]; + for l in lens { + starts.push(starts.last().unwrap() + l); + } + let total = *starts.last().unwrap(); + Genome { + transform_blocks: None, + sequence: vec![0u8; total as usize].into(), + n_genome: total, + n_genome_real: total, + n_chr_real: 4, + chr_start: starts, + chr_length: lens.to_vec(), + chr_name: ["chr1", "chrM", "chrX", "chrY"].map(String::from).to_vec(), + } + } + + fn rec(chr: &str, s: u64, e: u64, strand: char, gene: &str, tr: &str) -> GtfRecord { + let mut attributes = HashMap::new(); + attributes.insert("gene_id".to_string(), gene.to_string()); + attributes.insert("transcript_id".to_string(), tr.to_string()); + GtfRecord { + seqname: chr.to_string(), + feature: "exon".to_string(), + start: s, + end: e, + strand, + attributes, + } + } + + /// Synthetic annotation (GTF 1-based inclusive): + /// - G1 (+): T1a exons 1001-1200, 1501-1700, 2001-2200 (fully spliced); + /// T1b exons 1001-1700 (retains intron1), 2001-2200. + /// - G2 (-): antisense gene inside G1's intron 2: 1751-1950, single exon. + /// - G3 (+): 5001-5100, 5201-5300 (boundary tests). + /// - MT1 (+) on chrM: single-exon 101-600. + /// - PAR gene on chrX and its GENCODE `_PAR_Y` copy on chrY. + fn index() -> TranscriptomeIndex { + let recs = vec![ + rec("chr1", 1001, 1200, '+', "G1", "T1a"), + rec("chr1", 1501, 1700, '+', "G1", "T1a"), + rec("chr1", 2001, 2200, '+', "G1", "T1a"), + rec("chr1", 1001, 1700, '+', "G1", "T1b"), + rec("chr1", 2001, 2200, '+', "G1", "T1b"), + rec("chr1", 1751, 1950, '-', "G2", "T2"), + rec("chr1", 5001, 5100, '+', "G3", "T3"), + rec("chr1", 5201, 5300, '+', "G3", "T3"), + rec("chrM", 101, 600, '+', "MT1", "TM"), + rec("chrX", 101, 300, '+', "PAR1", "TP"), + rec("chrX", 501, 700, '+', "PAR1", "TP"), + rec("chrY", 101, 300, '+', "PAR1_PAR_Y", "TP_PAR_Y"), + rec("chrY", 501, 700, '+', "PAR1_PAR_Y", "TP_PAR_Y"), + ]; + TranscriptomeIndex::from_gtf_exons(&recs, &genome()).unwrap() + } + + fn gene(idx: &TranscriptomeIndex, id: &str) -> u32 { + idx.gene_ids.iter().position(|g| g == id).unwrap() as u32 + } + + /// SE transcript from 0-based chr-relative blocks `[s, e)`. + fn aln(g: &Genome, chr: usize, blocks: &[(u64, u64)], rev: bool) -> Transcript { + use noodles::sam::alignment::record::cigar::op::Kind as K; + let off = g.chr_start[chr]; + let mut cigar = Vec::new(); + let mut ex = Vec::new(); + let mut rpos = 0usize; + for (i, &(s, e)) in blocks.iter().enumerate() { + if i > 0 { + cigar.push(Op::new(K::Skip, (s - blocks[i - 1].1) as usize)); + } + cigar.push(Op::new(K::Match, (e - s) as usize)); + ex.push(Exon { + genome_start: off + s, + genome_end: off + e, + read_start: rpos, + read_end: rpos + (e - s) as usize, + i_frag: 0, + }); + rpos += (e - s) as usize; + } + Transcript { + chr_idx: chr, + genome_start: off + blocks[0].0, + genome_end: off + blocks.last().unwrap().1, + is_reverse: rev, + exons: ex, + cigar, + score: 0, + n_mismatch: 0, + n_gap: 0, + n_junction: (blocks.len() - 1) as u32, + junction_motifs: vec![], + junction_annotated: vec![], + } + } + + fn classify(t: &Transcript, idx: &TranscriptomeIndex) -> [ReadSplicing; 3] { + classify_read(&AlignBlocks::from_transcript(t), idx) + } + + #[test] + fn constitutive_exon_read_is_spliced() { + let g = genome(); + let idx = index(); + let g1 = gene(&idx, "G1"); + let r = classify(&aln(&g, 0, &[(2050, 2100)], false), &idx); + assert_eq!(r[0], ReadSplicing::Gene(g1, SpliceStatus::Spliced)); + assert_eq!(r[1], ReadSplicing::Gene(g1, SpliceStatus::Spliced)); + // Reverse library: a + read is antisense to G1, no feature. + assert_eq!(r[2], ReadSplicing::NoFeature); + } + + #[test] + fn junction_read_is_spliced() { + let g = genome(); + let idx = index(); + let g1 = gene(&idx, "G1"); + let r = classify(&aln(&g, 0, &[(1650, 1700), (2000, 2050)], false), &idx); + assert_eq!(r[1], ReadSplicing::Gene(g1, SpliceStatus::Spliced)); + } + + #[test] + fn retained_intron_read_is_ambiguous() { + // Inside intron1 of T1a, exonic in T1b: compatible with both. + let g = genome(); + let idx = index(); + let g1 = gene(&idx, "G1"); + let r = classify(&aln(&g, 0, &[(1300, 1350)], false), &idx); + assert_eq!(r[1], ReadSplicing::Gene(g1, SpliceStatus::Ambiguous)); + } + + #[test] + fn constitutive_intron_read_is_unspliced() { + // Intron2 of G1 (1700..2000), sense strand, clear of G2 (1750..1950). + let g = genome(); + let idx = index(); + let g1 = gene(&idx, "G1"); + let r = classify(&aln(&g, 0, &[(1955, 1990)], false), &idx); + assert_eq!(r[1], ReadSplicing::Gene(g1, SpliceStatus::Unspliced)); + // Exon/intron boundary read of a constitutive exon: unspliced. + let r = classify(&aln(&g, 0, &[(1680, 1730)], false), &idx); + assert_eq!(r[1], ReadSplicing::Gene(g1, SpliceStatus::Unspliced)); + } + + #[test] + fn antisense_overlapping_gene_is_split_by_strand() { + // A read in G2's exon, which lies in G1's intron 2 on the other strand. + let g = genome(); + let idx = index(); + let g1 = gene(&idx, "G1"); + let g2 = gene(&idx, "G2"); + // Read on - strand = sense for G2. + let r = classify(&aln(&g, 0, &[(1800, 1850)], true), &idx); + assert_eq!(r[0], ReadSplicing::MultiGene); // unstranded: G1 intron + G2 exon + assert_eq!(r[1], ReadSplicing::Gene(g2, SpliceStatus::Spliced)); + assert_eq!(r[2], ReadSplicing::Gene(g1, SpliceStatus::Unspliced)); + } + + #[test] + fn read_outside_transcripts_is_no_feature() { + let g = genome(); + let idx = index(); + let r = classify(&aln(&g, 0, &[(8000, 8050)], false), &idx); + assert_eq!(r, [ReadSplicing::NoFeature; 3]); + // Protruding past a transcript end: not contained, no feature. + let r = classify(&aln(&g, 0, &[(5280, 5330)], false), &idx); + assert_eq!(r, [ReadSplicing::NoFeature; 3]); + } + + #[test] + fn boundary_exon_reads() { + let g = genome(); + let idx = index(); + let g3 = gene(&idx, "G3"); + // First base of the first exon to within the exon. + let r = classify(&aln(&g, 0, &[(5000, 5050)], false), &idx); + assert_eq!(r[1], ReadSplicing::Gene(g3, SpliceStatus::Spliced)); + // Last bases of the last exon. + let r = classify(&aln(&g, 0, &[(5250, 5300)], false), &idx); + assert_eq!(r[1], ReadSplicing::Gene(g3, SpliceStatus::Spliced)); + } + + #[test] + fn single_exon_chrm_gene_is_spliced() { + let g = genome(); + let idx = index(); + let mt = gene(&idx, "MT1"); + let r = classify(&aln(&g, 1, &[(200, 250)], false), &idx); + assert_eq!(r[1], ReadSplicing::Gene(mt, SpliceStatus::Spliced)); + } + + #[test] + fn gencode_par_y_copies_are_separate_genes() { + let g = genome(); + let idx = index(); + let x = gene(&idx, "PAR1"); + let y = gene(&idx, "PAR1_PAR_Y"); + let rx = classify(&aln(&g, 2, &[(350, 400)], false), &idx); + let ry = classify(&aln(&g, 3, &[(350, 400)], false), &idx); + assert_eq!(rx[1], ReadSplicing::Gene(x, SpliceStatus::Unspliced)); + assert_eq!(ry[1], ReadSplicing::Gene(y, SpliceStatus::Unspliced)); + } + + #[test] + fn pair_blocks_are_pooled() { + // Mate 1 exonic in exon2, mate 2 intronic in intron2: ExonIntron in + // every G1 model -> unspliced. + let g = genome(); + let idx = index(); + let g1 = gene(&idx, "G1"); + let pair = PairedAlignment { + mate1_transcript: aln(&g, 0, &[(1550, 1600)], false), + mate2_transcript: aln(&g, 0, &[(1955, 1990)], true), + mate1_region: (0, 50), + mate2_region: (0, 35), + is_proper_pair: true, + insert_size: 440, + combined_wt_score: 0, + }; + let r = classify_read(&AlignBlocks::from_pair(&pair), &idx); + assert_eq!(r[1], ReadSplicing::Gene(g1, SpliceStatus::Unspliced)); + } + + #[test] + fn indels_do_not_split_blocks() { + use noodles::sam::alignment::record::cigar::op::Kind as K; + let g = genome(); + let mut t = aln(&g, 0, &[(2050, 2100)], false); + t.cigar = vec![ + Op::new(K::SoftClip, 3), + Op::new(K::Match, 20), + Op::new(K::Deletion, 2), + Op::new(K::Match, 10), + Op::new(K::Insertion, 1), + Op::new(K::Match, 18), + ]; + let b = AlignBlocks::from_transcript(&t); + assert_eq!(b.blocks, vec![(2050, 2099)]); + assert!(!b.has_junction); + } + + #[test] + fn collapse_matches_star_rules() { + let idx = index(); + let t1a = idx.tr_ids.iter().position(|t| t == "T1a").unwrap(); + let t1b = idx.tr_ids.iter().position(|t| t == "T1b").unwrap(); + let t2 = idx.tr_ids.iter().position(|t| t == "T2").unwrap(); + let g1 = gene(&idx, "G1"); + let span = AlignVsTranscript::ExonIntronSpan.type_bits(); + let cases: &[(&[(usize, u8)], ReadSplicing)] = &[ + (&[], ReadSplicing::NoFeature), + (&[(t1a, V_C), (t2, V_C)], ReadSplicing::MultiGene), + ( + &[(t1a, V_C), (t1b, V_C)], + ReadSplicing::Gene(g1, SpliceStatus::Spliced), + ), + ( + &[(t1a, V_I), (t1b, V_I)], + ReadSplicing::Gene(g1, SpliceStatus::Unspliced), + ), + ( + &[(t1a, V_EI)], + ReadSplicing::Gene(g1, SpliceStatus::Unspliced), + ), + ( + &[(t1a, span), (t1b, span)], + ReadSplicing::Gene(g1, SpliceStatus::Unspliced), + ), + ( + &[(t1a, V_I), (t1b, V_C)], + ReadSplicing::Gene(g1, SpliceStatus::Ambiguous), + ), + // STAR: a span in one model plus an only-exonic model is "spliced" + // (the span sets neither exonModel, intronModel nor mixedModel). + ( + &[(t1a, span), (t1b, V_C)], + ReadSplicing::Gene(g1, SpliceStatus::Spliced), + ), + ]; + for (types, want) in cases { + assert_eq!( + collapse_gene_category(types.iter().copied(), &idx), + *want, + "{types:?}" + ); + } + } + + #[test] + fn counts_table_and_summary() { + let g = genome(); + let idx = index(); + let counts = SplicingCounts::new(idx.gene_ids.len()); + // spliced (constitutive exon), unspliced (intron 2), ambiguous + // (retained intron), unmapped, multimapper, no feature. + counts.count_se_read(&[aln(&g, 0, &[(2050, 2100)], false)], &idx); + counts.count_se_read(&[aln(&g, 0, &[(1955, 1990)], false)], &idx); + counts.count_se_read(&[aln(&g, 0, &[(1300, 1350)], false)], &idx); + counts.count_se_read(&[], &idx); + let t = aln(&g, 0, &[(2050, 2100)], false); + counts.count_se_read(&[t.clone(), t], &idx); + counts.count_se_read(&[aln(&g, 0, &[(8000, 8050)], false)], &idx); + + let dir = tempfile::tempdir().unwrap(); + let tab = dir.path().join("t.tab"); + counts.write_table(&tab, &idx).unwrap(); + let tab = std::fs::read_to_string(tab).unwrap(); + let mut lines = tab.lines(); + assert_eq!( + lines.next().unwrap(), + "gene_id\tunstranded_spliced\tunstranded_unspliced\tunstranded_ambiguous\t\ + forward_spliced\tforward_unspliced\tforward_ambiguous\t\ + reverse_spliced\treverse_unspliced\treverse_ambiguous" + ); + assert_eq!(lines.next().unwrap(), "G1\t1\t1\t1\t1\t1\t1\t0\t0\t0"); + assert_eq!(lines.count(), idx.gene_ids.len() - 1); + + let sum = dir.path().join("s.tsv"); + counts.write_summary(&sum).unwrap(); + let sum = std::fs::read_to_string(sum).unwrap(); + let get = |k: &str| { + sum.lines() + .find(|l| l.split('\t').next() == Some(k)) + .unwrap() + .split('\t') + .skip(1) + .map(str::to_string) + .collect::>() + }; + assert_eq!(get("N_unmapped"), ["1", "1", "1"]); + assert_eq!(get("N_multimapping"), ["1", "1", "1"]); + assert_eq!(get("N_noFeature"), ["1", "1", "4"]); + assert_eq!(get("N_unspliced"), ["1", "1", "0"]); + assert_eq!(get("fraction_unspliced"), ["0.3333", "0.3333", "NA"]); + } + + #[test] + fn read_status_ignores_gene_boundaries() { + let g = genome(); + let idx = index(); + let st = |t: &Transcript| read_status(&AlignBlocks::from_transcript(t), &idx); + assert_eq!( + st(&aln(&g, 0, &[(2050, 2100)], false)), + Some(SpliceStatus::Spliced) + ); + assert_eq!( + st(&aln(&g, 0, &[(1955, 1990)], false)), + Some(SpliceStatus::Unspliced) + ); + assert_eq!( + st(&aln(&g, 0, &[(1300, 1350)], false)), + Some(SpliceStatus::Ambiguous) + ); + // Exonic in G2, intronic in G1 (other strand): both models count. + assert_eq!( + st(&aln(&g, 0, &[(1800, 1850)], true)), + Some(SpliceStatus::Ambiguous) + ); + assert_eq!(st(&aln(&g, 0, &[(8000, 8050)], false)), None); + } + + #[test] + fn sp_tag_on_records() { + use noodles::sam::alignment::record::data::field::Tag; + use noodles::sam::alignment::record_buf::data::field::Value; + let g = genome(); + let idx = index(); + let trs = vec![ + aln(&g, 0, &[(1955, 1990)], false), + aln(&g, 0, &[(8000, 8050)], false), + ]; + let mut recs = vec![noodles::sam::alignment::RecordBuf::default(); 2]; + tag_records_se(&mut recs, &trs, &idx); + assert_eq!( + recs[0].data().get(&Tag::new(b's', b'p')), + Some(&Value::Character(b'U')) + ); + assert_eq!(recs[1].data().get(&Tag::new(b's', b'p')), None); + let pair = PairedAlignment { + mate1_transcript: aln(&g, 0, &[(1550, 1600)], false), + mate2_transcript: aln(&g, 0, &[(2050, 2100)], true), + mate1_region: (0, 50), + mate2_region: (0, 50), + is_proper_pair: true, + insert_size: 550, + combined_wt_score: 0, + }; + let mut recs = vec![noodles::sam::alignment::RecordBuf::default(); 2]; + tag_records_pe(&mut recs, &[pair], &idx); + for r in &recs { + assert_eq!( + r.data().get(&Tag::new(b's', b'p')), + Some(&Value::Character(b'S')) + ); + } + } +} diff --git a/src/quant/transcriptome.rs b/src/quant/transcriptome.rs index 9444a0de..1283b906 100644 --- a/src/quant/transcriptome.rs +++ b/src/quant/transcriptome.rs @@ -95,6 +95,67 @@ impl QuantTranscriptomeSAMoutput { } } +/// `--quantTranscriptomeUnspliced` (rustar-aligner extension, not in STAR): +/// add one unspliced target per gene to the transcriptome so that pre-mRNA +/// reads have somewhere to go besides retained-intron isoforms. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum QuantTranscriptomeUnspliced { + /// STAR behaviour: annotated transcripts only. + #[default] + None, + /// `-I`: the union of the gene's annotated introns, each merged + /// interval extended by the flank on both sides (clipped to the gene body) + /// and merged again, concatenated in genome order (splici-style). + Intron, + /// `-I`: the whole gene body (first to last annotated base). + PreMRNA, +} + +impl std::str::FromStr for QuantTranscriptomeUnspliced { + type Err = String; + fn from_str(s: &str) -> Result { + match s { + "None" => Ok(Self::None), + "Intron" => Ok(Self::Intron), + "PreMRNA" => Ok(Self::PreMRNA), + _ => Err(format!( + "unknown --quantTranscriptomeUnspliced '{s}'; expected 'None', 'Intron' or 'PreMRNA'" + )), + } + } +} + +/// Unspliced intervals of one gene: `(chr, strand, merged [start, end))`. +pub type GeneIntervals = (usize, u8, Vec<(u64, u64)>); + +/// Suffix of the unspliced target names (`-I`, as in splici). +pub const UNSPLICED_SUFFIX: &str = "-I"; + +/// Marks a [`TranscriptomeIndex`] extended with unspliced targets: indices +/// `first..n_transcripts()` are the `-I` targets. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct UnsplicedTargets { + /// Index of the first unspliced target (= number of annotated transcripts). + pub first: usize, + /// How the targets were built. + pub mode: QuantTranscriptomeUnspliced, + /// Flank added around each merged intron (`Intron` mode). + pub flank: u64, +} + +/// Sort and merge overlapping or touching `[start, end)` intervals. +fn merge_intervals(mut iv: Vec<(u64, u64)>) -> Vec<(u64, u64)> { + iv.sort_unstable(); + let mut out: Vec<(u64, u64)> = Vec::with_capacity(iv.len()); + for (s, e) in iv { + match out.last_mut() { + Some(last) if s <= last.1 => last.1 = last.1.max(e), + _ => out.push((s, e)), + } + } + out +} + /// Per-transcript exon in absolute genome coordinates (0-based half-open), /// paired with the cumulative transcript-space length of all preceding exons /// (STAR's `exLenCum`). @@ -152,6 +213,9 @@ pub struct TranscriptomeIndex { pub tr_starts_sorted: Vec, /// Running max of `tr_end` along `tr_order` (STAR's `trEmax`). pub tr_end_max_sorted: Vec, + /// Set when unspliced `-I` targets were appended + /// (`--quantTranscriptomeUnspliced Intron|PreMRNA`). + pub unspliced: Option, } impl TranscriptomeIndex { @@ -357,6 +421,7 @@ impl TranscriptomeIndex { tr_order, tr_starts_sorted, tr_end_max_sorted, + unspliced: None, }) } @@ -377,6 +442,192 @@ impl TranscriptomeIndex { self.tr_ids.len() } + /// Whether transcript `tr` is an unspliced `-I` target. + pub fn is_unspliced_target(&self, tr: usize) -> bool { + self.unspliced.is_some_and(|u| tr >= u.first) + } + + /// Unspliced intervals of every gene, indexed by gene: `(chr, strand, + /// merged [start, end) intervals)`, or `None` for a gene with no interval + /// (a single-exon gene in `Intron` mode). A gene is taken on the + /// chromosome and strand of its first transcript; transcripts elsewhere + /// (e.g. pseudo-autosomal copies sharing a gene_id) are ignored here. + pub fn unspliced_intervals( + &self, + mode: QuantTranscriptomeUnspliced, + flank: u64, + ) -> Vec> { + let n_genes = self.gene_ids.len(); + let mut home: Vec> = vec![None; n_genes]; + let mut body: Vec<(u64, u64)> = vec![(u64::MAX, 0); n_genes]; + let mut introns: Vec> = vec![Vec::new(); n_genes]; + let n_annot = self.unspliced.map_or(self.n_transcripts(), |u| u.first); + for tr in 0..n_annot { + let g = self.tr_gene_idx[tr] as usize; + let key = (self.tr_chr_idx[tr], self.tr_strand[tr]); + if *home[g].get_or_insert(key) != key { + continue; + } + body[g].0 = body[g].0.min(self.tr_start[tr]); + body[g].1 = body[g].1.max(self.tr_end[tr]); + for w in self.tr_exons[tr].windows(2) { + if w[1].genome_start > w[0].genome_end { + introns[g].push((w[0].genome_end, w[1].genome_start)); + } + } + } + (0..n_genes) + .map(|g| { + let (chr, strand) = home[g]?; + let (bs, be) = body[g]; + let iv = match mode { + QuantTranscriptomeUnspliced::None => return None, + QuantTranscriptomeUnspliced::PreMRNA => vec![(bs, be)], + QuantTranscriptomeUnspliced::Intron => { + let merged = merge_intervals(std::mem::take(&mut introns[g])); + let flanked = merged + .into_iter() + .map(|(s, e)| (s.saturating_sub(flank).max(bs), (e + flank).min(be))) + .collect(); + merge_intervals(flanked) + } + }; + (!iv.is_empty()).then_some((chr, strand, iv)) + }) + .collect() + } + + /// A copy of this index with one unspliced target per gene appended after + /// the annotated transcripts, in gene order, named `-I`. The + /// annotated transcripts keep their indices and relative order. + #[must_use] + pub fn with_unspliced_targets(&self, mode: QuantTranscriptomeUnspliced, flank: u64) -> Self { + let mut idx = self.clone(); + let first = idx.n_transcripts(); + for (g, iv) in self + .unspliced_intervals(mode, flank) + .into_iter() + .enumerate() + { + let Some((chr, strand, iv)) = iv else { + continue; + }; + let mut cum = 0u32; + let exons: Vec = iv + .iter() + .map(|&(s, e)| { + let ex = TrExon { + genome_start: s, + genome_end: e, + ex_len_cum: cum, + }; + cum = cum.saturating_add((e - s) as u32); + ex + }) + .collect(); + idx.tr_ids + .push(format!("{}{UNSPLICED_SUFFIX}", self.gene_ids[g])); + idx.tr_chr_idx.push(chr); + idx.tr_strand.push(strand); + idx.tr_gene_idx.push(g as u32); + idx.tr_start.push(iv[0].0); + idx.tr_end.push(iv[iv.len() - 1].1); + idx.tr_exons.push(exons); + idx.tr_length.push(cum); + idx.tr_exi.push(0); + } + // Rebuild the sorted views. The sort is stable and annotated + // transcripts come first, so their relative order is unchanged. + let n_tr = idx.n_transcripts(); + let mut order: Vec = (0..n_tr).collect(); + order.sort_by(|&a, &b| { + idx.tr_start[a] + .cmp(&idx.tr_start[b]) + .then_with(|| idx.tr_end[a].cmp(&idx.tr_end[b])) + }); + let mut cum = 0u32; + for &i in &order { + idx.tr_exi[i] = cum; + cum = cum.saturating_add(idx.tr_exons[i].len() as u32); + } + idx.tr_starts_sorted = order.iter().map(|&i| idx.tr_start[i]).collect(); + let mut m = 0u64; + idx.tr_end_max_sorted = order + .iter() + .map(|&i| { + m = m.max(idx.tr_end[i]); + m + }) + .collect(); + idx.tr_order = order; + idx.unspliced = Some(UnsplicedTargets { first, mode, flank }); + idx + } + + /// Write `target_id, gene_id, gene_name, status, length` for every target + /// (annotated transcripts `spliced`, `-I` targets `unspliced`), + /// in BAM `@SQ` order: a tx2gene table for tximport / summarisation. + pub fn write_targets_tsv(&self, path: &Path) -> Result<(), Error> { + use std::fmt::Write as _; + let mut out = String::from("target_id\tgene_id\tgene_name\tstatus\tlength\n"); + for tr in 0..self.n_transcripts() { + let g = self.tr_gene_idx[tr] as usize; + let status = if self.is_unspliced_target(tr) { + "unspliced" + } else { + "spliced" + }; + let _ = writeln!( + out, + "{}\t{}\t{}\t{status}\t{}", + self.tr_ids[tr], + self.gene_ids.get(g).map_or("", String::as_str), + self.gene_names.get(g).map_or("", String::as_str), + self.tr_length[tr] + ); + } + std::fs::write(path, out).map_err(|e| Error::io(e, path)) + } + + /// Write the sequences of the unspliced targets (in transcript + /// orientation: reverse-complemented for `-` genes) as FASTA, 60 bases per + /// line. Salmon's alignment mode needs them next to the transcript FASTA. + pub fn write_unspliced_fasta(&self, path: &Path, genome: &Genome) -> Result<(), Error> { + let Some(u) = self.unspliced else { + return Ok(()); + }; + let f = std::fs::File::create(path).map_err(|e| Error::io(e, path))?; + let mut w = std::io::BufWriter::new(f); + let io = |e| Error::io(e, path); + const ACGTN: [u8; 5] = *b"ACGTN"; + for tr in u.first..self.n_transcripts() { + let mut seq: Vec = Vec::with_capacity(self.tr_length[tr] as usize); + for ex in &self.tr_exons[tr] { + for pos in ex.genome_start..ex.genome_end { + seq.push(ACGTN[(genome.sequence.base(pos as usize) as usize).min(4)]); + } + } + if self.tr_strand[tr] == 2 { + seq.reverse(); + for b in &mut seq { + *b = match *b { + b'A' => b'T', + b'C' => b'G', + b'G' => b'C', + b'T' => b'A', + x => x, + }; + } + } + writeln!(w, ">{}", self.tr_ids[tr]).map_err(io)?; + for chunk in seq.chunks(60) { + w.write_all(chunk).map_err(io)?; + w.write_all(b"\n").map_err(io)?; + } + } + w.flush().map_err(io) + } + /// Load from STAR-compatible index files in `dir`. /// /// Requires `transcriptInfo.tab`, `exonInfo.tab`, `geneInfo.tab` to be @@ -472,6 +723,7 @@ impl TranscriptomeIndex { tr_order, tr_starts_sorted, tr_end_max_sorted, + unspliced: None, }) } @@ -1150,6 +1402,12 @@ fn align_to_one_transcript( /// `read_bases_align_orientation` is the read in genome base encoding /// (A=0,C=1,G=2,T=3,N=4) already reversed/complemented when the alignment is /// on the reverse strand — STAR's `Read1[roStr==0 ? 0 : 2]`. +/// +/// `fragment_spliced` tells whether the read (SE) or either mate (PE) crosses +/// a splice junction. It is only read when `idx.unspliced` is set +/// (`--quantTranscriptomeUnspliced Intron|PreMRNA`): a spliced fragment is +/// processed RNA, so its projections onto `-I` targets are dropped. +#[allow(clippy::too_many_arguments)] pub fn filter_and_project( align: &Transcript, read_bases_align_orientation: &[u8], @@ -1158,6 +1416,7 @@ pub fn filter_and_project( lread: u32, mode: QuantTranscriptomeSAMoutput, params: &Parameters, + fragment_spliced: bool, ) -> Vec { if !mode.allow_indels() && align.n_gap > 0 { return Vec::new(); @@ -1172,7 +1431,12 @@ pub fn filter_and_project( } }; - align_to_transcripts(&align_for_projection, idx, lread) + let mut projected = align_to_transcripts(&align_for_projection, idx, lread); + if fragment_spliced && idx.unspliced.is_some() { + // A fragment crossing a junction is processed RNA: spliced targets only. + projected.retain(|p| !idx.is_unspliced_target(p.chr_idx)); + } + projected } fn has_soft_clip(align: &Transcript) -> bool { @@ -2375,6 +2639,7 @@ mod tests { 40, QuantTranscriptomeSAMoutput::BanSingleEndBanIndelsExtendSoftclip, ¶ms, + false, ); assert_eq!(results.len(), 0); } @@ -2402,6 +2667,7 @@ mod tests { 40, QuantTranscriptomeSAMoutput::BanSingleEnd, ¶ms, + false, ); assert_eq!(results.len(), 1); } @@ -2449,6 +2715,7 @@ mod tests { 44, QuantTranscriptomeSAMoutput::BanSingleEndBanIndelsExtendSoftclip, ¶ms, + false, ); assert_eq!(results.len(), 1); let r = &results[0]; @@ -2505,6 +2772,7 @@ mod tests { 44, QuantTranscriptomeSAMoutput::BanSingleEndBanIndelsExtendSoftclip, ¶ms, + false, ); // 4 extension mismatches > 2 → rejected assert_eq!(results.len(), 0); @@ -2535,6 +2803,7 @@ mod tests { 44, QuantTranscriptomeSAMoutput::BanSingleEnd, ¶ms, + false, ); assert_eq!(results.len(), 1); // CIGAR preserves the soft-clip. @@ -2607,4 +2876,226 @@ mod tests { c(&[(Kind::Match, 100)]) ); } + + /// G1 (chr1) has a retained-intron isoform; G2 (chr2) only has a cassette + /// exon and an alternative 5' splice site. + fn retained_intron_index() -> TranscriptomeIndex { + let genome = make_genome(); + let mut gtf = Vec::new(); + for (tr, exons) in [ + ("T1", &[(101, 200), (301, 400), (501, 600)][..]), + ("T1ri", &[(101, 400), (501, 600)][..]), + ( + "T1cas", + &[(101, 200), (241, 260), (301, 400), (501, 600)][..], + ), + ("T1alt", &[(101, 220), (301, 400), (501, 600)][..]), + ] { + for &(s, e) in exons { + gtf.push(make_exon("chr1", s, e, '+', "G1", tr)); + } + } + for (tr, exons) in [ + ("T2", &[(101, 200), (301, 400)][..]), + ("T2cas", &[(101, 200), (241, 260), (301, 400)][..]), + ("T2alt", &[(101, 220), (301, 400)][..]), + ] { + for &(s, e) in exons { + gtf.push(make_exon("chr2", s, e, '+', "G2", tr)); + } + } + TranscriptomeIndex::from_gtf_exons(>f, &genome).unwrap() + } + + #[test] + fn unspliced_mode_from_str() { + use std::str::FromStr; + for (v, m) in [ + ("None", QuantTranscriptomeUnspliced::None), + ("Intron", QuantTranscriptomeUnspliced::Intron), + ("PreMRNA", QuantTranscriptomeUnspliced::PreMRNA), + ] { + assert_eq!(QuantTranscriptomeUnspliced::from_str(v).unwrap(), m); + } + assert!(QuantTranscriptomeUnspliced::from_str("intron").is_err()); + let p = default_params(); + assert_eq!( + p.quant_transcriptome_unspliced, + QuantTranscriptomeUnspliced::None + ); + assert_eq!(p.quant_transcriptome_unspliced_flank, -1); + } + + #[test] + fn unspliced_intervals_intron_and_premrna() { + let idx = retained_intron_index(); + let g1 = idx.gene_ids.iter().position(|g| g == "G1").unwrap(); + let g2 = idx.gene_ids.iter().position(|g| g == "G2").unwrap(); + let iv = |mode, flank| idx.unspliced_intervals(mode, flank); + let intron0 = iv(QuantTranscriptomeUnspliced::Intron, 0); + // Union of every isoform's introns (the retained intron included). + assert_eq!(intron0[g1], Some((0, 1, vec![(200, 300), (400, 500)]))); + assert_eq!(intron0[g2], Some((1, 1, vec![(1200, 1300)]))); + let intron10 = iv(QuantTranscriptomeUnspliced::Intron, 10); + assert_eq!( + intron10[g1].as_ref().unwrap().2, + vec![(190, 310), (390, 510)] + ); + // Flanks are clipped to the gene body and merged when they meet. + let intron150 = iv(QuantTranscriptomeUnspliced::Intron, 150); + assert_eq!(intron150[g1].as_ref().unwrap().2, vec![(100, 600)]); + let pre = iv(QuantTranscriptomeUnspliced::PreMRNA, 0); + assert_eq!(pre[g1].as_ref().unwrap().2, vec![(100, 600)]); + assert_eq!(pre[g2].as_ref().unwrap().2, vec![(1100, 1400)]); + } + + #[test] + fn single_exon_gene_has_no_intron_target() { + let genome = make_genome(); + let gtf = vec![make_exon("chr1", 101, 300, '+', "G1", "T1")]; + let idx = TranscriptomeIndex::from_gtf_exons(>f, &genome).unwrap(); + let ext = idx.with_unspliced_targets(QuantTranscriptomeUnspliced::Intron, 10); + assert_eq!(ext.n_transcripts(), 1); + let ext = idx.with_unspliced_targets(QuantTranscriptomeUnspliced::PreMRNA, 10); + assert_eq!(ext.tr_ids, ["T1", "G1-I"]); + } + + #[test] + fn unspliced_targets_are_appended_in_gene_order() { + let idx = retained_intron_index(); + let ext = idx.with_unspliced_targets(QuantTranscriptomeUnspliced::Intron, 10); + let n = idx.n_transcripts(); + assert_eq!(ext.tr_ids[..n], idx.tr_ids[..]); + assert_eq!(ext.tr_ids[n..], ["G1-I", "G2-I"]); + assert_eq!(ext.unspliced.unwrap().first, n); + assert!(!ext.is_unspliced_target(n - 1) && ext.is_unspliced_target(n)); + assert_eq!(ext.tr_length[n], 240); + assert!(ext.tr_starts_sorted.windows(2).all(|w| w[0] <= w[1])); + // Annotated transcripts keep their relative order in the sorted view. + let annotated: Vec = ext.tr_order.iter().copied().filter(|&t| t < n).collect(); + assert_eq!(annotated, idx.tr_order); + } + + fn project_se_unspliced( + idx: &TranscriptomeIndex, + chr: usize, + start: u64, + end: u64, + spliced: bool, + ) -> Vec { + use cigar::op::{Kind, Op}; + let genome = make_genome(); + let len = (end - start) as usize; + let align = make_align( + chr, + false, + vec![(start, end, 0, len)], + vec![Op::new(Kind::Match, len)], + ); + let mut names: Vec = filter_and_project( + &align, + &vec![0u8; len], + &genome, + idx, + len as u32, + QuantTranscriptomeSAMoutput::BanSingleEnd, + &default_params(), + spliced, + ) + .iter() + .map(|p| idx.tr_ids[p.chr_idx].clone()) + .collect(); + names.sort(); + names + } + + #[test] + fn unspliced_target_projection() { + let idx = retained_intron_index(); + let ext = idx.with_unspliced_targets(QuantTranscriptomeUnspliced::Intron, 10); + // Inside the retained intron: retained-intron isoform AND G1-I. + assert_eq!(project_se_unspliced(&idx, 0, 270, 295, false), ["T1ri"]); + assert_eq!( + project_se_unspliced(&ext, 0, 270, 295, false), + ["G1-I", "T1ri"] + ); + // The fragment crosses a junction: spliced targets only. + assert_eq!(project_se_unspliced(&ext, 0, 270, 295, true), ["T1ri"]); + // Constitutive exon, away from the flanks: spliced targets only. + assert_eq!( + project_se_unspliced(&ext, 0, 520, 570, false), + ["T1", "T1alt", "T1cas", "T1ri"] + ); + // Intron 2 / exon 3 boundary: only an unspliced target can hold it, + // provided the flank covers the exonic part. + assert_eq!( + project_se_unspliced(&ext, 0, 480, 530, false), + Vec::::new() + ); + let ext49 = idx.with_unspliced_targets(QuantTranscriptomeUnspliced::Intron, 49); + assert_eq!(project_se_unspliced(&ext49, 0, 480, 530, false), ["G1-I"]); + // Cassette exon of G2: its inclusion isoform and G2-I. + assert_eq!( + project_se_unspliced(&ext, 1, 1241, 1259, false), + ["G2-I", "T2cas"] + ); + } + + #[test] + fn unspliced_target_on_minus_strand_is_flipped() { + use cigar::op::{Kind, Op}; + let genome = make_genome(); + let gtf = vec![ + make_exon("chr1", 101, 200, '-', "Gm", "Tm"), + make_exon("chr1", 301, 400, '-', "Gm", "Tm"), + ]; + let idx = TranscriptomeIndex::from_gtf_exons(>f, &genome) + .unwrap() + .with_unspliced_targets(QuantTranscriptomeUnspliced::Intron, 0); + let align = make_align( + 0, + false, + vec![(210, 260, 0, 50)], + vec![Op::new(Kind::Match, 50)], + ); + let p = align_to_transcripts(&align, &idx, 50); + assert_eq!(p.len(), 1); + assert_eq!(idx.tr_ids[p[0].chr_idx], "Gm-I"); + // Target [200, 300) of length 100 read 3' to 5': offset 100 - 60 = 40. + assert_eq!(p[0].genome_start, 40); + assert!(p[0].is_reverse); + } + + #[test] + fn targets_table_and_unspliced_fasta() { + let genome = make_genome(); + let gtf = vec![ + make_exon("chr1", 101, 200, '-', "Gm", "Tm"), + make_exon("chr1", 301, 400, '-', "Gm", "Tm"), + ]; + let mut ann = gtf.clone(); + for r in &mut ann { + r.attributes.insert("gene_name".into(), "GENEM".into()); + } + let idx = TranscriptomeIndex::from_gtf_exons(&ann, &genome) + .unwrap() + .with_unspliced_targets(QuantTranscriptomeUnspliced::Intron, 0); + let dir = tempfile::tempdir().unwrap(); + let tsv = dir.path().join("t.tsv"); + idx.write_targets_tsv(&tsv).unwrap(); + assert_eq!( + std::fs::read_to_string(&tsv).unwrap(), + "target_id\tgene_id\tgene_name\tstatus\tlength\n\ + Tm\tGm\tGENEM\tspliced\t200\n\ + Gm-I\tGm\tGENEM\tunspliced\t100\n" + ); + let fa = dir.path().join("u.fa"); + idx.write_unspliced_fasta(&fa, &genome).unwrap(); + let fa = std::fs::read_to_string(&fa).unwrap(); + let lines: Vec<&str> = fa.lines().collect(); + assert_eq!(lines[0], ">Gm-I"); + // The test genome is all A; the - strand target reads as T. + assert_eq!(lines[1..].concat(), "T".repeat(100)); + assert_eq!(lines[1].len(), 60); + } } diff --git a/src/solo/transcript3p.rs b/src/solo/transcript3p.rs index 5ba885da..0cb8957a 100644 --- a/src/solo/transcript3p.rs +++ b/src/solo/transcript3p.rs @@ -595,6 +595,7 @@ mod tests { tr_order: vec![0, 1, 2], tr_starts_sorted: vec![0; 3], tr_end_max_sorted: vec![1000; 3], + unspliced: None, } } } diff --git a/tests/bulk_unspliced.rs b/tests/bulk_unspliced.rs new file mode 100644 index 00000000..2eccc280 --- /dev/null +++ b/tests/bulk_unspliced.rs @@ -0,0 +1,622 @@ +//! Integration tests for the bulk total-RNA options (rustar-aligner +//! extensions, not in STAR): +//! +//! - `--quantGeneSplicing Yes` (per-gene spliced / unspliced / ambiguous); +//! - `--quantTranscriptomeUnspliced Intron|PreMRNA` (`-I` targets +//! in TranscriptomeSAM). +//! +//! A synthetic genome carries three genes, each with a fully spliced isoform +//! `A` and a retained-intron isoform `R` (exon 1 + intron 1 + exon 2 as one +//! exon). Reads are simulated from a known mixture of mature `A`, mature `R` +//! and unspliced pre-mRNA, as a reverse-stranded (dUTP) single-end library. + +use assert_cmd::cargo::cargo_bin_cmd; +use noodles::bam; +use std::collections::HashMap; +use std::fs; +use std::io::Write; +use std::path::{Path, PathBuf}; +use tempfile::TempDir; + +const READ_LEN: usize = 50; +const CHR_LEN: usize = 40_000; + +/// Gene layout relative to the gene start: exons `[0,300) [800,1000) +/// [1600,1900)`, introns `[300,800) [1000,1600)`. +const EXONS: [(usize, usize); 3] = [(0, 300), (800, 1000), (1600, 1900)]; +const GENE_LEN: usize = 1900; + +struct Gene { + name: &'static str, + start: usize, + reverse: bool, + /// Reads simulated from mature `A`, mature `R`, and pre-mRNA. + n_a: usize, + n_r: usize, + n_pre: usize, +} + +const GENES: [Gene; 3] = [ + Gene { + name: "G1", + start: 5_000, + reverse: false, + n_a: 300, + n_r: 0, + n_pre: 400, + }, + Gene { + name: "G2", + start: 15_000, + reverse: false, + n_a: 300, + n_r: 150, + n_pre: 400, + }, + Gene { + name: "G3", + start: 25_000, + reverse: true, + n_a: 300, + n_r: 0, + n_pre: 400, + }, +]; + +struct Lcg(u32); +impl Lcg { + fn next(&mut self) -> u32 { + self.0 = self.0.wrapping_mul(1_103_515_245).wrapping_add(12_345); + self.0 >> 16 + } +} + +fn rc(seq: &[u8]) -> Vec { + seq.iter() + .rev() + .map(|&b| match b { + b'A' => b'T', + b'T' => b'A', + b'C' => b'G', + b'G' => b'C', + _ => b, + }) + .collect() +} + +fn build_genome() -> Vec { + let mut rng = Lcg(424_242); + let mut g: Vec = (0..CHR_LEN) + .map(|_| b"ACGT"[(rng.next() & 3) as usize]) + .collect(); + // Plant canonical motifs at every intron: GT..AG on the transcribed + // strand, i.e. CT..AC in forward genome coordinates for a - gene. + for gene in &GENES { + for w in EXONS.windows(2) { + let (d, a) = (gene.start + w[0].1, gene.start + w[1].0); + let (left, right) = if gene.reverse { + (b"CT", b"AC") + } else { + (b"GT", b"AG") + }; + g[d..d + 2].copy_from_slice(left); + g[a - 2..a].copy_from_slice(right); + } + } + g +} + +/// Returns the FASTA, GTF and FASTQ paths, plus the gene-relative start of +/// every simulated pre-mRNA read. +fn write_inputs(dir: &Path, genome: &[u8]) -> (PathBuf, PathBuf, PathBuf, Vec) { + let fasta = dir.join("genome.fa"); + let mut f = fs::File::create(&fasta).unwrap(); + writeln!(f, ">chr1").unwrap(); + f.write_all(genome).unwrap(); + writeln!(f).unwrap(); + + let gtf = dir.join("ann.gtf"); + let mut f = fs::File::create(>f).unwrap(); + for gene in &GENES { + let strand = if gene.reverse { '-' } else { '+' }; + let a_exons = EXONS.to_vec(); + let r_exons = vec![(EXONS[0].0, EXONS[1].1), EXONS[2]]; + for (tr, exons) in [("A", a_exons), ("R", r_exons)] { + for (s, e) in exons { + writeln!( + f, + "chr1\tsim\texon\t{}\t{}\t.\t{strand}\t.\tgene_id \"{g}\"; transcript_id \"{g}_{tr}\";", + gene.start + s + 1, + gene.start + e, + g = gene.name + ) + .unwrap(); + } + } + } + + // Reads: fragments of the molecule in genome-forward orientation; a + // reverse-stranded read 1 is the reverse complement of the RNA, i.e. the + // reverse complement of the fragment for a + gene and the fragment itself + // for a - gene. + let fq = dir.join("reads.fq"); + let mut f = fs::File::create(&fq).unwrap(); + let mut rng = Lcg(7); + let mut pre_starts = Vec::new(); + for gene in &GENES { + let body = [(0, GENE_LEN)]; + let r_exons = [(EXONS[0].0, EXONS[1].1), EXONS[2]]; + for (kind, blocks, n) in [ + ("A", &EXONS[..], gene.n_a), + ("R", &r_exons[..], gene.n_r), + ("P", &body[..], gene.n_pre), + ] { + let seq: Vec = blocks + .iter() + .flat_map(|&(s, e)| genome[gene.start + s..gene.start + e].iter().copied()) + .collect(); + for i in 0..n { + let pos = rng.next() as usize % (seq.len() - READ_LEN + 1); + if kind == "P" { + pre_starts.push(pos); + } + let frag = &seq[pos..pos + READ_LEN]; + let read = if gene.reverse { + frag.to_vec() + } else { + rc(frag) + }; + writeln!(f, "@{}_{kind}_{i}", gene.name).unwrap(); + f.write_all(&read).unwrap(); + writeln!(f, "\n+\n{}", "I".repeat(READ_LEN)).unwrap(); + } + } + } + (fasta, gtf, fq, pre_starts) +} + +/// Build the index once per test and return (tmpdir, genome dir, gtf, +/// fastq, pre-mRNA read starts). +fn setup() -> (TempDir, PathBuf, PathBuf, PathBuf, Vec) { + let tmp = TempDir::new().unwrap(); + let genome = build_genome(); + let (fasta, gtf, fq, pre_starts) = write_inputs(tmp.path(), &genome); + let gdir = tmp.path().join("idx"); + fs::create_dir_all(&gdir).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .args(["--runMode", "genomeGenerate", "--genomeDir"]) + .arg(&gdir) + .arg("--genomeFastaFiles") + .arg(&fasta) + .args(["--genomeSAindexNbases", "8", "--sjdbOverhang", "49"]) + .arg("--sjdbGTFfile") + .arg(>f) + .arg("--outFileNamePrefix") + .arg(gdir.join("gen_")) + .assert() + .success(); + (tmp, gdir, gtf, fq, pre_starts) +} + +fn align(tmp: &TempDir, gdir: &Path, gtf: &Path, fq: &Path, name: &str, extra: &[&str]) -> PathBuf { + let out = tmp.path().join(name); + fs::create_dir_all(&out).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .arg("--genomeDir") + .arg(gdir) + .arg("--readFilesIn") + .arg(fq) + .arg("--sjdbGTFfile") + .arg(gtf) + .arg("--outFileNamePrefix") + .arg(format!("{}/", out.display())) + .args(extra) + .assert() + .success(); + out +} + +/// SAM body plus header lines that do not embed the command line. +fn sam_without_command_line(path: &Path) -> String { + fs::read_to_string(path) + .unwrap() + .lines() + .filter(|l| !l.starts_with("@PG") && !l.starts_with("@CO")) + .collect::>() + .join("\n") +} + +/// Every BAM record, in file order (the header carries the command line). +fn bam_records(path: &Path) -> Vec { + let mut reader = bam::io::Reader::new(fs::File::open(path).unwrap()); + reader.read_header().unwrap(); + reader + .records() + .map(|r| format!("{:?}", r.unwrap())) + .collect() +} + +/// `Log.final.out` minus the wall-clock lines. +fn log_final_without_times(path: &Path) -> String { + fs::read_to_string(path) + .unwrap() + .lines() + .filter(|l| { + !(l.contains("Started") || l.contains("Finished") || l.contains("Mapping speed")) + }) + .collect::>() + .join("\n") +} + +fn read(path: &Path) -> String { + fs::read_to_string(path).unwrap() +} + +/// Target names (in `@SQ` order) and, per read, the targets it was written +/// to in the transcriptome BAM. +/// `(target name, length)` in `@SQ` order, and read name -> target indices. +type TargetHits = (Vec<(String, usize)>, HashMap>); + +fn read_targets(path: &Path) -> TargetHits { + let mut reader = bam::io::Reader::new(fs::File::open(path).unwrap()); + let header = reader.read_header().unwrap(); + let refs: Vec<(String, usize)> = header + .reference_sequences() + .iter() + .map(|(k, r)| { + ( + String::from_utf8_lossy(k).to_string(), + usize::from(r.length()), + ) + }) + .collect(); + let mut per_read: HashMap> = HashMap::new(); + for rec in reader.records() { + let rec = rec.unwrap(); + let q = String::from_utf8_lossy(rec.name().unwrap()).to_string(); + let tid = rec.reference_sequence_id().unwrap().unwrap(); + per_read.entry(q).or_default().push(tid); + } + (refs, per_read) +} + +/// Reads per target estimated from the transcriptome BAM by a minimal +/// Salmon-like EM: a read compatible with targets `S` is shared in +/// proportion to `alpha_t / efflen_t`, `efflen_t = len_t - READ_LEN + 1`. +fn target_counts(path: &Path) -> HashMap { + let (refs, per_read) = read_targets(path); + let efflen: Vec = refs + .iter() + .map(|(_, l)| (l.saturating_sub(READ_LEN) + 1) as f64) + .collect(); + let classes: Vec> = per_read.into_values().collect(); + let n = refs.len(); + let mut alpha = vec![1.0 / n as f64; n]; + let mut counts = vec![0.0; n]; + for _ in 0..2000 { + counts.iter_mut().for_each(|c| *c = 0.0); + for tids in &classes { + let z: f64 = tids.iter().map(|&t| alpha[t] / efflen[t]).sum(); + if z > 0.0 { + for &t in tids { + counts[t] += alpha[t] / efflen[t] / z; + } + } + } + let total: f64 = counts.iter().sum(); + for (a, c) in alpha.iter_mut().zip(&counts) { + *a = c / total; + } + } + refs.into_iter().map(|(n, _)| n).zip(counts).collect() +} + +fn summary(path: &Path) -> HashMap> { + read(path) + .lines() + .skip(1) + .map(|l| { + let mut it = l.split('\t').map(str::to_string); + (it.next().unwrap(), it.collect()) + }) + .collect() +} + +#[test] +fn new_options_leave_existing_outputs_unchanged() { + let (tmp, gdir, gtf, fq, _) = setup(); + let base = align( + &tmp, + &gdir, + >f, + &fq, + "base", + &["--quantMode", "GeneCounts", "TranscriptomeSAM"], + ); + let splicing = align( + &tmp, + &gdir, + >f, + &fq, + "splicing", + &[ + "--quantMode", + "GeneCounts", + "TranscriptomeSAM", + "--quantGeneSplicing", + "Yes", + "--quantTranscriptomeUnspliced", + "None", + ], + ); + let unspliced = align( + &tmp, + &gdir, + >f, + &fq, + "unspliced", + &[ + "--quantMode", + "GeneCounts", + "TranscriptomeSAM", + "--quantTranscriptomeUnspliced", + "Intron", + ], + ); + + let base_sam = sam_without_command_line(&base.join("Aligned.out.sam")); + assert!(base_sam.lines().count() > 2000, "too few alignments"); + for other in [&splicing, &unspliced] { + assert_eq!( + base_sam, + sam_without_command_line(&other.join("Aligned.out.sam")) + ); + for f in ["ReadsPerGene.out.tab", "SJ.out.tab"] { + assert_eq!(read(&base.join(f)), read(&other.join(f)), "{f} differs"); + } + assert_eq!( + log_final_without_times(&base.join("Log.final.out")), + log_final_without_times(&other.join("Log.final.out")) + ); + } + // --quantGeneSplicing Yes and --quantTranscriptomeUnspliced None leave the + // transcriptome BAM untouched. + let base_bam = base.join("Aligned.toTranscriptome.out.bam"); + assert_eq!( + bam_records(&base_bam), + bam_records(&splicing.join("Aligned.toTranscriptome.out.bam")) + ); + // Intron adds the -I targets after the annotated transcripts. + let (base_refs, _) = read_targets(&base_bam); + let (refs, _) = read_targets(&unspliced.join("Aligned.toTranscriptome.out.bam")); + let names: Vec<&str> = refs.iter().map(|(n, _)| n.as_str()).collect(); + assert_eq!(refs[..base_refs.len()], base_refs[..]); + assert_eq!(names[base_refs.len()..], ["G1-I", "G2-I", "G3-I"]); + // New output files exist only when requested. + assert!(splicing.join("ReadsPerGeneSplicing.out.tab").exists()); + assert!(!base.join("ReadsPerGeneSplicing.out.tab").exists()); + assert!( + unspliced + .join("Aligned.toTranscriptome.targets.tsv") + .exists() + ); + assert!(!base.join("Aligned.toTranscriptome.targets.tsv").exists()); + assert!( + !unspliced + .join("Aligned.toTranscriptome.unspliced.fa") + .exists() + ); +} + +#[test] +fn simulated_total_rna_mixture() { + let (tmp, gdir, gtf, fq, pre_starts) = setup(); + let star = align( + &tmp, + &gdir, + >f, + &fq, + "star", + &[ + "--quantMode", + "TranscriptomeSAM", + "--quantGeneSplicing", + "Yes", + ], + ); + let intron = align( + &tmp, + &gdir, + >f, + &fq, + "intron", + &[ + "--quantMode", + "TranscriptomeSAM", + "--quantTranscriptomeUnspliced", + "Intron", + ], + ); + let premrna = align( + &tmp, + &gdir, + >f, + &fq, + "premrna", + &[ + "--quantMode", + "TranscriptomeSAM", + "--quantTranscriptomeUnspliced", + "PreMRNA", + ], + ); + + // --- Every simulated pre-mRNA read lying inside intron 1 (the intron + // that isoform R retains) is written to the gene's unspliced target, and + // no pre-mRNA read is written to a spliced target unless an isoform + // contains it. + let (refs, per_read) = read_targets(&intron.join("Aligned.toTranscriptome.out.bam")); + let (i1s, i1e) = (EXONS[0].1, EXONS[1].0); + let mut idx = 0; + for gene in &GENES { + for i in 0..gene.n_pre { + let p = pre_starts[idx]; + idx += 1; + if p >= i1s && p + READ_LEN <= i1e { + let targets = &per_read[&format!("{}_P_{i}", gene.name)]; + let names: Vec<&str> = targets.iter().map(|&t| refs[t].0.as_str()).collect(); + assert!( + names.contains(&format!("{}-I", gene.name).as_str()), + "{names:?}" + ); + assert!( + names.contains(&format!("{}_R", gene.name).as_str()), + "{names:?}" + ); + } + } + } + + // --- EM estimates: share of the mature reads given to the retained- + // intron isoform, and reads given to the spliced isoform A. + let counts: Vec> = [&star, &intron, &premrna] + .iter() + .map(|d| target_counts(&d.join("Aligned.toTranscriptome.out.bam"))) + .collect(); + let get = |c: &HashMap, t: String| c.get(&t).copied().unwrap_or(0.0); + let mut share_err = [0.0f64; 3]; + let mut a_err = [0.0f64; 3]; + for gene in &GENES { + let truth = gene.n_r as f64 / (gene.n_a + gene.n_r) as f64; + let mut line = format!("{}: RI share truth {truth:.3}", gene.name); + for (m, c) in counts.iter().enumerate() { + let a = get(c, format!("{}_A", gene.name)); + let r = get(c, format!("{}_R", gene.name)); + let u = get(c, format!("{}-I", gene.name)); + share_err[m] += (r / (a + r) - truth).abs(); + a_err[m] += (a - gene.n_a as f64).abs() / gene.n_a as f64; + line.push_str(&format!( + " | {} A {a:.0} R {r:.0} -I {u:.0} share {:.3}", + ["STAR", "Intron", "PreMRNA"][m], + r / (a + r) + )); + } + eprintln!( + "{line} (simulated A {}, R {}, pre-mRNA {})", + gene.n_a, gene.n_r, gene.n_pre + ); + } + eprintln!( + "summed |RI share error|: STAR {:.3}, Intron {:.3}, PreMRNA {:.3}", + share_err[0], share_err[1], share_err[2] + ); + eprintln!( + "summed relative error on A: STAR {:.3}, Intron {:.3}, PreMRNA {:.3}", + a_err[0], a_err[1], a_err[2] + ); + // Both modes stop the retained-intron isoform from absorbing pre-mRNA. + for m in [1, 2] { + assert!( + share_err[m] < share_err[0] / 2.0, + "RI share error {share_err:?}" + ); + } + // PreMRNA also improves the spliced isoform and recovers the pre-mRNA + // reads on -I. Intron (splici-style) cannot: the exonic part of + // pre-mRNA has no unspliced target to go to, so the spliced isoform + // over-counts it, as it does with splici. + assert!(a_err[2] < a_err[0], "A error {a_err:?}"); + for gene in &GENES { + let u = get(&counts[2], format!("{}-I", gene.name)); + assert!( + (u - gene.n_pre as f64).abs() < 0.15 * gene.n_pre as f64, + "{}: -I {u}", + gene.name + ); + } + + // --- GeneSplicing: the library is reverse-stranded, so the forward + // columns see nothing and the reverse columns see every gene read. + let s = summary(&star.join("ReadsPerGeneSplicing.summary.tsv")); + let n = |k: &str, col: usize| s[k][col].parse::().unwrap(); + let total_reads: usize = GENES.iter().map(|g| g.n_a + g.n_r + g.n_pre).sum(); + assert_eq!( + n("N_spliced", 1) + n("N_unspliced", 1) + n("N_ambiguous", 1), + 0 + ); + let assigned_rev = n("N_spliced", 2) + n("N_unspliced", 2) + n("N_ambiguous", 2); + assert!(assigned_rev as f64 > 0.95 * total_reads as f64); + // A read is unspliced when some model needs pre-mRNA and none is only + // exonic: here, exactly the pre-mRNA reads reaching more than 6 bases + // (STAR's minOverlapMinusOne) into intron 2, which no isoform retains. + let (i2s, i2e) = (EXONS[1].1, EXONS[2].0); + let expected_unspliced = pre_starts + .iter() + .filter(|&&p| (p + READ_LEN).min(i2e).saturating_sub(p.max(i2s)) > 6) + .count() as u64; + let got = n("N_unspliced", 2); + eprintln!("unspliced: expected {expected_unspliced}, got {got}"); + assert!(got.abs_diff(expected_unspliced) * 50 <= expected_unspliced); + assert!(n("N_ambiguous", 2) > 0); + let mature: usize = GENES.iter().map(|g| g.n_a).sum(); + assert!(n("N_spliced", 2) as usize >= mature * 95 / 100); +} + +#[test] +fn splicing_status_tag_in_genomic_sam() { + let (tmp, gdir, gtf, fq, pre_starts) = setup(); + let out = align( + &tmp, + &gdir, + >f, + &fq, + "sp", + &["--outSAMsplicingStatus", "Yes"], + ); + let base = align(&tmp, &gdir, >f, &fq, "nosp", &[]); + let sam = read(&out.join("Aligned.out.sam")); + let mut tag: HashMap = HashMap::new(); + let mut stripped = Vec::new(); + for line in sam.lines().filter(|l| !l.starts_with('@')) { + let f: Vec<&str> = line.split('\t').collect(); + if let Some(t) = f.iter().find(|x| x.starts_with("sp:A:")) { + tag.insert(f[0].to_string(), t[5..].to_string()); + } + stripped.push( + f.iter() + .filter(|x| !x.starts_with("sp:A:")) + .copied() + .collect::>() + .join("\t"), + ); + } + // Minus the tag, the records are the default ones. + let base_sam = read(&base.join("Aligned.out.sam")); + let base_records: Vec<&str> = base_sam.lines().filter(|l| !l.starts_with('@')).collect(); + assert_eq!(stripped, base_records); + // Pre-mRNA reads: inside exon 1 -> S, well inside the constitutive + // intron 2 -> U, well inside the retained intron 1 -> A. + let (i1s, i1e) = (EXONS[0].1, EXONS[1].0); + let (i2s, i2e) = (EXONS[1].1, EXONS[2].0); + let mut idx = 0; + let mut checked = 0; + for gene in &GENES { + for i in 0..gene.n_pre { + let p = pre_starts[idx]; + idx += 1; + let name = format!("{}_P_{i}", gene.name); + let want = if p >= i2s + 7 && p + READ_LEN + 7 <= i2e { + "U" + } else if p >= i1s + 7 && p + READ_LEN + 7 <= i1e { + "A" + } else if p + READ_LEN <= i1s { + "S" + } else { + continue; + }; + assert_eq!(tag[&name], want, "{name} at {p}"); + checked += 1; + } + } + assert!(checked > 500, "{checked}"); +}