Finding the genes that changed between two conditions takes one statistical test per gene. The MA and volcano plots show the results, and fold-change and q-value thresholds pick the candidates.
ChIP-seq showed where p53 binds near CDKN1A, but binding is not the same as transcription. Five Views predicted that CDKN1A expression should fall when p53 is broken. Testing that prediction, and finding the other genes p53 keeps switched on, starts with a differential expression experiment: RNA sequencing of a few samples per condition, a read count for each of about 20,000 protein-coding genes, and one statistical test per gene.
The figures use one simulated experiment: cells with working p53 (wildtype) against the same cells with TP53 knocked out, three samples of each to start, and 2,000 genes. Each gene gets a baseline expression level, spread over several orders of magnitude as in real tissue. 116 genes carry real names: the direct p53 targets in a census of such genes (Fischer, 2017), each planted with a fall in the knockout whose size we chose. Ten of them, CDKN1A, MDM2, BAX and seven others, fall clearly; the other 106 fall by only 0.2 to 0.6 log₂ units. Of the remaining genes, 10% are given a random change up or down. Read counts are drawn from a negative binomial distribution, the standard model for RNA-seq counts, with variance \(\mu + 0.08\,\mu^2\) for mean \(\mu\). Each gene is tested with limma's moderated t-test on log counts, and the p-values are adjusted with the Benjamini-Hochberg procedure. The code records each gene's planted change, so every call can be graded.
The test gives each gene an estimated log fold change and an adjusted p-value, or q-value. Two standard plots arrange that table.
The MA plot puts each gene's average abundance on the horizontal axis (A, the mean log count) and its log fold change on the vertical axis (M, the difference in mean log counts, knockout minus wildtype):
\[ \begin{gathered} M = \overline{\log_2 \text{knockout}} - \overline{\log_2 \text{wildtype}} \\ A = \tfrac{1}{2}\left(\overline{\log_2 \text{knockout}} + \overline{\log_2 \text{wildtype}}\right) \end{gathered} \]
The volcano plot puts the same log fold change on the horizontal axis and the strength of evidence, \(-\log_{10} q\), on the vertical. Hover or tap a gene in either plot to find it in the other, or tab to a plot and step through the genes from the smallest q with the arrow keys.
Figure 1. The MA plot and the volcano plot of the same 2,000 simulated genes. Genes that pass both thresholds are colored by direction: higher or lower in the knockout; gray genes fail at least one. Dashed lines mark the thresholds. CDKN1A is ringed and labeled; the other nine named p53 targets have a thin ring. Changing the number of replicates measures the same genes with more or fewer samples; New experiment draws new genes.
The MA plot shows a funnel: the spread of fold changes is widest at the left, among genes with few reads. A gene seen 3 times in one sample and 9 times in another has an apparent threefold change that means little, while the same ratio between 3,000 and 9,000 reads is hard to get by chance. So the volcano's arms rise only where a large change meets enough reads.
CDKN1A is the ringed gene on the left arm: a fall of 2.38 log₂ units (planted 2.6) with q = 3.5e-5. That is the prediction from Five Views, found by a test that did not know where to look.
More replicates narrow the funnel and lift the arms, since each group mean is estimated from more samples. With three samples per group, a gene's own variance estimate rests on four degrees of freedom and is itself very noisy. A plain t-test on these 2,000 genes, followed by the same correction, found 4.6 genes on average over 200 simulated experiments, and none at all in 50 of them. limma's moderated t-test (Smyth, 2004) borrows strength across genes. It fits a prior variance \(s_0^2(A)\) that follows abundance, the trend of limma-trend (Law et al., 2014), and estimates from the spread of all the genes' variances how much that prior should count, as \(d_0\) extra degrees of freedom. Each gene's variance \(s^2\), on \(d = 2n - 2\) degrees of freedom, is replaced by
\[ \tilde s^2 = \frac{d_0\, s_0^2(A) + d\, s^2}{d_0 + d}, \]
and t is computed with \(\tilde s^2\) on \(d_0 + d\) degrees of freedom. In the experiment shown, \(d_0\) = 18.9, so the prior counts for about five times as much as the gene's own four degrees of freedom. The estimate is noisy in its own right: over 200 experiments with three replicates its median was 19.6, and in 12 of them it was infinite, which gives every gene the trend's variance. The JavaScript here gives the same \(d_0\) and p-values as limma 3.68.5 run on the same counts, to at least ten significant digits. DESeq2 and edgeR model the counts directly as negative binomial and shrink each gene's dispersion toward a trend fitted against mean count (Love et al., 2014). DESeq2 can also shrink low-count fold changes toward zero in a separate step (lfcShrink), which narrows the funnel itself, and by default drops genes too faint to reach significance before correcting (independent filtering, Bourgon et al., 2010).
A candidate list usually requires both a minimum fold change and a maximum q-value. Figure 2 sweeps the fold-change cutoff, with and without the q filter, and grades each list against the planted changes, at three replicates and on the same genes measured with eight.
Figure 2. The simulated genes of Figure 1 graded against their planted changes as the |log₂ FC| cutoff rises, at q ≤ 0.05 (Figure 1's q). Solid curves use both filters, dashed curves the fold change alone; the vertical line is Figure 1's cutoff. Hover or tap to read a cutoff; click or use the arrow keys to move Figure 1's cutoff.
The fold-change filter alone admits many unchanged genes: 59 of the 209 on its list at three replicates, and 25.5% of the list on average over 200 experiments. Nearly all of them are low-count genes whose noisy estimates happened to land far from zero. The q filter keeps those out. At three replicates it is also the only filter doing any work: every gene with q ≤ 0.05 already has |log₂ FC| ≥ 1.08, since with so few samples only a large change reaches significance, so the solid curves stay flat until the cutoff passes that value.
At eight replicates the two filters part. The q filter finds 186 of the 221 changed genes, including real changes smaller than twofold, and requiring |log₂ FC| ≥ 1 as well cuts the list to 140. The solid power curve falls from the start while the unchanged share barely moves: the cut is a choice about effect size, not about error. A list of genes whose change exceeds a threshold, with the error still controlled, comes from testing against the threshold directly (TREAT, McCarthy and Smyth, 2009) rather than filtering afterwards.
With 2,000 tests, a p-value cutoff of 0.05 would let through about 84 unchanged genes by chance alone. Benjamini and Hochberg's procedure works on the ranked p-values instead. Sort the m p-values from smallest to largest, \(p_{(1)} \le \dots \le p_{(m)}\), find the largest rank k with
\[ p_{(k)} \le \frac{k}{m}\, q, \]
and call the k genes with the smallest p-values. Because it takes the largest such k, a gene can sit above the line and still be called, as long as a gene further down the ranking falls under it. Benjamini and Hochberg (1995) proved that for independent tests the false discovery rate, the expected unchanged share of the list, is at most q; their proof gives the tighter \((m_0/m)\,q\), where \(m_0\) is the number of unchanged genes. The histogram in Figure 3 shows where \(m_0\) can be seen: p-values of unchanged genes are spread evenly between 0 and 1, so they form a flat floor, while changed genes pile up near zero. Storey's estimate reads \(\pi_0 = m_0/m\) off the right half of the histogram, where nearly every gene is unchanged (Storey, 2002): here 0.900, against a true 0.890. It tends to run a little high, on the safe side, because small true changes leave some p-values in the right half too. With 84% of the genes unchanged, the bound at q = 0.05 is 4.2%.
Figure 3. Benjamini-Hochberg on the Figure 1 experiment. The histogram stacks the p-values of changed genes on those of unchanged genes (gray); the dashed line is Storey's estimate of the unchanged floor, \(\hat\pi_0 m / 20\) per bin. The rank plot shows each gene's p-value against its rank on log scales, with the line \(p = qk/m\); the genes ranked at or before the last crossing are called, colored changed or unchanged, and the vertical line marks that crossing, k. The strip reruns seeds 1 to 50 at the current replicates and q, one dot per run, with this experiment ringed.
The bound is a promise about the average over repeated experiments, not about one list. Over 200 simulated experiments at three replicates and q = 0.05, the unchanged share of the q list averaged 2.4%, under the bound of 4.2% (85% of genes unchanged on average). In 17 of the 200 runs it went above 5%, the worst reached 7.9%, and 26 lists had no unchanged gene at all. At eight replicates the average rose to 3.6%, closer to the bound.
With three replicates, that control costs power, because most real changes cannot beat the noise of three samples. The q filter found 32.0% of the truly changed genes on average, against 59.9% at eight replicates, and it misses nearly all of the 106 small falls among the p53 targets. At two replicates it found 4.3%, and it missed CDKN1A, a planted sixfold fall in a well-expressed gene, in 113 of 200 experiments.
Of the 116 genes in Fischer's census, 7 are on the q list of Figure 1. The other members fell too, by amounts too small to show one gene at a time. Gene-set enrichment asks instead whether a set's members sit lower in the ranking than random genes would. GSEA (Subramanian et al., 2005) ranks every gene by its test statistic, here the moderated t, and walks down the list: each member steps the running score up in proportion to its |t|, every other gene steps it down a little. The enrichment score is the walk's furthest point from zero.
Figure 4. Preranked GSEA of the 116 p53 targets of Fischer (2017, Table 1) on the Figure 1 experiment; the genes are real, their falls simulated. Genes are ranked by moderated t, most raised in the knockout on the left. Top: the running score of the set and of 30 random sets of its size; middle: one tick per member; bottom: each gene's t. The shaded stretch is the leading edge, the members past the lowest point. p comes from 10,000 random sets. Hover, tap or arrow along the ranking to read a gene.
The p-value compares the score with random sets of the same size, which treats genes as independent. They are here, but in real cells co-regulated genes move together and gene permutation overstates significance. Subramanian's method permutes the samples' labels instead, which keeps those correlations but needs more samples: three against three can be relabeled only 20 ways.
Now switch to the members that Figure 1 did not call. Without the 7 called members, the other 109 still score ES = −0.59, lower than any of the random sets. That is the case GSEA was invented for: Mootha et al. (2003) used it to find "modest but coordinate changes" in oxidative phosphorylation genes in diabetic muscle. Over 200 experiments, the set and its uncalled members both beat all 1,000 random sets in every run, even at two replicates, where on average fewer than one member reached q.
These tests treat each gene, or each named set, on its own. Clustered Heatmaps looks at the whole expression matrix at once, for genes that move together and samples that look alike.