edgeR for Differential Expression: A Practical Guide to Using edgeR and TMM Normalization

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

Key Takeaways

• edgeR employs a negative binomial generalized linear model framework, incorporating a quasi-likelihood approach for robust differential expression analysis of RNA-seq count data, particularly crucial for experiments with limited biological replicates.
• TMM (Trimmed Mean of M values) normalization is the default method in edgeR, designed to correct for RNA composition bias and library size differences, thereby enabling more reliable cross-sample expression comparisons.
• Dispersion estimation is a critical step, with edgeR using an empirical Bayes approach to borrow information across genes, stabilizing estimates and improving statistical power by accounting for biological variability beyond Poisson assumptions.
• The recommended edgeR workflow involves creating a DGEList object, filtering low-expression genes using filterByExpr, fitting a quasi-likelihood model with glmQLFit, and testing contrasts with glmQLFTest to identify differentially expressed genes.
• Reproducibility is paramount; analyses should be documented with version control (R, edgeR, and package versions via sessionInfo) and executable scripts, with results reported detailing normalization, dispersion estimation, and testing procedures.
• Careful quality control, including MDS plots for sample relationships and dispersion plots (plotBCV, plotQLDisp), is essential to identify potential issues like sample mislabeling or batch effects before statistical inference.

RNA sequencing has become the standard method for transcriptome-wide expression profiling, and differential expression analysis is the statistical core that converts raw read counts into biological insight. This guide addresses the specific problem of performing reliable differential expression analysis with the edgeR Bioconductor package, with particular attention to TMM normalization and generalized linear model fitting. The intended reader is a researcher, student, or laboratory professional who has generated or received RNA-seq count data and needs a defensible, reproducible path from count matrix to differentially expressed gene list. The guide assumes familiarity with R basics but does not assume prior experience with edgeR or with the statistical theory underlying count-based expression analysis.

The workflow described here follows the standard edgeR quasi-likelihood pipeline, which has been demonstrated in published computational workflows for RNA-seq experiments. A complete analysis of an RNA-seq experiment profiling epithelial cell subsets in the mouse mammary gland, using Rsubread for alignment and quantification and edgeR for statistical analysis, illustrates the practical steps of the pipeline from read alignment through pathway analysis. The quasi-likelihood functionality of edgeR is the recommended statistical approach within that workflow. This guide translates that pipeline into concrete decisions you can apply to your own dataset.

What edgeR Does and Where It Fits in the RNA-Seq Workflow

Differential expression analysis asks a simple question: which genes show statistically significant changes in expression between experimental conditions? The answer requires careful attention to data structure, normalization, dispersion estimation, and model specification. edgeR is an R package from the Bioconductor project that implements statistical methods specifically designed for count-based RNA-seq data. It operates on a matrix of read counts, where rows are genes and columns are samples, and it uses negative binomial models to account for the fact that RNA-seq counts exhibit more variability than a simple Poisson model would predict.

The position of edgeR within the broader RNA-seq workflow matters for practical decisions. Before edgeR can be used, raw sequencing reads must be aligned to a reference genome or transcriptome, and aligned reads must be summarized into counts per gene or transcript. These upstream steps are covered by tools such as Rsubread, which performs alignment and count quantification within the R environment. The edgeR analysis itself begins after count quantification, when you have a count matrix and a table describing your experimental design. Downstream of edgeR, differentially expressed gene lists feed into pathway analysis, enrichment testing, and biological interpretation.

The Bioconductor project provides the official package documentation, installation instructions, and workflow vignettes for edgeR. Bioconductor is the standard distribution channel for open-source genomic analysis software, and its packages are maintained with attention to reproducibility and version control. Installing edgeR through Bioconductor ensures that you receive a version compatible with your R installation and with other Bioconductor packages you may need for upstream or downstream analysis.

Core Principles of Count-Based Differential Expression

The Negative Binomial Model and Biological Variability

RNA-seq count data have two sources of variability. Technical variability arises from library preparation, sequencing depth, and other experimental steps. Biological variability arises from genuine differences between replicate samples, such as differences between individual animals, cell cultures, or tissue samples. A Poisson model can account for technical variability, but it assumes that the mean equals the variance. RNA-seq data typically show greater variance than the mean, a property called overdispersion. The negative binomial model addresses this by introducing a dispersion parameter that quantifies the extra variability beyond what a Poisson model would predict.

