Perform ClusterDE on a microbiome dataset
Christy Lee
lee700@g.ucla.eduYiwen Yang
yyang5@fredhutch.org27 August 2026
Source:vignettes/ClusterDE-microbiome.Rmd
ClusterDE-microbiome.RmdDownload data
Load lung cancer WGS microbiome profile processed in Seurat object, originally available at https://github.com/yge15/TCGA_Microbial_Content.
data(lung, package = "ClusterDE")Generate synthetic null data
Generate 20 null replicates using ClusterDE stable to eliminate randomness of null generation. Note that for microbiome data pass the meta column name that stores log sequencing depth to both formula and other_covariates.
RNGkind("L'Ecuyer-CMRG")
seed <- 123
null_data <- ClusterDE::constructNull(
lung,
data_type = "microbiome",
formula = "log_depth",
other_covariates = "log_depth",
n_rep = 20,
seed = seed
)
#> Registered S3 method overwritten by 'scDesign3':
#> method from
#> predict.gamlss gamlss
#> Input Data Construction Start
#> Input Data Construction End
#> Start Marginal Fitting
#> Marginal Fitting End
#> Start Copula Fitting
#> Convert Residuals to Multivariate Gaussian
#> Converting End
#> Copula group 1 starts
#> Copula Fitting End
#> Start Parameter Extraction
#> Parameter
#> Extraction End
#> Start Generate New Data
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> New Data Generating EndSample clustering
Perform Kmeans (K = 2) clustering on both real data and synthetic null data, the clustering method can be any method that is used in the original study.
DA analysis
Use any established DA methods to find differentially abundant taxa in both real and synthetic data. Using Wilcoxon, we can find 70 DA comparing the two clusters:
original_pval <- apply(Seurat::GetAssayData(lung, layer = "data"), 1, function(x) {
stats <- wilcox.test(x ~ clusters)
stats$p.value
})
original_pval <- p.adjust(original_pval, method = "BH")
print(paste0("Wilcox found ", sum(original_pval < 0.05), " DA"))
#> [1] "Wilcox found 70 DA"Next, apply the same DA method on all null replicates.
null_pvals <- lapply(seq_along(null_data), function(i) {
null_pval <- apply(null_data[[i]], 1, function(x) {
stats <- wilcox.test(x ~ null_clusters[[i]])
stats$p.value
})
p.adjust(null_pval, method = "BH")
})FDR control
Finally, compare original p-values and null p-values using ClusterDE FDR control, we observed 12 real DA:
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] BiocStyle_2.36.0
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 jsonlite_2.0.0
#> [3] magrittr_2.0.3 spatstat.utils_3.1-4
#> [5] farver_2.1.2 rmarkdown_2.29
#> [7] fs_1.6.6 ragg_1.4.0
#> [9] vctrs_0.6.5 ROCR_1.0-11
#> [11] memoise_2.0.1 spatstat.explore_3.4-3
#> [13] S4Arrays_1.8.1 htmltools_0.5.8.1
#> [15] SparseArray_1.8.0 sass_0.4.10
#> [17] sctransform_0.4.2 parallelly_1.45.0
#> [19] KernSmooth_2.23-26 bslib_0.9.0
#> [21] htmlwidgets_1.6.4 desc_1.4.3
#> [23] ica_1.0-3 plyr_1.8.9
#> [25] plotly_4.10.4 zoo_1.8-14
#> [27] cachem_1.1.0 igraph_2.1.4
#> [29] mime_0.13 lifecycle_1.0.4
#> [31] pkgconfig_2.0.3 Matrix_1.7-5
#> [33] R6_2.6.1 fastmap_1.2.0
#> [35] GenomeInfoDbData_1.2.14 MatrixGenerics_1.20.0
#> [37] fitdistrplus_1.2-2 future_1.58.0
#> [39] shiny_1.10.0 digest_0.6.37
#> [41] S4Vectors_0.46.0 patchwork_1.3.0
#> [43] Seurat_5.3.0 tensor_1.5
#> [45] RSpectra_0.16-2 irlba_2.3.5.1
#> [47] GenomicRanges_1.60.0 textshaping_1.0.1
#> [49] progressr_0.15.1 spatstat.sparse_3.1-0
#> [51] mgcv_1.9-3 httr_1.4.7
#> [53] polyclip_1.10-7 abind_1.4-8
#> [55] compiler_4.5.1 withr_3.0.2
#> [57] backports_1.5.0 fastDummies_1.7.5
#> [59] MASS_7.3-65 DelayedArray_0.34.1
#> [61] tools_4.5.1 lmtest_0.9-40
#> [63] httpuv_1.6.16 future.apply_1.20.0
#> [65] goftest_1.2-3 glue_1.8.0
#> [67] nlme_3.1-168 promises_1.3.3
#> [69] grid_4.5.1 checkmate_2.3.2
#> [71] Rtsne_0.17 cluster_2.1.8.1
#> [73] reshape2_1.4.4 generics_0.1.4
#> [75] gtable_0.3.6 spatstat.data_3.1-6
#> [77] tidyr_1.3.1 data.table_1.17.4
#> [79] XVector_0.48.0 sp_2.2-0
#> [81] BiocGenerics_0.54.0 spatstat.geom_3.4-1
#> [83] RcppAnnoy_0.0.22 ggrepel_0.9.6
#> [85] RANN_2.6.2 pillar_1.10.2
#> [87] stringr_1.5.1 spam_2.11-1
#> [89] RcppHNSW_0.6.0 later_1.4.2
#> [91] splines_4.5.1 dplyr_1.1.4
#> [93] lattice_0.22-9 survival_3.8-3
#> [95] deldir_2.0-4 gamlss.data_6.0-6
#> [97] tidyselect_1.2.1 SingleCellExperiment_1.30.1
#> [99] miniUI_0.1.2 pbapply_1.7-2
#> [101] knitr_1.50 gridExtra_2.3
#> [103] bookdown_0.43 IRanges_2.42.0
#> [105] SummarizedExperiment_1.38.1 scattermore_1.2
#> [107] stats4_4.5.1 xfun_0.52
#> [109] Biobase_2.68.0 matrixStats_1.5.0
#> [111] UCSC.utils_1.4.0 stringi_1.8.7
#> [113] ClusterDE_0.99.4 lazyeval_0.2.2
#> [115] yaml_2.3.10 evaluate_1.0.3
#> [117] codetools_0.2-20 tibble_3.3.0
#> [119] BiocManager_1.30.27 cli_3.6.5
#> [121] uwot_0.2.3 xtable_1.8-4
#> [123] reticulate_1.42.0 systemfonts_1.2.3
#> [125] jquerylib_0.1.4 GenomeInfoDb_1.44.0
#> [127] dichromat_2.0-0.1 Rcpp_1.0.14
#> [129] globals_0.18.0 spatstat.random_3.4-1
#> [131] scDesign3_1.7.6 png_0.1-8
#> [133] spatstat.univar_3.1-3 parallel_4.5.1
#> [135] pkgdown_2.0.9 ggplot2_3.5.2
#> [137] mclust_6.1.1 dotCall64_1.2
#> [139] bettermc_1.2.2.9000 listenv_0.9.1
#> [141] gamlss.dist_6.1-1 mvtnorm_1.3-3
#> [143] viridisLite_0.4.2 gamlss_5.4-22
#> [145] scales_1.4.0 ggridges_0.5.6
#> [147] crayon_1.5.3 SeuratObject_5.1.0
#> [149] purrr_1.0.4 rlang_1.1.6
#> [151] cowplot_1.1.3