AboutBlogMediaTags

Differential Expression Analysis: Guide with R & Python Code

Runcell Team,

What Is Differential Expression Analysis?

Differential expression (DE) analysis identifies which genes change their expression between conditions (e.g., treated vs. control). This guide walks you through the complete workflow for bulk RNA-seq and single-cell data, with ready-to-run code in R (DESeq2, edgeR, limma-voom) and Python (Scanpy).

The devil, as usual, is in the details — experimental design, normalization, dispersion modeling, multiple testing, and interpretation. This post gives you a clear, practical path through bulk RNA‑seq (and what changes for single‑cell), with concise definitions and copy-paste examples.

Updated 2026-10-19: section 6 now includes an executed check, on the public PBMC 3k dataset, of what Scanpy’s logfoldchanges column measures and how far it sits from a ratio of means. The R code on this page was not re-executed in this update.

TL;DR

  1. Plan the design (replicates, covariates, contrasts).
  2. Get counts (gene × sample matrix) and metadata (sample table).
  3. Filter low counts; normalize library size; model counts (NB).
  4. Test contrasts, control FDR, shrink log2FCs.
  5. QC and visualize (PCA, MA, volcano).
  6. Interpret with gene set enrichment.
  7. For single‑cell, prefer pseudobulk across biological replicates, or use cell‑level tests with care.

Key terms you’ll actually use

0) Experimental design first (don’t skip!)

1) From FASTQ to counts (what you need to know)

Everything on this page starts from a count matrix, which means it starts from libraries that were pooled at a target molarity, not a target mass. That conversion is where a surprising amount of uneven sequencing depth comes from — the ng/µL → nM step depends on a molecular-weight convention that different calculators do not agree on, to the tune of about 7%. If your library sizes came back lopsided, it is worth reading DNA molecular weight in Python before you blame the sequencer, and using the DNA concentration calculator  to redo the pooling arithmetic.

2) Bulk RNA‑seq with DESeq2 (R)

Below is a lean, reliable template using the Bioconductor airway dataset (glucocorticoid treatment). It shows filtering, modeling with a batch covariate, contrasts, LFC shrinkage, and common plots.

Install once (uncomment if needed):

# install.packages("BiocManager")
# BiocManager::install(c("DESeq2", "airway", "apeglm"))
# install.packages("ggplot2")
library(DESeq2)
library(airway)
library(ggplot2)
 
# 1) Load example counts + metadata
data(airway)
se <- airway
# Factors: cell line (~batch-like) and dex (treated vs untrt)
se$dex  <- relevel(se$dex, ref = "untrt")
 
# 2) Build DESeqDataSet, filter very low counts
dds <- DESeqDataSet(se, design = ~ cell + dex)
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep,]
 
# 3) Fit model (NB with empirical dispersion)
dds <- DESeq(dds)
 
# 4) Extract results for the contrast of interest
res <- results(dds, contrast = c("dex", "trt", "untrt"))  # treated vs control
res <- lfcShrink(dds, coef = "dex_trt_vs_untrt", type = "apeglm")  # effect-size shrinkage
 
# 5) QC & visualization
vsd <- vst(dds, blind = FALSE)
 
# PCA colored by condition
plotPCA(vsd, intgroup = "dex")
 
# MA plot of shrunk effects
plotMA(res, ylim = c(-4,4))
 
# Volcano with ggplot2
res_df <- as.data.frame(res)
res_df$neglog10padj <- -log10(res_df$padj)
res_df$hit <- with(res_df, !is.na(padj) & padj < 0.05 & abs(log2FoldChange) > 1)
 
ggplot(res_df, aes(x = log2FoldChange, y = neglog10padj)) +
  geom_point(alpha = 0.4) +
  geom_hline(yintercept = -log10(0.05), linetype = 2) +
  geom_vline(xintercept = c(-1, 1), linetype = 2) +
  labs(x = "log2 fold change (shrunk)", y = "-log10(FDR)", title = "Volcano plot")
 
# 6) Tidy significant table
sig <- subset(res_df, padj < 0.05)
sig <- sig[order(sig$padj), c("log2FoldChange","lfcSE","stat","pvalue","padj")]
head(sig, 10)

Notes on this workflow

3) Alternative bulk frameworks

