limma-voom for RNA-seq: How to Apply Linear Modeling to Count Data for Differential Expression

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

Key Takeaways

• The limma-voom method transforms RNA-seq count data into a format suitable for linear modeling by converting raw counts to log2-counts per million (log2-CPM) and estimating a mean-variance relationship to generate precision weights.
• These precision weights, inversely proportional to the predicted variance from the mean-variance trend, are incorporated into the limma empirical Bayes pipeline, allowing for robust variance moderation and differential expression analysis even with small biological replicates.
• The workflow involves importing count matrices, filtering low-expression genes (e.g., minimum 10 counts in smallest group size), applying TMM normalization for RNA composition differences, and then executing the voom transformation and limma linear modeling.
• Contrasts are defined to specify comparisons of interest, followed by empirical Bayes moderation (eBayes function) to compute moderated t-statistics, p-values, and adjusted p-values (e.g., Benjamini-Hochberg) for controlling the false discovery rate.
• While effective for designed experiments with moderate sample sizes, limma-voom can be sensitive to outliers in small-sample cases, and for large population studies with heterogeneous samples, robust methods like the Wilcoxon rank-sum test may offer superior false discovery rate control.
• Applications span diverse fields including single-cell pseudo-bulk analysis, neurological studies, cancer biomarker discovery, and respiratory disease research, demonstrating its versatility in identifying differentially expressed genes and gene sets.

RNA sequencing has become the standard approach for measuring gene expression across the transcriptome, and differential expression analysis is the statistical core of most RNA-seq studies. Analysts who learned differential expression with microarray platforms often find that the move to count-based RNA-seq data requires a different statistical framework. The limma-voom method bridges this gap by transforming RNA-seq count data into a form that can be analyzed with the linear modeling tools originally developed for microarrays. This article explains the voom transformation, the precision weights it generates, and the complete limma-voom workflow in R, from count matrix to ranked gene lists.

The limma-voom approach works by estimating the mean-variance relationship of log-counts, generating a precision weight for each observation, and entering those weights into the limma empirical Bayes analysis pipeline. This opens access for RNA-seq analysts to a large body of methodology developed for microarrays, including linear modeling, contrasts, and gene set testing. Simulation studies show that voom performs as well as or better than count-based RNA-seq methods even when the data are generated according to the assumptions of those earlier methods. Two case studies illustrate the use of linear modeling and gene set testing methods with voom.

The Statistical Problem with RNA-seq Count Data

RNA-seq experiments produce integer counts of sequencing reads that map to genomic features such as genes, transcripts, or exons. These counts are not normally distributed. For a given gene in a given sample, the count reflects both the true abundance of the transcript and the sequencing depth of that sample. A gene with low expression may produce zero counts in some samples and a handful of counts in others. A highly expressed gene may produce thousands of counts. This count nature creates a mean-variance relationship that differs fundamentally from the constant variance assumption of ordinary linear regression.

Microarray intensity data, by contrast, are continuous measurements with approximately normal errors after appropriate transformation. The limma package was developed for this setting and uses empirical Bayes moderation to borrow information across genes, stabilizing variance estimates even when the number of biological replicates is small. The challenge for RNA-seq analysis is that raw counts violate the assumptions of this linear modeling framework.

Several strategies have emerged for handling count data. Some methods model counts directly using the negative binomial distribution, as implemented in edgeR and DESeq2. Other methods transform the counts to make them approximately normal and then apply linear modeling tools. The voom method belongs to this second category. It converts counts to log-counts per million, estimates the mean-variance relationship, and assigns a precision weight to each observation. These weights are then used in the limma pipeline.

The choice of method matters for the final results. An extended review of RNA-seq differential expression methods evaluated six methods of mapping reads and nine methods of differential expression analysis using real RNA-seq data with qRT-PCR data as the reference standard. The results indicated that mapping methods have minimal impact on the final differential expression analysis when the data have an annotated reference genome. Among the differential expression identification methods, limma+voom, NOISeq, and DESeq2 produced the most consistent results. The review also found that a consensus among five differential expression identification methods guarantees a list of differentially expressed genes with great accuracy, indicating that combining different methods can produce more suitable results.

