What is Dr. Zubair Khalid's research focus?

Dr. Zubair Khalid specializes in molecular virology, mRNA vaccine development, and computational biology, with a focus on avian pathogens like IBDV and Avian Reovirus.

Where is Dr. Zubair Khalid currently working?

Dr. Zubair Khalid is a Postdoctoral Research Associate at the University of Maryland (UMD), specifically within the Department of Animal and Avian Sciences.

Differential Gene Expression Analysis with DESeq2: A Practical Guide

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

Introduction to Differential Gene Expression Analysis with DESeq2

What is Differential Gene Expression?

Differential gene expression (DGE) analysis is the statistical process of identifying genes whose transcript abundance changes significantly between two or more experimental conditions. The fundamental question is straightforward: which genes are transcribed at different levels in treated versus untreated cells, in diseased versus healthy tissue, or across a developmental time course? The answer, however, requires careful statistical modeling because RNA-seq data are count-based, overdispersed, and subject to substantial technical variation.

The biological rationale for DGE analysis rests on the central dogma: messenger RNA (mRNA) abundance is the primary readout of gene regulatory state. When a transcription factor such as p53 is activated, it binds to response elements and recruits RNA polymerase II, leading to increased transcription of target genes like CDKN1A (p21). Conversely, DNA methylation at promoter CpG islands can recruit methyl-CpG-binding domain proteins that compact chromatin and reduce transcription initiation, a mechanism detailed in DNA Methylation Decrease Gene Expression. DGE analysis quantifies these regulatory outcomes genome-wide, providing a snapshot of the transcriptional program under defined conditions.

RNA-seq measures gene expression by counting the number of sequencing reads that map to each gene's genomic coordinates. These counts are integer-valued and range from zero to hundreds of thousands. The dynamic range is enormous: a highly expressed gene like ACTB (beta-actin) may yield tens of thousands of counts per sample, while a low-abundance transcription factor may yield fewer than ten. This count-based nature distinguishes RNA-seq from microarray data, which produces continuous fluorescence intensities, and it dictates the statistical models used for analysis.

Why DESeq2?

DESeq2 is a Bioconductor R package developed by Michael Love, Simon Anders, and Wolfgang Huber for differential analysis of count-based expression data. It has become one of the most widely used tools for RNA-seq DGE analysis, alongside edgeR and limma-voom. DESeq2's popularity stems from several design features: it models count data with a negative binomial distribution, it implements a robust normalization procedure that accounts for library size and composition bias, and it applies empirical Bayes shrinkage to stabilize dispersion and fold change estimates, particularly for genes with low counts.

The package accepts a simple input—a count matrix and a metadata table—and produces a results table with log2 fold changes, p-values, and adjusted p-values for every gene. It handles complex experimental designs including multifactorial comparisons, interactions, and continuous covariates. For a broader overview of the DGE analysis workflow, see Differential Gene Expression Dge Analysis. This article provides a mechanistic guide to DESeq2: how it normalizes data, how it models count distributions, how it estimates dispersion, and how it performs statistical testing. Understanding these mechanisms is essential for interpreting results correctly and troubleshooting when analyses go wrong.

Input Data Requirements and Preparation

Count Matrix Generation

The starting point for DESeq2 is a matrix of integer counts, where rows represent genes and columns represent samples. Each entry $ c_{ij} $ is the number of sequencing reads that aligned to gene $ i $ in sample $ j $. This matrix can be generated from raw sequencing data through a standard pipeline: quality control, read alignment, and quantification.

Raw FASTQ files are first assessed for quality using tools like FastQC. Adapter sequences and low-quality bases are trimmed with Trimmomatic or cutadapt. The cleaned reads are then aligned to a reference genome or transcriptome. For model organisms with well-annotated genomes, splice-aware aligners such as STAR or HISAT2 are standard. STAR, for example, aligns reads by first building a genome index with --runMode genomeGenerate, then mapping reads with --runMode alignReads. Typical parameters include --outSAMtype BAM SortedByCoordinate and --quantMode GeneCounts, the latter producing per-gene counts directly.

