RNA-Seq Differential Expression: DESeq2, edgeR, and limma-voom Frameworks

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-05-10

Key Takeaways

• DESeq2, edgeR, and limma-voom are primary R-based frameworks for RNA-Seq differential expression analysis, each employing distinct statistical models and normalization strategies to address count data properties like overdispersion and library size variation.
• DESeq2 utilizes a negative binomial generalized linear model (GLM) with empirical Bayes shrinkage for dispersion estimation, making it robust for small sample sizes (n=2-3) and effective at handling outliers through automatic count replacement.
• edgeR also employs a negative binomial GLM but offers flexible dispersion estimation (common, trended, tagwise) and a computationally efficient trimmed mean of M-values (TMM) normalization, with its quasi-likelihood (QL) pipeline providing reliable inference for larger sample sizes.
• limma-voom transforms RNA-Seq counts to log2-counts per million (logCPM) and applies linear modeling with precision weights derived from variance modeling, offering high flexibility for complex experimental designs and benefiting from empirical Bayes moderation with moderate to large sample sizes (n>5).
• The selection of a framework should consider experimental design, sample size, and expected biological variability; for veterinary transcriptomics, DESeq2 is often preferred for small or highly variable cohorts, while limma-voom excels in complex multi-factor designs.

Introduction

RNA sequencing (RNA-Seq) has become the standard method for quantifying transcript abundance in biological samples [<a href="#ref-1">1</a>]. In veterinary research, RNA-Seq is employed to study host transcriptional responses to viral and bacterial pathogens, evaluate vaccine immunogenicity, and investigate tissue-specific gene expression in livestock and companion animals [<a href="#ref-2">2</a>]. The core analytical task in most RNA-Seq experiments is differential expression (DE) analysis: identifying genes whose expression levels change significantly between experimental conditions [<a href="#ref-3">3</a>]. Three widely used frameworks for DE analysis are DESeq2, edgeR, and limma-voom [<a href="#ref-4">4</a>]. Each framework employs distinct statistical models and normalization strategies to handle the unique properties of count-based RNA-Seq data, including overdispersion, library size variation, and heteroscedasticity [<a href="#ref-5">5</a>]. This article provides an exhaustive technical comparison of these three frameworks, focusing on their underlying algorithms, assumptions, and practical considerations for veterinary transcriptomics.

Overview of RNA-Seq Differential Expression Workflow

A typical RNA-Seq DE analysis begins with raw sequencing reads in FASTQ format [<a href="#ref-6">6</a>]. Reads are aligned to a reference genome or transcriptome using splice-aware aligners, producing BAM files [<a href="#ref-7">7</a>]. Gene-level quantification is performed by counting reads overlapping annotated exons, yielding a count matrix with rows representing genes and columns representing samples [<a href="#ref-8">8</a>]. This count matrix is the input for DE analysis. The counts are discrete, non-negative integers that exhibit overdispersion relative to a Poisson distribution due to biological variability [<a href="#ref-9">9</a>]. All three frameworks address overdispersion through parametric or non-parametric modeling of the mean-variance relationship [<a href="#ref-10">10</a>].

DESeq2 Framework

DESeq2 is an R package that models count data using a generalized linear model (GLM) with a negative binomial distribution [<a href="#ref-11">11</a>]. The negative binomial distribution has two parameters: the mean (mu) and the dispersion (alpha), where variance = mu + alpha * mu^2 [<a href="#ref-12">12</a>]. DESeq2 estimates dispersion for each gene by pooling information across genes using an empirical Bayes shrinkage procedure [<a href="#ref-13">13</a>]. This shrinkage stabilizes dispersion estimates, particularly for genes with low counts, improving statistical power [<a href="#ref-14">14</a>].

Normalization in DESeq2 is performed using the median-of-ratios method [<a href="#ref-15">15</a>]. For each gene, a size factor is calculated as the median of the ratios of observed counts to a geometric mean across samples [<a href="#ref-16">16</a>]. These size factors account for differences in library size and RNA composition [<a href="#ref-17">17</a>]. The DESeq2 model then fits a GLM with a log link function, and Wald tests or likelihood ratio tests are used to assess significance [<a href="#ref-18">18</a>]. DESeq2 also provides independent filtering to increase power by removing lowly expressed genes prior to multiple testing correction [<a href="#ref-19">19</a>].

A key advantage of DESeq2 is its robust handling of outliers through Cook's distance and automatic replacement of outlier counts [<a href="#ref-20">20</a>]. The framework is particularly suitable for experiments with small sample sizes (n=2-3 per group) due to its shrinkage estimation [<a href="#ref-21">21</a>]. In veterinary contexts, DESeq2 has been applied to study host transcriptomic responses to avian influenza virus infection in chicken lung tissue and to evaluate immune gene expression in porcine alveolar macrophages after bacterial challenge [<a href="#ref-22">22</a>].

edgeR Framework

edgeR also employs a negative binomial GLM but uses a different approach for dispersion estimation and normalization [<a href="#ref-23">23</a>]. edgeR estimates a common dispersion across all genes, a trended dispersion as a function of gene abundance, and a tagwise dispersion for each gene using an empirical Bayes strategy [<a href="#ref-24">24</a>]. The tagwise dispersion is shrunk toward the trended dispersion, with the degree of shrinkage controlled by the prior degrees of freedom [<a href="#ref-25">25</a>].

