Differential Gene Expression Analysis in R: A Practical Guide

By Dr. Zubair Khalid, DVM, MS, PhD ·

Differential Gene Expression Analysis in R: A Practical Guide

Introduction to Differential Gene Expression Analysis in R

What is Differential Gene Expression?

Differential gene expression (DGE) analysis is the systematic comparison of transcript abundance between two or more biological conditions to identify genes whose expression levels change significantly. The fundamental question is statistical: given the measured RNA or protein abundance for each gene across replicate samples from different conditions, which genes show differences larger than what would be expected from technical noise and biological variation alone?

The output of a DGE analysis is a ranked list of genes with associated effect sizes (typically log2 fold changes) and measures of statistical significance (p-values and adjusted p-values). This list serves as the entry point for downstream biological interpretation, including pathway enrichment, regulatory network inference, and biomarker discovery. The statistical challenge is substantial: a typical RNA-seq experiment quantifies 20,000–25,000 genes simultaneously, and the number of biological replicates is often small (3–5 per condition), making reliable variance estimation difficult.

DGE analysis differs fundamentally from simply comparing normalized expression values. The key distinction lies in the statistical framework: DGE methods model the relationship between mean expression and variance across the entire transcriptome, then use this relationship to assess whether observed differences for individual genes are credible. This is essential because genes with low expression have inherently higher sampling variability, and without proper modeling, lowly expressed genes would dominate the list of "significant" results.

Why R for DGE Analysis?

R has become the de facto standard environment for DGE analysis for several concrete reasons. First, the Bioconductor project provides a coordinated ecosystem of packages designed specifically for genomic data, with standardized data structures (e.g., SummarizedExperiment, DGEList) that enforce consistent practices. Second, the statistical rigor of R's base functions and the availability of advanced packages like DESeq2, edgeR, and limma mean that the most widely used and benchmarked methods are immediately accessible. Third, R's plotting capabilities, particularly through ggplot2, enable the publication-quality visualizations that are essential for both exploratory analysis and final reporting.

The typical workflow in R proceeds through distinct stages: data import and quality control, normalization, statistical modeling, hypothesis testing, and visualization. Each stage has dedicated tools, and the modular nature of R allows researchers to combine packages flexibly. This guide covers the three most widely used approaches: DESeq2 and edgeR for RNA-seq count data, and limma for both microarray data and RNA-seq data via the voom transformation.

Data Preparation and Quality Control

Importing Data into R

The starting point for DGE analysis is a count matrix: a rectangular table where rows represent genes (typically annotated with Ensembl or Entrez IDs) and columns represent individual samples. Each cell contains the number of sequencing reads that mapped to that gene in that sample. For microarray data, the analogous input is a matrix of fluorescence intensities or log2-transformed expression values.

Count matrices are typically generated by alignment and quantification tools such as STAR, HISAT2, or Salmon. The output should be imported into R using read.table() or read.csv(), with careful attention to row names (gene identifiers) and column names (sample identifiers). For example:

counts <- read.table("counts.txt", header = TRUE, row.names = 1)

The metadata table, containing information about each sample (condition, batch, sex, etc.), must be organized with rows matching the columns of the count matrix. The row names of the metadata must exactly match the column names of the count matrix; mismatches are a common source of errors. A typical metadata table includes columns for sample ID, condition (e.g., "treated" vs "control"), and any covariates to be modeled.

Quality Control Metrics

Before any statistical analysis, the data must pass through quality control at two levels: sample-level and gene-level. Sample-level QC aims to identify outlier samples, failed libraries, or sample swaps. Key metrics include:

  • Total read count per sample: Libraries with dramatically different total counts (e.g., 5-fold differences) may indicate technical problems.
  • Percentage of reads mapped to genes: Low mapping rates suggest contamination or alignment issues.
  • Percentage of reads in mitochondrial genes: High mitochondrial content (>20%) can indicate cell lysis or poor RNA quality.
  • RNA integrity number (RIN): For paired-end data, the RIN value from the Bioanalyzer should be recorded in the metadata.

Principal component analysis (PCA) on the most variable genes is the most informative sample-level QC visualization. Samples from the same condition should cluster together, and the first few principal components should separate known biological groups. If samples cluster by an unintended variable (e.g., sequencing batch) rather than the condition of interest, this indicates a batch effect that must be addressed in the statistical model.

