Differential Expression with DESeq2 on Your Own RNA-Seq Data

By the end of this guide you will have a table of genes ranked by evidence for differential expression between two sets of your RNA-seq samples. The table includes shrunken log2 fold changes and adjusted p-values. You will also have a set of diagnostic plots that tell you whether the model fit the data. To follow along you need a few things. Bring bulk RNA-seq FASTQ files or a count matrix someone else produced. You also need a sample table recording what distinguishes the samples, R 4.4 or newer with Bioconductor 3.20 or newer, and at least two samples per group.

The requirement for replicates is not arbitrary. DESeq2 estimates the variance of each gene from replicates, so with one sample per group there is no variance to estimate, and the package will refuse to give you dispersion estimates worth reading. For a single-person longitudinal profile, the replicates are usually timepoints: several draws under one condition compared against several draws under another. That is a legitimate design, but the inference it supports concerns within-person variation over time rather than a population, and you should say so when you write down what you found.

One caveat before any code. Gene expression differences are measurements rather than findings about your health. Transcript abundance in whole blood is dominated by cell composition, and RNA changes frequently fail to track protein changes in the same tissue ¹. If something in your results bears on a medical question, take it to a clinician rather than acting on it.

1. Produce counts you trust

Everything downstream depends on getting the input matrix right, so it is worth being precise about what DESeq2 expects. The package consumes a matrix of counts in which rows are genes, columns are samples, and each entry is the number of reads (or estimated reads) assigned to that gene in that sample. The entries must be counts rather than TPM, FPKM, or anything else already normalized for library size. The model is a negative binomial on integer counts and normalization happens inside the package.

For quantification we prefer Salmon followed by tximport. Salmon’s selective alignment handles multi-mapping across transcript isoforms properly, and it runs in a few minutes per sample on a laptop-class machine.

# index once, from a transcriptome FASTA plus decoys (genome) built with generateDecoyTranscriptome
salmon index -t gentrome.fa.gz -d decoys.txt -p 16 -i salmon_idx_gencode47 -k 31

# per sample
salmon quant -i salmon_idx_gencode47 -l A \
  -1 S01_R1.fastq.gz -2 S01_R2.fastq.gz \
  -p 16 --gcBias --seqBias --numBootstraps 30 \
  -o quant/S01

A few of those flags deserve comment. -l A asks Salmon to infer the library type; check the inferred value in quant/S01/lib_format_counts.json and confirm it is the same for every sample, because a mixed set points to a library prep or file pairing error. --gcBias and --seqBias cost little and correct fragment-level biases that would otherwise appear as apparent expression differences between batches.

With quantification done per transcript, the next step is to collapse transcripts to genes in R:

library(tximport); library(readr)
tx2gene <- read_tsv("tx2gene_gencode47.tsv", col_names = c("tx", "gene"))
files <- file.path("quant", samples$id, "quant.sf"); names(files) <- samples$id
txi <- tximport(files, type = "salmon", tx2gene = tx2gene, ignoreTxVersion = FALSE)

tximport returns counts along with an average transcript length per gene per sample, and DESeqDataSetFromTximport uses those lengths as a per-gene offset. That correction matters when isoform usage shifts between your groups, which happens often enough to be worth getting right. The alternative pipeline is STAR alignment followed by featureCounts -p --countReadPairs -s 2 -t exon -g gene_id. It is perfectly fine and gives you a BAM file you can inspect, at the cost of roughly an order of magnitude more CPU time and disk.

2. Build the DESeqDataSet

The design formula tells DESeq2 which variables explain the counts, and mistakes here are expensive because they are silent. Write the sample table first and read it carefully before you fit anything.

library(DESeq2)
samples <- data.frame(
  id        = c("S01","S02","S03","S04","S05","S06"),
  condition = factor(c("base","base","base","interv","interv","interv"),
                     levels = c("base","interv")),
  batch     = factor(c("A","A","B","A","B","B")),
  row.names = c("S01","S02","S03","S04","S05","S06")
)
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ batch + condition)

Notice that the factor levels are set explicitly. DESeq2 orders levels alphabetically by default, which silently makes “base” the numerator or the denominator depending on your naming, and a sign error in a fold change is the kind of thing nobody catches until much later. By convention the variable of interest goes last in the formula, because results() with no arguments reports the last coefficient.

One structural point about replicates. If you sequenced the same library on two lanes, sum those columns before modeling with collapseReplicates(dds, dds$sample_id, dds$run_id). Technical replicates are not additional degrees of freedom, and treating them as biological replicates inflates significance.