edgeR estimates dispersion from the data instead of assuming a fixed value. The dispersion estimation step is critical because it directly affects which genes are called differentially expressed. If dispersion is underestimated, the test will produce too many false positives. If dispersion is overestimated, the test will miss true differences. edgeR uses an empirical Bayes approach that borrows information across genes to stabilize dispersion estimates, which is especially important for experiments with few replicates.

TMM Normalization and Library Size Adjustment

Normalization in RNA-seq addresses the fact that different samples may have different sequencing depths and different RNA composition. A sample sequenced to twice the depth of another sample will have roughly twice the counts for most genes, but this does not mean those genes are differentially expressed. The simplest normalization divides each sample's counts by its total library size, producing counts per million. However, this approach can be biased when a small number of highly expressed genes differ dramatically between conditions, because those genes consume a large fraction of the total counts and distort the effective library size for all other genes.

TMM normalization, which stands for trimmed mean of M values, addresses this composition bias. The method selects a reference sample, computes log fold changes and absolute expression levels for each gene relative to that reference, trims the most extreme values, and computes a weighted average of the remaining log fold changes. The result is a scaling factor for each sample that adjusts the effective library size. edgeR applies TMM normalization by default through the calcNormFactors function, and the normalized library sizes are used in all subsequent statistical calculations.

The practical consequence of TMM normalization is that it makes expression comparisons between samples more reliable, particularly when the RNA composition of samples differs substantially. For example, in a plant experiment comparing roots treated with silicon to untreated roots, the treated samples may have large changes in a relatively small number of genes related to silicon response. Those genes could dominate the total read count and bias a simple total-count normalization. TMM normalization trims those extreme genes and produces scaling factors that better reflect the typical gene behavior across the transcriptome.

Dispersion Estimation and the Quasi-Likelihood Framework

edgeR offers two main statistical frameworks for testing differential expression: the classic likelihood ratio test and the quasi-likelihood F-test. The quasi-likelihood framework is generally recommended because it provides more robust error control, particularly for experiments with small numbers of replicates. The quasi-likelihood approach estimates a dispersion trend across genes, then estimates a per-gene dispersion that is squeezed toward the trend. The final test uses a quasi-likelihood F-test that accounts for the uncertainty in the dispersion estimates.

The quasi-likelihood pipeline in edgeR proceeds through several functions. estimateDisp estimates the dispersion parameters, glmQLFit fits a quasi-likelihood negative binomial generalized linear model to each gene, and glmQLFTest performs the quasi-likelihood F-test for specified contrasts. This pipeline is demonstrated in published workflows that use edgeR for complete RNA-seq analyses from read alignment through pathway analysis.

Preparing Your Data for edgeR

Required Input Files

The primary input to edgeR is a matrix of read counts with genes in rows and samples in columns. Each entry in the matrix is an integer count representing the number of sequencing reads that aligned to a particular gene in a particular sample. The count matrix can be generated by various alignment and quantification tools, including Rsubread within R, or by external tools such as STAR, HISAT2, or Salmon. The key requirement is that the counts are raw integer counts, not normalized values, not log-transformed values, and not counts per million.

The second required input is a sample information table, often called a targets file or experimental design table. This table has one row per sample and columns describing the experimental factors, such as treatment group, time point, tissue type, or batch. The sample information table must be aligned with the columns of the count matrix, meaning that the order of samples in the count matrix must match the order of rows in the sample information table. Mismatches between these files are a common source of errors that are difficult to detect after analysis.

Creating the DGEList Object

edgeR stores count data and sample information in a DGEList object. The DGEList function creates this object from the count matrix and the sample information. The counts are stored in the counts component, and the sample information is stored in the samples component. The DGEList object also has a genes component that can hold gene annotation information such as gene symbols, chromosome locations, or gene biotypes.

Creating the DGEList is the first step in every edgeR analysis. The function checks that the count matrix contains non-negative integers and that the sample information has the correct number of rows. After creating the DGEList, you should verify that the sample names in the count matrix match the sample identifiers in your experimental design table. This verification step is simple but essential, because a mismatch will produce results that are biologically meaningless.

Filtering Low-Expression Genes

Not all genes in the count matrix are informative. Genes with very low counts across all samples contribute little information and add noise to the dispersion estimation. Filtering these genes before analysis improves statistical power and reduces the multiple testing burden. The standard approach is to retain genes that have a minimum number of counts in a minimum number of samples. A common filter retains genes with at least one count per million in at least a certain number of samples, where the number of samples depends on the smallest group size in the experiment.

