Perform ClusterDE on a single-cell multiome dataset
Yihui Cen
ycen@fredhutch.org27 August 2026
Source:vignettes/ClusterDE-scmultiome.Rmd
ClusterDE-scmultiome.RmdDownload data
Here, we load single-cell multiome PBMC dataset done by 10X Genomics, originally from SeuratData. Specifically, we extract the CD4 Naive and CD4 TCM cell types and randomly sampled 500 cells for each cell type. To further save computational time, features with relatively low expression are filtered out. All the cells are paired in both the RNA and ATAC modalities.
Generate synthetic null data for the RNA data
We first check on the RNA data. Generate 20 null replicates using ClusterDE stable to eliminate randomness of null generation.
RNGkind("L'Ecuyer-CMRG")
seed <- 123
set.seed(seed)
null_data <- ClusterDE::constructNull(
obj = pbmc_rna,
n_rep = 20,
seed = seed
)
#> 123 genes have no more than 2 non-zero values; ignore fitting and return all 0s.
#> 26.8% of genes are used in correlation modelling.Next, we run the seurat pipeline on the null data
set.seed(seed)
null_seurat <- lapply(seq_along(null_data), function(i) {
obj <- Seurat::CreateSeuratObject(
counts = null_data[[i]],
assay = "RNA",
project = paste0("null_", i)
)
obj <- Seurat::NormalizeData(
obj,
normalization.method = "LogNormalize",
scale.factor = 10000,
verbose = FALSE
)
obj <- Seurat::FindVariableFeatures(
obj,
selection.method = "vst",
nfeatures = 3000,
verbose = FALSE
)
obj <- Seurat::ScaleData(
obj,
features = Seurat::VariableFeatures(obj),
verbose = FALSE
)
obj <- Seurat::RunPCA(
obj,
features = Seurat::VariableFeatures(obj),
npcs = 50,
reduction.name = "pca",
reduction.key = "PC_",
verbose = FALSE
)
obj <- Seurat::FindNeighbors(
obj,
reduction = "pca",
dims = 1:30,
verbose = FALSE
)
obj <- Seurat::FindClusters(
obj,
resolution = 0.4
)
obj
})
names(null_seurat) <- paste0("null_", seq_along(null_seurat))Find DEGs on RNA data
To find DEGs with ClusterDE, we first find DEGs with the original data and obtain the target scores.
original_levels <- levels(
droplevels(Seurat::Idents(pbmc_rna))
)
original_markers <- Seurat::FindMarkers(
object = pbmc_rna,
assay = "RNA",
ident.1 = original_levels[1],
ident.2 = original_levels[2],
test.use = "wilcox",
min.pct = 0,
logfc.threshold = 0,
only.pos = FALSE
)
#> Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
#> ℹ Please use the `layer` argument instead.
#> ℹ The deprecated feature was likely used in the Seurat package.
#> Please report the issue at <https://github.com/satijalab/seurat/issues>.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
#> ℹ Please use `rlang::check_installed()` instead.
#> ℹ The deprecated feature was likely used in the Seurat package.
#> Please report the issue at <https://github.com/satijalab/seurat/issues>.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
print(paste0("Wilcox found ", sum(p.adjust(original_markers$p_val, method = "BH") < 0.05), " DE"))
#> [1] "Wilcox found 676 DE"Then, we find DEGs for the synthetic null data
null_markers_list <- lapply(
null_seurat,
function(null_obj) {
null_levels <- levels(
droplevels(Seurat::Idents(null_obj))
)
null_markers <- Seurat::FindMarkers(
object = null_obj,
assay = "RNA",
ident.1 = null_levels[1],
ident.2 = null_levels[2],
test.use = "wilcox",
min.pct = 0,
logfc.threshold = 0,
only.pos = FALSE
)
null_markers$gene <- rownames(null_markers)
null_markers
}
)
names(null_markers_list) <- names(null_seurat)
nullScores_raw <- lapply(
null_markers_list,
function(null_markers) {
scores <- null_markers$p_val
names(scores) <- rownames(null_markers)
scores[
is.finite(scores)
]
}
)After that, we use ClusterDE to call DEG results on the RNA data.
original_markers$gene <- rownames(original_markers)
targetScores <- original_markers$p_val
names(targetScores) <- rownames(original_markers)
de <- ClusterDE::callDE(
targetScores = targetScores,
nullScores = nullScores_raw,
FDR = 0.05,
correct = TRUE,
nlogTrans = TRUE,
threshold = "DS"
)
print(paste0("ClusterDE found ", sum(de$record >= 0.5), " DE"))
#> [1] "ClusterDE found 96 DE"Generate synthetic null data for the ATAC data
Then, since the RNA data shows clear two clusters, we check on the ATAC data for differential accessibility analysis. Here, we generate 20 null replicates using ClusterDE stable to eliminate randomness of null generation.
null_data_atac <- ClusterDE::constructNull(
obj = pbmc_atac,
data_type = "scATAC",
n_rep = 20,
seed = seed
)
#> Loading required namespace: Signac
#> 100% of genes are used in correlation modelling.Next, we run the seurat pipeline on the null data
set.seed(seed)
null_signac <- lapply(seq_along(null_data_atac), function(i) {
atac_assay <- Signac::CreateChromatinAssay(
counts = null_data_atac[[i]],
sep = c(":", "-"),
min.cells = 0,
min.features = 0
)
obj <- Seurat::CreateSeuratObject(
counts = atac_assay,
assay = "ATAC",
project = paste0("null_atac_", i)
)
obj <- Signac::RunTFIDF(
obj
)
obj <- Signac::FindTopFeatures(
obj,
min.cutoff = "q0"
)
obj <- Signac::RunSVD(
obj,
assay = "ATAC",
n = 50,
reduction.name = "lsi",
reduction.key = "LSI_"
)
obj <- Seurat::FindNeighbors(
obj,
reduction = "lsi",
dims = 2:30,
verbose = FALSE
)
obj <- Seurat::FindClusters(
obj,
resolution = 0.7
)
obj
})
names(null_signac) <- paste0(
"null_atac_",
seq_along(null_signac)
)Find DE peaks on ATAC data
To find DE peaks with ClusterDE, we first find DE peaks with the original data and obtain the target scores.
original_levels_atac <- levels(
droplevels(Seurat::Idents(pbmc_atac))
)
original_markers_atac <- Seurat::FindMarkers(
object = pbmc_atac,
assay = "ATAC",
ident.1 = original_levels_atac[1],
ident.2 = original_levels_atac[2],
test.use = "wilcox",
min.pct = 0,
logfc.threshold = 0,
only.pos = FALSE
)
print(
paste0(
"Wilcox found ",
sum(
p.adjust(
original_markers_atac$p_val,
method = "BH"
) < 0.05
),
" DA peaks"
)
)
#> [1] "Wilcox found 262 DA peaks"Then, we find DE peaks for the synthetic null data
null_markers_atac_list <- lapply(
null_signac,
function(null_obj) {
null_levels <- levels(
droplevels(Seurat::Idents(null_obj))
)
null_markers <- Seurat::FindMarkers(
object = null_obj,
assay = "ATAC",
ident.1 = null_levels[1],
ident.2 = null_levels[2],
test.use = "wilcox",
min.pct = 0,
logfc.threshold = 0,
only.pos = FALSE
)
null_markers$peak <- rownames(null_markers)
null_markers
}
)
names(null_markers_atac_list) <- names(null_signac)
nullScores_atac_raw <- lapply(
null_markers_atac_list,
function(null_markers) {
scores <- null_markers$p_val
names(scores) <- rownames(null_markers)
scores[
is.finite(scores)
]
}
)After that, we use ClusterDE to call DEG results on the ATAC data.
original_markers_atac$peak <- rownames(
original_markers_atac
)
targetScores_atac <- original_markers_atac$p_val
names(targetScores_atac) <- rownames(
original_markers_atac
)
de_atac <- ClusterDE::callDE(
targetScores = targetScores_atac,
nullScores = nullScores_atac_raw,
FDR = 0.05,
nlogTrans = TRUE,
threshold = "DS"
)
print(
paste0(
"ClusterDE found ",
sum(de_atac$record >= 0.5),
" DA peaks"
)
)
#> [1] "ClusterDE found 44 DA peaks"Session information
sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS/LAPACK: FlexiBLAS OPENBLAS; LAPACK version 3.12.0
#>
#> Random number generation:
#> RNG: L'Ecuyer-CMRG
#> Normal: Inversion
#> Sample: Rejection
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/Los_Angeles
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] future_1.58.0 BiocStyle_2.36.0
#>
#> loaded via a namespace (and not attached):
#> [1] RcppAnnoy_0.0.22 splines_4.5.1 later_1.4.2
#> [4] bitops_1.0-9 tibble_3.3.0 polyclip_1.10-7
#> [7] fastDummies_1.7.5 lifecycle_1.0.4 globals_0.18.0
#> [10] lattice_0.22-9 MASS_7.3-65 backports_1.5.0
#> [13] magrittr_2.0.3 limma_3.64.1 plotly_4.10.4
#> [16] sass_0.4.10 rmarkdown_2.29 jquerylib_0.1.4
#> [19] yaml_2.3.10 httpuv_1.6.16 Seurat_5.3.0
#> [22] sctransform_0.4.2 spam_2.11-1 sp_2.2-0
#> [25] spatstat.sparse_3.1-0 reticulate_1.42.0 gld_2.6.7
#> [28] cowplot_1.1.3 pbapply_1.7-2 RColorBrewer_1.1-3
#> [31] abind_1.4-8 Rtsne_0.17 GenomicRanges_1.60.0
#> [34] purrr_1.0.4 presto_1.0.0 BiocGenerics_0.54.0
#> [37] kde1d_1.1.1 GenomeInfoDbData_1.2.14 IRanges_2.42.0
#> [40] S4Vectors_0.46.0 ggrepel_0.9.6 irlba_2.3.5.1
#> [43] listenv_0.9.1 spatstat.utils_3.1-4 goftest_1.2-3
#> [46] RSpectra_0.16-2 spatstat.random_3.4-1 fitdistrplus_1.2-2
#> [49] parallelly_1.45.0 pkgdown_2.0.9 codetools_0.2-20
#> [52] RcppRoll_0.3.1 tidyselect_1.2.1 UCSC.utils_1.4.0
#> [55] farver_2.1.2 randtoolbox_2.0.5 matrixStats_1.5.0
#> [58] stats4_4.5.1 spatstat.explore_3.4-3 jsonlite_2.0.0
#> [61] PairedData_1.1.1 e1071_1.7-16 progressr_0.15.1
#> [64] ggridges_0.5.6 survival_3.8-3 systemfonts_1.2.3
#> [67] bettermc_1.2.2.9000 tools_4.5.1 ragg_1.4.0
#> [70] ica_1.0-3 Rcpp_1.0.14 glue_1.8.0
#> [73] gridExtra_2.3 xfun_0.52 mvnfast_0.2.8
#> [76] GenomeInfoDb_1.44.0 dplyr_1.1.4 withr_3.0.2
#> [79] BiocManager_1.30.27 fastmap_1.2.0 digest_0.6.37
#> [82] R6_2.6.1 mime_0.13 textshaping_1.0.1
#> [85] scattermore_1.2 tensor_1.5 dichromat_2.0-0.1
#> [88] spatstat.data_3.1-6 tidyr_1.3.1 generics_0.1.4
#> [91] data.table_1.17.4 class_7.3-23 httr_1.4.7
#> [94] htmlwidgets_1.6.4 rngWELL_0.10-10 uwot_0.2.3
#> [97] pkgconfig_2.0.3 gtable_0.3.6 lmtest_0.9-40
#> [100] XVector_0.48.0 htmltools_0.5.8.1 dotCall64_1.2
#> [103] bookdown_0.43 SeuratObject_5.1.0 scales_1.4.0
#> [106] lmom_3.2 png_0.1-8 spatstat.univar_3.1-3
#> [109] knitr_1.50 Signac_1.14.0 reshape2_1.4.4
#> [112] checkmate_2.3.2 nlme_3.1-168 proxy_0.4-27
#> [115] cachem_1.1.0 zoo_1.8-14 stringr_1.5.1
#> [118] KernSmooth_2.23-26 parallel_4.5.1 miniUI_0.1.2
#> [121] desc_1.4.3 pillar_1.10.2 grid_4.5.1
#> [124] vctrs_0.6.5 RANN_2.6.2 promises_1.3.3
#> [127] xtable_1.8-4 cluster_2.1.8.1 gamlss.dist_6.1-1
#> [130] evaluate_1.0.3 mvtnorm_1.3-3 cli_3.6.5
#> [133] compiler_4.5.1 Rsamtools_2.24.0 rlang_1.1.6
#> [136] crayon_1.5.3 future.apply_1.20.0 plyr_1.8.9
#> [139] fs_1.6.6 stringi_1.8.7 BiocParallel_1.42.1
#> [142] viridisLite_0.4.2 deldir_2.0-4 assertthat_0.2.1
#> [145] Biostrings_2.76.0 rvinecopulib_0.7.3.1.0 lazyeval_0.2.2
#> [148] coop_0.6-3 spatstat.geom_3.4-1 Matrix_1.7-5
#> [151] RcppHNSW_0.6.0 patchwork_1.3.0 ggplot2_3.5.2
#> [154] statmod_1.5.0 shiny_1.10.0 ROCR_1.0-11
#> [157] igraph_2.1.4 memoise_2.0.1 ClusterDE_0.99.4
#> [160] bslib_0.9.0 fastmatch_1.1-6