What is Dr. Zubair Khalid's research focus?

Dr. Zubair Khalid specializes in molecular virology, mRNA vaccine development, and computational biology, with a focus on avian pathogens like IBDV and Avian Reovirus.

Where is Dr. Zubair Khalid currently working?

Dr. Zubair Khalid is a Postdoctoral Research Associate at the University of Maryland (UMD), specifically within the Department of Animal and Avian Sciences.

Differential Gene Expression (DGE) Analysis: Methods and Pitfalls

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-08-11

Introduction to Differential Gene Expression (DGE) Analysis

What is DGE analysis?

Differential gene expression (DGE) analysis is the computational and statistical framework for identifying genes whose transcript abundance changes significantly between two or more experimental conditions. At its core, DGE analysis asks a deceptively simple question: for each gene, is the observed difference in expression between conditions greater than what would be expected from technical noise and biological variation alone?

The fundamental unit of measurement in DGE analysis is the count of sequencing reads (in RNA-seq) or the fluorescence intensity (in microarrays) that maps to a given gene. These raw measurements are transformed, normalized, and modeled to estimate the magnitude of expression change—typically reported as a fold change—and the statistical confidence in that change, reported as a p-value or adjusted p-value. The output is a ranked list of genes, usually filtered by thresholds for fold change and statistical significance, that constitutes the "differentially expressed genes" (DEGs) for downstream interpretation.

DGE analysis is distinct from simply comparing expression levels at a single time point. It requires replication, appropriate statistical modeling of count distributions, and correction for the multiple testing problem that arises when tens of thousands of genes are interrogated simultaneously. The field has matured substantially since the early days of microarray analysis, with the current standard being count-based models that account for the discrete nature of RNA-seq data.

Applications in biology and medicine

DGE analysis underpins a vast range of biological discovery. In developmental biology, it identifies the transcriptional programs that drive cell fate decisions, such as the switch from pluripotency to differentiation. In oncology, it reveals the genes and pathways that distinguish tumor subtypes, predicts drug response, and identifies therapeutic targets. In immunology, it characterizes the activation states of immune cells in response to pathogens or cytokines. In toxicology, it provides mechanistic insight into how environmental chemicals alter cellular physiology.

A typical application is comparing gene expression between wild-type and knockout mice to identify the transcriptional targets of a transcription factor. Another is profiling patient tumors before and after treatment to identify resistance mechanisms. In all these cases, the quality of the biological conclusions depends entirely on the rigor of the statistical analysis. A poorly designed DGE experiment with inadequate replication or inappropriate normalization will produce misleading gene lists, and the downstream functional interpretation—whether Gene Ontology pathway enrichment or candidate biomarker validation—will inherit those errors.

Experimental Design and Data Preprocessing

Biological and technical replicates

The distinction between biological and technical replicates is the single most important design decision in DGE analysis. Technical replicates are repeated measurements of the same biological sample—for example, sequencing the same RNA library twice. They capture only instrument and protocol variability. Biological replicates are independent samples from separate organisms or cell culture wells—for example, RNA extracted from three different mice in the same treatment group. They capture the full range of biological variability, including differences between individuals, stochastic gene expression noise, and environmental variation.

Biological replicates are essential for DGE analysis because the statistical models used (discussed below) estimate biological variance from the replicate measurements. Technical replicates inflate the apparent precision of the experiment without providing any information about the true population variance. As a rule, technical replicates are useful for troubleshooting protocols but should not be used as a substitute for biological replication. The minimum number of biological replicates for a well-powered DGE experiment is generally considered to be three per condition, though this is a floor, not a target. More replicates (six to twelve) are needed when the expected effect sizes are small, when the biological system is highly variable, or when the goal is to detect subtle expression changes.

RNA-seq preprocessing: QC, alignment, quantification

RNA-seq preprocessing follows a well-established pipeline, and each step affects the final DGE results.

Quality control (QC): Raw FASTQ files are assessed with tools like FastQC. Metrics include per-base sequence quality (Phred scores), GC content distribution, adapter contamination, and overrepresented sequences. Low-quality bases (Phred score below 20) are trimmed, and adapter sequences are removed using tools like Trimmomatic or cutadapt. The Phred score Q20 corresponds to a 1 in 100 error rate; Q30 corresponds to 1 in 1000.