How the voom Transformation Works

The voom method proceeds through a series of steps that convert raw counts into precision-weighted log-expression values suitable for linear modeling.

From Counts to Log-Counts per Million

The first step is to convert raw counts to counts per million (CPM). For each gene in each sample, the raw count is divided by the total number of mapped reads in that sample and multiplied by one million. This normalization accounts for differences in sequencing depth between samples. The CPM values are then transformed to the log2 scale, producing log2-CPM values.

The choice of log base matters for interpretation. Log2-CPM values have a natural interpretation: a one-unit difference corresponds to a twofold change in expression. Most RNA-seq analyses use log2 transformations for this reason.

Estimating the Mean-Variance Relationship

The core innovation of voom is the explicit modeling of the mean-variance relationship. For genes with low expression, the variance of the log2-CPM values is high because small count differences produce large relative changes. For genes with high expression, the variance is lower and more stable. The voom method estimates this relationship by fitting a trend line through the variance of each gene as a function of its mean expression.

This mean-variance trend is the key to the precision weights. The trend captures the fact that low-expression genes are noisier and should receive less weight in the linear model, while high-expression genes are more reliable and should receive more weight.

Generating Precision Weights

For each observation, voom computes a precision weight that is inversely proportional to the predicted variance from the mean-variance trend. Observations with high predicted variance receive low weights. Observations with low predicted variance receive high weights. These weights are then entered into the limma linear modeling framework, which uses weighted least squares to fit the model.

The precision weights serve a role analogous to the variance-stabilizing transformation used in other methods. instead of transforming the data to make the variance constant, voom keeps the data on the log2-CPM scale and uses weights to account for the heteroscedasticity. This preserves the interpretability of the log2 scale while providing the statistical benefits of variance stabilization.

Empirical Bayes Moderation

After the weighted linear model is fitted, limma applies empirical Bayes moderation to the gene-wise variance estimates. This step borrows information across all genes to stabilize the variance estimates for individual genes. The moderation is particularly valuable when the number of biological replicates is small, because individual gene variance estimates based on only three or four replicates are noisy. The empirical Bayes approach shrinks these estimates toward a common value, reducing the number of false positives driven by genes with artificially low variance.

The output of the limma-voom pipeline includes moderated t-statistics, p-values, and adjusted p-values for each gene and each contrast of interest. The adjusted p-values control the false discovery rate across the set of tested genes.

The Complete limma-voom Workflow in R

The limma-voom workflow follows a standard sequence of steps from count matrix to differential expression results. The workflow is implemented in R using the Bioconductor packages edgeR and limma. The edgeR package is used for data import, filtering, and normalization, while the limma package provides the voom transformation and linear modeling.

Step 1: Import the Count Matrix

The starting point is a matrix of raw counts with genes in rows and samples in columns. This matrix can be generated by any read quantification tool, including featureCounts, HTSeq, or the quantification modules of alignment-based pipelines. The count matrix should contain integer counts for each gene in each sample.

The count matrix is imported into R and used to create a DGEList object from the edgeR package. The DGEList object stores the counts along with sample information such as group labels, batch identifiers, and other covariates.

Step 2: Filter Low-Expression Genes

Genes with very low counts across all samples provide little information and can interfere with the mean-variance estimation. The standard practice is to filter out genes that do not have a minimum level of expression in a minimum number of samples. A common filter keeps genes with at least 10 counts in at least the number of samples corresponding to the smallest group size.

Filtering reduces the number of tests performed in the multiple testing correction, which increases statistical power. It also removes genes whose counts are dominated by technical noise.

Step 3: Normalize the Counts

The voom method assumes that the library sizes, meaning the total number of mapped reads per sample, are accounted for in the CPM conversion. However, differences in RNA composition between samples can bias the analysis. The trimmed mean of M-values (TMM) normalization method from edgeR adjusts for these composition differences.