Alternatively, pseudo-alignment tools like Salmon or kallisto quantify transcript abundance without full alignment. These tools are substantially faster than splice-aware aligners and produce transcript-level counts that can be summarized to gene-level counts using the R package tximport. For example, Salmon quantifies transcripts using a dual-phase inference procedure that models sequence-specific biases and fragment length distributions, outputting estimated counts and transcript-per-million (TPM) values. The tximport function then aggregates transcript counts to gene counts, handling the mapping between transcript and gene identifiers.

Regardless of the quantification method, the final count matrix must contain raw integer counts, not normalized values like TPM or reads per kilobase per million (RPKM). DESeq2 performs its own normalization internally, and applying external normalization beforehand will distort the statistical model. The count matrix should be saved as a plain text file (tab-separated) or as an R object, with gene identifiers in the first column and sample names as column headers.

Sample Metadata and Experimental Design

The second required input is a metadata table, or column data, that describes each sample. This table must have one row per sample, with the sample names matching the column names of the count matrix. The metadata contains the experimental variables—condition, genotype, treatment, time point, batch—that will be used in the statistical model.

The experimental design is encoded in a design formula, which specifies how the variation in gene expression is attributed to the experimental variables. For a simple two-group comparison, the design formula is ~ condition, where condition is a factor with two levels, such as "control" and "treated". For more complex designs, the formula can include multiple factors and interactions, such as ~ genotype + treatment or ~ genotype + treatment + genotype:treatment.

The design formula determines what comparisons can be made and what variation is accounted for in the model. Including blocking factors like batch or sex in the design formula is critical because it reduces unwanted variation and increases statistical power. For example, if samples were processed in two batches, the design ~ batch + condition will estimate the condition effect while accounting for batch differences. Failing to include such factors can lead to spurious results if batch effects correlate with the condition of interest.

DESeq2 constructs a design matrix from the design formula using R's model matrix functionality. For a factor with $ k $ levels, the model matrix contains $ k-1 $ indicator variables, with one level serving as the reference or baseline. The choice of reference level affects the interpretation of coefficients but not the overall statistical conclusions. The reference level can be set using the relevel function in R.

Normalization and the DESeq2 Model

Median-of-Ratios Normalization

RNA-seq count data require normalization to account for two sources of technical variation: library size (total number of sequenced reads per sample) and composition bias (differences in the distribution of read counts across genes). DESeq2 implements the median-of-ratios method, which estimates a size factor for each sample that represents the relative sequencing depth.

The algorithm proceeds as follows. For each gene $ i $ and sample $ j $, the ratio of the observed count to the geometric mean of counts across all samples is computed:

\[ r_{ij} = \frac{c_{ij}}{\left( \prod_{k=1}^{n} c_{ik} \right)^{1/n}} \]

where $ n $ is the number of samples. The size factor $ s_j $ for sample $ j $ is the median of these ratios across all genes:

\[ s_j = \text{median}_{i} \left( r_{ij} \right) \]

The median is used rather than the mean because it is robust to the large number of genes with zero or very low counts, and to genes that are genuinely differentially expressed. If a gene is strongly upregulated in one sample, its ratio will be large, but the median across all genes remains stable because the majority of genes are not differentially expressed.

The size factor $ s_j $ represents the relative sequencing depth of sample $ j $. A sample with twice as many total reads as the median sample will have a size factor of approximately 2. The normalized count for gene $ i $ in sample $ j $ is then $ c_{ij} / s_j $. These normalized counts are used for visualization and quality control, but the statistical model operates on the raw counts with the size factors incorporated as offsets.

This normalization approach is compositionally robust. Consider a scenario where one gene is massively overexpressed in treated samples, consuming a large fraction of the sequencing reads. This would depress the counts of all other genes in those samples. The median-of-ratios method corrects for this because the size factor for the treated samples will be inflated, and dividing by this larger size factor restores the other genes to their true relative abundances.