Alignment: Cleaned reads are aligned to a reference genome or transcriptome. Spliced aligners such as STAR or HISAT2 are standard for mammalian genomes. STAR uses a seed-and-extend approach with an uncompressed suffix array index, achieving high speed at the cost of substantial memory usage (typically 30 GB for the human genome). Alignment rates above 80% for high-quality data are typical. Reads that map to multiple locations (multimappers) are usually discarded or assigned proportionally, as they are ambiguous.

Quantification: Gene-level counts are generated by counting the number of reads overlapping each gene's exons. Tools include featureCounts, HTSeq-count, and the pseudo-alignment tool Salmon. Salmon uses a different paradigm: it performs quasi-mapping to the transcriptome and estimates transcript abundances using an expectation-maximization algorithm, producing both transcript-level and gene-level estimates. The choice of quantification method affects the count matrix, but downstream DGE tools are generally robust to these differences if the counts are generated consistently.

Microarray preprocessing: normalization

Microarray analysis, while less common than RNA-seq, remains relevant for legacy datasets and certain clinical applications. The key preprocessing steps are background correction, normalization, and summarization. For Affymetrix arrays, the Robust Multi-array Average (RMA) method is the standard: it applies background correction, quantile normalization across arrays, and median-polish summarization of probe sets. For two-color arrays (e.g., Agilent), loess normalization is used to correct for intensity-dependent dye biases. The output is a continuous intensity value per gene, which is then log2-transformed for downstream analysis.

Statistical Models for DGE Analysis

Negative binomial distribution

RNA-seq count data are discrete and exhibit overdispersion: the variance across biological replicates is greater than the mean, violating the assumption of the Poisson distribution, where variance equals the mean. The negative binomial (NB) distribution models this overdispersion explicitly. For a gene with mean expression μ and dispersion parameter α, the variance is μ + αμ². The dispersion parameter captures the biological variability that cannot be explained by sampling noise alone.

The NB model is the foundation of the two most widely used DGE tools, DESeq2 and edgeR. Both estimate the dispersion for each gene, but they differ in how they handle the fact that dispersion is estimated from a small number of replicates. With only three replicates per condition, the dispersion estimate for a single gene is highly variable. Both tools address this by borrowing information across genes—a process called shrinkage or empirical Bayes moderation—where the gene-specific dispersion is pulled toward a global trend estimated from all genes with similar expression levels.

DESeq2 and edgeR

DESeq2 models count data with a generalized linear model (GLM) of the NB family. The key innovation is the use of a shrinkage estimator for both the dispersion and the log2 fold changes. Genes with low counts or high variance have their fold changes shrunk toward zero, which reduces false positives among weakly expressed genes. The DESeq2 model is:

log2(count_ij) = β0 + β1 * condition_j + offset(log(sizeFactor_j)) + error_ij

where the size factor accounts for library size differences (see Normalization below). The Wald test is used to assess whether the condition coefficient β1 is significantly different from zero.

edgeR also uses a NB GLM but employs a different dispersion estimation strategy. It uses the quantile-adjusted conditional maximum likelihood (qCML) method for experiments with a single factor, and the Cox-Reid profile-adjusted likelihood (CR) method for more complex designs. edgeR's empirical Bayes procedure, implemented in the estimateDisp function, moderates the dispersion estimates toward a common value or a trend. The test for differential expression is typically an exact test analogous to Fisher's exact test but adapted for the NB distribution, or a likelihood ratio test for multi-factor designs.

The practical differences between DESeq2 and edgeR are subtle. DESeq2 tends to be more conservative for genes with very low counts due to its fold change shrinkage. edgeR is slightly faster and may be more sensitive in some settings. Both tools are well-maintained and produce largely concordant results for well-powered experiments. A detailed comparison is available in the context of Differential Gene Expression Analysis Deseq2 and Differential Gene Expression Analysis in R.

limma-voom for RNA-seq

limma was originally developed for microarray analysis, where it models log2-transformed intensities with a linear model and uses empirical Bayes moderation of the gene-wise variances. The voom method extends limma to RNA-seq data by converting counts to log2-counts-per-million (log2-CPM) and estimating a mean-variance relationship. Each observation is assigned a precision weight based on this relationship, and the weighted linear model is fitted. The voom approach is computationally efficient and performs well, particularly for experiments with larger sample sizes. Its main advantage over DESeq2 and edgeR is speed and flexibility in handling complex experimental designs, including those with continuous covariates.