TMM normalization computes a scaling factor for each sample that corrects for the fact that a few highly expressed genes in one sample can suppress the apparent expression of all other genes in that sample. The scaling factors are applied to the library sizes before the CPM conversion.

Step 4: Estimate the Mean-Variance Relationship and Generate Weights

The voom function in limma takes the DGEList object, the design matrix, and optionally the normalization factors, and returns an EList object containing the log2-CPM values, the precision weights, and the design matrix. The function fits a trend to the mean-variance relationship and computes the precision weights for each observation.

The voom function also produces a plot showing the mean-variance trend. This plot is a useful diagnostic. The trend should be smooth and decreasing, with high variance at low expression levels and lower variance at high expression levels. An erratic or flat trend may indicate problems with the data or the normalization.

Step 5: Fit the Linear Model

The linear model is specified with a design matrix that encodes the experimental groups and any covariates to be adjusted for. The design matrix can be created with the model.matrix function in R. For a simple two-group comparison, the design matrix has one column for the intercept and one column for the group difference. For more complex designs, the design matrix can include multiple factors, interactions, and continuous covariates.

The lmFit function in limma fits the linear model for each gene using the precision weights from voom. The result is a fitted model object with coefficients and standard errors for each gene and each coefficient in the design matrix.

Step 6: Define Contrasts and Test

For designs with more than two groups, contrasts define the specific comparisons of interest. The contrasts.fit function in limma applies the contrasts to the fitted model. The eBayes function then performs the empirical Bayes moderation and computes moderated t-statistics, p-values, and adjusted p-values.

The topTable function extracts the results for a given contrast, sorted by p-value or by the log fold change. The output includes the log2 fold change, the average expression, the moderated t-statistic, the raw p-value, and the adjusted p-value.

Step 7: Explore and Interpret the Results

The results of the limma-voom analysis are a ranked list of genes for each contrast. The ranking can be based on the adjusted p-value, the log fold change, or a combination of both. The Glimma package provides interactive visualization of the results, allowing individual samples and genes to be examined in detail.

Gene set testing can be performed with the limma functions camera or roast, which test whether predefined gene sets are enriched for differential expression. These tests use the same linear modeling framework and precision weights as the gene-level analysis.

At a Glance: limma-voom Decision Table

Design Considerations for limma-voom Analysis

The performance of limma-voom depends on the experimental design and the characteristics of the data. Several factors should be considered before running the analysis.

Sample Size and Replication

The number of biological replicates is the most important factor in the power of a differential expression analysis. The empirical Bayes moderation in limma borrows information across genes, which helps when the number of replicates is small. However, the precision of the variance estimates and the reliability of the mean-variance trend both improve with more replicates.

A systematic benchmarking study of differential expression methods for circular RNAs compared 38 analysis pipelines and found that almost all methods improved when the number of replicates increased. The study also found that widely used methods such as DESeq2, edgeR, and limma-voom gave scarce results or unreliable predictions on data sets of typical size. Limma-voom achieved the most consistent performance throughout the different benchmark data sets and, along with SAMseq, reasonably balanced false discovery rate and recall rate.

Outliers and Robustness

The standard limma-voom approach is sensitive to outliers, particularly when the sample size is small. A study that robustified the voom approach using the minimum beta-divergence method found that limma with voom transformation is sensitive to outliers for small-sample cases. The robustified method improved performance over the competing methods in the presence of outliers.

For population-level RNA-seq studies with large sample sizes, a permutation analysis found that DESeq2 and edgeR had unexpectedly high false discovery rates, sometimes exceeding 20% when the target false discovery rate was 5%. The Wilcoxon rank-sum test was the most robust method in this setting. The authors recommended the Wilcoxon rank-sum test for population-level studies with large sample sizes. A subsequent response confirmed that the Wilcoxon rank-sum test remains the most robust method compared to DESeq2, edgeR, limma-voom, dearseq, and NOISeq in two-condition analyses after considering normalization and winsorization.

These findings indicate that limma-voom is well suited to designed experiments with moderate sample sizes and controlled conditions, but may not be the best choice for large population studies with heterogeneous samples and outliers.

Paired and Blocked Designs

