12 RNA-seq Quantification: Concepts

Concepts

12.1 The problem

RNA-seq measures how much of each transcript was present by sequencing fragments of it and counting. The counting is where the difficulty lives, because the thing you want to count is not the thing you observe.

You observe fragments. You want molecules. Between the two sit several transformations that do not preserve proportion: longer transcripts yield more fragments than short ones at the same abundance; library preparation and sequencing introduce composition-dependent biases; and a large share of fragments are consistent with more than one transcript. None of this is noise to be averaged away. It is systematic, and it has to be modeled.

A count is not a measurement of expression. It is a measurement of sequencing effort that landed on a feature, which depends on the feature’s abundance, its length, its sequence composition, and how much effort was spent on the library overall. Every method in this part exists to get from the second quantity back to the first.

12.2 What is the count a count of?

The first question is what you are counting, and there are two defensible answers that lead to different analyses.

Gene-level counting assigns each fragment to a gene and adds up. It is robust and it is what most differential expression is done on, because genes are what most biological questions are about.

Transcript-level quantification estimates the abundance of each isoform separately. It answers questions gene-level counting cannot — which isoform changed, whether a switch occurred at constant total expression — and it is much harder, because isoforms of a gene share most of their sequence.

Section 12.3 is why the second is hard, and it is also why the two are not related by simple addition.

12.3 Multi-mapping is the central problem

A fragment from a shared exon is consistent with every isoform containing that exon. A fragment from a conserved domain may be consistent with several paralogous genes. Discarding these fragments is not conservative — it systematically underestimates exactly the genes and isoforms that share sequence, which is a biased loss rather than a smaller sample.

The resolution used by every current method is to stop assigning fragments to transcripts at all, and instead ask which set of abundances would best explain the whole collection of fragments. A fragment compatible with three isoforms contributes to all three, weighted by how abundant each currently appears to be. Those weights change the abundance estimates, which change the weights, and the procedure is iterated to a fixed point.

This is expectation-maximization, and the important consequence is that a transcript-level estimate is not a count of anything. It is a parameter of a model fitted to all the fragments jointly. Estimates for isoforms that share most of their sequence are correlated and individually uncertain even when their sum is well determined — which is why gene-level conclusions are so much more stable than isoform-level ones, and why a change in a single isoform needs more evidence than its point estimate suggests.

12.4 Alignment-free quantification

If fragments are never assigned to positions, the alignment step is doing more work than the question requires. Knowing which transcripts a fragment is compatible with is enough to run the model above — the exact base-level correspondence is not used.

That observation is what alignment-free methods exploit. They index the transcriptome and determine compatibility directly, which is dramatically faster and uses far less memory than aligning to the genome and post-processing.

The trade is a real one and should be stated plainly. These methods quantify against a transcriptome, so they are bounded by the annotation: an unannotated transcript is not merely unquantified, it is invisible, and its fragments will be attributed to whichever annotated transcript they best match. They also produce no alignments, so they cannot support anything that needs them — variant calling, coverage inspection, novel junction discovery. For a well-annotated organism and a quantification question, that is an easy trade. For a poorly annotated one it may not be.

12.5 Normalization is not optional

Raw counts are not comparable, in two distinct ways that are often conflated.

Within a sample, a longer transcript produces more fragments at equal abundance, so counts cannot be compared between genes without a length correction. This is what TPM does.

Between samples, total sequencing depth differs, and — more subtly — so does composition. If a handful of genes are massively induced in one condition, they consume a larger share of the library, and every other gene’s count falls even at unchanged abundance. Scaling by total reads does not fix this; it propagates it.

The methods that work assume most genes are not differentially expressed and estimate a scaling factor from the bulk of the distribution rather than from the total. That assumption is usually reasonable and it is an assumption, and it fails in exactly the cases you might care about most — a global shift in transcription, or a comparison between very different tissues.

TPM is for looking, not for testing

Length-normalized units are the right thing for comparing genes within a sample and for eyeballing a value. They are the wrong input to a differential test, because normalizing away the count nature of the data destroys the information the test needs about how certain each measurement is. A gene with ten fragments and one with ten thousand can have the same TPM and very different reliability, and after the transform nothing downstream can tell.

12.6 Counts are not Gaussian

Differential expression asks whether a gene’s abundance differs between conditions by more than the variability within them. Doing that requires a model of the within-condition variability, and for counts the obvious choices are wrong.

Counts are non-negative integers, and their variance grows with their mean — a gene averaging 10 fragments does not vary the way one averaging 10,000 does. A Poisson model captures that but assumes variance equals mean, which holds for technical replicates and badly understates the variability between biological ones. The standard answer is the negative binomial, which adds a dispersion parameter for the extra biological variability.

That leaves the real difficulty: with three or four replicates you cannot estimate a dispersion per gene with any precision. The methods that work borrow strength across genes, fitting a trend of dispersion against expression and shrinking each gene’s estimate toward it. Shrinkage is what makes small RNA-seq experiments analysable at all, and it is worth understanding as an assumption rather than a detail: it presumes genes at similar expression behave similarly, which is usually true and is doing a great deal of work.

What the test is testing

A differential expression result is a statement about a gene’s mean, between the groups you defined, under a model of biological variability estimated with help from every other gene. It is not a statement that the gene matters, that the change is large, or that it would replicate in another cohort. With enough replicates, statistically significant changes of no biological consequence become easy to find, which is why an effect-size threshold alongside the adjusted p-value is a methodological choice worth making deliberately.

12.7 What to take forward

A count measures sequencing effort landing on a feature, not abundance; length, composition and library depth all sit in between.
Fragments compatible with several transcripts are the central problem, and discarding them biases against exactly the genes that share sequence.
Transcript-level estimates are model parameters, not counts, and isoform estimates are correlated and individually uncertain.
Alignment-free methods are bounded by the annotation and produce no alignments; for a well-annotated organism that is usually a good trade.
Normalization must handle length within a sample and composition between samples, and the composition methods assume most genes do not change.
Counts are over-dispersed, replicates are few, and shrinkage across genes is what makes the resulting tests usable.

Section 13.1 names what to run.