Introduction to GSABenchmark
Introduction
GSABenchmark is a package designed for benchmarking scRNA-seq gene set analysis methods. It provides both traditional and novel benchmark metrics as well as visualization tools. Currently, GSABenchmark supports 17 gene set analysis methods.
Installation
To install GSABenchmark, run the following commands in an R session:
Prerequisites
In addition to GSABenchmark, you need to install scRNAseq and scater for this tutorial.
Loading data
This tutorial uses an scRNA-seq human pancreas dataset with provided cell type annotations.
After loading the required packages, download the dataset using the
BaronPancreasData function from scRNAseq. The
dataset will be stored as a SingleCellExperiment
object.
GSABenchmark requires normalized and log-transformed data and PCA
coordinates. We will use the functions logNormCounts from
scuttle (loaded with scater) and runPCA from
scater for these steps:
scObj <- logNormCounts(scObj)
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Warning in .library_size_factors(assay(x, assay.type), ...): 'librarySizeFactors' is deprecated.
#> Use 'scrapper::centerSizeFactors' instead.
#> See help("Deprecated")
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Warning in .local(x, ...): 'normalizeCounts' is deprecated.
#> Use 'scrapper::normalizeCounts' instead.
#> See help("Deprecated")
scObj <- runPCA(scObj)
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'Creating gene sets
First, we will take a look at the cell types found in this dataset:
table(colData(scObj)$label)
#>
#> acinar activated_stellate alpha beta
#> 958 284 2326 2525
#> delta ductal endothelial epsilon
#> 601 1077 252 18
#> gamma macrophage mast quiescent_stellate
#> 255 55 25 173
#> schwann t_cell
#> 13 7We will create gene sets for a cell type found in the data, using cell type markers from PanglaoDB:
alphaMarkers <- c('GCG', 'TTR', 'PCSK2', 'FXYD5', 'LDB2', 'MAFB',
'CHGA', 'SCGB2A1', 'GLS', 'FAP', 'DPP4', 'GPR119',
'PAX6', 'NEUROD1', 'LOXL4', 'PLCE1', 'GC', 'KLHL41',
'FEV', 'PTGER3', 'RFX6', 'SMARCA1', 'PGR', 'IRX1',
'UCP2', 'RGS4', 'KCNK16', 'GLP1R', 'ARX', 'POU3F4',
'RESP18', 'PYY', 'SLC38A5', 'TM4SF4', 'CRYBA2', 'SH3GL2',
'PCSK1', 'PRRG2', 'IRX2', 'ALDH1A1','PEMT', 'SMIM24',
'F10', 'SCGN', 'SLC30A8')
geneSets <- setNames(list(alphaMarkers), 'alpha')Note: The gene set names cannot contain spaces or
underscores, and must all be found in the column of assessed class
identities. For this dataset, this column is named
label:
Running the methods
The supportedMethods function lists all the methods
currently available in GSABenchmark:
supportedMethods()
#> [1] "AddModuleScore" "AUCell" "CSOA" "GSVA"
#> [5] "JASMINE" "MDT" "MLM" "ORA"
#> [9] "Pagoda2" "PLAGE" "Singscore" "SiPSiC"
#> [13] "ssGSEA" "UCell" "UDT" "VAM"
#> [17] "Zscore"For the purposes of this tutorial, we will select three of these methods:
We now run the three methods for each gene set:
Running the benchmark
The benchmark is run in a single line of code. In this vignette, we
will only perform the correctness benchmark. To achieve this, we will
skip the speed and memory benchmarks by setting
runEFBenchmark to FALSE:
smr <- runBenchmark(scObj, 'label', geneSets, gsaMethods, runEFBenchmark=FALSE)
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Computing cosine distance matrix...
#> Computing silhouette for identity class: label...
#> Computing identity class-normalized silhouette...
#> Computing PCA-based distance matrix...
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Running class determination boundary benchmark...
#> Running MCC benchmark...
#> Constructing binary predictions...
#> Running global evaluation benchmark...Note: This tutorial illustrates using
GSABenchmark with Bioconductor’s native
SingleCellExperiment class. However,
GSABenchmark is also fully compatible with
Seurat, as runGSAMethods,
runBenchmark, runMethodShuffle and
runBenchmarkShuffle can all take a Seurat
object as their first argument.
Exploring the results
runBenchmark returns a list of lists of data frames. We
can further explore the structure of this object:
names(smr)
#> [1] "boundary" "MCC" "global" "predictions"
names(smr$boundary)
#> [1] "sensitivity" "specificity" "precision" "accuracy"
#> [5] "sizeProximity" "scoreCoverage" "avg" "metricSummary"
names(smr$MCC)
#> [1] "boundaryMCC" "directMCC"
names(smr$global)
#> [1] "AUROC" "PRAUC" "labRankAlignment" "silRankAlignment"
#> [5] "centrality" "labJaccard" "labCosine" "avg"
#> [9] "metricSummary"
names(smr$predictions)
#> [1] "alpha"We can access the main summaries of the results as follows:
smr$boundary$metricSummary
#> sensitivity specificity precision accuracy sizeProximity scoreCoverage
#> CSOA 0.9075666 0.9670030 0.9110919 0.9508694 0.9985584 0.8955639
#> PLAGE 0.9350817 0.9742111 0.9310788 0.9635897 0.9983982 0.5492890
#> Zscore 0.8383491 0.9593144 0.8847550 0.9264792 0.9804581 0.5068982
#> avg
#> CSOA 0.9384422
#> PLAGE 0.8919414
#> Zscore 0.8493757
smr$MCC
#> $boundaryMCC
#> alpha avg
#> PLAGE 0.9080720 0.9080720
#> CSOA 0.8756357 0.8756357
#> Zscore 0.8115518 0.8115518
#>
#> $directMCC
#> alpha avg
#> PLAGE 0.9081798 0.9081798
#> CSOA 0.8777333 0.8777333
#> Zscore 0.8128873 0.8128873
smr$global$metricSummary
#> AUROC PRAUC labRankAlignment silRankAlignment centrality
#> PLAGE 0.9921870 0.9789604 0.9934134 0.9804435 0.8691640
#> CSOA 0.9866448 0.9682932 0.9874876 0.9729358 0.9336223
#> Zscore 0.9718930 0.9356594 0.9761797 0.9640587 0.8642028
#> labJaccard labCosine avg
#> PLAGE 0.8745476 0.9330781 0.9459706
#> CSOA 0.8337283 0.9093276 0.9417199
#> Zscore 0.7558140 0.8612395 0.9041496To visualize the results, we will call the
allBenchmarkPlots function. This function plots each
summary result data frame:
Below, we display one of the resulting plots:
The output of runBenchmark can also be used to create
other plots.
We will visualize the aggregate ranks of the methods using the
aggregateRankPlot function:
Next, we will use the ratioRank function to visualize
the top results among all method-metric-gene set combinations as ranked
by the ratio of the maximum score and the mean:
To illustrate the similarity of the results obtained by different methods, we will create MDS plots for each gene set, as well as an aggregate MDS plot. We showcase the aggregate MDS plot below:
Another way to look at the similarity between the results of different methods is through correlation plots. Below, we display the aggregate correlation plot:
We can also display the Jaccard tile plots of the binary predictions made by each method. Here, we show the aggregate Jaccard tile plot:
Benchmarking a single method at different choices of gene loss and noise
GSABenchmark also offers the option of benchmarking a single method at different choices of:
- Gene loss (the fraction of genes removed from each gene signature)
- Noise (the fraction of retained genes in the final adjusted gene signatures after random genes have been added).
Each condition determined by a gene loss-noise combination can be run
in replicates (that is, different gene signatures will be constructed
with the same gene loss and noise parameters for each replicate). The
seeds argument controls the random seeds used in the
generation of gene sets. Its length determines the number of replicates.
To save execution time, we will use only one replicate:
scObj <- runMethodShuffle(scObj, 'label', geneSets, 'CSOA',
loss=1/3,
noise=0.5,
seeds=20,
averageReplicates=FALSE)
#> Shuffling genes: gene loss = 33.3%, noise = 50%, replicate = 1.
#> Warning in asMethod(object): sparse->dense coercion: allocating vector of size
#> 1.3 GiB
#> Running CSOA...
#> Computing percentile sets...
#> Warning in percentileSets(geneSetExp, percentile): 1 gene(s) had no top cells
#> at the indicated percentile. These are now excluded from the gene signature.
#> Generating cell set overlaps...
#> Processing overlaps...
#> 44 overlaps will be used in the calculation of CSOA scores.
#> Normalizing expression matrix by rows...
#> Computing per-cell scores for gene pairs...
#> Computing per-cell gene signature scores...Subsequently, the same benchmark as before can be run with
runBenchmarkShuffle, a thin wrapper around
runBenchmark:
smr <- runBenchmarkShuffle(scObj, 'label', geneSets, 'CSOA',
runEFBenchmark=FALSE)
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Computing cosine distance matrix...
#> Computing silhouette for identity class: label...
#> Computing identity class-normalized silhouette...
#> Computing PCA-based distance matrix...
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Found more than one class "package_version" in cache; using the first, from namespace 'SeuratObject'
#> Also defined by 'alabaster.base'
#> Running class determination boundary benchmark...
#> Running MCC benchmark...
#> Constructing binary predictions...
#> Running global evaluation benchmark...The output of runBenchmarkShuffle can be inspected as
before, e.g.:
Session information
sessionInfo()
#> R version 4.5.3 (2026-03-11)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 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.26.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] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] scRNAseq_2.25.0 scater_1.39.4
#> [3] ggplot2_4.0.2 scuttle_1.21.6
#> [5] SingleCellExperiment_1.33.2 SummarizedExperiment_1.41.1
#> [7] Biobase_2.71.0 GenomicRanges_1.63.2
#> [9] Seqinfo_1.1.0 IRanges_2.45.0
#> [11] S4Vectors_0.49.2 BiocGenerics_0.57.1
#> [13] generics_0.1.4 MatrixGenerics_1.23.0
#> [15] matrixStats_1.5.0 GSABenchmark_0.99.8
#> [17] BiocStyle_2.39.0
#>
#> loaded via a namespace (and not attached):
#> [1] igraph_2.3.0 graph_1.89.1
#> [3] ica_1.0-3 plotly_4.12.0
#> [5] rematch2_2.1.2 tidyselect_1.2.1
#> [7] bit_4.6.0 lattice_0.22-9
#> [9] rjson_0.2.23 brew_1.0-10
#> [11] blob_1.3.0 stringr_1.6.0
#> [13] lgr_0.5.2 S4Arrays_1.11.1
#> [15] RMTstat_0.3.1 parallel_4.5.3
#> [17] drat_0.2.5 png_0.1-9
#> [19] cli_3.6.6 qs2_0.1.7
#> [21] MLmetrics_1.1.3 ProtGenerics_1.43.0
#> [23] hammers_0.99.9 goftest_1.2-3
#> [25] textshape_1.7.5 BiocIO_1.21.0
#> [27] kernlab_0.9-33 purrr_1.2.2
#> [29] BiocNeighbors_2.5.4 uwot_0.2.4
#> [31] curl_7.0.0 mime_0.13
#> [33] evaluate_1.0.5 stringi_1.8.7
#> [35] XML_3.99-0.23 lubridate_1.9.5
#> [37] httpuv_1.6.17 AnnotationDbi_1.73.1
#> [39] paletteer_1.7.0 magrittr_2.0.5
#> [41] rappdirs_0.3.4 splines_4.5.3
#> [43] maketools_1.3.2 float_0.3-3
#> [45] ggraph_2.2.2 dplyr_1.2.1
#> [47] kerntools_1.2.1 sctransform_0.4.3
#> [49] ggbeeswarm_0.7.3 DBI_1.3.0
#> [51] HDF5Array_1.39.1 jquerylib_0.1.4
#> [53] withr_3.0.2 lmtest_0.9-40
#> [55] ggnewscale_0.5.2 rsparse_0.5.3
#> [57] GSEABase_1.73.0 tidygraph_1.3.1
#> [59] ggeasy_0.1.6 rtracklayer_1.71.3
#> [61] BiocManager_1.30.27 liver_1.28
#> [63] htmlwidgets_1.6.4 fs_2.1.0
#> [65] triebeard_0.4.1 ggrepel_0.9.8
#> [67] labeling_0.4.3 SparseArray_1.11.13
#> [69] cellranger_1.1.0 h5mread_1.3.3
#> [71] annotate_1.89.0 reticulate_1.46.0
#> [73] zoo_1.8-15 XVector_0.51.0
#> [75] knitr_1.51 SiPSiC_1.11.0
#> [77] UCSC.utils_1.7.1 dendsort_0.3.4
#> [79] RhpcBLASctl_0.23-42 timechange_0.4.0
#> [81] patchwork_1.3.2 grid_4.5.3
#> [83] data.table_1.18.2.1 rhdf5_2.55.16
#> [85] R.oo_1.27.1 RSpectra_0.16-2
#> [87] irlba_2.3.7 alabaster.schemas_1.11.0
#> [89] fastDummies_1.7.5 lazyeval_0.2.3
#> [91] yaml_2.3.12 survival_3.8-6
#> [93] SpatialExperiment_1.21.0 scattermore_1.2
#> [95] BiocVersion_3.23.1 crayon_1.5.3
#> [97] RcppAnnoy_0.0.23 RColorBrewer_1.1-3
#> [99] tidyr_1.3.2 progressr_0.19.0
#> [101] tweenr_2.0.3 later_1.4.8
#> [103] ggridges_0.5.7 codetools_0.2-20
#> [105] mltools_0.3.5 Seurat_5.4.0
#> [107] KEGGREST_1.51.1 sccore_1.0.7
#> [109] Rtsne_0.17 limma_3.67.3
#> [111] urltools_1.7.3.1 Rsamtools_2.27.2
#> [113] filelock_1.0.3 pkgconfig_2.0.3
#> [115] spatstat.univar_3.1-7 pagoda2_1.0.15
#> [117] text2vec_0.6.6 GenomicAlignments_1.47.0
#> [119] spatstat.sparse_3.1-0 alabaster.base_1.11.4
#> [121] viridisLite_0.4.3 xtable_1.8-8
#> [123] decoupleR_2.17.0 N2R_1.0.5
#> [125] plyr_1.8.9 httr_1.4.8
#> [127] tools_4.5.3 globals_0.19.1
#> [129] sys_3.4.3 SeuratObject_5.4.0
#> [131] beeswarm_0.4.0 nlme_3.1-169
#> [133] dbplyr_2.5.2 ExperimentHub_3.1.0
#> [135] digest_0.6.39 bookdown_0.46
#> [137] Matrix_1.7-5 farver_2.1.2
#> [139] tzdb_0.5.0 AnnotationFilter_1.35.0
#> [141] reshape2_1.4.5 SnowballC_0.7.1
#> [143] viridis_0.6.5 glue_1.8.1
#> [145] cachem_1.1.0 BiocFileCache_3.1.0
#> [147] polyclip_1.10-7 buildtools_1.0.0
#> [149] Biostrings_2.79.5 ggalluvial_0.12.6
#> [151] mlapi_0.1.1 parallelly_1.47.0
#> [153] statmod_1.5.1 escape_2.7.3
#> [155] RcppHNSW_0.6.0 ScaledMatrix_1.19.0
#> [157] pbapply_1.7-4 httr2_1.2.2
#> [159] spam_2.11-3 VAM_1.1.0
#> [161] henna_0.7.5 graphlayouts_1.2.3
#> [163] readxl_1.4.5 alabaster.se_1.11.0
#> [165] lsa_0.73.4 gridExtra_2.3
#> [167] shiny_1.13.0 GSVA_2.5.40
#> [169] R.utils_2.13.0 rhdf5filters_1.23.3
#> [171] RCurl_1.98-1.18 memoise_2.0.1
#> [173] alabaster.sce_1.11.0 rmarkdown_2.31
#> [175] fabR_2.1.1 scales_1.4.0
#> [177] R.methodsS3_1.8.2 future_1.70.0
#> [179] singscore_1.31.0 gypsum_1.7.0
#> [181] RANN_2.6.2 stringfish_0.19.0
#> [183] spatstat.data_3.1-9 cluster_2.1.8.2
#> [185] janitor_2.2.1 spatstat.utils_3.2-2
#> [187] hms_1.1.4 fitdistrplus_1.2-6
#> [189] cowplot_1.2.0 rlang_1.2.0
#> [191] CSOA_1.1.7 GenomeInfoDb_1.47.2
#> [193] Rook_1.2 DelayedMatrixStats_1.33.0
#> [195] sparseMatrixStats_1.23.0 jaccard_0.1.2
#> [197] dotCall64_1.2 ggforce_0.5.0
#> [199] mgcv_1.9-4 xfun_0.57
#> [201] alabaster.matrix_1.11.0 abind_1.4-8
#> [203] tibble_3.3.1 Rhdf5lib_1.99.6
#> [205] readr_2.2.0 bitops_1.0-9
#> [207] promises_1.5.0 RSQLite_2.4.6
#> [209] qvalue_2.43.0 DelayedArray_0.37.1
#> [211] compiler_4.5.3 forcats_1.0.1
#> [213] alabaster.ranges_1.11.0 distributional_0.7.0
#> [215] writexl_1.5.4 beachmat_2.27.5
#> [217] memuse_4.2-3 listenv_0.10.1
#> [219] Rcpp_1.1.1-1 edgeR_4.9.9
#> [221] AnnotationHub_4.1.0 BiocSingular_1.27.1
#> [223] tensor_1.5.1 usethis_3.2.1
#> [225] MASS_7.3-65 BiocParallel_1.45.0
#> [227] spatstat.random_3.4-5 R6_2.6.1
#> [229] fastmap_1.2.0 vipor_0.4.7
#> [231] ensembldb_2.35.0 ROCR_1.0-12
#> [233] rsvd_1.0.5 ggdist_3.3.3
#> [235] gtable_0.3.6 scLang_0.99.3
#> [237] KernSmooth_2.23-26 miniUI_0.1.2
#> [239] deldir_2.0-4 htmltools_0.5.9
#> [241] RcppParallel_5.1.11-2 bit64_4.8.0
#> [243] spatstat.explore_3.8-0 lifecycle_1.0.5
#> [245] S7_0.2.1-1 restfulr_0.0.16
#> [247] sass_0.4.10 vctrs_0.7.3
#> [249] spatstat.geom_3.7-3 snakecase_0.11.1
#> [251] haven_2.5.5 sp_2.2-1
#> [253] future.apply_1.20.2 bslib_0.10.0
#> [255] pillar_1.11.1 GenomicFeatures_1.63.2
#> [257] prismatic_1.1.2 magick_2.9.1
#> [259] locfit_1.5-9.12 otel_0.2.0
#> [261] jsonlite_2.0.0 abdiv_0.2.0
#> [263] cigarillo_1.1.0