limma-voom

The core innovation of the workflow. voom (variance modeling at the observational level) transforms raw count data to log-counts per million (log-cpm), but critically, it does not assume constant variance. Instead, it estimates the mean-variance trend of the log-cpm values across all genes. This trend is then used to calculate a precision weight for each individual observation.

• Low counts with high variance receive low weights.
• High counts with stable variance receive high weights.
• The output is a standard EList object ready for linear modeling.

Once the voom transformation is complete, the data is analyzed using the mature statistical engine of limma. This provides access to a suite of powerful tools unavailable in simpler count-based methods.

• Complex Design Matrices: Model factorial experiments, time courses, and paired samples with ease.
• Empirical Bayes Moderation: Borrows information across genes to stabilize variance estimates, a critical advantage when sample sizes are small.
• Quantitative Weights: The precision weights from voom are seamlessly integrated into the linear model fit, ensuring that noisy observations do not unduly influence results.

The limma-voom pipeline is designed to work with established normalization factors, not raw counts. It requires effective library sizes to be computed first.

• TMM Normalization: The default and recommended method, using calcNormFactors from the edgeR package to account for compositional biases between libraries.
• RLE Normalization: Relative Log Expression normalization from the DESeq2 package is also fully compatible.
• The workflow applies these factors to scale library sizes before computing log-cpm, ensuring that the transformation is performed on properly normalized data.

Genes with very low counts across all samples provide little statistical information and can destabilize the mean-variance estimation. The standard limma-voom pipeline includes a filtering step before analysis.

• A typical threshold is to retain genes with a minimum of 10-15 counts in a sufficient number of samples.
• Filtering is performed on raw counts, prior to TMM normalization and the voom transformation.
• This practice increases statistical power by reducing the burden of multiple testing correction and improves the accuracy of the empirical Bayes shrinkage.

The workflow provides specific diagnostic plots to verify the quality of the transformation and the model fit.

• Mean-Variance Trend Plot: Visualizes the square root of the standard deviation against average log-count. A successful voom transformation shows a smoothly decreasing trend that stabilizes at higher expression levels.
• voom With Quality Weights: An extension (voomWithQualityWeights) that estimates sample-level weights to down-weight outlier samples, increasing robustness in large cohort studies.
• SAF (Size-Adjusted Factors): The voom function outputs the normalized log-cpm matrix, precision weights, and the design matrix, all of which can be passed directly to standard limma functions like lmFit and eBayes.

While DESeq2 and edgeR model counts directly using a Negative Binomial distribution, limma-voom takes a different philosophical approach by transforming the data to approximate normality.

• Performance: limma-voom has been shown to control the false discovery rate well, especially in experiments with larger sample sizes.
• Speed: The linear modeling backend is computationally fast, even for very large datasets.
• Flexibility: Its primary advantage is the ability to analyze arbitrary experimental designs and incorporate quantitative weights, making it the method of choice for complex RNA-seq experiments beyond simple two-group comparisons.

The core innovation of the limma-voom workflow. Voom (variance modeling at the observational level) converts discrete, heteroskedastic RNA-seq counts into continuous log-counts per million (logCPM) values with associated precision weights. It estimates the mean-variance relationship of the log-counts through lowess smoothing, then assigns a weight to each observation inversely proportional to its variance. This allows the well-established limma linear modeling framework—originally designed for microarray data—to be applied to count-based RNA-seq data with full statistical rigor.

A critical output of the voom step that enables accurate modeling. Each gene-sample observation receives a precision weight that quantifies its reliability. Observations from genes with low counts receive low weights because their logCPM values are highly variable and unreliable. High-count observations receive high weights, anchoring the model fit. These weights are passed directly into the limma lmFit() function, effectively down-weighting noisy data points during linear model fitting and empirical Bayes moderation.

The statistical engine borrowed from limma that stabilizes variance estimates. In RNA-seq experiments with small sample sizes, gene-wise variance estimates are unstable. The empirical Bayes procedure shrinks each gene's estimated variance toward a common prior distribution estimated from all genes on the array. This moderated t-statistic prevents genes with tiny sample variances from being falsely ranked as highly significant, providing robust inference even when the number of biological replicates is limited.

The flexible statistical architecture that distinguishes limma-voom from purpose-built RNA-seq tools. The design matrix encodes the experimental setup, supporting complex designs beyond simple two-group comparisons:

• Multi-factor experiments with additive and interaction effects
• Continuous covariates like age or disease severity
• Time-course analyses with polynomial trends
• Paired or blocked designs with random effects This generality makes limma-voom the preferred choice for complex clinical and longitudinal studies.

A specialized limma function that handles technical replicates and repeated measures. The duplicateCorrelation() function estimates the intra-subject correlation when multiple samples come from the same biological unit. This correlation is then incorporated into the linear model, correctly accounting for the non-independence of observations. This is essential for experiments with repeated biopsies, technical replicates, or longitudinal sampling from the same patient.

The mechanism for extracting specific comparisons from a fitted model. A contrast matrix defines linear combinations of model coefficients to test. For example, in a two-factor experiment with treatment and time, a contrast can isolate the treatment effect averaged over time, or the treatment-by-time interaction. The eBayes() function then computes moderated t-statistics, p-values, and Benjamini-Hochberg adjusted p-values for each contrast, controlling the false discovery rate across thousands of genes.