How to Run DESeq2 for Differential Expression Analysis: A Step-by-Step Tutorial in R

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

Key Takeaways

• DESeq2 requires raw integer count matrices and corresponding sample metadata for differential expression analysis, utilizing the negative binomial distribution to model count data and estimate dispersion.
• The core analysis involves creating a DESeqDataSet object, running the DESeq() function for normalization and statistical testing (defaulting to the Wald test), and extracting results using the results() function.
• Quality control is critical, with diagnostic plots like plotDispEsts(), plotPCA(), and sample distance heatmaps essential for identifying outliers and assessing biological grouping.
• Batch effects and confounding factors must be addressed by incorporating them into the design formula (e.g., ~ batch + condition) to ensure accurate estimation of biological effects.
• Shrinkage of log2 fold change estimates, particularly using methods like apeglm via lfcShrink(), improves stability and interpretability for genes with low counts or high dispersion.
• Results are typically filtered by adjusted p-value (e.g., < 0.05) and log2 fold change (e.g., > 1) and then annotated and exported for downstream functional enrichment analysis using tools like clusterProfiler.

Differential expression analysis with DESeq2 is a core step in RNA sequencing projects that compares gene counts between experimental groups to identify transcripts with statistically supported changes. This tutorial provides a reproducible R workflow that takes a count matrix, builds a DESeq2 object, applies normalization and statistical testing, and extracts results tables suitable for downstream interpretation. The intended reader is a biology student, researcher, or laboratory professional who has raw or processed RNA-seq count data and needs a practical path from those counts to a defensible list of differentially expressed genes.

The workflow assumes you have R installed, have access to Bioconductor for package installation, and have a count matrix where rows are genes and columns are biological samples. The steps below follow the standard DESeq2 analysis pipeline documented by the Bioconductor project, which maintains the package and provides installation and usage guidance for reproducible genomic analysis [<a href="#ref-1">1</a>]. You will also need basic familiarity with the R console and with tabular data structures, skills covered by foundational computing lessons from The Carpentries [<a href="#ref-2">2</a>].

At a Glance

Understanding the Input Data Requirements

DESeq2 expects a matrix of raw integer counts, where each value represents the number of sequencing reads that aligned to a particular gene in a particular sample. The count matrix must contain non-negative integers because the statistical model in DESeq2 is based on the negative binomial distribution, which is designed for count data [<a href="#ref-3">3</a>]. Floating point values, normalized counts, or transcript-per-million values are not appropriate inputs for the core analysis because the model estimates dispersion from the relationship between the mean and variance of raw counts.

The sample metadata table is equally important. This table must have one row per sample and columns that describe experimental conditions, batch identifiers, or other covariates you intend to model. The column names in the count matrix must match the row names in the metadata table, and the design formula you supply to DESeq2 must reference column names that exist in the metadata. A common source of errors is a mismatch between sample names in the count matrix and the metadata file, so verify this alignment before building the DESeq2 object.

Public repositories such as the NCBI Sequence Read Archive and Gene Expression Omnibus provide access to raw and processed RNA-seq data that can be used for testing workflows or for reanalysis [<a href="#ref-4">4</a>]. The NCBI also maintains search systems and sequence resources that help researchers locate appropriate datasets for their biological questions [<a href="#ref-4">4</a>]. If you are learning the workflow, starting with a small published dataset can help you confirm that each step produces expected outputs before you apply the pipeline to your own data.

Installing and Loading DESeq2 from Bioconductor

DESeq2 is distributed through Bioconductor, a project that provides software for the analysis and comprehension of high-throughput genomic data [<a href="#ref-1">1</a>]. Bioconductor packages are installed using a dedicated function instead of the standard install.packages command used for CRAN packages. The Bioconductor website provides installation instructions and documentation for the DESeq2 package, including the release version appropriate for your R installation [<a href="#ref-1">1</a>].

The installation process requires an internet connection and may take several minutes because DESeq2 has multiple dependencies. After installation, load the package with the library function. You will also need to load other packages for data manipulation and visualization, such as tidyverse for data wrangling and ggplot2 for plotting, though these are not strictly required for the core DESeq2 analysis.

A typical session begins with these commands:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2")
library(DESeq2)