Negative Binomial Distribution

RNA-seq counts are integer-valued and exhibit overdispersion relative to the Poisson distribution. For a Poisson distribution, the variance equals the mean. In practice, RNA-seq data show variance substantially larger than the mean, particularly for highly expressed genes. This overdispersion arises from biological variability between replicates (different cells within a sample, different individuals) and technical variability beyond Poisson sampling.

DESeq2 models the count $ c_{ij} $ for gene $ i $ in sample $ j $ with a negative binomial (NB) distribution:

\[ c_{ij} \sim \text{NB}(\mu_{ij}, \alpha_i) \]

where $ \mu_{ij} $ is the mean and $ \alpha_i $ is the dispersion parameter for gene $ i $. The NB distribution has variance:

\[ \text{Var}(c_{ij}) = \mu_{ij} + \alpha_i \mu_{ij}^2 \]

The dispersion parameter $ \alpha_i $ quantifies the extra variance beyond the Poisson mean. When $ \alpha_i = 0 $, the NB reduces to the Poisson distribution. Typical RNA-seq dispersions range from 0.01 to 1, with highly variable genes having larger dispersions.

The mean $ \mu_{ij} $ is modeled as:

\[ \mu_{ij} = s_j q_{ij} \]

where $ s_j $ is the size factor for sample $ j $ and $ q_{ij} $ is the true expression level of gene $ i $ in sample $ j $, proportional to the concentration of mRNA in the sample. The size factor acts as an offset, accounting for library size differences.

The NB model is parameterized in terms of the log link function, so that:

\[ \log(\mu_{ij}) = \log(s_j) + \sum_{k} x_{jk} \beta_{ik} \]

where $ x_{jk} $ are the elements of the design matrix and $ \beta_{ik} $ are the coefficients for gene $ i $. The coefficients represent the log2 fold changes associated with each experimental variable. For a simple two-group comparison with design ~ condition, the model has two coefficients: the intercept $ \beta_{i0} $ (the log2 expression in the reference condition) and $ \beta_{i1} $ (the log2 fold change between conditions).

Dispersion Estimation and Shrinkage

Gene-Wise Dispersion

The dispersion parameter $ \alpha_i $ must be estimated from the data for each gene. The maximum likelihood estimate (MLE) of dispersion is obtained by maximizing the NB likelihood for each gene independently. For a gene with counts $ c_{i1}, c_{i2}, \ldots, c_{in} $ across $ n $ samples, the log-likelihood is:

\[ \ell(\alpha_i, \beta_i) = \sum_{j=1}^{n} \log \left[ \frac{\Gamma(c_{ij} + 1/\alpha_i)}{\Gamma(c_{ij} + 1) \Gamma(1/\alpha_i)} \left( \frac{1}{1 + \alpha_i \mu_{ij}} \right)^{1/\alpha_i} \left( \frac{\alpha_i \mu_{ij}}{1 + \alpha_i \mu_{ij}} \right)^{c_{ij}} \right] \]

The MLE is found by numerical optimization, typically using the BFGS algorithm. However, for genes with low counts, the MLE is unstable and can be severely underestimated or overestimated. A gene with counts of 0, 0, 0, 10 across four samples may have an MLE dispersion near zero, suggesting no overdispersion, when in fact the data are highly variable. Conversely, a gene with counts of 0, 0, 0, 1 may have an inflated MLE.

Trend Fitting

To stabilize dispersion estimates, DESeq2 fits a trend to the gene-wise dispersion estimates as a function of the mean normalized count. The trend captures the relationship that, on average, dispersion decreases with increasing expression. Highly expressed genes tend to have lower dispersion because the Poisson sampling noise is proportionally smaller and the biological variability is more consistent.