The linear modeling framework in limma naturally handles paired designs, blocked designs, and designs with technical covariates. A paired design, where each subject contributes a sample from two conditions, is encoded with a design matrix that includes a term for each subject. This approach accounts for the correlation between samples from the same subject and increases the power to detect differential expression.

An analysis of RNA-seq data from the Airway dataset demonstrated that correctly identifying paired experimental designs increased differential expression detection by 21 to 47%. The same analysis found that detecting technical covariates such as library type added 4 to 30% more differentially expressed genes. These results highlight the importance of specifying the design correctly.

Multi-Group and Factorial Designs

The contrast-based approach in limma supports comparisons among multiple groups and factorial designs. For a study with three or more groups, the design matrix includes a term for each group, and contrasts define the specific comparisons of interest. For a factorial design with two factors, the design matrix includes main effects and interactions, and contrasts can test for main effects, simple effects, or interaction effects.

The flexibility of the linear modeling framework is one of the main advantages of limma-voom over simpler methods that only support two-group comparisons.

Practical Implementation Steps

The following steps outline a complete limma-voom analysis from count matrix to results. The steps assume a basic familiarity with R and the Bioconductor package system.

Step 1: Install and Load Required Packages

The analysis requires the edgeR and limma packages from Bioconductor. The Glimma package is optional but useful for interactive visualization. The packages are installed with the BiocManager package and loaded with the library function.

Step 2: Prepare the Count Matrix and Sample Information

The count matrix should have genes in rows and samples in columns. The sample information should include the group labels and any covariates to be included in the model. The sample information is typically stored in a data frame with one row per sample.

Step 3: Create the DGEList Object

The DGEList function from edgeR creates the data object that stores the counts and sample information. The group labels are assigned to the samples at this stage.

Step 4: Filter and Normalize

Low-expression genes are filtered using the filterByExpr function or a manual filter based on counts per million. The calcNormFactors function computes the TMM normalization factors.

Step 5: Run voom and Fit the Model

The voom function transforms the counts and generates the precision weights. The lmFit function fits the linear model. The design matrix is specified at this stage.

Step 6: Apply Contrasts and Empirical Bayes Moderation

The contrasts.fit function applies the contrasts of interest. The eBayes function performs the empirical Bayes moderation and computes the test statistics.

Step 7: Extract and Save Results

The topTable function extracts the results for each contrast. The results are saved to a file for downstream analysis and reporting.

Records and Measurements for Quality Control

Quality control is an essential part of the limma-voom workflow. Several measurements and records should be maintained throughout the analysis.

Library Size and Mapping Statistics

The total number of mapped reads per sample is the library size. Samples with very low library sizes may need to be excluded or down-weighted. The mapping statistics from the alignment step, such as the percentage of reads mapped to genes, provide an indication of data quality.

Mean-Variance Trend Plot

The voom function produces a plot of the mean-variance trend. This plot should show a smooth decreasing trend. An irregular trend may indicate problems with normalization, filtering, or the presence of outliers.

Multidimensional Scaling Plot

A multidimensional scaling plot of the log2-CPM values shows the relationships between samples. Samples from the same group should cluster together. Samples that are distant from their group may be outliers or may indicate a batch effect.

Number of Differentially Expressed Genes

The number of differentially expressed genes at a given false discovery rate threshold provides a summary of the analysis. An unexpectedly high or low number may indicate problems with the data or the model specification.

Session Information

The R session information records the versions of R and all packages used in the analysis. This record is essential for reproducibility.

Common Failure Patterns in limma-voom Analysis

Several recurring problems can compromise a limma-voom analysis. Recognizing these patterns helps analysts diagnose and correct issues.

Failure to Filter Low-Expression Genes

Including genes with very low counts in the analysis increases the multiple testing burden and can distort the mean-variance trend. The voom transformation assumes that the mean-variance relationship is smooth, and genes with all-zero counts or near-zero counts violate this assumption.

Incorrect Design Matrix

The design matrix must correctly encode the experimental structure. A common error is treating a paired design as an unpaired design, which reduces power and can produce biased results. Another common error is including the group variable as a continuous covariate instead of a factor.