Normalization Methods and Their Impact

Library size normalization

The most basic normalization is for library size: the total number of reads sequenced for each sample. A sample with 30 million reads will have roughly twice the counts of a sample with 15 million reads, purely due to sequencing depth. The simplest correction is to divide each gene's count by the total library size and multiply by a constant (e.g., counts per million, CPM). However, this simple approach is inadequate when the RNA composition differs between samples.

Compositional biases

Compositional bias arises when a small number of highly expressed genes differ dramatically between conditions. For example, if a muscle-specific gene is massively upregulated in one condition, it consumes a larger fraction of the sequencing reads, causing all other genes to appear downregulated in that sample even if their absolute expression is unchanged. This is a fundamental problem with any method that normalizes by total read count.

The Trimmed Mean of M-values (TMM) method, implemented in edgeR, addresses this by computing a weighted trimmed mean of log2 fold changes between samples, excluding the most extreme genes. The scaling factor is calculated relative to a reference sample and applied to each library. DESeq2 uses the median-of-ratios method: for each gene, the geometric mean across samples is computed, and the ratio of each sample's count to this geometric mean is calculated. The median of these ratios is the size factor for that sample. Both methods are robust to compositional bias because they are based on the majority of genes, which are assumed to be non-differentially expressed.

Choosing a normalization method

The choice of normalization method can materially affect DGE results, particularly for genes with moderate expression changes. The table below summarizes the main options:

Method	Tool	Basis	Best for
TMM	edgeR	Trimmed mean of M-values	RNA-seq with compositional bias
RLE (median-of-ratios)	DESeq2	Median ratio to geometric mean	RNA-seq with compositional bias
CPM/TPM	Various	Counts per million / transcripts per million	Exploratory analysis, not DGE
Quantile	Microarray	Equalize distributions across samples	Microarray data
voom (with limma)	limma	Log2-CPM with precision weights	RNA-seq with complex designs

A critical point is that normalization should be performed within the DGE analysis framework, not as a separate step that produces "normalized counts" for downstream tools. The normalization factors are integrated into the statistical model, and the uncertainty in these factors is accounted for. Using pre-normalized data (e.g., TPM values) in DESeq2 or edgeR is incorrect because these tools expect raw integer counts.

Multiple Testing Correction and False Discovery Rate

Why correction is needed

A typical RNA-seq experiment measures expression for approximately 20,000 protein-coding genes. If we test each gene for differential expression at a significance threshold of α = 0.05, we expect 1,000 false positives by chance alone (20,000 × 0.05). Without correction, the list of "significant" genes would be dominated by noise. Multiple testing correction is therefore not optional; it is a mandatory step in DGE analysis.

FDR vs. family-wise error rate

The family-wise error rate (FWER) controls the probability of making at least one false positive across all tests. The Bonferroni correction, which multiplies each p-value by the number of tests, is the simplest FWER method but is extremely conservative for genome-wide data. With 20,000 genes, a gene must have a raw p-value below 2.5 × 10⁻⁶ to be significant at FWER = 0.05.

The false discovery rate (FDR), introduced by Benjamini and Hochberg in 1995, controls the expected proportion of false positives among the genes declared significant. If we declare 100 genes significant at FDR = 0.05, we expect 5 of them to be false positives. The Benjamini-Hochberg (BH) procedure ranks all p-values, and a gene with rank i is significant if its p-value is less than (i/m) × α, where m is the total number of tests. The BH procedure is less conservative than Bonferroni and is the standard for DGE analysis. The adjusted p-value (often called q-value) is the smallest FDR at which the gene would be declared significant.

Interpreting adjusted p-values

An adjusted p-value of 0.01 means that if you declare this gene significant, you can expect 1% of all genes declared significant at this threshold to be false positives. It does not mean there is a 1% chance that this specific gene is a false positive. This distinction is subtle but important. The FDR is a property of the set of significant genes, not of individual genes. In practice, researchers typically use an FDR threshold of 0.05 or 0.01, and they should report the adjusted p-values, not raw p-values, in their results.

Biological Interpretation of DGE Results

Gene ontology and pathway enrichment