If you would rather stay in Python than switch to R for the DESeq2 step, the whole model above has a native Python implementation — see the step-by-step PyDESeq2 tutorial, which covers count matrices, dispersion fitting, LFC shrinkage and plotting without an R install.

edgeR (R, quasi‑likelihood GLM; great for multifactor designs)

# install.packages("BiocManager"); BiocManager::install("edgeR")
library(edgeR)
 
counts_mat <- counts(dds)         # reuse counts from above for illustration
group <- colData(dds)$dex
batch <- colData(dds)$cell
 
y <- DGEList(counts = counts_mat)
keep <- filterByExpr(y, group = group)  # adaptive low-count filter
y <- y[keep,, keep.lib.sizes=FALSE]
y <- calcNormFactors(y, method = "TMM")
 
design <- model.matrix(~ batch + group)
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit, coef = "grouptrt")   # treated vs control
 
topTags(qlf, n = 10)

limma‑voom (R, fast and flexible; works exceptionally well with many samples)

# BiocManager::install(c("limma"))
library(limma)
 
v <- voom(counts_mat, design = design)  # uses the same 'design' from edgeR block
fit <- lmFit(v, design)
fit <- eBayes(fit)
topTable(fit, coef = "grouptrt", number = 10)

4) Going beyond a simple two‑group comparison

5) Common QC you should always look at

# Sample distance heatmap (VST) example
library(pheatmap)
d <- dist(t(assay(vsd)))
pheatmap(as.matrix(d), clustering_distance_rows = d, clustering_distance_cols = d)

QC plots are the part collaborators actually ask to see, and they travel badly as screenshots. If you are working in a notebook, the ipynb → HTML converter  turns the whole QC section into a single self-contained file with the figures inline — no kernel needed at the other end. For a fixed lab-record copy, ipynb → PDF  does the same without a local LaTeX install.

6) Single‑cell RNA‑seq: do it right

Two solid strategies

  1. Pseudobulk (recommended when you have replicates) Aggregate counts per sample × cell type (or cluster), then run bulk DE (DESeq2/edgeR/limma).

    • Captures between‑sample variability → valid p‑values.
    • Easy to add covariates (donor, batch).
    • Use the same steps as bulk, just do them per cell type.

    Sketch in R (using a SingleCellExperiment sce):

    # BiocManager::install(c("SingleCellExperiment","scuttle"))
library(SingleCellExperiment); library(scuttle)
# ids says how to group cells; here by sample and cell_type
pb <- aggregateAcrossCells(sce, ids = colData(sce)[, c("sample", "cell_type")])
# pb assays now contain summed counts; build DE model per cell_type or with interaction

  2. Cell‑level tests (when pseudobulk is impossible) Use methods like Wilcoxon rank‑sum (Scanpy), or model‑based (MAST). Control for library size and latent factors.

    Scanpy example (Python, cluster markers):

    import scanpy as sc
adata = sc.read_h5ad("my_data.h5ad")         # adata.X should be counts or use adata.layers['counts']
sc.pp.normalize_total(adata, target_sum=1e4) # for ranking/visualization
sc.pp.log1p(adata)
sc.tl.leiden(adata, resolution=0.5)
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")
sc.pl.rank_genes_groups(adata, n_genes=20)

    Condition DE within a cluster:

    cluster = "2"
ad = adata[adata.obs["leiden"] == cluster].copy()
sc.tl.rank_genes_groups(ad, groupby="condition", groups=["treated"], reference="control", method="wilcoxon")
sc.get.rank_genes_groups_df(ad, group="treated").head()

What logfoldchanges actually is in Scanpy (checked on a real dataset)

Added 2026-10-19. Everything in this subsection was executed on 10x Genomics’ public PBMC 3k dataset (sc.datasets.pbmc3k(), 2,700 cells) with scanpy 1.12.4, anndata 0.13.4, numpy 2.5.4 and pandas 3.0.6; the output blocks are pasted from that run, not typed.

The snippet above skips the steps between loading and clustering, so the run that produced the output below had to add them. This is the pipeline that was actually executed, with the one line that keeps a copy of the normalized-but-not-logged values for the comparison further down:

adata = sc.datasets.pbmc3k()                       # raw 10x counts, 2700 cells x 32738 genes
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
adata.layers["norm"] = adata.X.copy()              # normalized, not logged — used below
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
sc.pp.pca(adata, n_comps=30, mask_var="highly_variable")
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30, random_state=0)
sc.tl.leiden(adata, resolution=0.5, random_state=0, flavor="igraph", n_iterations=2)
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")

(flavor="igraph" because the default leidenalg flavor is not installed in that environment; it gave 8 clusters.) Asked for cluster 0, sc.get.rank_genes_groups_df returns this:

names    scores  logfoldchanges         pvals     pvals_adj
 CCL5 28.466553        5.961505 3.040742e-178 4.170073e-174
 NKG7 23.936359        4.918181 1.281658e-126 8.788332e-123
  B2M 19.863453        0.733833  8.431016e-88  3.854099e-84
 IL32 19.377666        2.669144  1.191286e-83  4.084324e-80
 GZMA 19.228132        4.025735  2.152405e-82  5.903617e-79

If you come from DESeq2, you will read logfoldchanges as log2(mean expression in the cluster / mean expression outside). It is not that number. Three candidate definitions were recomputed for those five genes from the same normalized matrix and compared with what Scanpy reports:

gene  scanpy_lfc  A_log2_ratio_of_means  B_scanpy_formula  C_diff_of_log_means_log2  mean_norm_in  mean_norm_out
CCL5       5.962                  3.689             5.962                     4.639     51.397999       3.984000
NKG7       4.918                  2.269             4.918                     3.746     42.792000       8.878000
 B2M       0.734                  0.649             0.734                     0.730    282.230011     180.044998
IL32       2.669                  1.489             2.669                     2.266     21.267000       7.578000
GZMA       4.026                  2.466             4.026                     2.265     10.218000       1.850000

max |scanpy_lfc - B| = 4.291534423828125e-06
max |scanpy_lfc - A| = 2.64912486076355
max |scanpy_lfc - C| = 1.7609531192757912

Column B reproduces Scanpy’s value to within 4.3 × 10⁻⁶ on all five genes. It is log2((expm1(mean of log1p values in group) + 1e-9) / (expm1(mean of log1p values outside) + 1e-9)): the mean is taken on the log scale first, then un-logged. Because expm1(mean(log)) is a geometric-style mean, it is pulled down by the zeros that dominate single-cell data, and the two groups rarely have the same fraction of zeros, so the ratio drifts away from a ratio of arithmetic means. For CCL5 that drift is 2.3 log2 units: Scanpy reports a 62-fold change where the ratio of mean normalized counts is 13-fold.

How often does it matter? Over the 5,690 genes with mean normalized count above 0.1 on both sides:

  B - A (Scanpy minus ratio-of-means), log2 units: median -0.077, 5th pct -0.388, 95th pct +0.238
  |B - A| > 1 log2 unit (i.e. > 2x) in 88 genes; sign differs in 642 genes
  sign differs AND ratio-of-means |LFC| > 0.5 (the two summaries disagree in direction by more than 1.4x): 43 genes
  five largest disagreements:
  gene  scanpy_lfc  ratio_of_means_lfc  mean_norm_in  mean_norm_out
NFE2L2      -0.108               2.808      2.195000          0.313
  NKG7       4.918               2.269     42.792000          8.878
  AIM2      -0.858               1.758      0.373000          0.110
   PF4      -0.706              -3.091      0.119000          1.017
  CCL5       5.962               3.689     51.397999          3.984

Most genes agree within ±0.4 log2 units, which is why this rarely gets noticed. The tail is the problem: NFE2L2 has a 7× higher mean normalized count in the cluster and Scanpy reports it as slightly down. The per-cell values show why:

NFE2L2: per-cell normalized counts, inside cluster 0 vs outside
  in  cells=  370  nonzero=  19 (5.1%)  mean=2.195  median=0.000  max=713.25  top5=[713.25, 9.890000343322754, 9.0, 8.640000343322754, 6.869999885559082]  top5 share of sum=92.1%  mean(log1p)=0.1071
  out cells= 2330  nonzero= 166 (7.1%)  mean=0.313  median=0.000  max=16.89  top5=[16.889999389648438, 16.420000076293945, 14.619999885559082, 13.329999923706055, 12.479999542236328]  top5 share of sum=10.1%  mean(log1p)=0.1150