The filterByExpr function in edgeR implements this filtering automatically. The function uses the experimental design to determine the minimum group size and computes a counts-per-million threshold based on the library sizes. Genes that do not meet the threshold are removed from the DGEList. The number of genes retained after filtering depends on the sequencing depth and the complexity of the transcriptome. For a typical mammalian RNA-seq experiment, filtering might retain 15,000 to 20,000 genes from an initial set of 30,000 to 60,000 annotated genes.

The edgeR Workflow Step by Step

Step 1: Load Packages and Read Data

The analysis begins by loading edgeR and any additional packages needed for data import and manipulation. The readDGE function can read count files directly, or you can read counts into a data frame and convert it to a DGEList. The sample information table should be read separately and checked for consistency with the count matrix.

library(edgeR)

## Read count matrix
counts <- read.delim("counts.txt", row.names = 1, header = TRUE)

## Read sample information
targets <- read.delim("targets.txt", stringsAsFactors = TRUE)

## Create DGEList
dge <- DGEList(counts = counts, group = targets$condition)

The group argument assigns each sample to a condition group. This grouping is used for the initial exploratory analysis and for the simple two-group comparison. For more complex designs with multiple factors, the design matrix is specified later in the workflow.

Step 2: Quality Control and Exploratory Analysis

Before fitting statistical models, you should examine the data for obvious problems. The plotMDS function produces a multidimensional scaling plot that shows the relationships between samples. Samples from the same biological condition should cluster together, and samples from different conditions should separate. Samples that cluster unexpectedly may indicate sample mislabeling, batch effects, or other technical problems.

The plotBCV function displays the biological coefficient of variation against the average log counts per million. This plot shows the dispersion trend and the per-gene dispersion estimates. The plot should show a decreasing trend, with higher dispersion at low expression levels and lower dispersion at high expression levels. An unusual pattern may indicate problems with the data or with the experimental design.

Quality control steps are recommended but not mandatory in RNA-seq analysis pipelines. Failing to check the characteristics of the dataset may lead to spurious results. Dedicated analysis pipelines such as SARTools include systematic quality control steps and diagnostic plots to help tune model parameters and prevent errors from misusing the statistical methods.

Step 3: TMM Normalization

TMM normalization is applied with the calcNormFactors function. This function computes a scaling factor for each sample and stores it in the DGEList object. The default method is TMM, and the default reference is the sample with the median library size.

dge <- calcNormFactors(dge)

After normalization, you can examine the normalization factors with dge$samples$norm.factors. The factors should be close to 1 for most samples. A normalization factor substantially different from 1, such as 0.5 or 2, indicates that the sample has a very different RNA composition from the reference. This is not necessarily an error, but it warrants investigation.

The normalized library sizes are the product of the raw library sizes and the normalization factors. These normalized library sizes are used in all subsequent statistical calculations, including the calculation of counts per million for visualization and the fitting of the negative binomial models.

Step 4: Design Matrix and Model Fitting

The design matrix specifies the experimental model that edgeR will fit. For a simple two-group comparison, the design matrix can be created with model.matrix(~ group). For more complex designs with multiple factors, interactions, or blocking variables, the design matrix must be constructed to reflect the experimental design.

design <- model.matrix(~ condition, data = targets)

The design matrix has one row per sample and one column per model parameter. The first column is typically the intercept, and subsequent columns represent the effects of the experimental factors. The interpretation of the model coefficients depends on the parameterization of the design matrix, so you should be familiar with the model you are fitting.

For experiments with multiple factors, such as a treatment factor and a time factor, the design matrix can include main effects and interactions. The DiCoExpress tool demonstrates the use of generalized linear models with automated contrast writing for multifactorial RNA-seq experiments. The tool uses edgeR for differential expression analysis and supports contrasts within generalized linear models, allowing users to test specific comparisons of interest.

Step 5: Dispersion Estimation

Dispersion estimation is performed with the estimateDisp function. This function estimates the common dispersion, the trended dispersion, and the tagwise dispersion. The common dispersion is a single value shared by all genes, the trended dispersion varies with expression level, and the tagwise dispersion is the per-gene estimate that is squeezed toward the trend.

dge <- estimateDisp(dge, design)

The estimateDisp function uses the design matrix to account for the experimental structure. The resulting dispersion estimates are stored in the DGEList object and are used by the quasi-likelihood fitting functions.

Step 6: Quasi-Likelihood Fitting and Testing

The quasi-likelihood pipeline fits a negative binomial generalized linear model to each gene using the glmQLFit function, then tests specified contrasts with the glmQLFTest function.

