How to Build and Read a Volcano Plot from RNA-seq Data

Oak
rna-seq deseq2 visualization differential-expression bioinformatics
A winged reptile specimen with wings spread in two symmetric plumes, speckled with glowing points, lit against a black backdrop.

By the end of this guide you will have a publication-quality volcano plot built from your own RNA-seq counts, plus the diagnostic plots you need to decide whether the volcano is telling you the truth. You need a gene-by-sample count matrix (from salmon quant summarized with tximport, or from featureCounts/STAR --quantMode GeneCounts), a sample metadata table with at least one grouping variable, and R with DESeq2, apeglm, ggplot2, and ggrepel installed. You also need biological replicates. A volcano plot is a display of a statistical test, and a test with one sample per group has no variance estimate and therefore no p-axis worth plotting. If your data are longitudinal measurements on a single person, your replicates are timepoints within a state (for example, eight samples across a baseline period versus eight across an intervention period), and you should read the section on within-person designs before you trust anything downstream.

1. Start from a count matrix and a DESeq2 result table

A volcano plot is a scatterplot with two coordinates per gene: effect size on the x-axis (log2 fold change between two conditions) and evidence against the null on the y-axis (usually −log10 of the adjusted p-value). Everything interesting about the plot is decided upstream, in how those two numbers are computed, so start there rather than at the plotting call.

Use raw integer counts, never TPM or FPKM, as input. The negative binomial models in DESeq2 and edgeR need counts to estimate the mean-variance relationship, and length-normalized values break that.

library(DESeq2)

dds <- DESeqDataSetFromMatrix(
  countData = counts,          # integer matrix, genes x samples
  colData   = meta,            # data.frame with rownames == colnames(counts)
  design    = ~ batch + condition
)

# Drop genes with essentially no signal. This is a speed and multiple-testing
# convenience, not a substitute for DESeq2's own independent filtering.
keep <- rowSums(counts(dds) >= 10) >= min(table(meta$condition))
dds  <- dds[keep, ]

dds$condition <- relevel(dds$condition, ref = "control")
dds <- DESeq(dds)
resultsNames(dds)   # e.g. "condition_treated_vs_control"

Put known technical structure in the design formula. Sequencing batch, library prep date, and sex all shift expression globally, and leaving them out inflates dispersion estimates and flattens the volcano. Sex is not a nuisance variable to be discarded automatically: a multi-omics study of human muscle after musculoskeletal injury found that the transcriptomic and proteomic response differed by sex, so pooling would have averaged away the actual biology.1 Model it, and if you have enough samples, test for an interaction rather than assuming one response.

2. Shrink the log fold changes before plotting them

The raw maximum-likelihood log2 fold change for a gene with 3 counts in one group and 0 in the other can be enormous and meaningless. If you plot those, the wings of your volcano fill with low-expression noise and the genes you care about are compressed near the center. The fix is moderated (shrunken) fold changes, which pull low-information estimates toward zero while leaving well-measured genes nearly untouched.

library(apeglm)

res <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")

Use type = "apeglm" when your contrast corresponds to a single coefficient in resultsNames(dds). For an arbitrary contrast between two levels not represented by a coefficient, either relevel and rerun nbinomWaldTest, or fall back to type = "ashr", which accepts a contrast argument. Do not use the deprecated normal shrinkage.

One consequence worth understanding: lfcShrink changes the x-coordinates but leaves the p-values from the Wald test in place. So a gene can sit at a modest shrunken fold change and still have a very small p-value, which is correct. The p-value reflects how confident we are that the effect is nonzero, and the shrunken LFC reflects our best estimate of how big it is.

3. Choose thresholds before you look at the plot

The two dashed lines everyone draws on a volcano plot (vertical at ±1 log2 fold change, horizontal at adjusted p < 0.05) are conventions, not inference. The vertical lines in particular are cosmetic: filtering on an estimate after testing against a null of zero does not control anything, and it preferentially keeps noisy genes whose fold change happens to be large.

If you care about an effect size floor, test for it directly. DESeq2 will compute p-values against a composite null of |LFC| ≤ θ:

res_thr <- results(dds,
  name          = "condition_treated_vs_control",
  lfcThreshold  = log2(1.5),
  altHypothesis = "greaterAbs",
  alpha         = 0.05
)

Now a gene that passes FDR 5% has statistical evidence that its fold change exceeds 1.5×, and a single horizontal line on the volcano is a complete decision rule. The lfcShrink function accepts the same lfcThreshold argument and will return s-values, which bound the expected rate of sign errors among the genes you call, a quantity that is usually closer to what you want than an FDR on “not exactly zero.”

The alternative, used in most published integrative analyses, is the conventional double filter: adjusted p < 0.05 and |log2FC| > 1, applied to nominate candidate genes for downstream pathway and network analysis. That is how hub genes are typically selected when transcriptomics is being intersected with metabolomics or proteomics.2 It is defensible as a candidate-generation heuristic as long as you do not present the fold-change cut as if it were a test.