Ignoring Composition Normalization

The TMM normalization is important for samples with different RNA compositions. Without normalization, samples with a few highly expressed genes can appear to have lower expression for all other genes. The voom function can use the normalization factors from edgeR, but only if they are computed and passed correctly.

Overlooking Outliers

The standard limma-voom approach is sensitive to outliers, particularly with small sample sizes. Outliers can be detected with multidimensional scaling plots or with diagnostic plots of the fitted model. Robust alternatives exist but are not part of the standard limma-voom workflow.

Misinterpreting Adjusted P-Values

The adjusted p-values from the Benjamini-Hochberg procedure control the false discovery rate, not the family-wise error rate. An adjusted p-value of 0.05 means that approximately 5% of the genes called significant are expected to be false positives. This interpretation differs from the traditional p-value interpretation.

Applications of limma-voom in Published Research

The limma-voom method has been applied across a wide range of biological and clinical research areas. The published applications illustrate the versatility of the method and the types of questions it can address.

Single-Cell and Tissue-Level Studies

A study of heart failure with preserved ejection fraction used single-nucleus RNA sequencing to characterize cell-specific gene expression patterns in human heart tissue. The study used limma-voom for differential expression analysis by cell type. After quality control, the study recovered 48,886 nuclei and identified 14 cell types. Cardiomyocytes and fibroblasts showed the most differentially expressed genes, while endothelial cells, pericytes, and macrophages had fewer. The study demonstrated the use of limma-voom in a pseudo-bulk analysis of single-cell data.

Brain and Neurological Studies

A study of alcohol use disorder performed RNA-seq analysis of long noncoding RNAs and messenger RNAs in 192 postmortem tissue samples from eight brain regions. Applying the limma-voom method, the study detected 57 long noncoding RNAs and 51 messenger RNAs with significant differential expression across at least one brain region. The study used weighted gene co-expression network analysis to identify modules associated with alcohol use disorder and constructed co-expression networks for each brain region.

Cancer Biomarker Discovery

A study of four hematologic malignancies mined public databases for RNA-seq data and used a standard analysis pipeline with DESeq2, limma+voom, and edgeR to identify differentially expressed genes. The study identified 2,136 highly ranked differentially expressed genes and performed gene enrichment, protein-protein interaction, survival, and tumor immune infiltration analyses. The study found that high expression of several genes and low expression of one gene were correlated with poor prognosis.

Respiratory Disease Studies

A study of COVID-19 used RNA-seq gene count data from healthy and treated cases and applied DESeq2, limma trend, and limma voom methods to identify differentially expressed genes. The study extracted 67 potential biomarkers and found that a support vector machine model achieved high classification accuracy for moderate differentially expressed genes. A related study compared the molecular signatures of COVID-19 and idiopathic pulmonary fibrosis using DESeq2, limma voom, and limma trend.

Equine Pregnancy Research

A study of early equine pregnancy analyzed conceptus and endometrial RNA-seq data using TMM normalization and voom-limma modeling. Differential expression was defined by an adjusted p-value threshold, and the B-statistic was used to prioritize high-confidence candidates. The study identified upregulated genes in the conceptus including steroidogenic enzymes, extracellular matrix components, protease regulators, and selective transporters.

Vaping and Smoking Studies

A study of oral epithelial cells from e-cigarette users, cigarette smokers, and non-users used covariate-adjusted limma-voom modeling with false discovery rate control. The study evaluated the extent to which exposure-specific dose metrics explained transcriptional changes and found that both vaping and smoking were associated with transcriptomic dysregulation relative to non-users.

Radiation Injury Research

A study of radiation-induced bone marrow injury used RNA-seq data from mouse hematopoietic stem cells and identified differentially expressed genes using two independent statistical frameworks. The study found that differentially expressed genes were primarily enriched in processes related to cytokine signaling, hematopoietic lineage regulation, immune response, and extracellular matrix remodeling.

Comparison with Other Differential Expression Methods