fit <- glmQLFit(dge, design)
qlf <- glmQLFTest(fit, contrast = c(-1, 1))

The contrast vector specifies the comparison being tested. For a two-group design with an intercept and a condition effect, the contrast c(-1, 1) tests the difference between the two conditions. The result is a test for each gene, with a log fold change, a log counts per million, a statistic, a p-value, and an adjusted p-value.

The topTags function extracts the top differentially expressed genes sorted by p-value or by absolute log fold change. The output includes the gene identifiers, the log fold change, the average log counts per million, the statistic, the p-value, and the adjusted p-value.

results <- topTags(qlf, n = Inf)

The adjusted p-values control the false discovery rate across all tested genes. The default adjustment method is the Benjamini-Hochberg method, which controls the expected proportion of false positives among the genes called significant.

Step 7: Extracting and Interpreting Results

The results table contains all the information needed to identify differentially expressed genes. The log fold change indicates the direction and magnitude of the expression change. A positive log fold change means the gene is upregulated in the second condition relative to the first, and a negative log fold change means the gene is downregulated.

The adjusted p-value indicates the statistical significance of the change after multiple testing correction. A common threshold is an adjusted p-value below 0.05, but the appropriate threshold depends on the research question and the expected number of true positives. Some analyses use a more stringent threshold, such as 0.01, while others combine the adjusted p-value threshold with a minimum absolute log fold change, such as 1 or 2, to focus on biologically meaningful changes.

The number of differentially expressed genes depends on the biological difference between conditions, the number of replicates, the sequencing depth, and the dispersion. Experiments with large biological differences and many replicates will identify more differentially expressed genes than experiments with small differences and few replicates.

At a Glance: edgeR Workflow Summary

Comparing edgeR with DESeq2

edgeR and DESeq2 are the two most widely used R packages for differential expression analysis of RNA-seq count data. Both packages use negative binomial models, but they differ in normalization methods, dispersion estimation, and testing procedures. Understanding these differences helps you choose the appropriate tool for your data and interpret results correctly.

DESeq2 uses a median-of-ratios normalization method, which computes a size factor for each sample based on the median ratio of each gene's count to the geometric mean across samples. edgeR uses TMM normalization, which computes scaling factors based on trimmed mean of M values. Both methods address composition bias, but they use different algorithms and may produce slightly different normalization factors for the same dataset.

The dispersion estimation also differs. DESeq2 fits a dispersion trend and shrinks per-gene dispersion estimates toward the trend using an empirical Bayes approach. edgeR uses a similar empirical Bayes approach but with a different implementation. The quasi-likelihood framework in edgeR provides an additional layer of robustness by accounting for uncertainty in the dispersion estimates.

A comparison of edgeR and DESeq2 using four real RNA-seq plant datasets found a large number of jointly identified differentially expressed genes between the two methods. However, the appropriate workflow depends on the research goal and the experimental design. Different approaches to statistical analysis and interpretation can be suggested depending on the dataset, and the goal is to minimize the number of falsely identified differentially expressed genes.

For most experiments, edgeR and DESeq2 will produce broadly similar results, with substantial overlap in the lists of differentially expressed genes. Differences are more likely at the margins, where genes have borderline significance or small fold changes. If the choice of method materially affects your biological conclusions, you should investigate the genes that differ between methods and consider whether the experimental design or data quality explains the discrepancy.

Handling Complex Experimental Designs

Multifactor Designs and Interactions

Many RNA-seq experiments involve more than one experimental factor. A plant experiment might compare two genotypes under two treatment conditions, producing a factorial design with four groups. A time course experiment might sample cells at multiple time points after a treatment. A clinical study might include patients with different disease subtypes and treatment responses.

edgeR handles these designs through the design matrix. The design matrix can include main effects for each factor and interaction terms that allow the effect of one factor to depend on the level of another factor. The model.matrix function constructs these designs from a formula, and the resulting design matrix is used in dispersion estimation and model fitting.

For a factorial design with two factors, the design matrix might include an intercept, a main effect for the first factor, a main effect for the second factor, and an interaction term. The interaction term tests whether the effect of the first factor differs across levels of the second factor. Testing the interaction requires a contrast that isolates the interaction effect.

The DiCoExpress tool demonstrates the use of contrasts within generalized linear models for multifactorial RNA-seq experiments. The tool automates contrast writing, allowing users to test specific comparisons of interest without manually specifying contrast vectors. This approach is particularly useful for complex designs with many possible comparisons.