4. Draw the plot

EnhancedVolcano produces a good default in one call, but building it yourself in ggplot2 takes fifteen lines and removes the mystery.

library(ggplot2); library(ggrepel); library(dplyr); library(tibble)

df <- as.data.frame(res) |>
  rownames_to_column("gene") |>
  filter(!is.na(padj)) |>                       # drop filtered-out genes, see step 5
  mutate(
    padj  = pmax(padj, .Machine$double.xmin),   # avoid -log10(0) = Inf
    y     = -log10(padj),
    class = case_when(
      padj < 0.05 & log2FoldChange >  1 ~ "up",
      padj < 0.05 & log2FoldChange < -1 ~ "down",
      TRUE                              ~ "ns"
    )
  )

top <- df |> filter(class != "ns") |> slice_min(padj, n = 15)

ggplot(df, aes(log2FoldChange, y, colour = class)) +
  geom_point(size = 0.8, alpha = 0.6) +
  scale_colour_manual(values = c(up = "#B2182B", down = "#2166AC", ns = "grey70")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", linewidth = 0.3) +
  geom_vline(xintercept = c(-1, 1),     linetype = "dashed", linewidth = 0.3) +
  geom_text_repel(data = top, aes(label = gene), size = 3,
                  max.overlaps = Inf, min.segment.length = 0) +
  coord_cartesian(xlim = c(-6, 6)) +
  labs(x = "log2 fold change (apeglm shrunken)",
       y = expression(-log[10]~adjusted~italic(p)),
       colour = NULL) +
  theme_bw(base_size = 11)

Three details matter. Clip the x-axis with coord_cartesian, not xlim, so that outlying points are pushed to the edge rather than silently dropped along with their labels. Cap the p-values before taking the log, because a handful of genes will return exactly 0 in double precision and Inf will either error or blow up your axis. Label a small, deliberate set of genes: fifteen labels are readable, sixty are a grey smear, and ggrepel will spend minutes trying.

For 20,000 points, overplotting hides density in the center. Add geom_point(alpha = 0.3) for the non-significant class, or rasterize that layer with ggrastr::rasterise() so the PDF stays under a megabyte.

5. Read the plot, including the parts that are missing

A healthy volcano is roughly symmetric about x = 0, with a dense low-y cloud narrowing as |LFC| grows, and it rises into two wings. The vertical axis is confidence, the horizontal axis is magnitude, and a gene in the upper corner has both. Genes high on the y-axis but near x = 0 are small effects measured precisely, which is common for very highly expressed genes and is not an artifact. Genes far out on x with low y are usually low-count genes, and after shrinkage there should be far fewer of them.

Several shapes signal a problem. A volcano shifted bodily to one side means your normalization has failed or a genuine global shift exists (a large fraction of the transcriptome changing, which violates the median-of-ratios assumption and calls for spike-ins or a housekeeping-based normalization). A volcano with essentially no height, where the maximum −log10(padj) is around 1.3, means you are underpowered or your groups are not separated. A volcano where thousands of genes reach padj < 1e-20 usually means an unmodeled batch that lines up with condition, or that you ran a per-cell test on single-cell data (see step 7).

Pay attention to what is not plotted. DESeq2’s independent filtering sets padj = NA for genes with low mean normalized counts, and results() also sets pvalue = NA for genes flagged as Cook’s distance outliers. Those genes vanish from the volcano entirely. Check sum(is.na(res$padj)) and metadata(res)$filterThreshold so you know how many genes were excluded and at what expression level. If a gene you expected is absent, that is usually why.

6. Keep the MA plot next to it

The MA plot has the same fold change but plots it against mean normalized expression on a log scale instead of significance. The two answer different questions. The volcano asks which genes have strong, confident effects. The MA plot asks whether the fold-change estimates behave sensibly across the dynamic range, which is a quality-control question you cannot answer from a volcano.

plotMA(res, ylim = c(-5, 5), alpha = 0.05)

A good MA plot is a horizontal band centered on zero, with the band widening at low counts and tightening at high counts, and significant points scattered above and below. A band that drifts off zero at high expression indicates a normalization problem. A trumpet flare of significant calls only at the extreme left indicates that shrinkage was not applied or failed. We look at the MA plot first and the volcano second, every time.

7. Single-cell data: aggregate to pseudobulk first

Volcano plots from single-cell RNA-seq are the most commonly misread version of this figure. If you run Seurat::FindMarkers() with the default Wilcoxon test across thousands of cells from a handful of donors, every gene with a modest difference reaches an adjusted p-value of 1e-100, because the test treats each cell as an independent replicate. Cells from one donor are not independent. The resulting volcano is tall, dense, and mostly reports how many cells you sequenced.

For comparisons between conditions or individuals, sum counts within each cell type and each sample, then run the same DESeq2 pipeline on that pseudobulk matrix:

library(Seurat)
pb <- AggregateExpression(obj, group.by = c("celltype", "sample"),
                          assays = "RNA", slot = "counts",
                          return.seurat = FALSE)$RNA
# then split by celltype and run DESeqDataSetFromMatrix per cell type

The volcano that comes out is shorter, has fewer points, and is a statement about your donors rather than your cells. Per-cell tests remain appropriate for their original purpose, identifying marker genes that distinguish clusters within one sample, where the question is descriptive.

8. The mass-spectrometry volcano is a different plot with the same axes

Proteomics volcano plots (from MaxQuant, DIA-NN, Spectronaut, or Perseus) use log2 ratios of normalized intensities on x and −log10 p on y, but the underlying statistics differ in three ways. Intensities are continuous, so you use a moderated t-test (limma::eBayes) rather than a negative binomial model. Missingness is structural: a protein absent from one condition because it fell below the detection limit is not the same as one missing at random, and imputation strategy (left-shifted Gaussian at 1.8 standard deviations below the sample mean is the Perseus default) visibly creates a cluster of points at extreme fold changes. And the significance boundary is often drawn as a curved line from an SAM-style permutation FDR with an s0 parameter, which is why proteomics volcanoes have hyperbolic rather than rectangular cutoffs.

If you are plotting transcript and protein volcanoes side by side, expect modest agreement in direction and poor agreement in magnitude. Integrated studies routinely find that only a minority of differentially expressed transcripts have a matching protein-level change, which is a real feature of post-transcriptional regulation rather than a technical failure.3 Cross-omics concordance is more convincing evidence than either layer alone, which is the logic behind using paired transcriptomic and metabolomic volcanoes to triangulate a pathway before committing to follow-up.4

Common problems

Your volcano is empty and nothing passes FDR. First check the PCA (plotPCA(vst(dds))) to see whether your conditions separate at all. If they do not, no plot will rescue the experiment. If they separate along PC2 while PC1 tracks batch, add batch to the design. If you have three versus three samples and a subtle perturbation, an empty volcano is the correct answer, and reporting it as such is better than dropping the FDR threshold to 0.1 after the fact.

Everything is significant. Look for a confound that aligns perfectly with condition, such as all treated samples sequenced in one run. Also check for a sample swap by correlating samples against each other on variance-stabilized counts.

Points at the top edge with Inf values. Adjusted p-values underflow to zero for very strong effects. Cap them as shown in step 4, and state in the legend that the axis is capped, so a reader does not interpret a flat row of points at the ceiling as a real cluster.

Gene labels are wrong or missing. Ensembl IDs with version suffixes (ENSG00000141510.16) will not join against an unversioned annotation table. Strip with sub("\\..*$", "", gene) before mapping, and expect a few percent of IDs to have no symbol.

Fold changes disagree with a heatmap. A heatmap of z-scored VST values shows relative patterns across samples, and the volcano shows a model coefficient adjusted for covariates in your design. When a design includes batch, the two will not match visually. Use limma::removeBatchEffect() on the VST matrix before heatmapping if you want them to correspond.

Finally, a note on what a volcano plot is for. It is a screening display that ranks candidates for further work: pathway enrichment, targeted assays, integration with other omics layers. It is not a result on its own, and a gene’s position on it says nothing about your health. If something in your own data prompts a health question, the next step is a conversation with a clinician who can order validated clinical assays, not a deeper stare at the plot.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Alexander R. Keeble, Allison M. Owen, Nicholas T. Thomas, et al. Multi-omics analysis reveals sex-specific etiology of human muscle weakness following musculoskeletal injury. BMC Medicine, 2026. https://doi.org/10.1186/s12916-026-04818-8 ↩

  2. Huan Liu, Guanming Qi, Shengrong Ouyang, et al. Integrated Multi‐Omics Analysis Identifies PDK4 and ACOT1 as Metabolic Hub Genes Associated With Myocardial Fibrosis in Diabetic Cardiomyopathy. Journal of Diabetes Research, 2026. https://doi.org/10.1155/jdr/1952672 ↩

  3. Die Dai, Fandie Dai, Jingchao Chen, et al. Integrated multi-omics reveal important roles of gut contents in intestinal ischemia–reperfusion induced injuries in rats. Communications Biology, 2022. https://doi.org/10.1038/s42003-022-03887-8 ↩

  4. Ziqi Cheng, Hua Zhu, Shi Feng, et al. Cross-Species Multi-Omics Analysis Reveals Myeloid-Driven Endothelial Oxidative Stress in Ischemic Stroke. Frontiers in Bioscience-Landmark, 2025. https://doi.org/10.31083/fbl37429 ↩