Single-cell assay for transposase-accessible chromatin sequencing (scATAC-seq) measures chromatin accessibility at single-cell resolution. Unlike scRNA-seq expression counts, scATAC-seq data have unique properties: * High Sparsity: Over 90% to 98% of values in peak-by-cell matrices are zeroes because each locus has at most 2 copies per cell. * Near-Binary Nature: Values primarily represent open (accessible) vs. closed chromatin states. * Peak Co-Accessibility: Genomic loci regulated together exhibit correlated accessibility across cells.
This tutorial shows how to evaluate synthetic scATAC-seq data using
scSimEval.
scSimEval includes an empirical reference and simulated
scATAC-seq peak accessibility dataset (example_scatac):
library(scSimEval)
# Load packaged scATAC-seq example data
data(example_scatac)
cat("Reference scATAC-seq dimensions:", dim(example_scatac$ref), "\n")
#> Reference scATAC-seq dimensions: 60 80
cat("Simulated scATAC-seq dimensions:", dim(example_scatac$sim), "\n")
#> Simulated scATAC-seq dimensions: 60 80
cat("Cell types:", levels(example_scatac$cell_types), "\n")
#> Cell types: TypeA TypeB
cat("Batches:", levels(example_scatac$batch_info), "\n")
#> Batches: Batch1 Batch2The matrices contain integer insertion counts across genomic peaks (\(60\text{ peaks} \times 80\text{ cells}\)).
We evaluate whether simulated peak insertion frequencies and cell
library sizes match real distributions using
evaluate_simulation_accuracy():
atac_res <- evaluate_simulation_accuracy(
ref_data = example_scatac$ref,
sim_data = example_scatac$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 peak accessibility and library size statistical distances
head(atac_res$metrics_summary_table, 8)
#> Category Property Metric Value
#> 1 Distributional Properties library_size MAD 2.0000000000
#> 2 Distributional Properties library_size KS 0.1750000000
#> 3 Distributional Properties library_size MAE 2.3375000000
#> 4 Distributional Properties library_size RMSE 2.4874685928
#> 5 Distributional Properties library_size OV 0.8548143838
#> 6 Distributional Properties library_size Bhattacharyya 0.0001055641
#> 7 Distributional Properties library_size Wasserstein 2.3375000000
#> 8 Distributional Properties library_size ECDF_DiffArea 0.0863425926In scATAC-seq, identifying cell types requires clustering cells based
on chromatin accessibility profiles. We test whether the simulator
preserves genuine biological separation using
evaluate_clustering_metrics():
atac_clust <- evaluate_clustering_metrics(
data = example_scatac$sim,
cell_types = example_scatac$cell_types,
ref_data = example_scatac$ref
)
cat("Average Silhouette Width (Simulated):", round(atac_clust$silhouette, 4), "\n")
#> Average Silhouette Width (Simulated): 0.0558
cat("Davies-Bouldin Index:", round(atac_clust$davies_bouldin, 4), "\n")
#> Davies-Bouldin Index: 3.6911
cat("Adjusted Rand Index (ARI):", round(atac_clust$ARI, 4), "\n")
#> Adjusted Rand Index (ARI): 0.7627
cat("Normalized Mutual Information (NMI):", round(atac_clust$NMI, 4), "\n")
#> Normalized Mutual Information (NMI): 0.721When simulating multi-sample experiments, batch variations can
confound biological signals. We evaluate batch mixing metrics using
evaluate_batch_metrics():
atac_batch <- evaluate_batch_metrics(
data = example_scatac$sim,
batch_info = example_scatac$batch_info,
cell_types = example_scatac$cell_types,
verbose = FALSE
)
cat("Batch Shannon Entropy:", round(atac_batch$shannon_entropy, 4), "\n")
#> Batch Shannon Entropy: 0.9671
cat("Principal Component Regression (PCR R2):", round(atac_batch$pcr_r2, 4), "\n")
#> Principal Component Regression (PCR R2): 0.0125
cat("Cross-Batch Transfer Accuracy:", round(atac_batch$cross_batch_accuracy, 4), "\n")
#> Cross-Batch Transfer Accuracy: 0.825In real biological cells, peaks that are close to each other or regulated by the same transcription factor complexes show correlated accessibility:
# Evaluate peak co-accessibility fidelity
coacc_res <- calc_peak_coaccessibility_fidelity(
ref_atac = example_scatac$ref,
sim_atac = example_scatac$sim
)
cat("Co-Accessibility RV Coefficient:", round(coacc_res$rv_coefficient, 4), "\n")
#> Co-Accessibility RV Coefficient: 0.6014
cat("Peak Co-Accessibility Pearson Correlation:", round(coacc_res$coaccessibility_pearson, 4), "\n")
#> Peak Co-Accessibility Pearson Correlation: 0.1128
cat("Peak Co-Accessibility Spearman Correlation:", round(coacc_res$coaccessibility_spearman, 4), "\n")
#> Peak Co-Accessibility Spearman Correlation: 0.0902
cat("Frobenius Matrix Distance:", round(coacc_res$frobenius_distance, 4), "\n")
#> Frobenius Matrix Distance: 0.0053Execute the comprehensive evaluation workflow in a single coordinated command:
master_atac <- evaluate_simulation_accuracy(
ref_data = example_scatac$ref,
sim_data = example_scatac$sim,
elapsed_time = 18.2,
memory_mb = 512.4,
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 tidy summary table
head(master_atac$metrics_summary_table)
#> Category Property Metric Value
#> 1 Distributional Properties library_size MAD 2.0000000000
#> 2 Distributional Properties library_size KS 0.1750000000
#> 3 Distributional Properties library_size MAE 2.3375000000
#> 4 Distributional Properties library_size RMSE 2.4874685928
#> 5 Distributional Properties library_size OV 0.8548143838
#> 6 Distributional Properties library_size Bhattacharyya 0.0001055641scSimEval effectively benchmarks scATAC-seq simulations,
validating peak insertion distributions, extreme sparsity, cell type
clustering, batch mixing, and chromatin peak co-accessibility without
requiring synthetic ground-truth annotations.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] scSimEval_0.99.7 rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 sys_3.4.3
#> [3] jsonlite_2.0.0 magrittr_2.0.5
#> [5] ggbeeswarm_0.7.3 SuppDists_1.1-9.9
#> [7] modeltools_0.2-25 farver_2.1.2
#> [9] nloptr_2.2.1 provenance_4.4
#> [11] vctrs_0.7.3 memoise_2.0.1
#> [13] minqa_1.2.8 htmltools_0.5.9
#> [15] S4Arrays_1.13.2 curl_8.0.0
#> [17] BiocNeighbors_2.7.3 SparseArray_1.13.4
#> [19] overlapping_2.5 sass_0.4.10
#> [21] bslib_0.12.0 testthat_3.3.2
#> [23] httr2_1.3.0 cachem_1.1.0
#> [25] buildtools_1.0.0 igraph_2.3.4
#> [27] lifecycle_1.0.5 pkgconfig_2.0.3
#> [29] rsvd_1.0.5 Matrix_1.7-6
#> [31] R6_2.6.1 fastmap_1.2.0
#> [33] rbibutils_2.4.1 MatrixGenerics_1.25.0
#> [35] clue_0.3-68 digest_0.6.39
#> [37] S4Vectors_0.51.10 scater_1.41.2
#> [39] irlba_2.4.1 IsoplotR_7.0
#> [41] GenomicRanges_1.65.4 RSQLite_3.53.3
#> [43] beachmat_2.29.3 philentropy_0.10.0
#> [45] filelock_1.0.3 labeling_0.4.3
#> [47] abind_1.4-8 compiler_4.6.1
#> [49] proxy_0.4-29 bit64_4.8.6
#> [51] withr_3.0.3 CellMixS_1.29.0
#> [53] S7_0.2.2 BiocParallel_1.47.0
#> [55] viridis_0.6.5 DBI_1.3.0
#> [57] MASS_7.3-66 DelayedArray_0.39.8
#> [59] bluster_1.23.1 tools_4.6.1
#> [61] vipor_0.4.7 otel_0.2.0
#> [63] beeswarm_0.4.0 prabclus_2.3-5
#> [65] nnet_7.3-21 glue_1.8.1
#> [67] nlme_3.1-171 grid_4.6.1
#> [69] waddR_1.27.0 cluster_2.1.8.3
#> [71] ade4_1.7-24 generics_0.1.4
#> [73] gtable_0.3.6 class_7.3-24
#> [75] tidyr_1.3.2 ScaledMatrix_1.21.0
#> [77] BiocSingular_1.29.2 XVector_0.53.0
#> [79] flexmix_2.3-21 BiocGenerics_0.59.12
#> [81] ggrepel_0.9.8 RANN_2.6.3
#> [83] pillar_1.11.1 limma_3.99.1
#> [85] clValid_0.7 robustbase_0.99-7
#> [87] splines_4.6.1 dplyr_1.2.1
#> [89] BiocFileCache_3.3.0 lattice_0.23-1
#> [91] bit_4.6.0 tidyselect_1.2.1
#> [93] SingleCellExperiment_1.35.2 locfit_1.5-9.12
#> [95] scuttle_1.23.2 maketools_1.3.2
#> [97] clusterSim_0.51-6 knitr_1.52
#> [99] gridExtra_2.3.1 reformulas_0.4.4
#> [101] IRanges_2.47.5 Seqinfo_1.3.2
#> [103] edgeR_4.99.6 SummarizedExperiment_1.43.0
#> [105] stats4_4.6.1 xfun_0.61
#> [107] Biobase_2.73.2 statmod_1.5.2
#> [109] diptest_0.77-2 randtests_1.0.2
#> [111] brio_1.1.5 arm_1.15-3
#> [113] matrixStats_1.5.0 DEoptimR_1.2-2
#> [115] yaml_2.3.12 boot_1.3-32
#> [117] eva_0.2.7 evaluate_1.0.5
#> [119] codetools_0.2-20 MLmetrics_1.1.3
#> [121] kSamples_1.2-12 kernlab_0.9-33
#> [123] tibble_3.3.1 cli_3.6.6
#> [125] RcppParallel_6.2.1 Rdpack_2.6.6
#> [127] jquerylib_0.1.4 Rcpp_1.1.2
#> [129] dbplyr_2.6.0 coda_0.19-4.1
#> [131] parallel_4.6.1 ggplot2_4.0.3
#> [133] blob_1.3.0 mclust_6.1.3
#> [135] lme4_2.0-6 viridisLite_0.4.3
#> [137] ggridges_0.5.7 scales_1.4.0
#> [139] e1071_1.7-17 purrr_1.2.2
#> [141] fpc_2.2-15 rlang_1.3.0
#> [143] cowplot_1.2.0