Volcano plot maker for differential expression

A volcano plot shows every gene (or protein, metabolite) at $x={\log }_2$ fold change and $y=-{\log }_{10}$ p-value, so large, reliable changes sit at the top left (down) and top right (up) and look like the plume of a volcano (Li 2012). A log₂ fold change of 1 is a 2-fold change; 0.585 is 1.5-fold. On the y axis, $p=0.05$ is at 1.30, $p=0.001$ at 3 and $p={10}^{-10}$ at 10.

Classification. A gene is up if ${\log }_2\mathrm{FC}\ge t$ and its p-value is below $\alpha$, down if ${\log }_2\mathrm{FC}\le -t$ and $p<\alpha$, and not significant otherwise. You choose whether $\alpha$ applies to the adjusted p-value (the false discovery rate, the standard for genome-wide data) or the raw p-value. The y axis always shows the raw p-value, as in DESeq2 and EnhancedVolcano; with FDR the horizontal line is drawn at the largest raw p whose adjusted p passes, so the line still separates significant from not.

Benjamini–Hochberg adjustment. With $m$ p-values sorted from smallest to largest, $p_{(1)}\le \cdots \le p_{(m)}$, the adjusted value of the $i$-th is $p_{\mathrm{adj}}(i)=\underset{j\ge i}{min}(p_{(j)}\cdot m/j)$, capped at 1. Calling every gene with $p_{\mathrm{adj}}<\alpha$ keeps the expected proportion of false discoveries among the hits at or below $\alpha$ (Benjamini and Hochberg 1995). This is R’s p.adjust(method = "BH") and the padj column of DESeq2. If your table already has an adjusted column (padj, FDR, q-value, adj.P.Val), the tool uses it.

By hand. Five p-values 0.005, 0.01, 0.03, 0.04, 0.5 ($m=5$): the raw products $p\cdot m/\text{rank}$ are 0.025, 0.025, 0.05, 0.05, 0.5; taking the running minimum from the largest down leaves them unchanged, so at $\mathrm{FDR}<0.05$ only the first two are significant, although four have $p<0.05$.

p-values reported as 0 (below the precision of the software) are drawn at the smallest non-zero p in the data divided by 10, so they stay on the plot.

Sources: Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B 1995;57:289–300, doi:10.1111/j.2517-6161.1995.tb02031.x. Li W. Volcano plots in analyzing differential expressions with mRNA microarrays. J Bioinform Comput Biol 2012;10:1231003, doi:10.1142/S0219720012310038.

The built-in table is synthetic: 1,000 genes with a log₂ fold change and a raw p-value, of which 40 were simulated as changed (24 up, including hypoxia-response genes such as VEGFA and CA9, and 16 down, proliferation genes such as MKI67). There is no adjusted column, so the tool computes Benjamini–Hochberg adjusted p-values.

1. Adjustment. The smallest p-value, DDIT4 (log₂FC 2.21, $p=3.42\times {10}^{-12}$), has rank 1, so $p_{\mathrm{adj}}=3.42\times {10}^{-12}\times 1000/1$ $=3.42\times {10}^{-9}$. VEGFA (log₂FC 3.06, $p=2.48\times {10}^{-10}$) ranks 5th: $2.48\times {10}^{-10}\times 1000/5=4.96\times {10}^{-8}$.
2. Thresholds. With $\mathrm{FDR}<0.05$ and $|{\log }_2\mathrm{FC}|\ge 1$ (2-fold), the plot shows 25 up and 16 down. 103 genes pass $\mathrm{FDR}<0.05$ on their own; the largest raw p among them is 0.00469, which is where the horizontal line is drawn ($-{\log }_{10}0.00469=2.33$). 230 genes have a raw $p<0.05$, which shows why the fold-change cut-off and the adjustment matter.
3. One false positive. Of the 25 up genes, 24 are the simulated changes and one (GENE0386, log₂FC 1.01, $p=0.00061$) is a null gene that crossed both lines by chance: 1 of 41 calls, within the 5% false discovery rate.
4. Looser threshold. At 1.5-fold ($|{\log }_2\mathrm{FC}|\ge 0.585$) the counts rise to 35 up and 27 down; at $\mathrm{FDR}<0.01$ and 2-fold they fall to 24 up and 15 down.