The trend is fitted using a smooth function, implemented as a gamma-family generalized additive model (GAM) or a local regression. The fitted trend provides a prior expectation for the dispersion of each gene based on its mean expression. For a gene with mean count of 100, the trend might predict a dispersion of 0.05; for a gene with mean count of 10,000, the trend might predict 0.01.

The trend fitting is performed on the log scale, and the relationship is approximately linear for most datasets. The strength of the trend depends on the biological variability in the dataset. In datasets with highly variable biological replicates, the trend is flatter; in datasets with tight replicates, the trend declines more steeply with mean expression.

Shrinkage to the Trend

The final dispersion estimate for each gene is a shrinkage of the gene-wise MLE toward the trend value. DESeq2 uses an empirical Bayes approach where the trend serves as the prior and the gene-wise estimate provides the likelihood. The posterior estimate is:

\[ \alpha_i^{\text{shrunken}} = \frac{d_0 \alpha_i^{\text{trend}} + d_i \alpha_i^{\text{MLE}}}{d_0 + d_i} \]

where $ d_0 $ is the prior degrees of freedom (a measure of the strength of the prior) and $ d_i $ is the residual degrees of freedom for gene $ i $, which equals the number of samples minus the number of model parameters. The prior degrees of freedom $ d_0 $ is estimated from the data by maximizing the likelihood of the observed dispersion estimates given the trend.

The shrinkage has a profound effect. For genes with high counts, the gene-wise MLE is well-determined and the shrinkage is minimal; the shrunken estimate is close to the MLE. For genes with low counts, the MLE is noisy and the shrinkage pulls the estimate strongly toward the trend. This reduces the number of false positives among low-count genes, which would otherwise show artificially low dispersion and hence artificially significant p-values.

The same shrinkage procedure is applied to the log2 fold changes. After fitting the model with the shrunken dispersions, DESeq2 applies empirical Bayes shrinkage to the fold change estimates. This is implemented in the lfcShrink function, which produces shrunken log2 fold changes that are more reliable for ranking genes, particularly for genes with low counts and high dispersion. The shrinkage uses a normal prior centered at zero, with the width determined by the distribution of observed fold changes.

Statistical Testing and Multiple Testing Correction

Wald Test

The Wald test is the default testing procedure in DESeq2 for comparing two conditions. After fitting the model and obtaining the maximum likelihood estimates of the coefficients $ \hat{\beta}_{ik} $, the Wald test evaluates whether a coefficient is significantly different from zero.

For a coefficient $ \hat{\beta}_{ik} $, the Wald statistic is:

\[ W = \frac{\hat{\beta}_{ik}}{\text{SE}(\hat{\beta}_{ik})} \]

where $ \text{SE}(\hat{\beta}_{ik}) $ is the standard error of the coefficient estimate, obtained from the observed Fisher information matrix. Under the null hypothesis that $ \beta_{ik} = 0 $, the Wald statistic follows a standard normal distribution. The p-value is computed as:

\[ p = 2 \left[ 1 - \Phi(|W|) \right] \]

where $ \Phi $ is the cumulative distribution function of the standard normal distribution.

The Wald test is used when the comparison is specified in the results call, such as results(dds, contrast = c("condition", "treated", "control")). The contrast specifies the numerator and denominator of the log2 fold change. The Wald test is appropriate for pairwise comparisons and for testing individual coefficients.

Likelihood Ratio Test

The likelihood ratio test (LRT) is used for comparing nested models, which is useful for multi-level factors or when testing the significance of a group of coefficients simultaneously. The LRT compares a full model (containing all terms) to a reduced model (with the terms of interest removed).

For example, consider a time course experiment with four time points (0h, 6h, 12h, 24h). The full model is ~ time, and the reduced model is ~ 1 (intercept only). The LRT tests whether any time point differs from the others. The test statistic is:

\[ D = 2 \left( \ell_{\text{full}} - \ell_{\text{reduced}} \right) \]