The BiocManager package handles the installation of Bioconductor packages and their dependencies. If you encounter installation errors, check that your R version is compatible with the current Bioconductor release. The Bioconductor website provides a version-specific installation command that matches your R version [<a href="#ref-1">1</a>].

Building the DESeqDataSet Object

The DESeqDataSet is the central object that stores the count matrix, the sample metadata, and the results of the statistical analysis. You create this object using the DESeqDataSetFromMatrix function, which takes three primary arguments: the count matrix, the column metadata, and the design formula.

The design formula specifies the experimental model. For a simple two-group comparison, the formula is ~ condition, where condition is the column in your metadata that identifies which samples belong to which experimental group. For experiments with batch effects or paired samples, the design formula can include additional terms, such as ~ batch + condition. The order of terms matters because the last term in the formula is the one tested for differential expression.

dds <- DESeqDataSetFromMatrix(
    countData = count_matrix,
    colData = sample_metadata,
    design = ~ condition
)

The count matrix must have gene identifiers as row names and sample identifiers as column names. The sample metadata must have sample identifiers as row names, and these must match the column names of the count matrix exactly. If the row names of the metadata do not match the column names of the count matrix, DESeq2 will return an error.

Before running the full analysis, you may want to filter out genes with very low counts across all samples. Genes with zero counts in every sample provide no information and can be removed to reduce the size of the object. A common filter is to keep genes that have at least 10 total counts across all samples, though the appropriate threshold depends on your sequencing depth and the number of samples. DESeq2 performs some internal filtering during the analysis, but pre-filtering can speed up computation and reduce memory usage.

Running the Differential Expression Test

The core analysis is performed with a single function call. The DESeq function runs the entire pipeline, including estimation of size factors, estimation of dispersion, fitting of the negative binomial model, and testing for differential expression [<a href="#ref-3">3</a>].

dds <- DESeq(dds)

This function applies the default Wald test, which compares the expression level of each gene between the groups specified in the design formula. The Wald test is appropriate for experiments with two groups or for comparing each level of a factor to a reference level. For experiments with more than two groups where you want to test whether any group differs from the others, you can use the likelihood ratio test by specifying test = "LRT" and reduced = ~ 1.

The DESeq function performs several steps automatically. It estimates size factors to account for differences in sequencing depth between samples, estimates dispersion parameters for each gene, and applies shrinkage to fold change estimates to improve stability for genes with low counts or high variability [<a href="#ref-3">3</a>]. The shrinkage of fold changes is a key feature of DESeq2 that distinguishes it from simpler approaches that use raw fold changes without accounting for uncertainty in the estimates.

The statistical model in DESeq2 assumes that count data follow a negative binomial distribution, which accommodates the overdispersion commonly observed in RNA-seq data [<a href="#ref-3">3</a>]. The dispersion parameter captures the variability that exceeds what would be expected from a Poisson distribution, and DESeq2 borrows information across genes to obtain stable dispersion estimates even when the number of biological replicates is small [<a href="#ref-3">3</a>].

Extracting and Interpreting Results

After running the DESeq function, you extract the results table using the results function. This function returns a data frame with columns for the base mean, log2 fold change, standard error, test statistic, p-value, and adjusted p-value for each gene.

res <- results(dds)

The base mean is the average of the normalized count values across all samples. The log2 fold change represents the magnitude and direction of the expression change, with positive values indicating higher expression in the treatment group relative to the control group. The p-value tests the null hypothesis that the gene is not differentially expressed, and the adjusted p-value corrects for multiple testing using the Benjamini-Hochberg procedure.

The default results table compares the last level of the condition factor to the first level. You can specify the comparison explicitly using the contrast argument:

res <- results(dds, contrast = c("condition", "treated", "control"))

This command compares the treated group to the control group, where treated is the numerator and control is the denominator. The log2 fold change is positive for genes upregulated in the treated group and negative for genes downregulated in the treated group.

For experiments with more than two groups, you can extract results for each pairwise comparison using the contrast argument. You can also use the results function with the alpha argument to set the significance level for the adjusted p-value, which affects the summary statistics reported by the summary function.

summary(res)

The summary function provides a count of genes that are significantly differentially expressed at the chosen threshold, broken down by direction of change. This summary helps you assess whether the analysis produced a reasonable number of candidate genes before proceeding to downstream analysis.

Applying Shrinkage to Fold Change Estimates