Blocking Factors and Batch Effects

Batch effects are technical sources of variation that affect multiple samples processed together. Samples processed in the same batch may share systematic differences in library preparation, sequencing runs, or other technical steps. If batch effects are not accounted for in the analysis, they can inflate the variability and obscure true biological differences, or they can create spurious differences that reflect technical artifacts instead of biology.

edgeR can account for batch effects by including the batch as a factor in the design matrix. This approach, called blocking, adjusts for the batch effect while testing the biological factors of interest. The design matrix includes a main effect for the batch factor, and the differential expression test evaluates the biological factor while holding the batch effect constant.

The SARTools pipeline supports designs with a blocking factor such as a batch effect or sample pairing. The pipeline is based on DESeq2 and edgeR and includes systematic quality control steps to check the characteristics of the dataset. The diagnostic plots help tune the model parameters and prevent errors from misusing the statistical methods.

Contrasts for Specific Comparisons

Contrasts allow you to test specific comparisons that are not directly represented by individual model coefficients. In a factorial design, the interaction effect is tested with a contrast that combines coefficients. In a time course experiment, you might compare the expression at each time point to the baseline time point, requiring multiple contrasts.

The makeContrasts function in edgeR constructs contrast vectors from a named list of comparisons. The contrasts are then passed to glmQLFTest for testing. Multiple contrasts can be tested simultaneously, and the results are combined into a single table with one row per gene per contrast.

The automated contrast writing in DiCoExpress simplifies this process for users who are not familiar with the details of contrast construction. The tool generates contrasts based on the experimental design and the comparisons of interest, reducing the risk of errors in contrast specification.

Quality Control and Diagnostic Checks

Sample-Level Quality Control

Before fitting statistical models, you should examine the data for sample-level problems. The multidimensional scaling plot from plotMDS shows the relationships between samples based on the top differentially expressed genes. Samples from the same condition should cluster together, and the separation between conditions should be visible. Samples that cluster with the wrong condition may be mislabeled, contaminated, or otherwise problematic.

The library sizes, which are the total number of reads assigned to genes for each sample, should be similar across samples within an experiment. A sample with a much smaller library size than the others may have failed during library preparation or sequencing. A sample with an unusually large library size may have been sequenced more deeply or may contain contamination from another sample.

The plotMDS function can also be used to examine the effect of normalization. Comparing the MDS plots before and after normalization can show whether normalization improves the separation between conditions or introduces unexpected patterns.

Gene-Level Quality Control

Gene-level quality control focuses on the behavior of individual genes across samples. The biological coefficient of variation plot from plotBCV shows the relationship between dispersion and expression level. The plot should show a decreasing trend, with higher dispersion at low expression levels. Genes with unusually high dispersion at high expression levels may be outliers or may represent genes with genuinely variable expression.

The plotQLDisp function displays the quasi-likelihood dispersion estimates against the abundance. This plot is produced after glmQLFit and shows the trend and the per-gene estimates. The plot helps identify genes with unusually high dispersion that may warrant further investigation.

Diagnostic Plots from Analysis Pipelines

Dedicated analysis pipelines such as SARTools include systematic quality control steps and diagnostic plots that help tune model parameters. The SARTools pipeline produces an HTML report that displays diagnostic plots for quality control and model hypothesis checking. The report also keeps track of the whole analysis process, parameter values, and versions of the R packages used.

The diagnostic plots in SARTools include the MDS plot, the dispersion plot, and plots of the normalization factors. The pipeline also checks the model assumptions and provides warnings when the assumptions appear to be violated. These checks are important because failing to check the characteristics of the dataset may lead to spurious results.

Common Failure Patterns and How to Avoid Them

Sample Mislabeling and Design Mismatches

The most common and most damaging error in differential expression analysis is a mismatch between the count matrix and the sample information. If the sample order in the count matrix does not match the sample order in the targets file, every comparison in the analysis will be wrong. This error is difficult to detect after the analysis because the results will look plausible, with differentially expressed genes that reflect the mislabeled groups instead of the true biological conditions.

The best defense is to verify the sample order before analysis. Check that the column names of the count matrix match the sample identifiers in the targets file. Use the MDS plot to confirm that samples cluster by condition as expected. If a sample clusters with the wrong condition, investigate before proceeding.

Inappropriate Normalization

TMM normalization is appropriate for most RNA-seq experiments, but it can produce unexpected results in certain situations. If the normalization factors are extreme, such as below 0.5 or above 2, the sample may have a very different RNA composition from the reference. This can happen when a sample is contaminated, when a sample has a very different distribution of expression levels, or when the reference sample is not representative.