The choice of differential expression method can affect the final results. Several studies have compared limma-voom with other methods under different conditions.

Consistency with DESeq2 and edgeR

The extended review of RNA-seq differential expression methods found that limma+voom, NOISeq, and DESeq2 had the most consistent results when compared with qRT-PCR data. The review also found that a consensus among five methods produced a list of differentially expressed genes with great accuracy. This finding suggests that combining methods can improve the reliability of the results.

Performance on Circular RNA Data

The benchmarking study of circular RNA differential expression methods found that limma-voom achieved the most consistent performance throughout different benchmark data sets. The study also found that almost all methods improved when the number of replicates increased. The study concluded that circular RNA expression studies require careful design, choice of method, and method configuration.

Performance on Population Samples

The permutation analysis of human population RNA-seq samples found that limma-voom, along with DESeq2 and edgeR, failed to control the false discovery rate in some settings. The Wilcoxon rank-sum test was the most robust method. This finding is relevant for studies with large sample sizes and heterogeneous populations.

Cross-Platform Data Integration

A study of combined RNA-seq experiments from different platforms proposed a procedure that combines popular methods with a data-driven simulation to calibrate the levels of significance for multiple testing. The procedure was demonstrated with edgeR, DESeq2, and limma+voom. The study noted that combining data from different experimental platforms creates a complex data distribution that causes problems in statistical modeling and multiple testing.

Limitations and Interpretation Constraints

The limma-voom method has several limitations that analysts should understand before applying it.

Sensitivity to Outliers

The standard limma-voom approach is sensitive to outliers, particularly with small sample sizes. The precision weights are computed from the mean-variance trend, and outliers can distort the trend and the resulting weights. Robustified versions of voom exist but are not part of the standard workflow.

Assumptions of the Linear Model

The linear modeling framework assumes that the log2-CPM values are approximately normally distributed after weighting. This assumption is reasonable for moderately and highly expressed genes but may be questionable for genes with very low expression. The filtering step removes the lowest-expression genes, which mitigates this concern.

Interpretation of Fold Changes

The log2 fold change from the linear model is an estimate of the difference in log2-CPM between conditions. This estimate is not the same as the fold change in raw counts, because the TMM normalization and the precision weights affect the estimate. The log2 fold change should be interpreted as the change in normalized expression, not the change in raw counts.

Multiple Testing Correction

The adjusted p-values control the false discovery rate across all tested genes. The number of tests is the number of genes that pass the filtering step. Increasing the stringency of the filtering reduces the number of tests and increases the power to detect differential expression among the remaining genes.

Generalizability of Results

The results of a limma-voom analysis are specific to the samples and conditions in the study. The differentially expressed genes identified in one study may not replicate in another study with different samples, conditions, or platforms. Validation with an independent method such as qRT-PCR is recommended for key findings.

Reproducibility and Reporting Standards

Reproducibility is a central concern in bioinformatics analysis. The limma-voom workflow should be documented and reported in a way that allows others to reproduce the results.

Version Control and Session Information

The R session information records the versions of R and all packages used in the analysis. This information should be saved and reported. The sessionInfo function in R produces this record.

Pipeline Documentation

The analysis steps should be documented in a script or notebook that can be rerun. The script should include the filtering thresholds, the normalization method, the design matrix specification, and the contrasts tested. The Bioconductor project provides documentation and workflows for reproducible genomic analysis.

Containerization and Workflow Tools

Containerization tools such as Docker and Singularity can package the analysis environment, including the operating system, R version, and package versions. Workflow tools such as nf-core provide standardized pipelines for RNA-seq analysis that include quality control, alignment, quantification, and differential expression. The nf-core documentation describes the usage and configuration of these pipelines.

Training and Education

The EMBL-EBI Training program provides learning pathways for bioinformatics, including data-resource training and practical analysis education. The Galaxy Training Network offers accessible workflow training and analysis tutorials. The Carpentries provides foundational computing, data, shell, Git, and programming training. These resources can help analysts develop the skills needed for reproducible RNA-seq analysis.

Professional Escalation Criteria