3. Filter low-count genes

Genes with almost no reads carry no information about differential expression, and every one you keep costs you power once p-values are adjusted for multiple testing. The minimal filter from the DESeq2 vignette is a single line:

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

We use something stricter, for a concrete reason. A gene with 40 reads in one sample and zero in eleven others passes rowSums >= 10 and will frequently come back with a small p-value driven entirely by that single sample. Regulatory-grade frameworks for transcriptomics require a gene to exceed a counts-per-million threshold in at least about three quarters of the samples of at least one experimental group before it enters the test ². That group-wise prevalence rule is what we implement:

cpm  <- edgeR::cpm(counts(dds))
grp  <- dds$condition
keep <- sapply(levels(grp), function(g) rowMeans(cpm[, grp == g, drop = FALSE] > 1) >= 0.75)
dds  <- dds[rowSums(keep) >= 1, ]
nrow(dds)   # expect roughly 12,000-16,000 for a human tissue or blood library

Whatever filter you choose, it must depend only on the counts and never on the outcome you are testing, or the p-values are no longer valid. Prevalence-and-abundance filters are independent of the group labels in the relevant sense, so they are safe to use here.

4. Understand what normalization DESeq2 performs

Libraries differ in depth and composition, and DESeq2 corrects for this in a specific way that is useful to understand before you start interpreting size factors. The package does not divide by total reads. Instead it computes a size factor per sample by the median-of-ratios method: for each gene it takes the geometric mean of that gene’s counts across all samples, divides each sample’s count by that reference to form a ratio, and takes the median ratio over genes as that sample’s size factor. Using the median is what makes the estimate resistant to a handful of very highly expressed genes, which total-count scaling is not.

dds <- estimateSizeFactors(dds)
sizeFactors(dds)

Expect values within roughly 0.5 to 2. A size factor of 4 or 0.25 tells you that one library is far deeper or far shallower than the rest, and it is worth checking whether the cause is depth or composition. This is where whole-blood RNA-seq bites: hemoglobin transcripts can consume a large fraction of reads, and if that fraction varies between draws, the median-of-ratios estimate absorbs part of the difference while the effective depth for the remaining genes still varies substantially. Globin depletion at library prep, or removing HB* genes before computing size factors, both help.

One practical exception: if some genes are zero in every sample after filtering, use estimateSizeFactors(dds, type = "poscounts"), which computes the geometric mean over nonzero entries only.

5. Estimate dispersion and fit the model

With normalization settled, the model itself can be fit, and the key quantity is dispersion. Read counts for a gene across replicates are overdispersed relative to Poisson, because biological variability adds a term proportional to the square of the mean. DESeq2 therefore fits a negative binomial generalized linear model with a log link, in which each gene carries its own dispersion parameter.

With three or four replicates, a per-gene dispersion estimate on its own is terrible. DESeq2 handles this by fitting a smooth trend of dispersion against mean expression across all genes and shrinking each gene’s estimate toward that trend, with the amount of shrinkage set by how far the gene sits from the curve and how many residual degrees of freedom the design leaves. Genes whose dispersion lies far above the trend are not shrunk, which keeps the test conservative for genuinely noisy genes.

dds <- DESeq(dds)          # size factors, dispersions, fit, and Wald test
plotDispEsts(dds)

DESeq() is a wrapper around estimateSizeFactors, estimateDispersions, and nbinomWaldTest. Look at the dispersion plot before you look at any individual gene. A healthy plot shows a cloud of black per-gene estimates, a red trend line falling from roughly 0.1-1 at low counts toward 0.01-0.05 at high counts, and blue final estimates pulled toward that line with a scatter of un-shrunk outliers above it. A trend that is flat and high, or estimates that do not follow the curve at all, usually means the design formula is missing a variable that explains a lot of variance.

6. Choose the Wald test or the likelihood ratio test

DESeq2 offers two hypothesis tests, and which one you want depends on the question. The Wald statistic is the point of the default test: for each gene DESeq2 divides the fitted log2 fold change by its standard error and compares the ratio to a standard normal distribution, giving a p-value for the null hypothesis that the coefficient is zero. It tests one coefficient, or one contrast of coefficients, at a time. It returns an effect size with an interval, which is what you want for a two-group comparison.

resultsNames(dds)
res <- results(dds, contrast = c("condition", "interv", "base"), alpha = 0.05)
summary(res)