The default results table contains log2 fold changes that have been shrunk toward zero using the empirical Bayes procedure implemented in DESeq2 [<a href="#ref-3">3</a>]. This shrinkage improves the stability and interpretability of fold change estimates, particularly for genes with low counts or high dispersion [<a href="#ref-3">3</a>]. The shrinkage is more pronounced for genes with high dispersion and less pronounced for genes with low dispersion and high counts.

You can access the shrunken fold changes using the lfcShrink function, which produces a results table with shrunken log2 fold changes and associated standard errors. The lfcShrink function requires you to specify the contrast or the coefficient of interest:

res_shrunk <- lfcShrink(dds, contrast = c("condition", "treated", "control"), type = "apeglm")

The apeglm type uses the adaptive shrinkage estimator, which provides more accurate fold change estimates for genes with low counts. Other options include ashr and normal, each with different assumptions about the prior distribution of fold changes. The choice of shrinkage method can affect the ranking of genes, so it is worth comparing the results from different methods for your dataset.

Shrunken fold changes are particularly useful for ranking genes and for downstream analyses such as gene set enrichment, where the magnitude of the fold change matters in addition to the statistical significance. The unshrunken fold changes from the results function are appropriate for reporting the observed effect size, while the shrunken fold changes are more suitable for ranking and for analyses that weight genes by effect size.

Quality Control and Diagnostic Checks

Quality control is an essential part of the differential expression workflow. DESeq2 provides several diagnostic plots that help you assess whether the data meet the assumptions of the statistical model and whether the experimental groups separate as expected [<a href="#ref-5">5</a>]. Dedicated pipelines such as SARTools emphasize the importance of systematic quality control steps to prevent spurious results from misusing the statistical methods [<a href="#ref-5">5</a>].

The dispersion plot is one of the first diagnostics to examine. This plot shows the dispersion estimates for each gene against the mean normalized count. The black points are the gene-wise dispersion estimates, the red points are the fitted dispersion values, and the blue points are the final dispersion estimates used in the analysis. If the blue points deviate substantially from the red curve, the dispersion fit may be poor, and the analysis may produce unreliable results.

plotDispEsts(dds)

The PCA plot is another useful diagnostic. This plot projects the samples onto the first two principal components of the variance-stabilized or regularized log-transformed data. Samples from the same experimental group should cluster together, and samples from different groups should separate along the principal components. If the samples do not separate according to the condition of interest, the experimental design may not have sufficient power, or there may be a dominant batch effect that needs to be modeled.

vsd <- vst(dds, blind = TRUE)
plotPCA(vsd, intgroup = "condition")

The sample distance heatmap provides a complementary view of sample relationships. This heatmap shows the Euclidean distance between samples based on the transformed data, with samples that are more similar to each other appearing closer together. The heatmap can reveal outliers that cluster far from their group, which may indicate sample quality issues or labeling errors.

sampleDists <- dist(t(assay(vsd)))
library(pheatmap)
pheatmap(as.matrix(sampleDists))

The MA plot is a standard visualization for differential expression results. This plot shows the log2 fold change against the mean normalized count for each gene, with significantly differentially expressed genes highlighted in red. The MA plot helps you assess whether the fold changes are consistent across the range of expression levels and whether the significant genes are concentrated in a particular expression range.

plotMA(res)

Handling Batch Effects and Confounding Factors

Batch effects are systematic technical variations that affect all samples processed in a particular batch, such as samples sequenced on different runs or prepared on different days. If batch effects are present and not accounted for in the model, they can inflate the number of false positives or obscure true biological differences.

The design formula is the primary tool for handling batch effects in DESeq2. By including the batch variable in the design formula, the model estimates the effect of the batch and adjusts the condition comparison accordingly:

dds <- DESeqDataSetFromMatrix(
    countData = count_matrix,
    colData = sample_metadata,
    design = ~ batch + condition
)

The order of terms in the design formula matters. The condition of interest should be the last term in the formula, and the batch or blocking factors should come before it. This ordering ensures that the differential expression test evaluates the condition effect after accounting for the batch effect.

Confounding occurs when the batch variable is completely confounded with the condition of interest, meaning that all samples in one condition come from one batch and all samples in the other condition come from another batch. In this situation, it is impossible to distinguish the biological effect from the batch effect, and no statistical method can separate them. The only solution is to design the experiment so that conditions are balanced across batches.