Analysts should seek additional expertise when certain conditions arise in a limma-voom analysis.

When to Consult a Bioinformatics Specialist

A bioinformatics specialist should be consulted when the experimental design is complex, such as with multiple factors, interactions, or batch effects that are difficult to model. A specialist should also be consulted when the data have unusual characteristics, such as a high proportion of outliers, strong batch effects, or a large number of samples with very low library sizes.

When to Consider Alternative Methods

Alternative methods should be considered when the assumptions of limma-voom are violated. For population-level studies with large sample sizes and heterogeneous samples, the Wilcoxon rank-sum test may be more robust. For data with many outliers, robustified versions of voom or other robust methods may be appropriate.

When to Validate with Independent Methods

Key findings should be validated with an independent method such as qRT-PCR or with an alternative differential expression method. The consensus approach, where multiple methods are applied and the intersection of results is used, can increase confidence in the gene list.

When to Revisit the Experimental Design

If the analysis produces unexpected results, such as very few or very many differentially expressed genes, the experimental design should be revisited. The sample size may be insufficient, the groups may not be well separated, or the covariates may not be correctly specified.

Frequently Asked Questions

What is the difference between voom and limma-trend?

The voom method transforms count data to log2-CPM values and generates precision weights based on the mean-variance relationship. The limma-trend method fits a trend to the mean-variance relationship and uses this trend to moderate the variance estimates without generating per-observation weights. Both methods allow the use of limma linear modeling tools for RNA-seq data, but voom is generally preferred for count data because the precision weights account for the heteroscedasticity more directly.

How many biological replicates are needed for a limma-voom analysis?

The number of biological replicates needed depends on the effect size, the variability between samples, and the desired statistical power. The empirical Bayes moderation in limma helps with small sample sizes, but the reliability of the variance estimates and the mean-variance trend improves with more replicates. A benchmarking study of circular RNA methods found that almost all methods improved when the number of replicates increased. Three biological replicates per group is a common minimum, but more replicates are recommended when the expected effect size is small or the variability is high.

Can limma-voom be used for single-cell RNA-seq data?

Limma-voom can be used for single-cell RNA-seq data after the data are aggregated to the pseudo-bulk level. In a pseudo-bulk analysis, the counts from all cells of a given cell type in a given sample are summed, and the resulting pseudo-bulk counts are analyzed with standard bulk RNA-seq methods. A study of heart failure with preserved ejection fraction used limma-voom for differential expression analysis by cell type after pooling nuclei from multiple patients. The study demonstrated that limma-voom can be applied to pseudo-bulk single-cell data.

What is the role of TMM normalization in the limma-voom workflow?

TMM normalization computes scaling factors that correct for differences in RNA composition between samples. Without this correction, a few highly expressed genes in one sample can suppress the apparent expression of all other genes in that sample. The normalization factors are applied to the library sizes before the CPM conversion in the voom transformation. The voom function can use the normalization factors from edgeR if they are computed and passed correctly.

How should the results of a limma-voom analysis be reported?

The results should include the number of samples per group, the filtering criteria, the normalization method, the design matrix specification, the contrasts tested, and the thresholds used for significance. The session information should be saved and reported to document the versions of R and all packages used. The results table should include the log2 fold change, the average expression, the moderated t-statistic, the raw p-value, and the adjusted p-value for each gene.

What is the difference between a raw p-value and an adjusted p-value?

The raw p-value is the probability of observing a test statistic as extreme as the one observed, assuming the null hypothesis is true. The adjusted p-value accounts for the multiple testing problem by controlling the false discovery rate across all tested genes. An adjusted p-value of 0.05 means that approximately 5% of the genes called significant are expected to be false positives. The adjusted p-value should be used for deciding which genes are differentially expressed.

Can limma-voom handle paired or blocked experimental designs?

Yes, the linear modeling framework in limma naturally handles paired designs, blocked designs, and designs with technical covariates. A paired design is encoded with a design matrix that includes a term for each subject. This approach accounts for the correlation between samples from the same subject and increases the power to detect differential expression. An analysis of the Airway dataset found that correctly identifying paired experimental designs increased differential expression detection by 21 to 47%.

