Introduction
In single-cell RNA sequencing (scRNA-seq), evaluating the quality of simulated count data is important before testing analysis tools such as clustering, batch correction, or differential expression algorithms.
This tutorial shows how to evaluate synthetic scRNA-seq data
step-by-step using scSimEval without
needing any artificial ground truth.
1. Loading the Package and Example Dataset
scSimEval includes pre-packaged real and simulated
scRNA-seq count matrices (example_scrna):
library(scSimEval)
# Load packaged scRNA-seq example data
data(example_scrna)
# Examine dimensions and structure
cat("Reference scRNA-seq dimensions:", dim(example_scrna$ref), "\n")
#> Reference scRNA-seq dimensions: 60 80
cat("Simulated scRNA-seq dimensions:", dim(example_scrna$sim), "\n")
#> Simulated scRNA-seq dimensions: 60 80
cat("Cell types represented:", levels(example_scrna$cell_types), "\n")
#> Cell types represented: TypeA TypeB
cat("Batches represented:", levels(example_scrna$batch_info), "\n")
#> Batches represented: Batch1 Batch2The dataset consists of: * ref: Feature-by-cell matrix
of real experimental scRNA-seq read counts (60 genes
80 cells). * sim: Feature-by-cell matrix of simulated
counts generated by a simulation method. * cell_types:
Factor of cell type labels. * batch_info: Factor of
technical batch identifiers.
2. Step-by-Step Category Evaluation
2.1 Distribution Properties & Summary Measures (Categories 1 and 2)
We evaluate how closely the simulated gene expression distributions
match the real reference data using
evaluate_simulation_accuracy():
unimodal_res <- evaluate_simulation_accuracy(
ref_data = example_scrna$ref,
sim_data = example_scrna$sim,
compute_bivariate = FALSE,
verbose = FALSE
)
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
# Inspect statistical distance metrics for library size
lib_summary <- subset(unimodal_res$metrics_summary_table, Property == "library_size")
head(lib_summary, 6)
#> Category Property Metric Value
#> 1 Distributional Properties library_size MAD 1.500000e+01
#> 2 Distributional Properties library_size KS 2.125000e-01
#> 3 Distributional Properties library_size MAE 1.506250e+01
#> 4 Distributional Properties library_size RMSE 1.571504e+01
#> 5 Distributional Properties library_size OV 8.370258e-01
#> 6 Distributional Properties library_size Bhattacharyya 1.188839e-042.2 Cell Clustering and Concordance (Category 3)
We evaluate whether the simulated data preserves distinct cell types
in unsupervised and supervised clustering using
evaluate_clustering_metrics():
clust_res <- evaluate_clustering_metrics(
data = example_scrna$sim,
cell_types = example_scrna$cell_types,
ref_data = example_scrna$ref
)
# Preview clustering concordance metrics
cat("Average Silhouette Width (Simulated):", round(clust_res$silhouette, 4), "\n")
#> Average Silhouette Width (Simulated): 0.0337
cat("Davies-Bouldin Index:", round(clust_res$davies_bouldin, 4), "\n")
#> Davies-Bouldin Index: 4.9116
cat("Adjusted Rand Index (ARI):", round(clust_res$ARI, 4), "\n")
#> Adjusted Rand Index (ARI): 0.0968
cat("Normalized Mutual Information (NMI):", round(clust_res$NMI, 4), "\n")
#> Normalized Mutual Information (NMI): 0.11382.3 Batch Effects and Confounder Mixing (Category 4)
We test whether technical batch differences are properly mixed
without masking biological differences using
evaluate_batch_metrics():
batch_res <- evaluate_batch_metrics(
data = example_scrna$sim,
batch_info = example_scrna$batch_info,
cell_types = example_scrna$cell_types,
verbose = FALSE
)
# Preview batch mixing metrics
cat("Batch Shannon Entropy:", round(batch_res$shannon_entropy, 4), "\n")
#> Batch Shannon Entropy: 0.9796
cat("Principal Component Regression (PCR R2):", round(batch_res$pcr_r2, 4), "\n")
#> Principal Component Regression (PCR R2): 0.0121
cat("Cross-Batch Transfer Accuracy:", round(batch_res$cross_batch_accuracy, 4), "\n")
#> Cross-Batch Transfer Accuracy: 0.61252.4 Differential Expression (DEG) Fidelity (Category 5)
We assess marker gene effect sizes and top marker overlap using
evaluate_deg_fidelity():
deg_res <- evaluate_deg_fidelity(
ref_data = example_scrna$ref,
sim_data = example_scrna$sim,
ref_celltypes = example_scrna$cell_types,
sim_celltypes = example_scrna$cell_types
)
# Display DEG summary metrics
head(deg_res$deg_summary_table)
#> Framework Metric Value
#> 1 Simpipe (Duo et al., 2024) DEG_Ratio 0.0000000
#> 2 Simpipe (Duo et al., 2024) PValue_Uniformity_Chisq 8.0000000
#> 3 Simpipe (Duo et al., 2024) Distribution_Score 1.0000000
#> 4 Simpipe (Duo et al., 2024) Classifier_Accuracy 0.5625000
#> 5 Simpipe (Duo et al., 2024) Classifier_Macro_F1 0.5607843
#> 6 Simpipe (Duo et al., 2024) Classifier_Macro_Recall 0.5793651
#> Direction
#> 1 Target = 1.0 (Optimal balance)
#> 2 Lower is better (closer to Uniform)
#> 3 Higher is better (1.0 = Uniform null)
#> 4 Higher is better (Group predictability)
#> 5 Higher is better (Balanced F1)
#> 6 Higher is better (Sensitivity)2.5 Trajectory and Lineage Dynamics (Category 6)
We infer developmental pseudotime and lineage branching directly from
scRNA-seq expression counts using
evaluate_trajectory_metrics():
traj_res <- evaluate_trajectory_metrics(
ref_data = example_scrna$ref,
sim_data = example_scrna$sim,
cell_types_ref = example_scrna$cell_types,
cell_types_sim = example_scrna$cell_types
)
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
cat("Pseudotime Spearman Correlation:", round(traj_res$pseudotime_correlation, 4), "\n")
#> Pseudotime Spearman Correlation: 1
cat("Lineage Tree Branch Height RMSE:", round(as.numeric(traj_res$tree_height_rmse), 4), "\n")
#> Lineage Tree Branch Height RMSE: 1.86723. Master Unimodal Pipeline:
evaluate_simulation_accuracy
Instead of running individual steps, execute the entire unimodal evaluation suite in a single coordinated command:
eval_res <- evaluate_simulation_accuracy(
ref_data = example_scrna$ref,
sim_data = example_scrna$sim,
elapsed_time = 12.5,
memory_mb = 340.2,
compute_bivariate = FALSE,
verbose = FALSE
)
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
#> Metric: 'bhattacharyya' with unit: 'log'; comparing: 2 vectors
# Preview tidy evaluation summary
head(eval_res$metrics_summary_table)
#> Category Property Metric Value
#> 1 Distributional Properties library_size MAD 1.500000e+01
#> 2 Distributional Properties library_size KS 2.125000e-01
#> 3 Distributional Properties library_size MAE 1.506250e+01
#> 4 Distributional Properties library_size RMSE 1.571504e+01
#> 5 Distributional Properties library_size OV 8.370258e-01
#> 6 Distributional Properties library_size Bhattacharyya 1.188839e-044. Diagnostic Visualizations
scSimEval includes built-in functions for
high-resolution diagnostic plotting:
# 14-Panel Comparative Distribution QC Plot
plot_distribution_qc(
ref_data = example_scrna$ref,
sim_data = example_scrna$sim,
title = "scRNA-seq Simulation Quality Control"
)5. Conclusion
This demonstration illustrates how scSimEval
comprehensively assesses scRNA-seq simulations across distributions,
zero patterns, clustering, batch integration, marker genes, and
trajectory dynamics without requiring artificial ground-truth
labels.
sessionInfo()
#> R version 4.6.1 (2026-06-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#> LAPACK version 3.12.1
#>
#> locale:
#> [1] LC_COLLATE=English_India.utf8 LC_CTYPE=English_India.utf8
#> [3] LC_MONETARY=English_India.utf8 LC_NUMERIC=C
#> [5] LC_TIME=English_India.utf8
#>
#> time zone: Asia/Calcutta
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] scSimEval_0.99.7
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 jsonlite_2.0.0
#> [3] magrittr_2.0.5 ggbeeswarm_0.7.3
#> [5] SuppDists_1.1-9.9 modeltools_0.2-24
#> [7] farver_2.1.2 nloptr_2.2.1
#> [9] rmarkdown_2.32 provenance_4.4
#> [11] fs_2.1.0 ragg_1.5.2
#> [13] vctrs_0.7.3 memoise_2.0.1
#> [15] minqa_1.2.8 htmltools_0.5.9
#> [17] S4Arrays_1.13.0 curl_8.0.0
#> [19] BiocNeighbors_2.7.3 SparseArray_1.13.2
#> [21] overlapping_2.5 sass_0.4.10
#> [23] bslib_0.12.0 htmlwidgets_1.6.4
#> [25] desc_1.4.3 testthat_3.3.2
#> [27] httr2_1.3.0 cachem_1.1.0
#> [29] igraph_2.3.3 lifecycle_1.0.5
#> [31] pkgconfig_2.0.3 rsvd_1.0.5
#> [33] Matrix_1.7-6 R6_2.6.1
#> [35] fastmap_1.2.0 clue_0.3-68
#> [37] rbibutils_2.4.1 MatrixGenerics_1.25.0
#> [39] digest_0.6.39 S4Vectors_0.51.10
#> [41] scater_1.41.2 irlba_2.3.7
#> [43] IsoplotR_7.0 textshaping_1.0.5
#> [45] GenomicRanges_1.65.4 RSQLite_3.53.3
#> [47] beachmat_2.29.3 philentropy_0.10.0
#> [49] filelock_1.0.3 abind_1.4-8
#> [51] compiler_4.6.1 proxy_0.4-29
#> [53] bit64_4.8.6 withr_3.0.3
#> [55] CellMixS_1.29.0 S7_0.2.2
#> [57] BiocParallel_1.47.0 viridis_0.6.5
#> [59] DBI_1.3.0 MASS_7.3-66
#> [61] DelayedArray_0.39.8 bluster_1.23.1
#> [63] tools_4.6.1 vipor_0.4.7
#> [65] otel_0.2.0 beeswarm_0.4.0
#> [67] prabclus_2.3-5 nnet_7.3-21
#> [69] glue_1.8.1 nlme_3.1-171
#> [71] grid_4.6.1 waddR_1.27.0
#> [73] cluster_2.1.8.3 ade4_1.7-24
#> [75] generics_0.1.4 gtable_0.3.6
#> [77] class_7.3-24 tidyr_1.3.2
#> [79] ScaledMatrix_1.21.0 BiocSingular_1.29.1
#> [81] XVector_0.52.0 flexmix_2.3-21
#> [83] BiocGenerics_0.59.12 ggrepel_0.9.8
#> [85] RANN_2.6.3 pillar_1.11.1
#> [87] limma_3.99.0 clValid_0.7
#> [89] robustbase_0.99-7 splines_4.6.1
#> [91] dplyr_1.2.1 BiocFileCache_3.3.0
#> [93] lattice_0.22-9 bit_4.6.0
#> [95] tidyselect_1.2.1 SingleCellExperiment_1.35.2
#> [97] locfit_1.5-9.12 scuttle_1.23.2
#> [99] clusterSim_0.51-6 knitr_1.52
#> [101] gridExtra_2.3.1 reformulas_0.4.4
#> [103] IRanges_2.47.5 Seqinfo_1.3.2
#> [105] edgeR_4.99.6 SummarizedExperiment_1.43.0
#> [107] stats4_4.6.1 xfun_0.60
#> [109] Biobase_2.73.2 statmod_1.5.2
#> [111] diptest_0.77-2 randtests_1.0.2
#> [113] brio_1.1.5 arm_1.15-3
#> [115] matrixStats_1.5.0 DEoptimR_1.2-2
#> [117] yaml_2.3.12 boot_1.3-32
#> [119] kSamples_1.2-12 eva_0.2.7
#> [121] evaluate_1.0.5 codetools_0.2-20
#> [123] MLmetrics_1.1.3 kernlab_0.9-33
#> [125] tibble_3.3.1 cli_3.6.6
#> [127] RcppParallel_6.2.1 systemfonts_1.3.2
#> [129] Rdpack_2.6.6 jquerylib_0.1.4
#> [131] dichromat_2.0-1 Rcpp_1.1.2
#> [133] dbplyr_2.6.0 coda_0.19-4.1
#> [135] parallel_4.6.1 pkgdown_2.2.1
#> [137] ggplot2_4.0.3 blob_1.3.0
#> [139] mclust_6.1.3 lme4_2.0-6
#> [141] viridisLite_0.4.3 ggridges_0.5.7
#> [143] scales_1.4.0 e1071_1.7-17
#> [145] purrr_1.2.2 fpc_2.2-15
#> [147] rlang_1.3.0 cowplot_1.2.0