The choice of reference sample can affect the normalization factors. edgeR selects the reference sample with the median library size by default, but you can specify a different reference. If the default reference produces extreme normalization factors, try a different reference or investigate the sample composition.

Overdispersion and Underdispersion

The negative binomial model assumes that the dispersion is positive and that the data are overdispersed relative to a Poisson model. If the dispersion estimates are very small, the data may be underdispersed, which can happen with technical replicates or with data that have been processed in a way that removes biological variability. If the dispersion estimates are very large, the data may have more variability than the model can accommodate, which can happen with poor-quality samples or with hidden batch effects.

The dispersion plot helps identify these problems. If the dispersion estimates are unusually small or large, investigate the data quality and consider whether the experimental design is appropriate.

Multiple Testing Burden

Testing thousands of genes simultaneously creates a multiple testing problem. Without correction, many genes will appear significant by chance alone. edgeR uses the Benjamini-Hochberg method to control the false discovery rate, but the adjusted p-values depend on the number of tests and the distribution of p-values.

The number of differentially expressed genes should be interpreted in the context of the multiple testing correction. A small number of significant genes may reflect a weak biological effect, a small number of replicates, or high variability. A very large number of significant genes may reflect a strong biological effect or may indicate problems with the data or the model.

Reproducibility and Reporting

Version Control and Session Information

Reproducibility requires that the analysis can be repeated with the same results. The version of R, the version of edgeR, and the versions of all other packages used in the analysis should be recorded. The sessionInfo function in R prints the versions of R and all loaded packages, and this information should be saved with the analysis results.

The Bioconductor project maintains versioned releases of its packages, and the release version should be recorded. The BiocManager::version function reports the Bioconductor version. The combination of R version, Bioconductor version, and package versions determines the exact behavior of the analysis.

Scripts and Workflow Documentation

The analysis should be documented in a script that can be rerun from the raw count data. The script should include comments explaining each step and the rationale for parameter choices. The script should be saved with the analysis results so that the analysis can be reproduced or modified.

Workflow management systems such as nf-core provide standardized pipelines for RNA-seq analysis, including quality control, alignment, and quantification. These pipelines are designed for reproducibility and can be configured for different experimental designs. The nf-core documentation describes the usage, configuration, and customization of these pipelines.

Reporting Results in Publications

The methods section of a publication should describe the differential expression analysis in sufficient detail that a reader can reproduce it. The description should include the version of edgeR, the normalization method, the dispersion estimation approach, the statistical test, and the thresholds for significance. The number of differentially expressed genes and the criteria for calling significance should be reported.

The European Bioinformatics Institute provides training resources for bioinformatics analysis, including practical analysis education and data-resource training. These resources can help researchers develop the skills needed to perform and report differential expression analysis correctly.

Limitations and Interpretation Caveats

Statistical Significance Does Not Imply Biological Importance

A gene with a statistically significant adjusted p-value may have a small fold change that is not biologically meaningful. Conversely, a gene with a large fold change may not reach statistical significance if the variability is high or the number of replicates is small. The interpretation of differential expression results should consider both the statistical significance and the magnitude of the change.

The log fold change reported by edgeR is the log2 of the ratio of normalized expression levels between conditions. A log fold change of 1 corresponds to a 2-fold change, and a log fold change of 2 corresponds to a 4-fold change. The threshold for biological importance depends on the biological context and the expected magnitude of changes for the system under study.

Replicate Number and Statistical Power

The number of biological replicates is the most important factor determining the statistical power of a differential expression analysis. With few replicates, the dispersion estimates are uncertain, and the test has limited power to detect small or moderate changes. With more replicates, the dispersion estimates are more precise, and the test can detect smaller changes.

The recommended minimum number of replicates is three per condition, but more replicates are better. The DiCoExpress tutorial used a dataset with three replicates for each condition, which is sufficient for a basic analysis but may not be sufficient for detecting subtle changes. The relationship between replicate number and statistical power should be considered when designing experiments and interpreting results.

Composition Bias and Normalization Assumptions

TMM normalization assumes that most genes are not differentially expressed between conditions. If a large fraction of the transcriptome changes between conditions, the normalization may be biased. This situation can arise in comparisons of very different cell types or tissues, or in experiments with strong global responses such as heat shock or immune activation.

