Nothing
Pseudobulk differential expression with DESeq2
In presto: Fast Functions for Differential Expression Using Wilcoxon and AUC
Why 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. It is fetched with
load_ircolitis_cd8(), a small build-time helper (see
vignettes/ircolitis.R in the package sources), not part of presto.
library(presto)
library(Matrix)
library(dplyr)
## load_ircolitis_cd8() is defined in vignettes/ircolitis.R (build-time
## only), not exported by presto.
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 4
The 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 110
We 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 5
Step 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-28
Result 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 ICOS
The 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-32
summarize_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-21
Step 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-09
Top 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_sample and
present_in_min_samples are 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 (in mode = "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.1.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
Try the presto package in your browser
Any scripts or data that you put into this service are public.
presto documentation built on Sept. 30, 2026, 5:13 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Why 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. It is fetched with
load_ircolitis_cd8(), a small build-time helper (see
vignettes/ircolitis.R in the package sources), not part of presto.
library(presto) library(Matrix) library(dplyr) ## load_ircolitis_cd8() is defined in vignettes/ircolitis.R (build-time ## only), not exported by presto. 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 4
The 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 110
We 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 5
Step 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-28
Result 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 ICOS
The 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-32
summarize_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-21
Step 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-09
Top 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.1.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
Try the presto package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.