For experiments with paired samples, such as before-and-after treatment designs, the design formula can include the subject or patient identifier as a blocking factor:

dds <- DESeqDataSetFromMatrix(
    countData = count_matrix,
    colData = sample_metadata,
    design = ~ subject + condition
)

This paired design accounts for the baseline differences between subjects and increases the power to detect condition effects.

Filtering and Annotating Significant Genes

After extracting the results table, you will typically filter for genes that meet your significance thresholds. A common approach is to select genes with an adjusted p-value below 0.05 and an absolute log2 fold change above a threshold such as 1, which corresponds to a twofold change. The appropriate thresholds depend on your biological question and the expected magnitude of effects in your system.

sig_genes <- res[!is.na(res$padj) & res$padj < 0.05 & abs(res$log2FoldChange) > 1, ]

The filtering step removes genes with missing adjusted p-values, which occur when the gene was filtered out during the analysis due to low counts or extreme dispersion. The number of significant genes depends on the biological difference between the groups, the number of biological replicates, and the sequencing depth.

Published RNA-seq studies commonly report differential expression thresholds. For example, a study of intramuscular fat deposition in cattle used DESeq2 with a threshold of p-value less than 0.05 and absolute log2 fold change greater than 0.5 to identify significantly regulated genes [<a href="#ref-6">6</a>]. A study of formaldehyde treatment effects on chicken embryo liver microRNAs used a threshold of adjusted p-value less than 0.05 and absolute log2 fold change greater than 1 [<a href="#ref-7">7</a>]. These examples illustrate that the choice of thresholds varies by study and should be justified in the context of the biological question.

Once you have a list of significant genes, you will typically annotate the results with gene symbols, descriptions, and other biological information. The annotation can be added using a mapping file that links your gene identifiers to gene symbols and functional descriptions. The Bioconductor project provides annotation packages for many organisms that facilitate this mapping [<a href="#ref-1">1</a>].

library(AnnotationDbi)
library(org.Hs.eg.db)
res$symbol <- mapIds(org.Hs.eg.db,
                     keys = rownames(res),
                     column = "SYMBOL",
                     keytype = "ENSEMBL",
                     multiVals = "first")

The annotation step is essential for interpreting the biological meaning of the results and for preparing the gene list for downstream functional enrichment analysis.

Exporting Results for Downstream Analysis

The results table can be exported to a CSV or TSV file for sharing with collaborators or for use in other software. The write.csv or write.table functions in base R provide the simplest export options:

write.csv(as.data.frame(res), file = "deseq2_results.csv")

For downstream functional enrichment analysis, you will typically export the full results table with all genes, beyond the significant ones. Many enrichment tools accept a ranked gene list, where the ranking is based on the test statistic or the shrunken log2 fold change. The clusterProfiler package, which is available through Bioconductor, provides functions for over-representation analysis and functional class scoring that accept DESeq2 results as input [<a href="#ref-8">8</a>].

Functional enrichment analysis helps identify which biological pathways or gene ontology terms are overrepresented among the differentially expressed genes [<a href="#ref-8">8</a>]. This analysis uses knowledge databases such as Gene Ontology, KEGG, and Reactome to map genes to functional categories [<a href="#ref-8">8</a>]. The clusterProfiler package can access the latest KEGG knowledge database directly, and it provides visualization options including dot plots, enrichment maps, and ridge plots [<a href="#ref-8">8</a>].

The results of the differential expression analysis can also be used for more specialized analyses. For example, the shrunken fold changes can be used as input for gene set enrichment analysis, which tests whether predefined gene sets are coordinately upregulated or downregulated. The choice of downstream analysis depends on the biological question and the availability of appropriate gene set databases for your organism.

Reproducibility and Documentation

Reproducibility is a central concern in bioinformatics analysis. The same dataset analyzed with different parameter settings or different software versions can produce different results, so it is important to document the analysis environment and all parameter choices.

The SARTools pipeline emphasizes the importance of keeping track of the whole analysis process, including parameter values and versions of the R packages used [<a href="#ref-5">5</a>]. The HTML report generated by SARTools displays diagnostic plots for quality control and model hypothesis checking, and it records the parameter values and package versions [<a href="#ref-5">5</a>]. This level of documentation is essential for reproducing the analysis and for meeting the reporting requirements of many journals.