The normalization factors should be examined for extreme values, and the results should be interpreted with caution if the normalization assumptions are violated. Alternative normalization methods may be considered, but the choice of normalization method should be justified based on the data characteristics.

Gene-Specific Outliers and Missing Observations

edgeR can show weak performance against gene-specific outliers, where a single sample has an unusually high or low count for a particular gene. These outliers can inflate the dispersion estimate for the gene and reduce the power to detect differential expression. A proposed pre-processing approach combines outlier detection and missing imputation to boost the performance of edgeR Robust. The approach uses leave-one-out cross-validation for outlier detection and random forest-based imputation for missing observations, and it outperformed conventional edgeR Robust in both simulation and real data analysis.

If your data contain gene-specific outliers, consider whether the outliers reflect genuine biology or technical artifacts. Outliers caused by alignment errors, mapping ambiguities, or other technical issues should be investigated and potentially removed. Outliers that reflect genuine biological variability should be retained, but the results should be interpreted with caution.

Practical Implementation Steps

Step 1: Verify Data Integrity

Before starting the analysis, verify that the count matrix contains non-negative integers and that the sample names match the experimental design. Check the library sizes for each sample and investigate any samples with unusually small or large libraries. Confirm that the number of samples in the count matrix matches the number of rows in the targets file.

Step 2: Create the DGEList and Filter Genes

Create the DGEList object and filter low-expression genes using filterByExpr. Examine the number of genes retained after filtering and compare it to the total number of genes in the count matrix. If the filtering removes an unexpectedly large fraction of genes, investigate whether the sequencing depth is adequate or whether the annotation is appropriate.

Step 3: Perform Quality Control

Generate the MDS plot and the BCV plot. Examine the MDS plot for sample clustering by condition and investigate any samples that cluster unexpectedly. Examine the BCV plot for the dispersion trend and investigate any unusual patterns.

Step 4: Normalize and Estimate Dispersion

Apply TMM normalization with calcNormFactors and examine the normalization factors. Estimate dispersion with estimateDisp and examine the dispersion plot. If the dispersion estimates are unusual, investigate the data quality and the experimental design.

Step 5: Fit the Model and Test Contrasts

Construct the design matrix and fit the quasi-likelihood model with glmQLFit. Specify the contrasts of interest and test with glmQLFTest. Extract the results with topTags and examine the number of differentially expressed genes at your chosen thresholds.

Step 6: Validate and Interpret Results

Examine the top differentially expressed genes and verify that they make biological sense. Check whether known marker genes for the conditions under study appear in the results. Consider whether the direction and magnitude of the changes are consistent with the biology of the system.

Step 7: Document and Report

Save the analysis script, the session information, and the results. Record the versions of R, Bioconductor, and edgeR. Document the parameter choices and the rationale for those choices. Prepare the methods description for publication or reporting.

Records and Measurements to Maintain

The following records should be maintained for each differential expression analysis:

Professional Escalation Criteria

Some situations require consultation with a bioinformatics specialist or statistician before proceeding with the analysis. These situations include:

• Samples that cluster unexpectedly on the MDS plot, which may indicate sample mislabeling, contamination, or batch effects that are not accounted for in the design
• Normalization factors that are extreme, which may indicate composition bias or sample quality problems
• Dispersion estimates that are unusually large or small, which may indicate model misspecification or data quality problems
• Results that are highly sensitive to the choice of normalization method, dispersion estimation approach, or statistical test
• Experiments with very few replicates, where the statistical power is limited and the results may be unreliable
• Complex experimental designs with multiple factors, interactions, or blocking variables, where the design matrix specification requires statistical expertise

The Galaxy Training Network provides accessible workflow training and analysis tutorials that can help researchers develop the skills needed to address these situations. The training resources are backed by expert-reviewed content and are designed to empower researchers to perform and interpret their own analyses.

Frequently Asked Questions

What is the difference between TMM normalization and simple counts per million normalization?

TMM normalization computes a scaling factor for each sample based on the trimmed mean of M values, which accounts for composition bias caused by a small number of highly expressed genes that differ between conditions. Simple counts per million normalization divides each sample's counts by its total library size, which does not account for composition bias. TMM normalization is the default in edgeR and is generally preferred for differential expression analysis because it produces more reliable comparisons between samples with different RNA compositions.

How many biological replicates do I need for edgeR analysis?

The recommended minimum is three biological replicates per condition, but more replicates provide greater statistical power and more precise dispersion estimates. The number of replicates needed depends on the biological variability of the system, the magnitude of the expression changes, and the sequencing depth. Experiments with high variability or small expected changes require more replicates to achieve adequate power.