One cell in the cluster carries a normalized count of 713 — about 88% of the cluster’s total of roughly 812 (370 cells × mean 2.195); the top five cells together are 92.1% — so the arithmetic mean is 7× higher inside, while the gene is detected in a smaller fraction of cells inside (5.1% vs 7.1%) and the mean of the log values is slightly lower inside. The arithmetic mean sees the outlier; the mean of logs sees the detection rate. For 43 genes the two summaries disagree in sign while the ratio-of-means version is larger than 0.5 log2 units. These are differences between two descriptive summaries of the same matrix; nothing here tests which genes are significantly different, and none of it involves replicates, covariates or shrinkage the way the DESeq2 example above does.

What to do with that:

names    scores  logfoldchanges         pvals     pvals_adj  pct_nz_group  pct_nz_reference
 CCL5 28.466553        5.961505 3.040742e-178 4.170073e-174      0.986486          0.213305
 NKG7 23.936359        4.918181 1.281658e-126 8.788332e-123      0.913514          0.205150
  B2M 19.863453        0.733833  8.431016e-88  3.854099e-84      1.000000          0.999571
 IL32 19.377666        2.669144  1.191286e-83  4.084324e-80      0.948649          0.491416
 GZMA 19.228132        4.025735  2.152405e-82  5.903617e-79      0.732432          0.115880

The formula is in Scanpy’s source (scanpy/tools/_rank_genes_groups.py, the expm1 branch) and the check above is the only evidence this page offers for it; the run is deterministic (random_state=0), so the numbers should reproduce exactly with the same package versions.

Tip: If you have multiple donors per condition, pseudobulk generally yields the most reliable inferences for scRNA‑seq DE because it respects biological replication.

7) Multiple testing, thresholds, and effect sizes

8) Gene set enrichment (quick start)

Once you have a ranked list, run GSEA/ORA against pathways (GO/KEGG/MSigDB). In R, fgsea works well:

# BiocManager::install("fgsea")
library(fgsea)
# ranks: named vector of gene scores (e.g., signed -log10 p; or log2FC)
ranks <- res_df$stat; names(ranks) <- rownames(res_df)
# pathways <- your list of gene sets as a named list
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize=500)
fgseaRes[order(fgseaRes$pval),][1:10, c("pathway","NES","padj")]

9) Frequent pitfalls & how to dodge them

10) Handy snippets (copy/paste)

DESeq2 custom contrast (multi‑group):

# Suppose condition has levels: ctrl, drugA, drugB
res_AB <- results(dds, contrast = c("condition", "drugA", "drugB"))

DESeq2 interaction test (does treatment effect differ by genotype?):

# design = ~ batch + genotype * treatment
# Test the interaction term directly (naming may vary by factor levels)
results(dds, name = "genotypeKO.treatmenttreated")

tximport for Salmon/Kallisto → gene‑level counts:

# BiocManager::install(c("tximport","readr"))
library(tximport); library(readr)
samples <- read_csv("samples.csv") # columns: sample, condition, path_to_quant
files <- setNames(samples$path_to_quant, samples$sample)
tx2gene <- read_csv("tx2gene.csv") # columns: transcript, gene
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ batch + condition)
dds <- DESeq(dds)

edgeR user‑defined contrasts (limma-style):

library(limma); library(edgeR)
design <- model.matrix(~ 0 + condition)         # no intercept
colnames(design) <- levels(condition)
y <- estimateDisp(DGEList(counts), design)
fit <- glmQLFit(y, design)
contr <- makeContrasts(DrugA_vs_Ctrl = DrugA - Ctrl, levels = design)
qlf <- glmQLFTest(fit, contrast = contr[,"DrugA_vs_Ctrl"])
topTags(qlf)

11) A short checklist you can tape next to your monitor

Closing thoughts

Differential expression analysis is conceptually simple but statistically rich. If you (1) plan your design, (2) model with the correct likelihood (NB for counts), (3) control FDR, and (4) keep an eye on effect sizes and QC plots, you will produce results that are both statistically defensible and biologically interpretable. The code above is a practical starting point you can adapt to bulk and single‑cell projects alike.

The design formula is the part worth spending your time on. Everything downstream — shrinkage, FDR, enrichment — is comparatively mechanical, and a formula that omits a real covariate cannot be rescued by any of it.

© Runcell.RSS