Perform ClusterDE on a pure cell-line dataset
Dongyuan Song
Department of Genetics & Genome Sciences, UConn HealthBioinformatics IDP, University of California, Los Angelesdongyuansong@ucla.edu
27 August 2026
Source:vignettes/ClusterDE-pure-cellline.Rmd
ClusterDE-pure-cellline.RmdDownload data
We download the cell line data set H2228. The original data is from Tian et al., Nature Methods 2019 as the gold standard for benchmarking the accuracy of clustering. Since the data is from the pure cell line, it should not have cell types, and, of course, between cell type DE genes. Note: it does not mean that there are no variations within one cell line; of course there are, e.g., cell cycle, total UMI, etc. However, these variations do not represent discrete cell groups, and essentially it means you should not use your obtained clusters to explain the variation.
data(cellline, package = "ClusterDE")Run the regular Seurat pipeline
We perform the default Seurat clustering. Please note that ClusterDE is designed for 1 vs 1 comparison; therefore, we set the resolution as 0.2 here to obtain two clusters for illustration purpose.
RNGkind("L'Ecuyer-CMRG")
seed <- 123
set.seed(seed)
cellline <- Seurat::UpdateSeuratObject(cellline)
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in originalexp
#> Validating object structure for Assay 'originalexp'
#> Object representation is consistent with the most current Seurat version
cellline <- Seurat::NormalizeData(cellline)
#> Warning: The `slot` argument of `SetAssayData()` 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: 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.
cellline <- Seurat::FindVariableFeatures(cellline)
cellline <- Seurat::ScaleData(cellline)
#> Centering and scaling data matrix
cellline <- Seurat::RunPCA(cellline)
#> PC_ 1
#> Positive: RPS14, RPL18AP3, RPL36, RPS23, RPL28, LRRC75A-AS1, AC079250.1, FTH1, ZFAS1, EEF2
#> RPL7P9, RPL13A, RPS3AP26, EEF1A1P13, RPS16, RPS23P8, RPL13AP5, RPL29, FTH1P10, RPL13AP25
#> SNHG5, FTH1P8, RPL4, RPS3AP6, AC064799.1, C6orf48, FTH1P7, C1orf56, RPL7AP6, TMSB4X
#> Negative: PSMB2, PSMA7, U2AF1, NUDC, RBM8A, CALM1, BUB3, CLIC1, U2AF1L5, XRCC5
#> VPS29, RBM8B, CACYBP, RPA3, SSBP1, PSMC5, MRPL47, PSMD8, BRIX1, CNIH4
#> PCMT1, PSMD13, CYC1, PRDX2, SEPT7, S100A11, VDAC3, PSME2P2, ZWINT, HMGB1
#> PC_ 2
#> Positive: NACA, RPL7AP6, SKP1, UBA52, BTF3, SSR2, RPL7A, ARPC3, RPL9P9, PPIA
#> PSMD4, EIF1, RPL10, LGALS3BP, RPL10P16, SNRPB2, RPL10P9, S100A11, PPIB, ANXA5
#> EEF2, PSME1, SSBP1, SSR4, RPL7P9, COPE, BSG, MGST1, VPS28, COPS6
#> Negative: SIVA1, HNRNPAB, RPL39L, DEK, CDCA5, TMPO, FAM111A, ASF1B, CENPK, ESCO2
#> BRCA1, H2AFV, RAD51AP1, MT-RNR2, ORC6, CENPX, SNRNP25, FBXO5, RRM1, DIAPH3
#> USP1, CDCA4, TMEM106C, PGP, LSM4, C21orf58, CENPN, BRI3BP, SGO1, CHAF1A
#> PC_ 3
#> Positive: RPL13AP5, RPL13AP25, AC024293.1, RPL29, RPSAP19, RPS5, RPL18, RPS3AP26, RPS3AP6, RPL15
#> RPL28, DRAP1, RPS11, RPL9P9, RPL13AP7, RPS19, DCBLD2, FXYD5, FEN1, SLBP
#> RPS15, COTL1, RPSA, FLNA, RPL7AP6, RPL36, C1orf21, CPA4, ORC6, RPS16
#> Negative: SMIM22, TSPAN13, ST14, PERP, CRB3, MT-CO1, SERINC2, ATP1B1, CDH1, F11R
#> B2M, MT-RNR2, SPINT1, NMB, PLA2G16, SPDEF, CD55, ADGRF1, TSPAN1, LIMA1
#> ERBB3, ERO1A, ASS1, CDA, ALCAM, SYNGR2, MT-CO2, CDH3, C3, LSR
#> PC_ 4
#> Positive: IFNGR1, NAMPT, NAP1L1, CPD, LMAN1, CALR, ITGA2, NAMPTP1, ITM2B, RRM1
#> C3, RHOBTB3, CTHRC1, EEF2, HSD17B11, C1S, IFI16, SMC2, CPE, EPHX1
#> DST, HLA-DMB, NUCB2, MT-ND6, TMEM45A, BRCA1, CDK5RAP2, HINT1, C1R, FAM111A
#> Negative: S100A16, TMA7, TIMM8B, PFN1, SLIRP, GPX1, LAMC2, POLR2L, MRPL52, RPS19
#> CDH1, TOMM40, ATP5MD, HSPE1, NAA10, GPX1P1, RPS16, MRPL12, MCRIP2, PDCD5
#> RPL18, PLEC, S100A13, RPL36AL, LAD1, MGLL, BOLA2B, MISP, MRPL36, SEC61G
#> PC_ 5
#> Positive: BTG1, PPP1R15A, EIF1, JUN, CEBPG, H3F3B, TMEM132A, C6orf48, HIST2H4B, SGK1
#> KPNA4, PMEPA1, KLF6, CDKN1A, WARS, PEA15, GARS, MAP1LC3B, SNHG12, SERTAD1
#> LAMC2, EPB41L4A-AS1, NAP1L1, SNHG5, KLF10, SLC7A5, KIF5B, ATP2B1, EIF5, ABL2
#> Negative: HINT1, COX5B, TXN, SOD1, NDUFA4, NDUFS6, ATP5PO, ATP5MC3, S100A10, ATP6V0E1
#> CYB5A, SLIRP, NDUFB4, ATP5MC1, HSPE1, TXNP6, POMP, POLR2L, RPS14, HSPE1P4
#> RPS15, CBR1, NDUFB3, HSPE1P3, ATP5PD, COX7B, ADGRF1, PSMB9, AC079250.1, NDUFAB1
cellline <- Seurat::FindNeighbors(cellline)
#> Computing nearest neighbor graph
#> Computing SNN
cellline <- Seurat::FindClusters(cellline, resolution = 0.2)
#> Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
#>
#> Number of nodes: 758
#> Number of edges: 24895
#>
#> Running Louvain algorithm...
#> Maximum modularity in 10 random starts: 0.8257
#> Number of communities: 2
#> Elapsed time: 0 seconds
cellline <- Seurat::RunUMAP(cellline, dims = 1:10)
#> Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
#> To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
#> This message will be shown once per session
#> 15:29:02 UMAP embedding parameters a = 0.9922 b = 1.112
#> 15:29:02 Read 758 rows and found 10 numeric columns
#> 15:29:02 Using Annoy for neighbor search, n_neighbors = 30
#> 15:29:02 Building Annoy index with metric = cosine, n_trees = 50
#> 0% 10 20 30 40 50 60 70 80 90 100%
#> [----|----|----|----|----|----|----|----|----|----|
#> **************************************************|
#> 15:29:02 Writing NN index file to temp file /loc/scratch/2560765/RtmpgdLPJH/filedb7725b3f3997
#> 15:29:02 Searching Annoy index using 1 thread, search_k = 3000
#> 15:29:02 Annoy recall = 100%
#> 15:29:03 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
#> 15:29:03 Initializing from normalized Laplacian + noise (using RSpectra)
#> 15:29:03 Commencing optimization for 500 epochs, with 27914 positive edges
#> 15:29:03 Using rng type: pcg
#> 15:29:04 Optimization finished
Seurat::DimPlot(cellline, reduction = "umap") + ggplot2::ggtitle("Clustering result")
From the UMAP, the two clusters seem to be dubious. Although we do not expect the existence of cell types, when we perform Seurat DE test between the two clusters, and we get DE genes with FDR < 0.05. It means that the double-dipping introduces a huge number of false discoveries.
original_deg <- Seurat::FindMarkers(
cellline,
ident.1 = 0,
ident.2 = 1,
min.pct = 0,
logfc.threshold = 0
)
#> 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.
original_deg <- original_deg[original_deg$avg_log2FC > 0,]
original_pval <- original_deg$p_val_adj
names(original_pval) <- rownames(original_deg)
print(paste0("Seurat found ", sum(original_deg$p_val_adj < 0.05), " DEG"))
#> [1] "Seurat found 109 DEG"Find DEGs using ClusterDE
We generate the synthetic null replicate from the target data using
ClusterDE::constructNULL and perform the same preprocessing
and clustering pipeline as the target cell line data.
set.seed(seed)
null_data <- ClusterDE::constructNull(cellline, data_type = "cellline", seed = seed)
#> 1 genes have no more than 2 non-zero values; ignore fitting and return all 0s.
#> 98.5% of genes are used in correlation modelling.
null_obj <- Seurat::CreateSeuratObject(null_data)
null_obj <- Seurat::NormalizeData(null_obj)
#> Normalizing layer: counts
null_obj <- Seurat::FindVariableFeatures(null_obj)
#> Finding variable features for layer counts
null_obj <- Seurat::ScaleData(null_obj)
#> Centering and scaling data matrix
null_obj <- Seurat::RunPCA(null_obj)
#> PC_ 1
#> Positive: VPS29, PSMC5, MRPL47, PSMB2, CACYBP, CNIH4, ATP5PD, SEPT7, PSMA7, SKP1
#> BUB3, ATP5PO, U2AF1, NUCKS1, MGST1, NUDC, RPA3, PSMB9, COX5A, CLIC1
#> RBM8B, NDUFAB1, CALM1, VDAC3, NDUFB4, XRCC5, PSME2, PHB, RBM8A, NOL7
#> Negative: RPL28, RPL13AP5, LRRC75A-AS1, RPL13AP25, EEF2, RPS16, RPS3AP6, RPL13A, TMSB4X, ZFAS1
#> FTH1, RPS11, C9orf3, AHNAK2, RPL13AP7, RPL10, SNHG5, RPL4, AC064799.1, UGCG
#> RPS5, RPL10P16, FTH1P8, PLEC, RPL7AP6, FTH1P7, RPL9P9, RPL10P9, RPS6, C6orf48
#> PC_ 2
#> Positive: NACA, RPL7AP6, RPL9P9, SKP1, RPL7A, SSR2, PSMD4, BTF3, EIF3H, RPL10P16
#> RPL10, UBA52, LGALS3BP, EEF2, ARPC3, RPL10P9, EIF1, SNRPB2, ZFAS1, RPL4
#> SQSTM1, APEX1, BSG, C6orf48, ADH5, COPE, LRRC75A-AS1, PSME1, MGST1, CTSD
#> Negative: SIVA1, DEK, RPL39L, CENPX, HNRNPAB, H2AFV, TMPO, CENPK, SAC3D1, CDCA5
#> USP1, ESCO2, CDCA4, DNAJC9, FBXO5, CENPN, CHCHD2, FAM111A, LSM4, ASF1B
#> SUZ12, PGP, CMC2, BRCA1, CARHSP1, DTYMK, RRM1, CKLF, WDR34, ZWINT
#> PC_ 3
#> Positive: MT-CO1, MTCO1P12, SMIM22, PERP, TSPAN13, ADGRF1, MT-CO2, B2M, ST14, SPINT1
#> NMB, LIMA1, MT-CYB, CPD, GPX3, SERINC2, ASS1, CRB3, VAMP8, GCHFR
#> C3, PPIC, SPDEF, F11R, LINC00467, PLA2G16, SMIM14, EPS8, TMED10, CD55
#> Negative: RPSAP19, RPL18, RPL13AP25, RPL15, RPL13AP5, RPS11, RPSA, RPS5, CFL1, FXYD5
#> RPL13AP7, DRAP1, RPS3AP6, RPL9P9, RPS19, SLBP, NAP1L1, RPL28, DCBLD2, PDCD5
#> PEBP1, FLNA, RPS3, RPS16, NME1, RPL10P16, HNRNPD, BASP1, GLIPR1, CACYBP
#> PC_ 4
#> Positive: IFNGR1, NAMPT, CPD, NAP1L1, BRCA1, RFC4, RRM1, C1S, CTHRC1, C1R
#> NAMPTP1, RHOBTB3, CARHSP1, DTL, MCM10, HSD17B4, CASP4, AC092683.1, NCOA7, SYNE2
#> RMI1, ITGA2, EZH2, ITM2B, CDT1, MCM6, C3, UGP2, EPHX1, POLD3
#> Negative: S100A16, SLIRP, RPL36AL, CYC1, TOMM40, TSSC4, PFN1, POLR2L, H3F3B, DYNLL1
#> NDUFB9, RPS16, RPS19, SNRPF, TMA7, COX7B, MRPL21, BOLA2B, HSPE1, AHNAK
#> PDCD5, MGLL, HSPE1P4, PFDN2, RPL18, NDUFS6, NAA10, HSPE1P3, PPP1R14B, ATP5MC3
#> PC_ 5
#> Positive: BTG1, KIF14, C6orf48, H3F3B, EIF1, KIF5B, CDCA2, SRGAP2B, GARS, KPNA4
#> RAD21, NKAP, PPP1R15A, KIF4A, SPAG9, MAP1LC3B, SMOX, SARS, TCEA1P2, SGK1
#> CEBPG, CDKN1A, SDC4, CDCA8, RPL10P16, TCEAL3, MARS, MARCKS, C1orf56, ARID4A
#> Negative: ATP6V0E1, TXN, NDUFA4, POMP, HSPE1P3, ATP5MPL, SOD1, ATP5MC3, NDUFS6, COX5B
#> NDUFB3, SH3BGRL3, NDUFB4, S100A10, ATP5PO, TRMT112, NDUFAB1, POLR2L, HSPE1, HINT1
#> ATP5MC1, HSPE1P4, CBR1, MRPL17, SEC61G, TXNP6, SLIRP, ATP5PD, COTL1, COX4I1
null_obj <- Seurat::FindNeighbors(null_obj)
#> Computing nearest neighbor graph
#> Computing SNN
null_obj <- Seurat::FindClusters(null_obj, resolution = 0.3)
#> Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
#>
#> Number of nodes: 758
#> Number of edges: 30040
#>
#> Running Louvain algorithm...
#> Maximum modularity in 10 random starts: 0.7174
#> Number of communities: 2
#> Elapsed time: 0 secondsWe then perform Wilcoxon DE test on the null data and compare it with the p-values of target genes. We can observe that after ClusterDE FDR control, no genes pass the test which means we don’t detect any DEG in this cell line data.
null_deg <- Seurat::FindMarkers(
null_obj,
ident.1 = 0,
ident.2 = 1,
min.pct = 0,
logfc.threshold = 0
)
null_deg <- null_deg[null_deg$avg_log2FC > 0,]
null_pval <- null_deg$p_val_adj
names(null_pval) <- rownames(null_deg)
deg <- ClusterDE::callDE(original_pval, null_pval)
print(paste0("ClusterDE found ", sum(deg$record >= 0.5), " DEG"))
#> [1] "ClusterDE found 0 DEG"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] RColorBrewer_1.1-3 jsonlite_2.0.0 magrittr_2.0.3
#> [4] spatstat.utils_3.1-4 farver_2.1.2 rmarkdown_2.29
#> [7] fs_1.6.6 ragg_1.4.0 vctrs_0.6.5
#> [10] ROCR_1.0-11 memoise_2.0.1 spatstat.explore_3.4-3
#> [13] htmltools_0.5.8.1 sass_0.4.10 sctransform_0.4.2
#> [16] parallelly_1.45.0 KernSmooth_2.23-26 bslib_0.9.0
#> [19] htmlwidgets_1.6.4 desc_1.4.3 ica_1.0-3
#> [22] plyr_1.8.9 plotly_4.10.4 zoo_1.8-14
#> [25] cachem_1.1.0 igraph_2.1.4 mime_0.13
#> [28] lifecycle_1.0.4 pkgconfig_2.0.3 Matrix_1.7-5
#> [31] R6_2.6.1 fastmap_1.2.0 fitdistrplus_1.2-2
#> [34] shiny_1.10.0 digest_0.6.37 patchwork_1.3.0
#> [37] Seurat_5.3.0 tensor_1.5 RSpectra_0.16-2
#> [40] irlba_2.3.5.1 kde1d_1.1.1 textshaping_1.0.1
#> [43] labeling_0.4.3 progressr_0.15.1 spatstat.sparse_3.1-0
#> [46] coop_0.6-3 httr_1.4.7 polyclip_1.10-7
#> [49] abind_1.4-8 compiler_4.5.1 withr_3.0.2
#> [52] backports_1.5.0 fastDummies_1.7.5 MASS_7.3-65
#> [55] tools_4.5.1 lmtest_0.9-40 httpuv_1.6.16
#> [58] future.apply_1.20.0 goftest_1.2-3 glue_1.8.0
#> [61] nlme_3.1-168 promises_1.3.3 grid_4.5.1
#> [64] checkmate_2.3.2 Rtsne_0.17 cluster_2.1.8.1
#> [67] reshape2_1.4.4 generics_0.1.4 gtable_0.3.6
#> [70] spatstat.data_3.1-6 tidyr_1.3.1 data.table_1.17.4
#> [73] sp_2.2-0 spatstat.geom_3.4-1 RcppAnnoy_0.0.22
#> [76] ggrepel_0.9.6 RANN_2.6.2 pillar_1.10.2
#> [79] stringr_1.5.1 spam_2.11-1 RcppHNSW_0.6.0
#> [82] limma_3.64.1 later_1.4.2 splines_4.5.1
#> [85] dplyr_1.1.4 lattice_0.22-9 survival_3.8-3
#> [88] deldir_2.0-4 tidyselect_1.2.1 miniUI_0.1.2
#> [91] pbapply_1.7-2 knitr_1.50 gridExtra_2.3
#> [94] bookdown_0.43 scattermore_1.2 xfun_0.52
#> [97] statmod_1.5.0 matrixStats_1.5.0 ClusterDE_0.99.4
#> [100] stringi_1.8.7 lazyeval_0.2.2 yaml_2.3.10
#> [103] evaluate_1.0.3 codetools_0.2-20 tibble_3.3.0
#> [106] BiocManager_1.30.27 cli_3.6.5 uwot_0.2.3
#> [109] xtable_1.8-4 reticulate_1.42.0 systemfonts_1.2.3
#> [112] jquerylib_0.1.4 dichromat_2.0-0.1 Rcpp_1.0.14
#> [115] globals_0.18.0 spatstat.random_3.4-1 png_0.1-8
#> [118] rngWELL_0.10-10 spatstat.univar_3.1-3 parallel_4.5.1
#> [121] randtoolbox_2.0.5 assertthat_0.2.1 pkgdown_2.0.9
#> [124] ggplot2_3.5.2 presto_1.0.0 mvnfast_0.2.8
#> [127] dotCall64_1.2 bettermc_1.2.2.9000 gamlss.dist_6.1-1
#> [130] listenv_0.9.1 viridisLite_0.4.2 scales_1.4.0
#> [133] ggridges_0.5.6 SeuratObject_5.1.0 purrr_1.0.4
#> [136] crayon_1.5.3 rlang_1.1.6 rvinecopulib_0.7.3.1.0
#> [139] cowplot_1.1.3