What is the quasi-likelihood approach in edgeR and why is it recommended?

The quasi-likelihood approach fits a negative binomial generalized linear model to each gene and uses a quasi-likelihood F-test for differential expression. This approach accounts for uncertainty in the dispersion estimates and provides more robust error control than the classic likelihood ratio test, particularly for experiments with small numbers of replicates. The quasi-likelihood pipeline uses glmQLFit and glmQLFTest and is the recommended approach in published edgeR workflows.

How do I choose the threshold for calling genes differentially expressed?

The choice of threshold depends on the research question and the expected number of true positives. A common approach is to use an adjusted p-value below 0.05, which controls the false discovery rate at 5 percent. Some analyses use a more stringent threshold, such as 0.01, or combine the adjusted p-value threshold with a minimum absolute log fold change to focus on biologically meaningful changes. The threshold should be specified before the analysis and justified in the methods description.

What should I do if my samples do not cluster by condition on the MDS plot?

Samples that do not cluster by condition may indicate sample mislabeling, contamination, batch effects, or a weak biological effect. Investigate the sample information and the library preparation records to identify potential problems. Consider whether batch effects should be included in the design matrix. If the samples still do not cluster after accounting for known technical factors, the biological effect may be weak or the data quality may be inadequate for the intended comparison.

Can edgeR handle experiments with more than two conditions?

Yes, edgeR can handle experiments with any number of conditions through the design matrix and contrasts. For a single factor with multiple levels, the design matrix includes an intercept and indicator variables for each condition relative to a reference condition. Contrasts can be used to test specific comparisons between any pair of conditions or combinations of conditions. The DiCoExpress tool demonstrates the use of contrasts within generalized linear models for multifactorial experiments.

How do edgeR results compare with DESeq2 results?

edgeR and DESeq2 generally produce broadly similar results, with substantial overlap in the lists of differentially expressed genes. A comparison using four real plant RNA-seq datasets found a large number of jointly identified differentially expressed genes between the two methods. Differences are more likely at the margins, where genes have borderline significance or small fold changes. The choice of method should be based on the experimental design and the specific features of each package.

What is the role of quality control in differential expression analysis?

Quality control is essential for identifying problems that can produce spurious results. The quality control steps include examining sample clustering, library sizes, normalization factors, and dispersion estimates. Dedicated analysis pipelines such as SARTools include systematic quality control steps and diagnostic plots that help tune model parameters and prevent errors from misusing the statistical methods. Failing to check the characteristics of the dataset may lead to spurious results.

Related Bioinformatics Guides

• Proteomics Data Analysis in R: A Practical Workflow for Differential Expression and Visualization
• RNA Sequencing Data Analysis: From Raw Reads to Differential Expression
• RNA-Seq Differential Expression: DESeq2, edgeR, and limma-voom Frameworks
• RNA-Seq Normalization Methods: TPM, RPKM, and Beyond
• RNA-Seq Data Analysis in Galaxy: A User-Friendly Platform

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.
• DiCoExpress: a tool to process multifactorial RNAseq experiments from quality controls to co-expression analysis through differential analysis based on contrasts inside GLM models.. Plant methods, 2020.
• Protocol for integrating and interpreting multi-omics data combining unsupervised and supervised data integrating approaches.. 2026.
• Galaxy single-cell & spatial omics community update: Navigating new frontiers in 2025.. 2025.
• A cell-intrinsic glucocorticoid biosynthesis and sensing circuit maintains a homeostatic Th17 cell state.. 2026.
• scPASU: A computational protocol for quantifying polyadenylation site usage and alternative polyadenylation from 3' scRNA-seq data.. 2026.
• SigFuge (Tutorial). 2013.
• useR2013!: High-throughput sequence analysis with R and Bioconductor. 2012.
• Boosting edgeR (robust) by dealing with missing observations and gene-specific outliers in RNA-Seq profiles and its application to explore biomarker genes for diagnosis and therapies of ovarian cancer.. Genomics, 2024.
• Current best practices in single-cell RNA-seq analysis: a tutorial. Molecular Systems Biology, 2019.
• From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Research, 2016.
• SARTools: A DESeq2- and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data. bioRxiv, 2015.
• Tutorial: RNA-seq differential expression & pathway analysis with Sailfish, DESeq2, GAGE, and Pathview. 2015.
• Application of edgeR and DESeq2 methods in plant experiments based on RNA-seq technology. Biometrical Letters, 2022.

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