Introduction to GSABenchmark
March 4, 2025
Package
GSABenchmark 0.99.6
Contents
1 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.
2 Installation
To install GSABenchmark, run the following commands in an R session:
if (!require("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("GSABenchmark")3 Prerequisites
In addition to GSABenchmark, you need to install scRNAseq and scater for this tutorial.
4 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.
library(GSABenchmark)
library(scater)
library(scRNAseq)
scObj <- BaronPancreasData('human')
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'5 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:
setdiff(names(geneSets), colData(scObj)[['label']])
#> character(0)6 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:
gsaMethods <- c('CSOA', 'PLAGE', 'Zscore')We now run the three methods for each gene set:
scObj <- runGSAMethods(scObj, 'label', geneSets, gsaMethods)
#> Running CSOA...
#> Running PLAGE...
#> Running Zscore...7 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.
8 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.9804436 0.8691637
#> CSOA 0.9866448 0.9682932 0.9874876 0.9729359 0.9336222
#> Zscore 0.9718930 0.9356594 0.9761797 0.9640588 0.8642024
#> labJaccard labCosine avg
#> PLAGE 0.8745476 0.9330781 0.9459705
#> 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:
plots <- allBenchmarkPlots(smr)
length(plots)
#> [1] 19Below, we display one of the resulting plots:
plots[[3]]
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:
aggregateRankPlot(smr)
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:
ratioPlot(smr, labelSize=2.5)
#> Cannot return 25 top ratios. Will retun only 15.
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:
plots <- mdsPlots(scObj, smr)
plots[[length(geneSets) + 1]]
Another way to look at the similarity between the results of different methods is through correlation plots. Below, we display the aggregate correlation plot:
plots <- corrPlots(scObj, smr)
plots[[length(geneSets) + 1]]
We can also display the Jaccard tile plots of the binary predictions made by each method. Here, we show the aggregate Jaccard tile plot:
plots <- predJaccardPlots(smr$predictions)
plots[[length(geneSets) + 1]]
9 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.:
smr$boundary$metricSummary
#> sensitivity specificity precision accuracy sizeProximity
#> CSOA 0.9075666 0.9670030 0.9110919 0.9508694 0.9985584
#> CSOA_33.3_50_1 0.9092863 0.9575525 0.8886555 0.9444509 0.9913503
#> scoreCoverage avg
#> CSOA 0.8955639 0.9384422
#> CSOA_33.3_50_1 0.8401581 0.9219089Session information
sessionInfo()
#> R version 4.6.0 alpha (2026-04-05 r89794)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.23-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [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: America/New_York
#> 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.5
#> [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.1 BiocGenerics_0.57.0
#> [13] generics_0.1.4 MatrixGenerics_1.23.0
#> [15] matrixStats_1.5.0 GSABenchmark_0.99.6
#> [17] BiocStyle_2.39.0
#>
#> loaded via a namespace (and not attached):
#> [1] SpatialExperiment_1.21.0 R.methodsS3_1.8.2 dichromat_2.0-0.1
#> [4] GSEABase_1.73.0 goftest_1.2-3 Biostrings_2.79.5
#> [7] HDF5Array_1.39.1 vctrs_0.7.2 spatstat.random_3.4-5
#> [10] digest_0.6.39 png_0.1-9 gypsum_1.7.0
#> [13] ggrepel_0.9.8 deldir_2.0-4 parallelly_1.46.1
#> [16] alabaster.sce_1.11.0 magick_2.9.1 MASS_7.3-65
#> [19] reshape2_1.4.5 httpuv_1.6.17 qvalue_2.43.0
#> [22] withr_3.0.2 escape_2.7.3 xfun_0.57
#> [25] survival_3.8-6 memoise_2.0.1 liver_1.28
#> [28] ggbeeswarm_0.7.3 janitor_2.2.1 zoo_1.8-15
#> [31] pbapply_1.7-4 ggdist_3.3.3 R.oo_1.27.1
#> [34] rematch2_2.1.2 KEGGREST_1.51.1 promises_1.5.0
#> [37] otel_0.2.0 httr_1.4.8 restfulr_0.0.16
#> [40] globals_0.19.1 fitdistrplus_1.2-6 rhdf5filters_1.23.3
#> [43] mltools_0.3.5 stringfish_0.18.0 rhdf5_2.55.16
#> [46] UCSC.utils_1.7.1 miniUI_0.1.2 ggalluvial_0.12.6
#> [49] curl_7.0.0 ScaledMatrix_1.19.0 ggraph_2.2.2
#> [52] h5mread_1.3.3 polyclip_1.10-7 ExperimentHub_3.1.0
#> [55] lgr_0.5.2 SparseArray_1.11.13 xtable_1.8-8
#> [58] stringr_1.6.0 evaluate_1.0.5 S4Arrays_1.11.1
#> [61] BiocFileCache_3.1.0 hms_1.1.4 bookdown_0.46
#> [64] irlba_2.3.7 filelock_1.0.3 ROCR_1.0-12
#> [67] qs2_0.1.7 reticulate_1.45.0 float_0.3-3
#> [70] readxl_1.4.5 spatstat.data_3.1-9 magrittr_2.0.5
#> [73] lmtest_0.9-40 snakecase_0.11.1 readr_2.2.0
#> [76] later_1.4.8 viridis_0.6.5 lattice_0.22-9
#> [79] VAM_1.1.0 spatstat.geom_3.7-3 future.apply_1.20.2
#> [82] triebeard_0.4.1 scattermore_1.2 XML_3.99-0.23
#> [85] cowplot_1.2.0 RcppAnnoy_0.0.23 pillar_1.11.1
#> [88] nlme_3.1-169 singscore_1.31.0 compiler_4.6.0
#> [91] beachmat_2.27.4 RSpectra_0.16-2 stringi_1.8.7
#> [94] tensor_1.5.1 lubridate_1.9.5 GenomicAlignments_1.47.0
#> [97] plyr_1.8.9 drat_0.2.5 BiocIO_1.21.0
#> [100] crayon_1.5.3 abind_1.4-8 rsparse_0.5.3
#> [103] locfit_1.5-9.12 haven_2.5.5 sp_2.2-1
#> [106] graphlayouts_1.2.3 bit_4.6.0 dplyr_1.2.1
#> [109] codetools_0.2-20 BiocSingular_1.27.1 bslib_0.10.0
#> [112] alabaster.ranges_1.11.0 textshape_1.7.5 paletteer_1.7.0
#> [115] plotly_4.12.0 mime_0.13 splines_4.6.0
#> [118] Rcpp_1.1.1 fastDummies_1.7.5 dbplyr_2.5.2
#> [121] henna_0.7.5 sparseMatrixStats_1.23.0 prismatic_1.1.2
#> [124] cellranger_1.1.0 pagoda2_1.0.15 brew_1.0-10
#> [127] N2R_1.0.5 knitr_1.51 blob_1.3.0
#> [130] BiocVersion_3.23.1 AnnotationFilter_1.35.0 fs_2.0.1
#> [133] listenv_0.10.1 SiPSiC_1.11.0 GSVA_2.5.24
#> [136] tibble_3.3.1 Matrix_1.7-5 statmod_1.5.1
#> [139] tzdb_0.5.0 tweenr_2.0.3 pkgconfig_2.0.3
#> [142] tools_4.6.0 Rook_1.2 cachem_1.1.0
#> [145] cigarillo_1.1.0 RhpcBLASctl_0.23-42 RSQLite_2.4.6
#> [148] viridisLite_0.4.3 DBI_1.3.0 fastmap_1.2.0
#> [151] rmarkdown_2.31 scales_1.4.0 grid_4.6.0
#> [154] usethis_3.2.1 ica_1.0-3 Seurat_5.4.0
#> [157] Rsamtools_2.27.1 AnnotationHub_4.1.0 abdiv_0.2.0
#> [160] sass_0.4.10 patchwork_1.3.2 BiocManager_1.30.27
#> [163] dotCall64_1.2 graph_1.89.1 alabaster.schemas_1.11.0
#> [166] RANN_2.6.2 farver_2.1.2 mgcv_1.9-4
#> [169] tidygraph_1.3.1 yaml_2.3.12 rtracklayer_1.71.3
#> [172] cli_3.6.5 purrr_1.2.1 writexl_1.5.4
#> [175] lifecycle_1.0.5 uwot_0.2.4 kernlab_0.9-33
#> [178] BiocParallel_1.45.0 annotate_1.89.0 timechange_0.4.0
#> [181] gtable_0.3.6 rjson_0.2.23 ggridges_0.5.7
#> [184] progressr_0.19.0 limma_3.67.0 parallel_4.6.0
#> [187] SnowballC_0.7.1 edgeR_4.9.5 jsonlite_2.0.0
#> [190] bitops_1.0-9 RcppHNSW_0.6.0 CSOA_1.1.6
#> [193] bit64_4.6.0-1 hammers_0.99.9 Rtsne_0.17
#> [196] alabaster.matrix_1.11.0 BiocNeighbors_2.5.4 spatstat.utils_3.2-2
#> [199] SeuratObject_5.3.0 urltools_1.7.3.1 RcppParallel_5.1.11-2
#> [202] alabaster.se_1.11.0 jquerylib_0.1.4 fabR_2.1.1
#> [205] mlapi_0.1.1 spatstat.univar_3.1-7 jaccard_0.1.2
#> [208] distributional_0.7.0 R.utils_2.13.0 alabaster.base_1.11.4
#> [211] lazyeval_0.2.3 shiny_1.13.0 text2vec_0.6.6
#> [214] htmltools_0.5.9 sctransform_0.4.3 rappdirs_0.3.4
#> [217] tinytex_0.59 ensembldb_2.35.0 glue_1.8.0
#> [220] spam_2.11-3 httr2_1.2.2 XVector_0.51.0
#> [223] RCurl_1.98-1.18 decoupleR_2.17.0 RMTstat_0.3.1
#> [226] kerntools_1.2.1 gridExtra_2.3 sccore_1.0.7
#> [229] igraph_2.2.3 R6_2.6.1 tidyr_1.3.2
#> [232] labeling_0.4.3 forcats_1.0.1 GenomicFeatures_1.63.2
#> [235] scLang_0.99.3 cluster_2.1.8.2 Rhdf5lib_1.33.6
#> [238] memuse_4.2-3 GenomeInfoDb_1.47.2 ProtGenerics_1.43.0
#> [241] vipor_0.4.7 DelayedArray_0.37.1 tidyselect_1.2.1
#> [244] ggforce_0.5.0 dendsort_0.3.4 AnnotationDbi_1.73.1
#> [247] future_1.70.0 rsvd_1.0.5 KernSmooth_2.23-26
#> [250] S7_0.2.1 data.table_1.18.2.1 ggeasy_0.1.6
#> [253] htmlwidgets_1.6.4 RColorBrewer_1.1-3 rlang_1.2.0
#> [256] spatstat.sparse_3.1-0 lsa_0.73.4 spatstat.explore_3.8-0
#> [259] ggnewscale_0.5.2 MLmetrics_1.1.3 beeswarm_0.4.0