For your own analysis, you should record the following information:

• The version of R and Bioconductor used
• The version of DESeq2 and all other packages used
• The exact commands used to generate each result
• The thresholds applied for significance filtering
• The reference genome and annotation version used for read alignment
• The alignment and quantification tools used to generate the count matrix

The Carpentries lessons provide foundational training in reproducible computing practices, including version control with Git and the use of scripts to document analysis steps [<a href="#ref-2">2</a>]. Adopting these practices from the start of your analysis will save time and prevent errors when you need to revisit or share your work.

Workflow management systems can also support reproducibility. The nf-core project provides community-developed pipelines with standardized usage and configuration documentation [<a href="#ref-9">9</a>]. These pipelines follow best practices for reproducible analysis and can be adapted for RNA-seq differential expression workflows [<a href="#ref-9">9</a>]. The Galaxy Training Network offers accessible workflow training that emphasizes reproducibility and provides tutorials for RNA-seq analysis [<a href="#ref-10">10</a>].

Common Failure Patterns and Troubleshooting

Several common problems arise when running DESeq2 for the first time. Understanding these failure patterns can help you diagnose issues quickly and avoid wasting time on incorrect analyses.

The most common error is a mismatch between the column names of the count matrix and the row names of the sample metadata. DESeq2 will return an error message indicating that the row names of the metadata do not match the column names of the count data. This error is usually caused by inconsistent sample naming between the two files. Check that the sample identifiers are identical in both files, including any prefixes or suffixes.

Another common issue is the presence of non-integer values in the count matrix. DESeq2 expects raw integer counts, and the DESeqDataSetFromMatrix function will round non-integer values to the nearest integer with a warning. If your count matrix contains normalized values or transcript-per-million values, the analysis will produce incorrect results because the dispersion estimation assumes raw counts. Re-run the quantification step to obtain raw counts before proceeding.

A third common problem is a design formula that references a column that does not exist in the metadata. This error occurs when the column name in the design formula does not match any column name in the colData. Check the column names of your metadata using the colnames function and verify that the design formula uses the exact column names.

A fourth issue is the presence of samples with very low total counts. These samples may have failed library preparation or sequencing, and they can distort the size factor estimation. The PCA plot and sample distance heatmap will typically reveal such samples as outliers. If a sample has an extremely low total count, consider whether it should be removed from the analysis or whether the experimental design can accommodate it.

A fifth issue is the complete separation of groups in the PCA plot that does not match the condition of interest. This pattern may indicate a dominant batch effect or a technical artifact. If the samples separate by sequencing run or preparation date instead of by condition, you should include the batch variable in the design formula.

Limitations of DESeq2 and Alternative Approaches

DESeq2 is a well-established method for bulk RNA-seq differential expression analysis, but it has limitations that you should understand before applying it to your data. The method assumes that the count data follow a negative binomial distribution, which is a reasonable assumption for most bulk RNA-seq datasets but may not hold for all data types [<a href="#ref-3">3</a>].

For single-cell RNA-seq data, the assumptions of DESeq2 may be less appropriate because single-cell data have different characteristics, including higher dropout rates and greater technical variability [<a href="#ref-11">11</a>]. Best practices for single-cell RNA-seq analysis recommend specialized methods for normalization and differential expression that account for these characteristics [<a href="#ref-11">11</a>]. The DESeq2 package can be applied to pseudobulk counts from single-cell data, where counts are aggregated across cells within each sample, but the analysis should be designed with the single-cell data structure in mind.

Recent work has proposed alternative approaches to differential expression analysis that do not rely on distributional assumptions. One approach uses weighted averaging of transcript counts based on measured noise variances, which reduces both false-positive and false-negative findings without requiring parametrization of data distributions [<a href="#ref-12">12</a>]. This approach is closely related to statistics of cluster-randomized experiments and may produce more consistent differential gene expression estimates [<a href="#ref-12">12</a>]. While DESeq2 remains a standard tool, you should be aware that alternative methods exist and may be more appropriate for certain data types.