Gene-level QC involves filtering out genes with very low counts across all samples. These genes carry little information and their inclusion inflates the multiple testing burden. A common filter retains genes with at least 10 counts in at least 3 samples (for a typical experiment with 3 replicates per condition). More formally, DESeq2 uses the independent filtering approach, which automatically determines an optimal filter threshold based on the relationship between mean counts and statistical power.

Filtering and Normalization

Normalization is the process of adjusting raw counts to make samples comparable. The necessity arises because sequencing depth varies between samples: a sample with 30 million reads will have systematically higher counts than one with 20 million reads, even if the underlying expression is identical. The simplest normalization is counts per million (CPM), which divides each count by the total library size and multiplies by one million.

However, CPM is inadequate for DGE analysis because it does not account for compositional differences between samples. If a few highly expressed genes are upregulated in one condition, they consume a larger fraction of the sequencing budget, causing all other genes to appear downregulated. The Trimmed Mean of M-values (TMM) method, implemented in edgeR, and the median-of-ratios method, implemented in DESeq2, both address this issue by computing scaling factors based on the majority of genes that are not differentially expressed.

For DESeq2, the normalization is performed internally during the DESeq() function call. For edgeR, the calcNormFactors() function computes TMM scaling factors. For limma with voom, the voom() function incorporates library size normalization internally. The key point is that normalized counts are used only for visualization and exploratory analysis; the statistical model for DGE uses the raw counts with the scaling factors incorporated into the model.

Statistical Models for Differential Expression

Negative Binomial Model

RNA-seq count data are integer-valued and exhibit overdispersion: the variance across biological replicates is greater than the mean, violating the assumption of the Poisson distribution. The negative binomial (NB) distribution addresses this by introducing a dispersion parameter that models the extra variability.

For gene g in sample j, the count K<sub>gj</sub> is modeled as NB with mean μ<sub>gj</sub> and dispersion φ<sub>g</sub>. The mean is modeled as:

μ<sub>gj</sub> = q<sub>gj</sub> × s<sub>j</sub>

where s<sub>j</sub> is the size factor for sample j (accounting for library size and composition) and q<sub>gj</sub> is the expression level proportional to the true abundance. The variance is:

Var(K<sub>gj</sub>) = μ<sub>gj</sub> + φ<sub>g</sub> × μ<sub>gj</sub>²

The dispersion parameter φ<sub>g</sub> captures biological variability beyond technical noise. With only 3–5 replicates, estimating φ<sub>g</sub> independently for each gene is unreliable. Both DESeq2 and edgeR solve this by sharing information across genes: they model the dispersion as a function of the mean expression and shrink individual gene dispersions toward this trend. This empirical Bayes approach is the key innovation that makes DGE analysis feasible with small sample sizes.

Limma and Empirical Bayes

Limma was originally developed for microarray data, where expression measurements are continuous and approximately normally distributed after log2 transformation. The core model is a linear regression for each gene:

y<sub>g</sub> = Xβ<sub>g</sub> + ε<sub>g</sub>

where y<sub>g</sub> is the vector of log2 expression values for gene g, X is the design matrix, β<sub>g</sub> is the vector of coefficients, and ε<sub>g</sub> is the error term.

The innovation of limma is the empirical Bayes moderation of the gene-wise variance. Rather than using each gene's variance estimate directly, limma borrows information across all genes to shrink individual variances toward a common prior. This results in moderated t-statistics that follow a t-distribution with augmented degrees of freedom. The effect is particularly beneficial for experiments with few replicates, where individual variance estimates are noisy.

For RNA-seq data, the voom function transforms count data to log2-counts per million (log2-CPM) and estimates the mean-variance relationship, then assigns a precision weight to each observation. The weighted linear model is then fitted with limma's empirical Bayes machinery. This approach has been shown to perform comparably to DESeq2 and edgeR while offering greater flexibility in experimental design.

Design Matrices and Contrasts

