Pseudobulk differential expression with DESeq2
pseudobulk.RmdWhy pseudobulk?
Single cells from the same donor are not statistically independent: they share genetic background, batch effects, and ambient RNA. Treating each cell as an independent observation in a regression inflates the effective sample size and distorts p-values.
A common, principled remedy is pseudobulk: sum each
donor’s cells within a group (e.g. cluster) to produce a per-donor count
vector, then run a count-based regression like
DESeq2::DESeq over donors. The unit of replication becomes
the donor, not the cell.
presto provides:
-
collapse_counts()— fast collapse of a cell-by-gene count matrix to pseudobulks defined by one or more metadata columns. -
pseudobulk_deseq2()— run DESeq2 across the resulting pseudobulks in a one-vs-all, pairwise, or within-group design. -
top_markers_dds()— extract top features per group from the result.
Demo data
We use the colon-tissue CD8 T-cell dataset from Thomas et al. Nat. Med. 2024 (GSE206299). 27 donors, 9 cell clusters, 25,341 cells.
library(presto)
library(Matrix)
library(dplyr)
d <- load_ircolitis_cd8(verbose = FALSE)
dim(d$counts)
#> [1] 28165 25341
table(d$obs$donor, d$obs$cluster)[1:5, ]
#>
#> 1 2 3 4 5 6 7 8 9
#> MC_1 323 110 46 100 80 51 61 7 5
#> MC_2 580 229 113 154 179 115 78 27 5
#> MC_9 698 449 143 174 230 119 91 32 10
#> SIC_100 428 373 119 152 150 164 59 25 10
#> SIC_109 172 85 33 68 45 21 170 6 4The crosstab confirms most donors contribute cells to most clusters — good news for pseudobulk, since clusters with few represented donors will have low statistical power.
Step 1: collapse cells to pseudobulks
collapse_counts() sums the columns of a count matrix
according to one or more metadata columns. Here we collapse by
(donor, cluster), drop any (donor, cluster) pair with fewer
than 10 cells, and keep the cell count N so we can inspect
what came through.
md <- data.frame(
donor = as.character(d$obs$donor),
cluster = as.character(d$obs$cluster),
case = as.character(d$obs$case),
stringsAsFactors = FALSE
)
cc <- collapse_counts(
d$counts, md, c("donor", "cluster", "case"),
min_cells_per_group = 10,
keep_n = TRUE
)
dim(cc$counts_mat)
#> [1] 28165 204
head(cc$meta_data)
#> donor cluster case N
#> 1 MC_1 5 Control 80
#> 2 MC_1 3 Control 46
#> 3 MC_1 4 Control 100
#> 4 MC_1 1 Control 323
#> 5 MC_1 6 Control 51
#> 6 MC_1 2 Control 110We collapse by (donor, cluster, case) instead of just
(donor, cluster) because case is constant
within donor — including it now makes it available for the
within-cluster contrast in Step 4.
cc$counts_mat is now a 28,165-gene × 204-pseudobulk
matrix. Each column is one (donor, cluster) sample;
cc$meta_data is the matching data.frame of
pseudobulk identities.
table(cc$meta_data$cluster)
#>
#> 1 2 3 4 5 6 7 8 9
#> 27 27 26 27 25 23 26 18 5Step 2: one-vs-all DESeq2
pseudobulk_deseq2() with
mode = "one_vs_all" tests each cluster against the union of
the others. With collapse_background = FALSE we keep each
background pseudobulk as its own observation; setting it to
TRUE collapses the background pseudobulks per donor first,
which is preferable when donors are unbalanced across clusters.
The first variable in the formula is the contrast variable.
meta_data should contain only the pseudobulk identifying
variables used in the formula.
md_dge <- cc$meta_data[, "cluster", drop = FALSE]
res <- pseudobulk_deseq2(
~cluster, md_dge, cc$counts_mat,
verbose = FALSE,
min_counts_per_sample = 10,
present_in_min_samples = 5,
collapse_background = FALSE,
mode = "one_vs_all"
)
head(res)
#> group feature baseMean log2FoldChange lfcSE stat pvalue padj
#> 1 1 ENSG00000128951|DUT 42.72 1.35 0.108 12.5 6.67e-36 6.69e-32
#> 2 1 ENSG00000105486|LIG1 4.18 1.69 0.138 12.2 2.49e-34 1.25e-30
#> 3 1 ENSG00000134291|TMEM106C 13.14 1.67 0.138 12.1 7.54e-34 2.52e-30
#> 4 1 ENSG00000188486|H2AFX 10.49 1.86 0.155 12.0 2.92e-33 7.33e-30
#> 5 1 ENSG00000120802|TMPO 14.95 1.34 0.113 11.9 1.13e-32 2.27e-29
#> 6 1 ENSG00000163535|SGO2 2.85 2.17 0.186 11.7 1.50e-31 2.50e-28Result columns come straight from DESeq2::results(),
prepended with the group identifier:
| column | description |
|---|---|
group |
foreground cluster |
feature |
gene identifier |
baseMean |
mean of normalized counts |
log2FoldChange |
effect size, group vs background |
lfcSE |
standard error of log2FoldChange
|
stat |
Wald statistic |
pvalue |
nominal p-value |
padj |
Benjamini-Hochberg adjusted p-value |
Step 3: top markers per cluster
top_markers_dds() filters by padj_max and
lfc_min and returns a wide table of the top features per
group:
top5 <- top_markers_dds(res, n = 5, padj_max = 1e-4, lfc_min = 1)
top5 %>%
dplyr::mutate(across(-rank, ~ sub(".*\\|", "", .x)))
#> # A tibble: 5 × 9
#> rank `1` `2` `3` `4` `5` `6` `7` `8`
#> <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 1 DUT CCL4L2 KLRG1 SDF2L1 IL7R LGALS1 LINC00996 ICA1
#> 2 2 LIG1 CCL4 CST7 HSPA5 NT5E <NA> CD247 NFKBID
#> 3 3 TMEM106C CCL3L1 DTHD1 TNFRSF9 FLT3LG <NA> MCTP2 MIR155HG
#> 4 4 H2AFX <NA> CD81 TSPAN17 ANXA1 <NA> IL2RB SOD1
#> 5 5 TMPO <NA> SH2D1A HYOU1 MID1IP1 <NA> PTPN12 ICOSThe signal is biologically reasonable: cluster 1 is dominated by
S-phase / DNA-replication genes (DUT, LIG1,
H2AFX, TMPO), cluster 2 by effector chemokines
(CCL4, CCL3L1).
Step 4: pairwise cluster contrasts
mode = "pairwise" tests each cluster against each other
cluster individually, returning a long-form table with
group1 and group2 columns. This is more
conservative than one-vs-all because a true marker has to differentiate
the cluster from every other cluster, not just the average
background.
The cost grows quickly — N clusters means N × (N−1) directional DESeq2 fits — so it’s worth subsetting to the levels you care about first. Here we contrast clusters 1, 2, and 3:
keep <- cc$meta_data$cluster %in% c("1", "2", "3")
md_p <- cc$meta_data[keep, "cluster", drop = FALSE]
mat_p <- cc$counts_mat[, keep]
res_p <- pseudobulk_deseq2(
~cluster, md_p, mat_p,
verbose = FALSE,
min_counts_per_sample = 10,
present_in_min_samples = 5,
mode = "pairwise"
)
head(res_p)
#> group1 group2 feature baseMean log2FoldChange lfcSE stat pvalue padj
#> 1 1 3 ENSG00000117632|STMN1 142.2 3.55 0.264 13.5 3.03e-41 7.42e-38
#> 2 1 3 ENSG00000176890|TYMS 48.3 5.06 0.382 13.2 5.09e-40 8.31e-37
#> 3 1 2 ENSG00000117632|STMN1 201.2 3.78 0.289 13.1 5.56e-39 5.52e-35
#> 4 1 2 ENSG00000166508|MCM7 30.7 2.60 0.206 12.6 2.03e-36 1.01e-32
#> 5 1 2 ENSG00000128951|DUT 142.1 1.65 0.131 12.5 4.28e-36 1.42e-32
#> 6 1 2 ENSG00000276043|UHRF1 12.7 4.96 0.400 12.4 2.82e-35 7.01e-32summarize_dge_pairs() collapses the directional pairs to
one row per (group, gene). Pass "min" to keep each gene’s
worst comparison (the most conservative effect, useful for
high-confidence markers) or "max" for the best.
summarize_dge_pairs(res_p, "min") %>%
head(10) %>%
dplyr::mutate(feature = sub(".*\\|", "", feature))
#> group feature baseMean log2FoldChange lfcSE stat pvalue padj
#> <char> <char> <num> <num> <num> <num> <num> <num>
#> 1: 1 STMN1 201.22 3.78 0.289 13.1 5.56e-39 5.52e-35
#> 2: 1 TYMS 70.03 4.38 0.381 11.5 1.56e-30 1.94e-27
#> 3: 1 TUBA1B 282.72 2.49 0.222 11.3 2.31e-29 2.55e-26
#> 4: 1 PCNA 49.39 2.16 0.198 10.9 6.87e-28 5.68e-25
#> 5: 1 TUBB 309.64 1.75 0.163 10.7 7.74e-27 5.91e-24
#> 6: 1 MCM7 22.69 2.34 0.219 10.7 1.29e-26 6.30e-24
#> 7: 1 UHRF1 9.14 4.73 0.451 10.5 1.11e-25 4.74e-23
#> 8: 1 TMEM106C 34.53 1.88 0.180 10.4 2.55e-25 1.04e-22
#> 9: 1 DUT 102.42 1.49 0.145 10.2 1.27e-24 4.45e-22
#> 10: 1 ASF1B 13.35 3.48 0.343 10.2 3.25e-24 1.06e-21Step 5: within-cluster Case vs Control
A more interesting biological question than “what defines each
cluster” is “for each cluster, which genes change between irColitis
cases and healthy controls?”. mode = "within" answers this.
The first variable in the formula is the split variable
(the cluster), and the second is the contrast variable
(case). Within each level of the split variable, DESeq2
fits the inner formula (here ~case).
md_w <- cc$meta_data[, c("cluster", "case")]
res_w <- pseudobulk_deseq2(
~cluster + case, md_w, cc$counts_mat,
verbose = FALSE,
min_counts_per_sample = 10,
present_in_min_samples = 5,
mode = "within"
)
head(res_w)
#> group feature baseMean log2FoldChange lfcSE stat pvalue padj
#> 1 1 ENSG00000158517|NCF1 73.62 1.45 0.194 7.50 6.58e-14 3.38e-11
#> 2 1 ENSG00000171867|PRNP 30.67 1.80 0.242 7.45 9.11e-14 4.45e-11
#> 3 1 ENSG00000162496|DHRS3 8.52 2.62 0.386 6.77 1.29e-11 4.51e-09
#> 4 1 ENSG00000185101|ANO9 10.32 1.76 0.261 6.75 1.52e-11 4.63e-09
#> 5 1 ENSG00000265972|TXNIP 202.85 1.34 0.199 6.74 1.60e-11 4.74e-09
#> 6 1 ENSG00000268804|LINC02132 9.84 2.78 0.414 6.72 1.78e-11 4.99e-09Top genes upregulated in cases per cluster:
top_markers_dds(res_w, n = 5, padj_max = 1e-4, lfc_min = 1) %>%
dplyr::mutate(across(-rank, ~ sub(".*\\|", "", .x)))
#> # A tibble: 5 × 9
#> rank `1` `2` `3` `4` `5` `6` `7` `8`
#> <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 1 NCF1 SPINK2 LAIR1 DHRS7 SPINK2 GNPTAB LDLRAD4 RPS12
#> 2 2 PRNP SORBS3 <NA> CA10 <NA> DHRS7 <NA> RPL30
#> 3 3 DHRS3 KIFC3 <NA> NCF1 <NA> TXNIP <NA> EEF1B2
#> 4 4 ANO9 GNPTAB <NA> KIFC3 <NA> FCER1G <NA> MT-ND3
#> 5 5 TXNIP CALHM6 <NA> MYBL1 <NA> TGFBR1 <NA> <NA>Cluster-1 cases up-regulate NCF1, PRNP,
TXNIP — interferon-response and oxidative-stress genes
consistent with inflamed tissue.
Tips
-
Filter low-count genes early.
min_counts_per_sampleandpresent_in_min_samplesare applied before fitting; tighter filters speed up the DESeq2 step substantially. -
Always include batch / channel covariates when
possible. Add them to the formula after the contrast variable
(e.g.
~cluster + channel) if the design supports it. -
collapse_background = TRUE(inmode = "one_vs_all") can stabilize results when donor representation is uneven across clusters. -
mode = "within"supports both 2-level and ordinal contrasts. Two-level contrasts (factor or character) get a Wald test on the level-vs-reference coefficient. Three or more levels are integer-encoded and treated as an ordinal trend.
Session info
sessionInfo()
#> R version 4.5.2 (2025-10-31)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.6.1
#>
#> Matrix products: default
#> BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] en_US/en_US/en_US/C/en_US/en_US
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] dplyr_1.2.1 Matrix_1.7-4 presto_1.0.0 knitr_1.51 ggplot2_4.0.2
#>
#> loaded via a namespace (and not attached):
#> [1] utf8_1.2.6 generics_0.1.4 tidyr_1.3.2
#> [4] SparseArray_1.10.2 DESeq2_1.50.2 lattice_0.22-7
#> [7] magrittr_2.0.5 evaluate_1.0.5 grid_4.5.2
#> [10] RColorBrewer_1.1-3 purrr_1.2.1 scales_1.4.0
#> [13] codetools_0.2-20 abind_1.4-8 cli_3.6.5
#> [16] rlang_1.1.7 XVector_0.50.0 Biobase_2.70.0
#> [19] DelayedArray_0.36.0 withr_3.0.2 otel_0.2.0
#> [22] S4Arrays_1.10.0 tools_4.5.2 parallel_4.5.2
#> [25] BiocParallel_1.44.0 Rhdf5lib_1.32.0 locfit_1.5-9.12
#> [28] SummarizedExperiment_1.40.0 BiocGenerics_0.56.0 vctrs_0.7.2
#> [31] R6_2.6.1 matrixStats_1.5.0 stats4_4.5.2
#> [34] lifecycle_1.0.5 rhdf5_2.54.1 Seqinfo_1.0.0
#> [37] S4Vectors_0.48.0 IRanges_2.44.0 pkgconfig_2.0.3
#> [40] pillar_1.11.1 gtable_0.3.6 data.table_1.18.2.1
#> [43] glue_1.8.0 Rcpp_1.1.1 xfun_0.56
#> [46] tibble_3.3.1 GenomicRanges_1.62.1 tidyselect_1.2.1
#> [49] MatrixGenerics_1.22.0 dichromat_2.0-0.1 rhdf5filters_1.22.0
#> [52] farver_2.1.2 compiler_4.5.2 S7_0.2.1