The choice of differential expression method can affect the results, and different methods may produce inconsistent findings even for high-quality samples [<a href="#ref-12">12</a>]. This inconsistency raises concerns about the assumptions behind widely accepted data analysis approaches [<a href="#ref-12">12</a>]. For this reason, it is good practice to validate key findings with an independent method or with experimental validation.

Interpreting Results in a Biological Context

The output of DESeq2 is a statistical list of genes with associated fold changes and p-values. The biological interpretation of this list requires additional analysis and domain knowledge. A gene that is statistically significant may not be biologically meaningful, and a gene that is not statistically significant may still be important if the effect is small or the statistical power is limited.

The magnitude of the fold change should be interpreted in the context of the biological system. A log2 fold change of 1 indicates a twofold change, which may be biologically meaningful for some genes but not for others. The biological relevance of a fold change depends on the gene function, the cell type, and the experimental context.

Functional enrichment analysis is a standard approach for interpreting the biological meaning of differential expression results. This analysis tests whether the differentially expressed genes are enriched for particular gene ontology terms, KEGG pathways, or Reactome pathways [<a href="#ref-8">8</a>]. The results of the enrichment analysis can identify the biological processes that are most affected by the experimental condition.

A study of intramuscular fat deposition in cattle used DESeq2 to identify differentially expressed genes and then performed gene ontology and KEGG pathway analysis to link these genes to lipid metabolism and intramuscular fat-associated pathways such as PPAR signaling, fatty acid metabolism, and PI3K-Akt signaling [<a href="#ref-6">6</a>]. This example illustrates how differential expression results can be connected to biological mechanisms through pathway analysis.

A study of formaldehyde treatment effects on chicken embryo liver microRNAs used DESeq2 to identify differentially expressed microRNAs and then predicted target genes and performed pathway enrichment analysis [<a href="#ref-7">7</a>]. The pathway analysis revealed effects on signaling pathways, post-translational protein modification, immune system, and oxidative stress pathways [<a href="#ref-7">7</a>]. This example shows how differential expression analysis can be extended to non-coding RNAs and integrated with target prediction.

Reporting Results for Publication

When reporting DESeq2 results in a publication, you should provide sufficient detail for readers to evaluate the analysis and reproduce the results. The methods section should describe the following:

• The version of DESeq2 and R used
• The reference genome and annotation used for read alignment
• The alignment and quantification tools used to generate the count matrix
• The design formula used for the analysis
• The significance thresholds applied
• The number of biological replicates per condition
• The total number of genes tested and the number of significant genes

The results section should report the number of differentially expressed genes and the direction of change. A common format is to state that X genes were significantly upregulated and Y genes were significantly downregulated in the treatment group compared to the control group, based on an adjusted p-value threshold and a fold change threshold.

The supplementary materials should include the full results table with all genes and their statistics. This table allows readers to examine the results for specific genes of interest and to perform their own downstream analyses.

The SARTools pipeline provides a template for reporting that includes diagnostic plots and parameter tracking [<a href="#ref-5">5</a>]. Adopting a similar approach in your own analysis will improve the transparency and reproducibility of your reporting.

Professional Escalation Criteria

Some analysis situations require consultation with a bioinformatics specialist or a statistician. You should escalate the analysis if you encounter any of the following situations:

• The dispersion plot shows a poor fit, with the final dispersion estimates deviating substantially from the fitted curve
• The PCA plot shows no separation between experimental groups, suggesting that the experimental design may not have sufficient power
• The PCA plot shows separation that does not match the condition of interest, suggesting a dominant batch effect or technical artifact
• The number of significant genes is unexpectedly high or low, suggesting a problem with the data or the analysis
• The results are highly sensitive to small changes in the analysis parameters, suggesting instability in the estimates
• The experimental design has complete confounding between the condition and a technical variable, making it impossible to separate biological and technical effects

A bioinformatics specialist can help with more complex experimental designs, with the integration of multiple data types, or with the application of alternative statistical methods when DESeq2 assumptions are not met.

Frequently Asked Questions

What is the difference between raw counts and normalized counts in DESeq2?

Raw counts are the integer numbers of sequencing reads that aligned to each gene in each sample. DESeq2 uses raw counts as input because the statistical model is based on the negative binomial distribution, which is appropriate for count data [<a href="#ref-3">3</a>]. Normalized counts are raw counts adjusted for sequencing depth and other technical factors. DESeq2 estimates size factors to account for differences in sequencing depth between samples, and these size factors are used to normalize the counts for visualization and for the statistical analysis. You should not input normalized counts into DESeq2 because the dispersion estimation assumes raw counts.