The design matrix encodes the experimental design and is the bridge between the biological question and the statistical model. Each row corresponds to a sample, and each column corresponds to a coefficient in the model. For a simple two-group comparison (treated vs. control), the design matrix can be parameterized in two ways:

  1. Treatment-contrast parameterization: The first column is the intercept (all 1s), and the second column is an indicator variable (0 for control, 1 for treated). The coefficient for the second column directly estimates the log2 fold change between groups.
  1. Group-means parameterization: Two columns, each an indicator for one group. The coefficients estimate the mean expression in each group, and the contrast treated - control estimates the log2 fold change.

For more complex designs (e.g., time course, factorial designs, or paired samples), the design matrix must include all relevant covariates. The key principle is that the design matrix must account for all sources of variation that could confound the comparison of interest. Failing to include a known batch variable in the design matrix is a common and serious error.

Contrasts are linear combinations of coefficients that define the comparisons of interest. For example, in a factorial design with two factors (treatment and time), the interaction contrast tests whether the treatment effect differs between time points. The contrasts.fit() function in limma and the contrast argument in DESeq2's results() function implement these comparisons.

Running Differential Expression Analysis with DESeq2

Creating DESeqDataSet

The DESeq2 workflow begins with the construction of a DESeqDataSet object from the count matrix and metadata. The DESeqDataSetFromMatrix() function requires three arguments: the count matrix, the column data (metadata), and the design formula.

library(DESeq2)

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData = metadata,
  design = ~ condition
)

The design formula specifies the model to be fitted. The simplest design is ~ condition, where condition is a column in the metadata. For paired samples or batch effects, the design would include additional terms, e.g., ~ batch + condition. The condition of interest should be the last term in the formula.

Before running the analysis, the reference level of the condition factor should be set explicitly to ensure the correct comparison:

dds$condition <- relevel(dds$condition, ref = "control")

Running the DESeq Workflow

The core analysis is executed with a single function call:

dds <- DESeq(dds)

This function performs, in sequence: estimation of size factors (median-of-ratios normalization), estimation of gene-wise dispersions, fitting of the dispersion-mean trend, shrinkage of dispersions toward the trend, and fitting of the NB GLM for each gene. The default Wald test is used for hypothesis testing, but a likelihood ratio test can be specified with test = "LRT" for comparing nested models.

The dispersion estimation step is worth understanding in detail. For each gene, DESeq2 first estimates a raw dispersion using maximum likelihood. It then fits a curve through these estimates as a function of the mean normalized count. Finally, each gene's dispersion is shrunk toward the curve, with the amount of shrinkage determined by the number of replicates and the gene's expression level. Genes with more replicates and higher expression receive less shrinkage because their raw estimates are more reliable.

Extracting and Shrinking Results

The results are extracted with the results() function:

res <- results(dds, contrast = c("condition", "treated", "control"))

The resulting results object is a DataFrame with columns for baseMean (mean normalized count across all samples), log2FoldChange, lfcSE (standard error of the log2 fold change), stat (Wald statistic), pvalue, and padj (Benjamini-Hochberg adjusted p-value).

A crucial step is the application of log2 fold change shrinkage. The raw log2 fold changes from the GLM are not shrunk, and for genes with low counts or high dispersion, these estimates can be unrealistically large. The lfcShrink() function applies an adaptive shrinkage estimator (apeglm or normal) that produces more accurate effect size estimates:

res_shrunk <- lfcShrink(dds, contrast = c("condition", "treated", "control"), type = "apeglm")

The shrunk log2 fold changes are preferred for ranking genes and for downstream visualization, as they reduce the noise in effect size estimates while preserving the statistical significance from the Wald test.

Running Differential Expression Analysis with edgeR

Creating DGEList and Normalization

The edgeR workflow uses the DGEList object. The count matrix and metadata are combined as follows:

library(edgeR)

dge <- DGEList(counts = counts, group = metadata$condition)

The DGEList function also calculates library sizes and CPM values. The next step is filtering low-count genes. edgeR provides the filterByExpr() function, which automatically determines a filter threshold based on the experimental design:

keep <- filterByExpr(dge, design = model.matrix(~ metadata$condition))
dge <- dge[keep, , keep.lib.sizes = FALSE]

Normalization is performed with the TMM method:

dge <- calcNormFactors(dge, method = "TMM")

TMM computes scaling factors by comparing each sample to a reference sample, using the weighted mean of log ratios of gene expression after trimming the most extreme values. The scaling factors are typically close to 1, with values between 0.7 and 1.3 indicating reasonable library composition.