Set alpha to whatever false discovery rate you intend to report. The argument is not cosmetic. results() performs independent filtering, choosing a mean-normalized-count threshold that maximizes the number of rejections at that alpha, and genes below the threshold get padj = NA. If you pass the default 0.1 and then apply a 0.05 cutoff yourself, the filtering was tuned for the wrong target.

Use the likelihood ratio test instead when you want to ask whether a variable with more than two levels matters at all, for example whether expression changes across five timepoints in any pattern:

dds_t <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ batch + timepoint)
dds_t <- DESeq(dds_t, test = "LRT", reduced = ~ batch)
res_t <- results(dds_t, alpha = 0.05)

The LRT compares deviance between the full and reduced models against a chi-square distribution with degrees of freedom equal to the difference in parameter count. It yields one p-value per gene for the whole variable and no single interpretable fold change, so pair it with a per-contrast Wald test or with clustering of the significant genes.

7. Shrink the fold changes

Raw maximum likelihood log2 fold changes are unreliable for low-count genes, and the size of the problem is easy to see. A gene with 12 reads in one group and 2 in the other gets an LFC near 2.6 with an enormous standard error. Shrinkage replaces these with estimates pulled toward zero in proportion to their noise, and we would always use it when ranking or plotting.

resLFC <- lfcShrink(dds, coef = "condition_interv_vs_base", type = "apeglm")
plotMA(resLFC, ylim = c(-4, 4))

apeglm is the default recommendation and works well at small sample sizes. Keep in mind that shrinkage changes the effect sizes but leaves the Wald p-values untouched, so ranking by padj and ranking by shrunken log2FoldChange give different orders. Our practice is to rank by padj, filter by abs(log2FoldChange) computed from the shrunken values, and never filter on the unshrunken ones. If you need a formal test against a nonzero threshold, use results(dds, lfcThreshold = 1, altHypothesis = "greaterAbs"), which tests the composite null properly rather than relying on post-hoc filtering.

8. Run the diagnostics that catch real errors

Three plots find most of the problems that occur in practice, and they are quick to produce. The first is a PCA on variance-stabilized data, which shows you what dominates the variation between samples:

vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = c("condition", "batch"))

If PC1 separates by batch rather than condition, your design needs the batch term and your power is lower than you hoped. The second diagnostic is the p-value histogram from results(..., independentFiltering = FALSE). A healthy result is roughly uniform between 0.2 and 1 with a spike near zero; a histogram rising toward 1, or humped in the middle, indicates a misspecified model or correlated samples. The third is the dispersion plot from step 5.

It also pays to check mcols(res)$maxCooks behavior. DESeq2 flags genes whose fit is dominated by a single sample using Cook’s distance and sets their p-value to NA when there are at least three replicates per group. Genes with padj = NA therefore fall into three buckets. They may have all counts zero, be outlier-driven, or be filtered by independent filtering. Running summary(res) tells you how many are in each.

9. Handle batch and other covariates

Batch effects are the most common nuisance variable in this kind of data, and the design formula is the right place to deal with them. Adding + batch lets the GLM absorb an additive batch offset on the log scale, which is the right move whenever batch is not confounded with condition. That approach fails when the batch effect differs between conditions, which happens in multi-factor designs, and a method has been developed specifically to correct multi-factorial RNA-seq designs where batch interacts with the factors of interest rather than shifting everything uniformly ³. If your sample table has a batch column that is perfectly correlated with condition, no software can separate the two, and the only defensible analysis is to report that the comparison is confounded.

One thing to avoid: never run a batch-correction tool on the counts and then feed the corrected values to DESeq2. The model needs raw counts with the batch structure declared in the formula. Corrected matrices belong in visualizations, produced with limma::removeBatchEffect on the VST output.

10. If you would rather work in Python

The same model is available outside R. PyDESeq2 reimplements the DESeq2 method for bulk RNA-seq and reproduces the R package’s results closely, while offering an anndata-based interface that fits into a scanpy or scikit-learn workflow ⁴.

from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

dds = DeseqDataSet(counts=counts_df, metadata=meta_df,
                   design="~batch + condition", refit_cooks=True, n_cpus=8)
dds.deseq2()
st = DeseqStats(dds, contrast=["condition", "interv", "base"], alpha=0.05, n_cpus=8)
st.summary()
st.lfc_shrink(coeff="condition[T.interv]")

Two differences are worth knowing. Counts go in samples-by-genes orientation, the transpose of the R convention. And the R package remains ahead on peripheral features, including DESeqDataSetFromTximport with length offsets and some of the shrinkage estimators, so if your quantification came from Salmon we still run the main analysis in R and move the results to Python afterward.

A related question is DESeq2 versus edgeR. Both fit negative binomial GLMs and both are defensible choices. A systematic comparison across studies found that the two tools differ in sensitivity and cross-study consistency in ways that depend on the dataset, so the choice is not neutral and agreement between them is a useful signal ⁵. We run DESeq2 as primary because of independent filtering and apeglm shrinkage, then rerun with edgeR::glmQLFTest on the same filtered matrix and treat the intersection as the confident set. If you want a graphical route for exploration before committing to a pipeline, integrated interfaces built on these same engines exist and will get you a first pass without writing R ⁶.

Common problems

Most failed analyses trace back to a short list of causes, so it helps to know what they look like.

The most frequent error is passing normalized values into DESeqDataSetFromMatrix. If your matrix carries decimals from TPM, or has already been CPM-scaled, the negative binomial model is wrong and the p-values are meaningless. Estimated counts from Salmon or kallisto do contain decimals and are fine, because tximport and DESeq2 handle them.

The second is a design formula that omits a variable driving most of the variance. In longitudinal single-subject data, the usual culprits are draw time of day, cell composition, and library prep batch. Whole-blood signal in particular is largely cell-composition signal, and multi-omic studies of blood routinely find that the differentially expressed genes track immune cell programs rather than anything tissue-specific ⁷. Estimating cell fractions from the expression matrix and including them as covariates changes the result list substantially.

The third is too few replicates. With two samples per group, dispersion shrinkage is doing nearly all the work and the result list is dominated by whichever genes happen to agree across the two pairs. We would not report an FDR from n=2 without an independent confirmation, and three or more per group is the practical floor.

The fourth is treating padj = NA as “not significant” when it means “not tested.” Count how many genes were dropped and why before interpreting the fraction of the transcriptome that appears to have changed.

The fifth is over-reading the gene list. A significant transcript change is a change in RNA abundance in the sampled tissue at the sampled time. Paired transcriptome and proteome measurements on the same cells frequently disagree on direction and magnitude for a large share of genes ¹, and integrating across layers is what turns an expression list into a claim about biology ⁸. Where a result touches a clinical question, that interpretation belongs with a physician who can see the rest of your record.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Letícia B. C. Penariol, Carolina H. Thomé, Patrícia A. Tozetti, et al. What Do the Transcriptome and Proteome of Menstrual Blood-Derived Mesenchymal Stem Cells Tell Us about Endometriosis?. International Journal of Molecular Sciences, 2022. https://doi.org/10.3390/ijms231911515 ↩ ↩²

2. M.C. Verheijen, T.W. Gant, W. Tong, et al. R-ODAF: Omics data analysis framework for regulatory application. Toxicology Letters, 2021. https://doi.org/10.1016/s0378-4274(21)00539-7 ↩

3. Julien Roy, Adrian S. Monthony, Davoud Torkamaneh. DESeq2-MultiBatch: batch correction for multi-factorial RNA-seq experiments. Genome, 2026. https://doi.org/10.1139/gen-2025-0049 ↩

4. Boris Muzellec, Maria Teleńczuk, Vincent Cabeli, et al. PyDESeq2: a python package for bulk RNA-seq differential expression analysis. Bioinformatics, 2023. https://doi.org/10.1093/bioinformatics/btad547 ↩

5. Mostafa Rezapour. Tool choice matters: Evaluating edgeR vs. DESeq2 for sensitivity, robustness, and cross-study performance. PLOS One, 2026. https://doi.org/10.1371/journal.pone.0353788 ↩

6. Brandon Monier, Adam McDermaid, Cankun Wang, et al. IRIS-EDA: An integrated RNA-Seq interpretation system for gene expression data analysis. PLOS Computational Biology, 2019. https://doi.org/10.1371/journal.pcbi.1006792 ↩

7. Ze Wang, Yi Huang, Ziyu Guo, et al. Interferon-Linked Lipid and Bile Acid Imbalance Uncovered in Ankylosing Spondylitis in a Sibling-Controlled Multi-Omics Study. International Journal of Molecular Sciences, 2025. https://doi.org/10.3390/ijms26167919 ↩

8. Sarah Kate Nyquist, Laasya Devi Annepureddy, Kristija Sejane, et al. Integrated ‘omics analysis reveals human milk oligosaccharide biosynthesis programs in human lactocytes. iScience, 2025. https://doi.org/10.1016/j.isci.2025.113269 ↩