Once a list of DEGs is obtained, the next step is to determine whether these genes share biological functions. Gene Ontology (GO) enrichment analysis tests whether genes annotated with a particular GO term (e.g., "DNA repair" or "mitochondrial translation") are overrepresented in the DEG list compared to the background set of all measured genes. The standard statistical test is the hypergeometric distribution or Fisher's exact test, applied to a 2×2 contingency table for each GO term. As with DGE analysis itself, multiple testing correction is essential because thousands of GO terms are tested simultaneously.

Pathway enrichment analysis extends this concept to curated pathway databases such as KEGG, Reactome, and WikiPathways. Tools like clusterProfiler, Enrichr, and DAVID implement these analyses. A more sophisticated approach is gene set enrichment analysis (GSEA), which does not require a thresholded DEG list. Instead, GSEA ranks all genes by their differential expression statistic and tests whether genes in a given pathway are enriched at the top or bottom of the ranking. This approach is more sensitive for detecting coordinated changes in pathways where individual genes may not reach significance. For practical guidance on running these analyses, see Gene Ontology Analysis Online and Gene Ontology Analysis Tool.

Visualizing DGE results

The two most common visualizations in DGE analysis are the volcano plot and the heatmap.

A volcano plot displays the log2 fold change on the x-axis and the negative log10 of the adjusted p-value on the y-axis. Genes with large fold changes and high significance appear in the upper left and upper right corners. Horizontal and vertical lines indicate the chosen thresholds (e.g., |log2FC| > 1 and adjusted p < 0.05). The volcano plot provides a rapid overview of the distribution of effects and the relationship between effect size and significance.

A heatmap shows the expression values of the top DEGs across all samples, typically with rows representing genes and columns representing samples. The values are usually z-scored (subtract the mean across samples, divide by the standard deviation) so that the pattern of up- and downregulation is visible. Hierarchical clustering of rows and columns reveals co-expressed gene modules and sample groupings. Heatmaps are useful for quality control (checking that replicates cluster together) and for identifying expression patterns that distinguish conditions.

Caveats in interpretation

DGE analysis identifies statistical associations, not causal relationships. A gene that is differentially expressed between conditions may be a driver of the phenotype, a downstream consequence, or a bystander. Functional validation—through knockdown, overexpression, or rescue experiments—is required to establish causality. Furthermore, mRNA levels do not always correlate with protein levels due to post-transcriptional regulation, including microRNA-mediated degradation and differential translation rates. The relationship between DNA methylation and gene expression is similarly indirect: __MASK_9 is a common pattern at promoters, but MASK_10__ can occur at gene bodies, illustrating that the direction of the effect depends on genomic context.

Common Pitfalls and Best Practices in DGE Analysis

Replicate number and power

The most common failure in DGE analysis is insufficient biological replication. With two replicates per condition, the dispersion cannot be estimated reliably, and the statistical test has almost no power to detect anything but the most extreme fold changes. Three replicates is the accepted minimum, but this provides limited power for detecting small effects. A power analysis should be performed before the experiment, using pilot data or published datasets to estimate the expected dispersion and effect sizes. As a rough guide, detecting a 2-fold change at FDR = 0.05 with 80% power typically requires 5–10 biological replicates per condition for mammalian systems, depending on the variability of the tissue or cell type.

Handling batch effects

