Implementationstooling reviewed 2026-08
13 RNA-seq Quantification: Implementations
Everything on this page is generated from tools/rnaseq.json. If something here is wrong or out of date, the fix belongs in that file, not in this prose.
For what a count is a count of, and why counts are not Gaussian, see Section 12.1.
13.1 Start here
A bulk RNA-seq analysis is two decisions, not one: how to get from reads to a count matrix, and how to test it. The cards below are the default answer to each, plus the spliced aligner from Section 7.3 for the route that goes through alignments.
Fits the negative binomial model described in Section 12.6, and its useful behavior is what it does when a gene has too few replicates to estimate its own dispersion: it shrinks toward a trend fitted across all genes. The same shrinkage applies to fold changes, which is why its effect sizes for low-count genes are usable rather than merely large. (Love et al. 2014)
Expects raw counts, not normalized values, because the normalization is part of the model. Feeding it TPMs or FPKMs produces numbers and invalidates the test.
Quantifies against the transcriptome without producing alignments, which is both much faster and a different question from the one an aligner answers. Its distinguishing feature is bias modeling — sequence, position and GC — estimated from the data rather than assumed away, which is what makes its estimates comparable across libraries prepared differently. (Patro et al. 2017)
It quantifies only what is in the index, so a transcript missing from your annotation is invisible rather than unmapped. Selective alignment against a decoy-aware index matters for the same reason: without decoys, reads from unannotated regions are forced onto whichever transcript fits best.
The default spliced aligner for short-read RNA-seq. Its manual recommends two-pass mode for the most sensitive novel-junction discovery: the first pass collects junctions, which are inserted into the index so the second pass can align across them. Upstream is careful about what this buys — not many more junctions, but more reads mapped to the ones found. (Dobin et al. 2013)
Memory-hungry: STAR's README asks for at least 16 GB of RAM for a mammalian genome and recommends 32 GB, which rules it out on small machines regardless of its merits.
13.2 Getting to a count matrix
Two routes, and the choice is about what else you need rather than about accuracy.
The alignment-free route indexes the transcriptome and quantifies directly. It is faster by a wide margin and it is bounded by your annotation — nothing outside it can be seen. The align-then-count route produces alignments you can inspect, call variants from, or use to find novel junctions, and then counts reads per gene, which discards isoform information entirely.
If you need both quantification and alignments, note that the alignment-based estimator below works from a BAM you already have, rather than requiring a third pass.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| kallisto | active | C++ |
The tool that established alignment-free quantification, using pseudoalignment against a transcriptome de Bruijn graph. Still fast, still maintained and still a defensible choice; it is listed below Salmon because it models fewer of the biases described in Section 12.5, not because of its speed or its accuracy on clean data. (Bray et al. 2016) |
| RSEM | maintained | C++ |
The alignment-based expectation-maximization approach the faster tools were built to approximate. It still produces excellent estimates and it works from an existing BAM, which matters when you need quantification and alignments from the same run rather than two passes over the data. (Li and Dewey 2011) |
| Rsubread | active | RC |
The counting half of the align-then-count route: given a BAM and an annotation, assign reads to features. This is the R packaging of the same engine distributed as the featureCounts command in Subread, and it is the one with a current, verifiable release channel. Choose this route when you need gene-level counts from alignments you already have. (Liao et al. 2014) |
| tximport | active | R |
The join between transcript-level quantification and gene-level testing, and it is listed because the naive version of that join is wrong. Summing transcript estimates to a gene ignores that different transcripts have different effective lengths, so the correct summary carries an offset that the downstream model then uses. Skipping this step biases gene-level results in a way nothing downstream will flag. (Soneson et al. 2016) |
WarningDo not skip the summarization step
Going from transcript-level estimates to gene-level counts is not addition. Transcripts of a gene have different effective lengths, so the correct summary carries an offset the downstream model uses. The tool for this is in the table above, and omitting it biases gene-level results silently — nothing downstream will warn you.
13.3 Testing for differential expression
The default is the card in Section 13.1 above; the alternatives are below. All three implement Section 12.6 competently and agree far more often than the discourse suggests, so choose on the shape of your experiment rather than on published comparisons. Two of them model counts directly and differ mainly in how conservative their small-sample test is. The third transforms counts into a linear model framework, which expresses complex designs — covariates, batch terms, interactions — more naturally than either.
They all want raw counts. Feeding a differential test length-normalized values produces output and invalidates it, for the reason in Section 12.5.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| edgeR | active | R |
The other negative binomial implementation, developed alongside DESeq2 rather than after it, and reaching very similar conclusions on the same data. Its quasi-likelihood test is the more conservative of the two on small experiments, which is a reason to prefer it rather than a reason to run both and pick. (Robinson et al. 2010) |
| limma | active | R |
Takes the other route of Section 12.6: rather than modeling counts directly, voom transforms them and estimates a precision weight per observation, which puts the data into the linear model framework limma was built for. That framework is the reason to choose it — complex designs with covariates, batch terms and interactions are expressed more naturally here than in either count model. (Ritchie et al. 2015; Law et al. 2014) |
NotePick before you look
Running two or three of these and reporting the one with the most hits is a multiple comparisons problem you have concealed from yourself rather than solved. Choose on the design, run one, and if you do compare tools, say so.
13.4 What is not here
Single-cell quantification is a different part. The unit is a cell rather than a sample, the counts are sparse in a way bulk methods do not model, and the questions asked of the matrix afterwards are largely not differential expression. Its tooling barely overlaps with this page.
Novel transcript assembly is also absent. Everything above takes an annotation as given, and discovering transcripts that are not in one is a distinct problem with distinct tools.
13.5 Language coverage: R, and why
This is the first part of the book where R is not a secondary option but the place the work actually happens. Quantification is command-line C++, as everything performance-bound in this book is. Everything after it — the count matrix, the model, the shrinkage, the design formula — is Bioconductor, and the Python equivalents are genuinely less complete for this specific task rather than merely less popular.
The reason is historical and worth knowing: the statistical machinery here grew out of microarray analysis, which was an R field a decade before RNA-seq existed, and the linear-model and empirical-Bayes tooling was already mature when counts arrived. Some of it was adapted rather than rebuilt.
NoteWhat this means in practice
A bulk RNA-seq analysis is usually a shell pipeline that stops at a count matrix and an R session that starts there. Planning for that boundary is easier than discovering it after you have written half the analysis somewhere else.