How many biological replicates do I need for DESeq2?

DESeq2 can analyze experiments with as few as two biological replicates per condition, but the statistical power increases with the number of replicates [<a href="#ref-3">3</a>]. The shrinkage estimation of dispersion and fold changes in DESeq2 is designed to provide stable estimates even with small replicate numbers [<a href="#ref-3">3</a>]. For most experiments, three biological replicates per condition is a common minimum, and more replicates are recommended when the expected effect size is small or the biological variability is high. The appropriate number of replicates depends on the biological system and the magnitude of the expected differences.

What design formula should I use for my experiment?

The design formula should include the condition of interest as the last term and any blocking factors or batch variables before it. For a simple two-group comparison, use ~ condition. For experiments with batch effects, use ~ batch + condition. For paired samples, use ~ subject + condition. The condition of interest must be the last term in the formula because DESeq2 tests the last term for differential expression. If you have multiple factors of interest, you may need to use a more complex design or perform separate analyses.

How do I choose the significance threshold for adjusted p-value and fold change?

The choice of thresholds depends on your biological question and the expected magnitude of effects. A common threshold is an adjusted p-value below 0.05 and an absolute log2 fold change above 1, which corresponds to a twofold change. Published studies have used various thresholds, including adjusted p-value below 0.05 with absolute log2 fold change above 1 [<a href="#ref-7">7</a>] and p-value below 0.05 with absolute log2 fold change above 0.5 [<a href="#ref-6">6</a>]. The thresholds should be justified in the context of your study and should be specified before the analysis to avoid bias.

What is the difference between the Wald test and the likelihood ratio test in DESeq2?

The Wald test is the default test in DESeq2 and is used to compare the expression level of each gene between two groups or between each level of a factor and a reference level [<a href="#ref-3">3</a>]. The likelihood ratio test is used when you want to test whether any level of a factor differs from the others, such as in experiments with more than two groups. The likelihood ratio test compares a full model to a reduced model and tests whether the additional terms in the full model significantly improve the fit. The choice of test depends on the experimental design and the specific question being asked.

Why are my log2 fold changes different between the results function and the lfcShrink function?

The results function returns log2 fold changes that have been shrunk using the empirical Bayes procedure in DESeq2, but the shrinkage is applied to the dispersion estimates and the fold changes are then calculated from the fitted model [<a href="#ref-3">3</a>]. The lfcShrink function applies additional shrinkage to the fold changes using a specified method, such as apeglm or ashr. The additional shrinkage in lfcShrink produces fold changes that are more conservative for genes with low counts or high dispersion, which improves the ranking of genes for downstream analysis. The unshrunken fold changes from the results function are appropriate for reporting the observed effect size, while the shrunken fold changes are more suitable for ranking and for analyses that weight genes by effect size.

Can I use DESeq2 for single-cell RNA-seq data?

DESeq2 is designed for bulk RNA-seq data, and its assumptions may not be appropriate for single-cell RNA-seq data, which have different characteristics including higher dropout rates and greater technical variability [<a href="#ref-11">11</a>]. Best practices for single-cell RNA-seq analysis recommend specialized methods for normalization and differential expression [<a href="#ref-11">11</a>]. However, DESeq2 can be applied to pseudobulk counts from single-cell data, where counts are aggregated across cells within each sample. This approach treats each sample as a biological replicate and applies the bulk RNA-seq analysis framework to the aggregated counts.

What should I do if my samples do not separate by condition in the PCA plot?

If the PCA plot shows no separation between experimental groups, the experimental design may not have sufficient power, or there may be a dominant batch effect. First, check whether the samples separate by any technical variable such as sequencing run or preparation date. If so, include the batch variable in the design formula. If the samples still do not separate, the biological difference between the groups may be small relative to the technical variability, and you may need more biological replicates or a more sensitive experimental design. Consider consulting a bioinformatics specialist if the lack of separation persists.

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 Data Analysis Workflow: From Raw Reads to Insights
• RNA-Seq Differential Expression: DESeq2, edgeR, and limma-voom Frameworks
• Metabolomics Data Analysis in R: A Practical Workflow

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

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