Normalization in edgeR is performed using the trimmed mean of M-values (TMM) method [<a href="#ref-26">26</a>]. TMM computes a scaling factor for each sample by trimming the most extreme log-fold changes and absolute expression levels, then averaging the remaining values [<a href="#ref-27">27</a>]. This approach is robust to a small number of highly differentially expressed genes [<a href="#ref-28">28</a>]. edgeR offers exact tests for simple two-group comparisons and likelihood ratio tests or quasi-likelihood F-tests for more complex designs [<a href="#ref-29">29</a>]. The quasi-likelihood (QL) pipeline in edgeR models the gene-level variability with an additional dispersion parameter, providing more reliable inference when sample sizes are larger [<a href="#ref-30">30</a>].

edgeR is computationally efficient and scales well to large datasets [<a href="#ref-31">31</a>]. It has been used extensively in veterinary transcriptomics, for example in characterizing the mammary gland transcriptome of dairy cows during mastitis and in profiling the immune response of dogs with leishmaniasis [<a href="#ref-32">32</a>].

limma-voom Framework

limma-voom takes a different approach by transforming count data to log2-counts per million (logCPM) and then applying linear modeling methods originally developed for microarrays [<a href="#ref-33">33</a>]. The "voom" (variance modeling at the observational level) function estimates the mean-variance relationship of the logCPM values and assigns a precision weight to each observation [<a href="#ref-34">34</a>]. These weights are incorporated into the limma linear model, which uses empirical Bayes moderation to shrink the gene-wise variances toward a common value [<a href="#ref-35">35</a>].

Normalization in limma-voom typically involves TMM normalization applied to the counts before voom transformation, followed by quantile normalization of the logCPM values [<a href="#ref-36">36</a>]. The linear model framework allows flexible design matrices for complex experimental designs, including factorial designs and batch effects [<a href="#ref-37">37</a>]. Hypothesis testing is performed using moderated t-statistics and F-statistics [<a href="#ref-38">38</a>].

limma-voom is particularly powerful when sample sizes are moderate to large (n>5 per group) because the empirical Bayes variance moderation benefits from many degrees of freedom [<a href="#ref-39">39</a>]. It is also well suited for experiments with continuous covariates or time-course data [<a href="#ref-40">40</a>]. In veterinary research, limma-voom has been applied to study gene expression changes in feline oral squamous cell carcinoma and to analyze the transcriptomic response of equine synovial cells to inflammatory stimuli [<a href="#ref-41">41</a>].

Comparative Summary

The following table summarizes key features of the three frameworks.

Workflow Decision Diagram

The following Mermaid diagram illustrates a decision tree for selecting an appropriate DE framework based on experimental parameters.

flowchart TD
 A["RNA-Seq Count Matrix"] --> B{"Sample size per group?"}
 B -->|"n <= 3"| C["DESeq2"]
 B -->|"n = 4-5"| D{"Expected dispersion?"}
 D -->|"High biological variability"| C
 D -->|"Low variability"| E["edgeR QL or limma-voom"]
 B -->|"n >= 6"| F{"Design complexity?"}
 F -->|"Simple two-group"| G["edgeR or DESeq2"]
 F -->|"Multi-factor or continuous"| H["limma-voom"]
 C --> I["Perform DE analysis"]
 E --> I
 G --> I
 H --> I
 I --> J["Output: gene list, logFC, p-value, FDR"]

Practical Considerations for Veterinary Transcriptomics

Veterinary RNA-Seq experiments often involve outbred populations with high genetic variability, which can increase biological noise [<a href="#ref-42">42</a>]. In such cases, DESeq2's robust dispersion shrinkage and outlier handling are advantageous [<a href="#ref-43">43</a>]. For studies with limited sample sizes due to ethical or cost constraints, DESeq2 is generally recommended [<a href="#ref-44">44</a>]. When sample sizes are larger and the experimental design includes multiple factors (e.g., breed, treatment, time), limma-voom offers greater flexibility [<a href="#ref-45">45</a>]. edgeR provides a good balance between speed and accuracy and is suitable for routine DE analysis in well-replicated designs [<a href="#ref-46">46</a>].

All three frameworks require careful quality control of input counts, including removal of lowly expressed genes and assessment of sample clustering via principal component analysis [<a href="#ref-47">47</a>]. Batch effects should be modeled explicitly in the design matrix or corrected using methods such as ComBat-seq [<a href="#ref-48">48</a>]. For veterinary species with incomplete genome annotations, transcript-level quantification using pseudoalignment tools (e.g., Salmon, kallisto) can improve accuracy, and the resulting counts can be used with any of the three frameworks [<a href="#ref-49">49</a>].

Conclusion

DESeq2, edgeR, and limma-voom are robust and well-validated frameworks for RNA-Seq differential expression analysis. Each has distinct strengths: DESeq2 excels with small sample sizes and outlier sensitivity; edgeR offers computational efficiency and a mature quasi-likelihood pipeline; limma-voom provides flexibility for complex designs and leverages the power of empirical Bayes moderation. The choice of framework should be guided by experimental design, sample size, and biological variability. In veterinary transcriptomics, these tools enable the discovery of biomarkers, elucidation of host-pathogen interactions, and evaluation of therapeutic interventions, ultimately supporting animal health and disease management.

References

Related Clinical & Scientific Guides

• A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data
• Computational Immunology: Modeling the Immune System
• How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices