Overview
Differential expression (DE) analysis is the statistical core of transcriptomics, determining which genes show meaningful expression changes between experimental conditions. Whether comparing treated versus untreated cells, tumor versus normal tissue, or time-course samples, DE analysis transforms raw count data into biological insight. The challenge lies in distinguishing true biological signals from technical noise while accounting for the multiple testing burden inherent in genome-wide measurements. Modern DE methods use sophisticated statistical models that have been benchmarked extensively on real and simulated datasets.
Methods
DE analysis typically starts with a count matrix of reads per gene per sample. Normalization methods (TMM, RLE, or quantile normalization) adjust for library size and compositional biases. Popular tools include DESeq2, which models counts with a negative binomial distribution and uses shrinkage estimation for dispersion; edgeR, which uses empirical Bayes methods; and limma-voom, which applies linear modeling to log-transformed counts transformed with precision weights. For non-model organisms or experiments without [a reference genome, tools like Sailfish or Kallisto perform alignment-free quantification. Results are summarized as log2 fold changes and adjusted p-values (Benjamini-Hochberg correction). Principal component analysis (PCA) and heatmaps provide global views of expression patterns.
Applications
DE analysis is central to virtually every transcriptomics study. It identifies biomarkers for disease diagnosis and prognosis, reveals drug mechanisms of action, and characterizes cellular responses to environmental stimuli. In clinical settings, DE analysis of patient biopsies can stratify cancers for targeted therapy. The method is closely related to RT-PCR validation experiments, which confirm candidate genes, and builds upon DNA microarrays and gene expression techniques. DE results also power gene set enrichment analysis (GSEA) to identify affected pathways and functional categories.
Practical Protocol
A standard DE analysis in R starts by reading the count matrix and sample metadata into a DESeqDataSet object: dds <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~ condition). The design formula specifies the experimental model; a simple two-group comparison uses ~ condition, while paired or batch-corrected designs use ~ batch + condition or ~ patient + condition. Run the full DE pipeline with dds <- DESeq(dds), which performs normalization, dispersion estimation, and Wald testing. Extract results with results(dds, contrast = c("condition", "treated", "control")), which returns a table with columns for baseMean, log2FoldChange, lfcSE, stat, pvalue, and padj. The padj column contains Benjamini-Hochberg adjusted p-values accounting for multiple testing across thousands of genes. An alternative workflow with edgeR uses DGEList, calcNormFactors, estimateDisp, and exactTest or glmFit/glmLRT for complex designs. Interpret log2 fold changes: a value of 1 means a 2-fold increase in the treated condition, while -1 means a 2-fold decrease. Shrinkage of log2 fold changes with lfcShrink(dds, coef = 2, type = "apeglm") improves ranking and removes noise from lowly expressed genes. Filter results using thresholds such as padj < 0.05 and |log2FoldChange| > 1. Visualize with plotMA for the mean-average plot, plotPCA on variance-stabilized data for sample clustering, and pheatmap on the top 50 most significant genes. The EnhancedVolcano package produces publication-ready volcano plots highlighting significantly upregulated and downregulated genes.