Batch effects are systematic technical variations that affect all samples processed together—for example, samples sequenced on different days, in different flow cells, or by different technicians. If batches are confounded with the biological condition of interest (e.g., all treated samples in batch 1, all controls in batch 2), it is impossible to distinguish biological from technical effects. The solution is experimental: randomize samples across batches so that each batch contains a mix of conditions. If batch effects are known but not confounded, they can be included as covariates in the statistical model (e.g., in DESeq2's design formula: __MASK_2__). If batch effects are unknown, tools like RUVseq or SVA can estimate surrogate variables from the data, but these methods require careful validation.

Over-interpreting small fold changes

A gene with a log2 fold change of 0.3 (a 1.23-fold change) can be statistically significant with sufficient replication, but its biological relevance is questionable. Small fold changes may reflect true biology, but they are more susceptible to technical artifacts and are harder to validate experimentally. A common practice is to apply a fold change threshold (e.g., |log2FC| > 1, corresponding to a 2-fold change) in addition to the FDR threshold. However, this practice is debated: some argue that biologically important genes can have small fold changes, and that filtering by fold change discards information. The safest approach is to report all genes that pass the statistical threshold and to discuss the fold change distribution explicitly.

Misapplying statistical models

A frequent error is using a tool designed for one data type on another. For example, running DESeq2 on TPM-normalized data, or running limma on raw counts without voom transformation, produces invalid results. Another error is using a paired test when the data are unpaired, or vice versa. The design formula must reflect the actual experimental structure: if samples are paired (e.g., tumor and normal tissue from the same patient), the design should include the patient as a blocking factor. Ignoring pairing reduces power; incorrectly assuming pairing when samples are independent inflates false positives.

Frequently Asked Questions

What is the difference between DESeq2 and edgeR?

Both DESeq2 and edgeR use negative binomial models to test for differential expression in RNA-seq count data. The main differences are in dispersion estimation and fold change handling. DESeq2 uses a shrinkage estimator for both dispersion and log2 fold changes, which makes it more conservative for genes with low counts. edgeR uses quantile-adjusted conditional maximum likelihood for dispersion estimation and offers both exact tests and likelihood ratio tests. In practice, the results are largely concordant for well-powered experiments. The choice between them is often a matter of preference, though DESeq2's fold change shrinkage can be advantageous for ranking genes.

How many biological replicates are needed for DGE analysis?

Three biological replicates per condition is the absolute minimum, but this provides limited statistical power. For most experiments, 5–10 replicates per condition are recommended. The exact number depends on the biological variability of the system, the magnitude of the expected effect, and the desired statistical power. A power analysis using pilot data is the best way to determine the required sample size.

What is a volcano plot in DGE analysis?

A volcano plot is a scatter plot that displays the log2 fold change on the x-axis and the negative log10 of the adjusted p-value on the y-axis for every gene. Genes in the upper corners are both highly significant and have large fold changes. The plot provides a quick visual summary of the DGE results and helps identify genes that pass both statistical and effect size thresholds.

Why do I get different results from different DGE tools?

Different DGE tools make different statistical assumptions and use different normalization and dispersion estimation methods. DESeq2 shrinks fold changes toward zero, edgeR uses a different dispersion prior, and limma-voom models the mean-variance relationship differently. These differences are most pronounced for genes with low counts and for experiments with small sample sizes. For well-powered experiments with adequate replication, the overlap between tools is typically high. If results diverge substantially, it is worth investigating whether the data meet the assumptions of the tools being used.

What is the false discovery rate (FDR) in DGE analysis?

The FDR is the expected proportion of false positives among the genes declared significant. If you declare 100 genes significant at FDR = 0.05, you expect 5 of them to be false positives. The Benjamini-Hochberg procedure is the standard method for controlling the FDR in DGE analysis. The adjusted p-value (q-value) for each gene is the minimum FDR at which that gene would be declared significant.

How do I handle batch effects in DGE analysis?

The best approach is to prevent batch effects through experimental design: randomize samples across batches and include batch as a covariate in the statistical model. If batch effects are present but not confounded with the condition of interest, adding + batch to the design formula in DESeq2 or edgeR will account for them. If batch effects are unknown, methods like RUVseq or surrogate variable analysis (SVA) can estimate them from the data, but these require careful validation.

What is the difference between fold change and statistical significance?

Fold change measures the magnitude of the expression difference between conditions (e.g., a 2-fold change means expression is doubled). Statistical significance measures the confidence that the observed difference is not due to chance. A gene can have a large fold change but be statistically insignificant (if variance is high or replication is low), or a small fold change but be highly significant (if variance is low and replication is high). Both metrics are important: fold change indicates biological relevance, while significance indicates reliability.

Key Takeaways

DGE analysis identifies genes whose expression changes between conditions using statistical models that account for biological variability and technical noise.
Biological replicates are essential; three per condition is the minimum, but 5–10 are typically needed for adequate power.
Negative binomial models (DESeq2, edgeR) and limma-voom are the standard statistical frameworks for RNA-seq DGE analysis.
Normalization (TMM, median-of-ratios) corrects for library size and compositional biases and must be integrated into the statistical model.
Multiple testing correction using the Benjamini-Hochberg FDR procedure is mandatory when testing thousands of genes simultaneously.
Batch effects must be addressed through experimental design (randomization) and statistical modeling (including batch as a covariate).
DGE results identify statistical associations, not causal relationships; functional validation and pathway analysis are required for biological interpretation.