Dispersion Estimation

edgeR estimates three types of dispersion: common, trended, and tagwise (gene-wise). The common dispersion is a single value shared by all genes, the trended dispersion varies as a function of the mean expression, and the tagwise dispersion is the gene-specific value shrunk toward the trend.

dge <- estimateDisp(dge, design = design)

This function estimates the common dispersion, fits the trend, and computes tagwise dispersions using an empirical Bayes approach. The output includes the common.dispersion, trended.dispersion, and tagwise.dispersion elements. The biological coefficient of variation (BCV), which is the square root of the dispersion, is typically between 0.1 and 0.4 for well-controlled experiments.

Testing and Extracting Results

edgeR offers several testing approaches. The most commonly used is the exact test, analogous to Fisher's exact test but adapted for the negative binomial distribution:

et <- exactTest(dge, pair = c("control", "treated"))
topTags(et, n = 10)

For more complex designs, the quasi-likelihood F-test is recommended:

fit <- glmQLFit(dge, design)
qlf <- glmQLFTest(fit, contrast = c(-1, 1))
topTags(qlf, n = 10)

The glmQLFit function fits a negative binomial GLM with quasi-likelihood dispersion estimation, which provides robust inference when the dispersion estimates are uncertain. The topTags() function extracts the results table with columns for logFC, logCPM, F (or LR), PValue, and FDR.

Running Differential Expression Analysis with limma

Limma for Microarray

For microarray data, limma operates directly on log2-transformed expression values. The workflow is:

library(limma)

# Assume expr is a matrix of log2 intensities
design <- model.matrix(~ metadata$condition)
fit <- lmFit(expr, design)
fit <- eBayes(fit)
results <- topTable(fit, coef = 2, number = Inf)

The lmFit() function fits a linear model for each gene. The eBayes() function applies empirical Bayes moderation to the gene-wise variances, shrinking them toward a common value. The topTable() function extracts the results with moderated t-statistics, p-values, and adjusted p-values.

Voom for RNA-seq

For RNA-seq count data, limma's voom transformation converts counts to log2-CPM and models the mean-variance relationship:

design <- model.matrix(~ metadata$condition)
v <- voom(dge, design, plot = TRUE)
fit <- lmFit(v, design)
fit <- eBayes(fit)
results <- topTable(fit, coef = 2, number = Inf)

The voom() function takes a DGEList object (after filtering and TMM normalization) and produces an EList object containing log2-CPM values and precision weights. The precision weights account for the fact that the variance of log2-CPM is not constant across the expression range: lowly expressed genes have higher variance. The subsequent linear modeling uses these weights, providing a flexible framework that handles both simple and complex designs.

Contrasts and Hypothesis Testing

Limma's contrast machinery is particularly powerful for complex designs. After fitting the linear model, contrasts are defined and tested:

contrast_matrix <- makeContrasts(
  treated_vs_control = treated - control,
  levels = design
)
fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)
results <- topTable(fit2, coef = "treated_vs_control", number = Inf)

This approach allows testing multiple comparisons simultaneously and is essential for factorial designs, time courses, and other designs where the comparisons of interest are not simply individual coefficients.

Interpreting and Visualizing Results

Understanding the Results Table

The results table from any DGE analysis contains the same essential columns. The baseMean (DESeq2) or logCPM (edgeR) indicates the average expression level, which is important for interpreting the biological relevance of a change. The log2FoldChange or logFC is the effect size: a value of 1 means a 2-fold increase, -1 means a 2-fold decrease. The standard error (lfcSE) quantifies the uncertainty in the effect size estimate.

The p-value tests the null hypothesis that the true log2 fold change is zero. The adjusted p-value (padj or FDR) corrects for multiple testing using the Benjamini-Hochberg procedure, which controls the false discovery rate (FDR). An FDR of 0.05 means that approximately 5% of the genes called significant are expected to be false positives.

A common threshold for calling differentially expressed genes is padj < 0.05 and |log2FoldChange| > 1 (corresponding to a 2-fold change). However, these thresholds are arbitrary and should be adjusted based on the biological context. For example, a genome-wide screen might use a more stringent FDR of 0.01, while a focused candidate gene study might accept an FDR of 0.1.