What should I do if my limma-voom results do not match the results from DESeq2 or edgeR?

Differences between methods are expected because each method makes different assumptions and uses different statistical approaches. The extended review of RNA-seq differential expression methods found that limma+voom, NOISeq, and DESeq2 had the most consistent results when compared with qRT-PCR data. The review also found that a consensus among five methods produced a list of differentially expressed genes with great accuracy. If the results differ substantially, the data should be examined for outliers, the normalization should be checked, and the design matrix should be reviewed.

Related Bioinformatics Guides

• RNA-Seq Differential Expression: DESeq2, edgeR, and limma-voom Frameworks
• RNA-Seq Data Analysis Workflow: From Raw Reads to Insights
• RNA Sequencing Data Analysis: From Raw Reads to Differential Expression
• Proteomics Data Analysis in R: A Practical Workflow for Differential Expression and Visualization
• RNA-Seq Databases: Accessing and Using Public RNA-Seq Data

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

References and Further Reading

• NCBI Data Resources. National Center for Biotechnology Information.
• EMBL-EBI Training. European Bioinformatics Institute.
• Bioconductor. Bioconductor Project.
• Galaxy Training Network. Galaxy Project.
• nf-core Documentation. nf-core.
• The Carpentries Lessons. The Carpentries.
• Single cell transcriptomic analyses of human heart failure with preserved ejection fraction.. bioRxiv : the preprint server for biology, 2025.
• RNA-Seq differential expression analysis: An extended review and a software tool.. PloS one, 2017.
• Robust identification of differentially expressed genes from RNA-seq data.. Genomics, 2020.
• Systematic benchmarking of statistical methods to assess differential expression of circular RNAs.. Briefings in bioinformatics, 2023.
• Brain lncRNA-mRNA co-expression regulatory networks and alcohol use disorder.. Genomics, 2024.
• Gene expression analysis of combined RNA-seq experiments using a receiver operating characteristic calibrated procedure.. Computational biology and chemistry, 2021.
• New Clues to Prognostic Biomarkers of Four Hematological Malignancies.. Journal of Cancer, 2022.
• Exaggerated false positives by popular differential expression methods when analyzing human population samples.. Genome biology, 2022.
• Identification of Radiation-Induced Injury Pathways and Hub Genes from RNA-Seq Data Based on Integrative Bioinformatics Approach.. 2026.
• Confidence: a web app for cross-platform differential gene expression analysis, gene scoring, and enrichment analysis.. 2026.
• Matched embryo-endometrium RNA-seq reveals coordinated but asymmetric transcriptomic reprogramming at the onset of early equine pregnancy.. 2026.
• ARIA: Adaptive Reasoning for Integrated Analysis - An LLM-Powered Framework for Autonomous Transcriptome Analysis with Decision-Aware Workflow Orchestration. 2026.
• BRPtools: An AutoML-Powered web platform for multiclass disease prediction from bulk blood RNA-seq data.. 2026.
• Multidimensional exposure architecture shapes vaping-associated transcriptomic dysregulation in oral epithelium.. 2026.
• Substantial unannotated noncoding transcripts in tumors may transcriptionally regulate cancer-related genes.. 2026.
• voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology, 2014.
• RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Research, 2016.
• Response to "Neglecting normalization impact in semi-synthetic RNA-seq data simulation generates artificial false positives" and "Winsorization greatly reduces false positives by popular differential expression methods when analyzing human population samples". Genome Biology, 2024.
• Integrated COVID-19 Predictor: Differential expression analysis to reveal potential biomarkers and prediction of coronavirus using RNA-Seq profile data. Comput. Biol. Medicine, 2022.
• Comprehensive RNA-Seq analysis of molecular signatures between COVID-19 and Idiopathic Pulmonary Fibrosis. Journal of high school science, 2023.
• DEGnet: Identifying Differentially Expressed Genes Using Deep Neural Network from RNA-Seq Datasets. Lecture Notes in Computer Science Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics, 2019.

This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.