where $ \ell_{\text{full}} $ and $ \ell_{\text{reduced}} $ are the maximized log-likelihoods of the full and reduced models, respectively. Under the null hypothesis, $ D $ follows a chi-squared distribution with degrees of freedom equal to the difference in the number of parameters between the models. For the time course example with four levels, the difference is 3 degrees of freedom.

The LRT is specified in DESeq2 by providing the reduced design to the results function: results(dds, test = "LRT", reduced = ~ 1). The LRT is more appropriate than the Wald test when testing multiple coefficients jointly, such as the effect of a factor with more than two levels or the interaction between two factors.

Multiple Testing Correction

When testing thousands of genes simultaneously, the probability of false positives is substantial. If 20,000 genes are tested at a significance threshold of $ \alpha = 0.05 $, approximately 1,000 genes would be expected to show significant p-values by chance alone, even if no genes are truly differentially expressed.

DESeq2 controls the false discovery rate (FDR) using the Benjamini-Hochberg (BH) procedure. The BH procedure works as follows:

Sort all p-values in ascending order: $ p_{(1)} \le p_{(2)} \le \ldots \le p_{(m)} $, where $ m $ is the total number of genes tested.
For each rank $ i $, compute the adjusted p-value: \[ p_{\text{adj}(i)} = \min_{k \ge i} \left( \frac{m}{k} p_{(k)} \right) \]
A gene is declared significant if its adjusted p-value is below the chosen threshold (typically 0.05 or 0.01).

The BH procedure ensures that the expected proportion of false positives among the rejected hypotheses is controlled at the specified level. This is a less stringent criterion than the Bonferroni correction (which controls the family-wise error rate) and is more appropriate for exploratory genomic studies where the goal is to identify a candidate list of genes for further investigation.

DESeq2 reports the BH-adjusted p-values in the padj column of the results table. These adjusted p-values should be used for significance filtering, not the raw p-values.

Interpreting DESeq2 Results

Results Table Columns

The results function in DESeq2 returns a DataFrame with one row per gene and the following columns:

baseMean: The mean of normalized counts across all samples, averaged over the size factor. This is a measure of overall expression level.
log2FoldChange: The log2 of the fold change between the conditions compared. A value of 1 means a 2-fold increase; -1 means a 2-fold decrease.
lfcSE: The standard error of the log2 fold change estimate.
stat: The Wald statistic (for Wald test) or the likelihood ratio statistic (for LRT).
pvalue: The raw p-value from the statistical test.
padj: The Benjamini-Hochberg adjusted p-value.

The baseMean is useful for filtering low-expression genes. Genes with very low baseMean (e.g., below 10) have unreliable fold change estimates and are often filtered out before downstream analysis. The log2FoldChange is the effect size, and its sign indicates the direction of change: positive means higher expression in the numerator condition (e.g., treated), negative means higher expression in the denominator condition (e.g., control).

Significance Thresholds

The standard significance threshold for calling differentially expressed genes is padj < 0.05. However, the choice of threshold depends on the biological question and the downstream validation strategy. For exploratory analyses where the goal is to generate hypotheses, a more lenient threshold (e.g., padj < 0.1) may be appropriate. For studies where the goal is to identify a small set of high-confidence genes for experimental validation, a stricter threshold (e.g., padj < 0.01 and |log2FoldChange| > 1) is recommended.

The fold change threshold is equally important. A gene with a log2 fold change of 0.1 (a 7% increase) may be statistically significant with a large sample size, but biologically irrelevant. Conversely, a gene with a log2 fold change of 3 (an 8-fold increase) may be biologically important even if the p-value is marginal due to high variability. The biological relevance of the fold change depends on the gene and the context. For transcription factors and signaling molecules, even 2-fold changes can be functionally significant; for housekeeping genes, 2-fold changes are unusual and likely pathological.

Visualization of Results

Several visualization approaches are standard for DESeq2 results. The MA plot displays the log2 fold change (y-axis) against the mean normalized count (x-axis). Each gene is a point, and significant genes are colored. The plot shows the distribution of fold changes and the relationship between expression level and effect size. In DESeq2, the plotMA function generates this plot directly from the results object.

The volcano plot displays the log2 fold change (x-axis) against the negative log10 of the p-value (y-axis). This plot emphasizes both statistical significance and effect size, with significant genes appearing in the upper left (downregulated) and upper right (upregulated) corners. The volcano plot is useful for identifying genes with large fold changes and high significance.

Heatmaps of the top differentially expressed genes, with sample annotations, provide a visual summary of the expression patterns. The pheatmap or ComplexHeatmap R packages are commonly used. The heatmap should use variance-stabilized or regularized-log transformed counts, not raw counts, to avoid the dominance of highly expressed genes. The vst and rlog functions in DESeq2 provide these transformations.

For downstream functional interpretation, the list of significant genes can be subjected to Gene Ontology Analysis Online or Gene Ontology Pathway Enrichment to identify overrepresented biological processes and pathways. These analyses are typically performed with tools like clusterProfiler or DAVID, and they provide biological context for the observed expression changes.

Common Pitfalls and Troubleshooting

Design Formula Mistakes

The design formula is the most common source of errors in DESeq2 analyses. A common mistake is using a design formula that does not include all relevant sources of variation, such as batch or sex. If samples were processed in two batches and the batch is not included in the design, the batch effect will be confounded with the condition effect, potentially producing false positives or masking true differences.

Another mistake is using a continuous variable as a factor without specifying it as such. If a time variable is encoded as numeric (0, 6, 12, 24), DESeq2 will treat it as a continuous covariate and model a linear trend. If the relationship is nonlinear, this will miss the true pattern. The variable should be converted to a factor using factor(time) in the metadata.

A related issue is the reference level. The default reference level for a factor is alphabetical, which may not be the desired baseline. For example, if the condition factor has levels "control" and "treated", the reference will be "control" (alphabetically first). If the levels are "mutant" and "wildtype", the reference will be "mutant". The reference level should be set explicitly using relevel(metadata$condition, ref = "control") to ensure the correct comparison.

Handling Outliers

Outliers are samples that show extreme expression values for a gene, far outside the range of other samples in the same condition. Outliers can arise from technical issues (sample contamination, mapping errors) or biological causes (a single cell with an extreme expression state). DESeq2 has built-in outlier detection: for each gene, it computes the Cook's distance for each sample, which measures the influence of that sample on the model fit. Samples with Cook's distance greater than the 0.99 quantile of the F distribution are flagged as outliers.

By default, DESeq2 replaces the count for an outlier sample with the predicted value from the model, and adds a flag in the results table. This prevents the outlier from dominating the fold change estimate. However, if more than half of the samples for a gene are flagged as outliers, the gene is removed from the results table (the p-value and padj are set to NA).

If outliers are prevalent, it is worth investigating the cause. Check the alignment metrics for the flagged samples, examine the PCA plot for sample clustering, and consider whether the outlier reflects a genuine biological phenomenon. If the outliers are technical, the samples may need to be removed or the experiment repeated. If they are biological, they may indicate a subpopulation of cells with distinct expression states, which could be analyzed separately.

Overinterpretation of Small Fold Changes

A common error is to interpret any statistically significant gene as biologically important, regardless of the fold change magnitude. With sufficient sample size, DESeq2 can detect fold changes as small as 10-20% as statistically significant. These small changes may be reproducible and statistically robust, but they may not be biologically meaningful.

The biological relevance of a fold change depends on the gene and the context. For example, a 1.5-fold change in a metabolic enzyme may have minimal impact on flux through a pathway, while a 1.5-fold change in a transcription factor could have substantial downstream effects. The fold change should be interpreted in the context of the gene's function, the pathway it participates in, and the expected magnitude of change based on the experimental perturbation.

The shrunken log2 fold changes from lfcShrink are more reliable for ranking genes than the raw estimates, because the shrinkage reduces the exaggerated fold changes that occur for low-count genes. When filtering genes for downstream analysis, consider using both the adjusted p-value and an absolute fold change threshold, such as |log2FoldChange| > 1 (2-fold change). This dual filtering reduces the number of biologically irrelevant hits.

Practical Workflow Summary

Step-by-Step Analysis Pipeline

A typical DESeq2 analysis follows a standardized workflow. The steps are enumerated below, with the key R commands provided.

Step 1: Prepare the count matrix and metadata. The count matrix should have genes as rows and samples as columns. The metadata should have samples as rows and experimental variables as columns. Ensure that the sample names match exactly between the two tables.

Step 2: Create the DESeqDataSet object. The DESeqDataSetFromMatrix function combines the count matrix, metadata, and design formula into a single object:

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                              colData = metadata,
                              design = ~ condition)

Step 3: Pre-filter low-count genes. Remove genes with very low total counts to reduce the number of tests and improve the dispersion estimation:

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

Step 4: Run the DESeq2 pipeline. The DESeq function performs normalization, dispersion estimation, and fitting in one step:

dds <- DESeq(dds)

Step 5: Extract results. The results function extracts the results table for the comparison of interest:

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

For shrunken fold changes:

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

Step 6: Filter significant genes. Apply thresholds for adjusted p-value and fold change:

sig_genes <- subset(res_shrunk, padj < 0.05 & abs(log2FoldChange) > 1)

Step 7: Visualize results. Generate MA and volcano plots, and heatmaps of the top genes.

Step 8: Export results. Write the results table to a file for downstream analysis:

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

Key R Commands

The essential DESeq2 commands are summarized in the table below:

Function	Purpose
`DESeqDataSetFromMatrix`	Create the DESeqDataSet object from count matrix and metadata
`DESeq`	Run the full analysis pipeline (normalization, dispersion, fitting)
`results`	Extract the results table for a specified comparison
`lfcShrink`	Apply empirical Bayes shrinkage to log2 fold changes
`plotMA`	Generate an MA plot of log2 fold change vs. mean expression
`vst` / `rlog`	Apply variance-stabilizing transformation for visualization
`plotPCA`	Generate a PCA plot of sample relationships

For a more detailed walkthrough of the R implementation, see Differential Gene Expression Analysis in R. The complete workflow from count matrix to significant gene list typically takes less than an hour for a standard dataset of 20,000 genes and 12 samples.

Frequently Asked Questions

What is the difference between DESeq2 and edgeR?

Both DESeq2 and edgeR use the negative binomial distribution to model RNA-seq counts, but they differ in normalization and dispersion estimation. edgeR uses the trimmed mean of M-values (TMM) normalization, while DESeq2 uses median-of-ratios. Both methods are compositionally robust, but they can produce slightly different size factors for the same dataset. For dispersion, edgeR uses a weighted likelihood approach with a trend, while DESeq2 uses empirical Bayes shrinkage with a trend. In practice, the two methods produce highly concordant results for most datasets. DESeq2 is often preferred for its ease of use and built-in visualization functions, while edgeR offers more flexibility in some advanced designs. The choice between them is often based on personal preference and the specific features of the dataset.

How many biological replicates are needed for DESeq2?

DESeq2 requires at least two biological replicates per condition to estimate dispersion. With two replicates, the dispersion estimates are very noisy, and the analysis will have low power to detect small fold changes. Three replicates per condition is the minimum recommended for reliable results, and four to six replicates are preferable for detecting subtle changes. The number of replicates needed depends on the biological variability of the system and the magnitude of the expected fold changes. Highly variable systems, such as primary cells from different individuals, require more replicates than tightly controlled systems, such as cell lines. Increasing the number of replicates improves power more than increasing sequencing depth beyond approximately 20-30 million reads per sample.

What does the DESeq2 size factor represent?

The size factor for a sample represents its relative sequencing depth, adjusted for composition bias. A size factor of 1.5 means that the sample has approximately 1.5 times the sequencing depth of the median sample, after accounting for the fact that some genes may be differentially expressed. The size factor is used as an offset in the negative binomial model, so that the mean expression for a gene is proportional to the size factor times the true expression level. The size factors are estimated by the median-of-ratios method, which is robust to the presence of differentially expressed genes. The size factors can be accessed with the sizeFactors function and should be inspected to ensure that no sample has an extreme value, which would indicate a problem with library preparation or sequencing.

Why are p-values adjusted in DESeq2?

P-values are adjusted to control the false discovery rate (FDR) when testing thousands of genes simultaneously. If 20,000 genes are tested at a raw p-value threshold of 0.05, approximately 1,000 genes would be expected to show significant p-values by chance alone. The Benjamini-Hochberg adjustment controls the expected proportion of false positives among the genes declared significant. An adjusted p-value of 0.05 means that approximately 5% of the genes declared significant are expected to be false positives. This is a more appropriate criterion for exploratory genomic studies than the Bonferroni correction, which controls the family-wise error rate and is overly conservative for thousands of tests. The adjusted p-values should always be used for significance filtering, not the raw p-values.

Can DESeq2 be used for non-model organisms?

Yes, DESeq2 can be used for any organism with RNA-seq data, provided that a reference genome or transcriptome is available for read alignment and quantification. The method does not require any organism-specific information; it operates on the count matrix and metadata. For non-model organisms without a reference genome, de novo transcriptome assembly can be performed with tools like Trinity, and the resulting transcripts can be quantified with Salmon or kallisto. The count matrix can then be analyzed with DESeq2. The main challenge for non-model organisms is the quality of the transcriptome assembly and the annotation of gene functions for downstream interpretation. For functional analysis, Gene Ontology Analysis Tool can be used if Gene Ontology annotations are available for the organism.

What is the effect of outliers on DESeq2 results?

Outliers can have a substantial effect on DESeq2 results if not handled properly. A single outlier sample with an extreme count for a gene can inflate the estimated dispersion and reduce the significance of that gene, or it can produce a spurious fold change if the outlier is in one condition. DESeq2 detects outliers using Cook's distance and replaces the outlier count with the predicted value from the model, preventing the outlier from dominating the analysis. However, if the outlier is not detected (e.g., because it is not extreme enough), it can still affect the results. It is good practice to examine the PCA plot and the Cook's distances for the samples, and to investigate any samples that cluster separately from their biological replicates. If outliers are technical, the samples should be removed; if they are biological, they may indicate a subpopulation that should be analyzed separately.

How do I choose the design formula in DESeq2?

The design formula should include all sources of variation that are relevant to the biological question and that could affect gene expression. At a minimum, it should include the primary variable of interest (e.g., condition, treatment, genotype). It should also include any blocking factors that account for known sources of technical variation, such as batch, sequencing lane, or sex. The design formula should not include variables that are not relevant or that have too few samples per level. For example, if there are only two batches and one batch is confounded with the condition, the batch cannot be included because it is perfectly confounded with the condition. The design formula should be specified before running the DESeq function, and it can be tested with different formulations to assess the sensitivity of the results. A common approach is to start with a simple design (~ condition) and then add blocking factors to see if they improve the analysis.

Key Takeaways

DESeq2 models RNA-seq count data with a negative binomial distribution, accounting for overdispersion that the Poisson distribution cannot capture.
The median-of-ratios normalization estimates size factors that correct for library size and composition bias, enabling accurate comparisons across samples.
Empirical Bayes shrinkage stabilizes dispersion and fold change estimates, particularly for low-count genes, reducing false positives.
The Wald test is used for pairwise comparisons, while the likelihood ratio test is used for multi-level factors and nested model comparisons.
The Benjamini-Hochberg procedure controls the false discovery rate and should always be used for significance filtering.
The design formula must include all relevant sources of variation, including blocking factors like batch, to avoid confounding and spurious results.
Biological interpretation requires considering both statistical significance (adjusted p-value) and effect size (fold change), as small but significant changes may not be biologically meaningful.