Visualization Techniques

The MA plot (also called a Bland-Altman plot) displays the log2 fold change (M) versus the mean expression (A). Each point represents a gene, with significant genes highlighted. The plot should show a symmetric distribution around zero for non-significant genes, with significant genes appearing as points far from zero, particularly at higher mean expression where the variance is lower.

The volcano plot displays log2 fold change on the x-axis and -log10(p-value) on the y-axis. This visualization emphasizes both the effect size and the statistical significance, with significant genes appearing in the upper left (downregulated) and upper right (upregulated) corners.

Heatmaps of the most variable or most significant genes provide a sample-level view of the expression patterns. The pheatmap package or ComplexHeatmap can be used, with rows scaled to z-scores to make genes comparable. Hierarchical clustering of samples should recapitulate the biological conditions.

Gene Annotation and Enrichment

The results table contains gene identifiers (e.g., Ensembl IDs) that must be mapped to gene symbols and functional annotations. The AnnotationDbi and org.Hs.eg.db packages provide this mapping. Once annotated, the list of significant genes can be tested for enrichment of Gene Ontology (GO) terms or KEGG pathways using packages like clusterProfiler or topGO. For a detailed overview of the available tools and their applications, see the Gene Ontology Analysis Tool and Gene Ontology Pathway Enrichment resources. The Gene Ontology Analysis Online guide provides practical guidance for web-based tools.

Common Pitfalls and Best Practices

Batch Effects and Confounding

The most serious error in DGE analysis is failing to account for batch effects. If samples from different conditions are processed in different batches (e.g., different sequencing runs, different RNA extraction dates), the batch effect is confounded with the condition effect, and any observed differences could be due to technical rather than biological variation.

The solution is to include batch as a covariate in the design formula: ~ batch + condition. This is only possible if each batch contains samples from both conditions. If the design is confounded (all treated samples in one batch, all controls in another), no statistical correction can recover the true biological signal. The experiment must be redesigned.

The ComBat function from the sva package can be used to remove batch effects when they are known but not included in the design. However, this should be done with caution, as it can remove biological signal if the batch variable correlates with the condition of interest.

Multiple Testing Correction

With 20,000 genes tested simultaneously, the expected number of false positives at p < 0.05 is 1,000. Multiple testing correction is therefore essential. The Benjamini-Hochberg procedure, which controls the FDR, is the standard approach and is implemented in all major packages. The more conservative Bonferroni correction controls the family-wise error rate but is often too stringent for exploratory analyses.

A related issue is the interpretation of p-values near the significance threshold. A gene with an adjusted p-value of 0.051 is not meaningfully different from one with 0.049, and the binary classification of "significant" vs. "not significant" should not be over-interpreted. The ranked list of genes with their effect sizes and p-values is more informative than the binary call.

Reproducibility and Reporting

DGE analysis involves many decisions—filtering thresholds, normalization methods, statistical tests, significance cutoffs—that can substantially affect the results. Best practices include:

  1. Record the version numbers of R and all packages used.
  2. Set a random seed if any stochastic steps are involved.
  3. Document all filtering and normalization steps in the analysis script.
  4. Report the number of genes tested and the number called significant at each threshold.
  5. Provide the full results table as a supplementary file, not just the significant genes.

The sessionInfo() function in R provides a complete record of the software environment. The renv package can be used to create reproducible project environments.

Frequently Asked Questions

What is the best R package for differential gene expression analysis?

There is no single "best" package; DESeq2, edgeR, and limma-voom all perform well and produce largely concordant results on typical datasets. The choice depends on the experimental design and personal preference. DESeq2 is particularly well-suited for designs with multiple factors and provides robust dispersion estimation. edgeR offers a wider range of statistical tests, including quasi-likelihood methods. limma-voom is the most flexible for complex designs and provides the most comprehensive framework for contrasts. For most analyses, any of the three will give reliable results, and the choice can be based on familiarity and the specific features needed.

How do I choose between DESeq2 and edgeR?

Both packages use the negative binomial model and empirical Bayes shrinkage, and their results are typically highly correlated. DESeq2 uses a different normalization method (median-of-ratios vs. TMM) and a different dispersion shrinkage approach. In practice, the choice often comes down to the experimental design: DESeq2 handles designs with multiple factors and continuous covariates more naturally, while edgeR's quasi-likelihood tests are more robust when the dispersion estimates are uncertain. If you are unsure, run both and compare the results; genes that are significant in both are the most reliable.

What is the difference between normalized counts and differential expression results?

Normalized counts are the adjusted expression values that account for library size and composition differences between samples. They are suitable for visualization and exploratory analysis but do not provide any statistical inference. Differential expression results are the output of the statistical model: for each gene, they provide an effect size (log2 fold change), a measure of uncertainty (standard error), and a p-value testing whether the effect is significantly different from zero. The differential expression analysis uses the normalized counts as input but applies the statistical model to determine which differences are credible given the biological variability.

How do I correct for multiple testing in differential expression analysis?

The standard approach is the Benjamini-Hochberg procedure, which controls the false discovery rate (FDR). This is implemented automatically in DESeq2 (as padj), edgeR (as FDR), and limma (as adj.P.Val). The FDR is the expected proportion of false positives among the genes called significant. An FDR threshold of 0.05 is commonly used, meaning that approximately 5% of the significant genes are expected to be false positives. For more stringent control, the Bonferroni correction can be applied, but this is rarely recommended for exploratory analyses.

What is a design formula in DESeq2 and how do I set it up?

The design formula specifies the statistical model that DESeq2 fits for each gene. It is written in R formula syntax, with the tilde (~) separating the response (implicitly the expression) from the predictors. For a simple two-group comparison, the design is ~ condition. For a paired design, it is ~ patient + condition. For a factorial design with two factors, it is ~ factor1 + factor2 + factor1:factor2 (the last term is the interaction). The condition of interest should be the last term in the formula, and the reference level of the factor should be set explicitly with relevel().

How do I handle batch effects in differential expression analysis?

The preferred approach is to include the batch variable in the design formula, e.g., ~ batch + condition. This requires that each batch contains samples from all conditions. If the batch effect is not known but is suspected (e.g., from PCA analysis), the ComBat function from the sva package can be used to remove it. However, this should be done with caution, as it can remove biological signal if the batch variable correlates with the condition. The best approach is to design the experiment to avoid confounding and to include known technical variables in the design.

What is the voom transformation in limma?

The voom transformation converts RNA-seq count data to log2-counts per million (log2-CPM) and models the mean-variance relationship. For each gene, the log2-CPM values are computed, and the variance is estimated as a function of the mean. A precision weight is then assigned to each observation, reflecting the reliability of that measurement. The weighted linear model is fitted using limma's empirical Bayes machinery. This approach allows limma, originally designed for microarray data, to be applied to RNA-seq data while accounting for the heteroscedasticity inherent in count data.

Why are my p-values all NA in DESeq2 results?

This typically occurs when the design formula is misspecified or when there is no residual degrees of freedom. Common causes include: (1) the design includes a factor with only one level, (2) the design is saturated (more coefficients than samples), (3) the condition factor has a level with no replicates, or (4) the contrast specified in results() does not match the levels of the factor. Check the design formula, ensure that each group has at least two replicates, and verify that the contrast is correctly specified. The resultsNames(dds) function lists the available coefficients and can help diagnose the issue.

Key Takeaways

  • Differential gene expression analysis in R relies on negative binomial models (DESeq2, edgeR) or linear models with empirical Bayes moderation (limma), all of which share information across genes to handle small sample sizes.
  • Proper experimental design, including adequate replication (at least 3 biological replicates per condition) and balanced batch assignment, is more important than any statistical refinement.
  • Normalization (TMM in edgeR, median-of-ratios in DESeq2) corrects for library size and composition differences and is a prerequisite for reliable cross-sample comparisons.
  • The design matrix must include all known sources of variation, including batch effects, to avoid confounding and false positive results.
  • Multiple testing correction using the Benjamini-Hochberg FDR procedure is essential when testing thousands of genes simultaneously.
  • Log2 fold change shrinkage (apeglm in DESeq2, quasi-likelihood in edgeR) produces more accurate effect size estimates, particularly for lowly expressed genes.
  • Reproducibility requires documenting all analysis decisions, recording software versions, and providing full results tables, not just lists of significant genes.

Related Clinical & Scientific Guides