Testing escape. Processing and differences for Seurat and SingleCellExperiment objects

[cathalgking]

escape is a great package that I would like to use more often. I have just noticed some technical issues, detailed below.

Processing

Issue with the enrichment method used.

I successfully execute escape::runEscape() with method = "ssGSEA" and method = "UCell", for example:

# ssGSEA
song_seu <- runEscape(song_seu, 
                           method = "ssGSEA",
                           gene.sets = all_geneset, 
                           groups = 5000, 
                           min.size = 0,
                           new.assay.name = "escape.ssGSEA")

> song_seu
An object of class Seurat 
21881 features across 13322 samples within 3 assays 
Active assay: RNA (21877 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 2 other assays present: escape.ssGSEA, escape.UCell
 3 dimensional reductions calculated: pca, tsne, umap

However, I get an error when I change method = "AUC"

song_seu <- runEscape(song_seu, 
+                            method = "AUC",
+                            gene.sets = all_geneset, 
+                            groups = 5000, 
+                            min.size = 0,
+                            new.assay.name = "escape.AUC")
[1] "Using sets of 5000 cells. Running 3 times."
Error in as.data.frame.default(a) : 
  cannot coerce class ‘c("gg", "ggplot")’ to a data.frame

Visualisations

geyserEnrichment

I can plot geyserEnrichment on my SCE object but not on my Seurat object due to the gene.set argument. The error is:

Error in cnts[gene.set, ] : subscript out of bounds

Plot with:

# seu
geyserEnrichment(song_seu, 
                 assay = "escape.ssGSEA",
                 gene.set = "BMP_geneset", 
                 group.by = "malignancy")
# sce
geyserEnrichment(song_sce, 
                 assay = "escape.ssGSEA",
                 gene.set = "BMP_geneset", 
                 group.by = "malignancy")

I notice that the _ is converted to a - in the gene set name BMP_geneset.

> str(song_seu@assays$escape.ssGSEA)
Formal class 'Assay' [package "SeuratObject"] with 8 slots
  ..@ counts       : num[0 , 0 ] 
  ..@ data         : num [1:2, 1:13322] 5150 -1800 5785 -3084 5832 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "PUFA-GENESET" "BMP-geneset"
  .. .. ..$ : chr [1:13322] "AUG_PB1" "AUG_PB1.1" "AUG_PB1.2" "AUG_PB1.3" ...
  ..@ scale.data   : num[0 , 0 ] 
  ..@ assay.orig   : NULL
  ..@ var.features : logi(0) 
  ..@ meta.features:'data.frame':	2 obs. of  0 variables
  ..@ misc         : Named list()
  ..@ key          : chr "escapessgsea_"

So plotting with the - instead does work. When I assign my gene set as BMPgeneset then there is no issue.

#
geyserEnrichment(song_seu, 
                 assay = "escape.ssGSEA",
                 gene.set = "BMP-geneset", 
                 group.by = "malignancy")

FeaturePlot

This plot will work only if one method of enrichment has been calculated (on a Seurat object). It does not work if for example 2 other assays present: escape.ssGSEA, escape.UCell.
If only one enrichment method is calculated then it does plot but still has trouble locating it within the Seurat object with a returned warning for example:

FeaturePlot(song_seu, 
+             "BMPgeneset", 
+             reduction = "umap") + 
+   scale_color_gradientn(colors = colorblind_vector) + 
+   theme(plot.title = element_blank())
Error in `FeaturePlot()`:
! None of the requested features were found: BMPgeneset in slot  data
Run `rlang::last_trace()` to see where the error occurred.
Warning messages:
1: Found the following features in more than one assay, excluding the default. We will not include these in the final data frame: BMPgeneset 
2: The following requested variables were not found: BMPgeneset 

How can I plot the same feature plot for a SCE object?

The same issue seems to happen for escape::heatmapEnrichment().

> song_seu@version
[1] ‘5.0.2’

And

> sessionInfo()
R version 4.4.0 (2024-04-24)
Platform: aarch64-apple-darwin20
Running under: macOS 15.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Adelaide
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggplot2_3.5.1               UCell_2.10.1                AUCell_1.28.0               cellula_0.0.900             SingleCellExperiment_1.28.1 SummarizedExperiment_1.36.0
 [7] Biobase_2.66.0              GenomicRanges_1.58.0        GenomeInfoDb_1.42.3         IRanges_2.40.1              S4Vectors_0.44.0            BiocGenerics_0.52.0        
[13] MatrixGenerics_1.18.1       matrixStats_1.4.1           Seurat_5.1.0                SeuratObject_5.0.2          sp_2.1-4                    escape_2.2.3               

loaded via a namespace (and not attached):
  [1] GSVA_2.0.7                spatstat.sparse_3.1-0     httr_1.4.7                RColorBrewer_1.1-3        tools_4.4.0               sctransform_0.4.1         R6_2.6.1                 
  [8] ResidualMatrix_1.16.0     HDF5Array_1.34.0          lazyeval_0.2.2            uwot_0.2.2                ggdist_3.3.3              rhdf5filters_1.18.1       withr_3.0.2              
 [15] gridExtra_2.3             progressr_0.15.1          cli_3.6.5                 spatstat.explore_3.4-3    fastDummies_1.7.5         labeling_0.4.3            spatstat.data_3.1-6      
 [22] ggridges_0.5.6            pbapply_1.7-2             R.utils_2.12.3            scater_1.34.1             parallelly_1.44.0         limma_3.62.2              rstudioapi_0.17.1        
 [29] RSQLite_2.3.7             generics_0.1.4            ica_1.0-3                 spatstat.random_3.4-1     dplyr_1.1.4               distributional_0.5.0      Matrix_1.7-3             
 [36] ggbeeswarm_0.7.2          abind_1.4-8               R.methodsS3_1.8.2         lifecycle_1.0.4           yaml_2.3.10               edgeR_4.4.2               rhdf5_2.50.2             
 [43] SparseArray_1.6.2         Rtsne_0.17                grid_4.4.0                blob_1.2.4                promises_1.3.2            dqrng_0.4.1               crayon_1.5.3             
 [50] miniUI_0.1.1.1            lattice_0.22-6            beachmat_2.22.0           msigdbr_7.5.1             cowplot_1.1.3             annotate_1.84.0           KEGGREST_1.46.0          
 [57] magick_2.8.6              pillar_1.10.2             knitr_1.50                metapod_1.14.0            rjson_0.2.23              future.apply_1.11.3       codetools_0.2-20         
 [64] leiden_0.4.3.1            glue_1.8.0                spatstat.univar_3.1-3     data.table_1.17.4         vctrs_0.6.5               png_0.1-8                 spam_2.11-1              
 [71] gtable_0.3.6              cachem_1.1.0              xfun_0.52                 S4Arrays_1.6.0            mime_0.13                 rsconnect_1.4.1           survival_3.8-3           
 [78] statmod_1.5.0             bluster_1.16.0            fitdistrplus_1.2-2        ROCR_1.0-11               nlme_3.1-168              bit64_4.6.0-1             RcppAnnoy_0.0.22         
 [85] irlba_2.3.5.1             vipor_0.4.7               KernSmooth_2.23-24        colorspace_2.1-1          DBI_1.2.3                 tidyselect_1.2.1          bit_4.6.0                
 [92] compiler_4.4.0            graph_1.84.1              BiocNeighbors_2.0.1       DelayedArray_0.32.0       plotly_4.10.4             scales_1.3.0              lmtest_0.9-40            
 [99] stringr_1.5.1             SpatialExperiment_1.16.0  digest_0.6.37             goftest_1.2-3             spatstat.utils_3.1-4      rmarkdown_2.29            XVector_0.46.0           
[106] htmltools_0.5.8.1         pkgconfig_2.0.3           sparseMatrixStats_1.18.0  fastmap_1.2.0             rlang_1.1.6               htmlwidgets_1.6.4         UCSC.utils_1.2.0         
[113] shiny_1.10.0              DelayedMatrixStats_1.28.1 farver_2.1.2              zoo_1.8-12                jsonlite_2.0.0            BiocParallel_1.40.2       R.oo_1.27.1              
[120] BiocSingular_1.22.0       magrittr_2.0.3            scuttle_1.16.0            GenomeInfoDbData_1.2.13   dotCall64_1.2             patchwork_1.3.0           Rhdf5lib_1.28.0          
[127] munsell_0.5.1             Rcpp_1.0.14               babelgene_22.9            viridis_0.6.5             reticulate_1.42.0         stringi_1.8.4             zlibbioc_1.52.0          
[134] MASS_7.3-65               plyr_1.8.9                parallel_4.4.0            listenv_0.9.1             ggrepel_0.9.6             deldir_2.0-4              Biostrings_2.74.1        
[141] splines_4.4.0             tensor_1.5                locfit_1.5-9.12           igraph_2.1.4              spatstat.geom_3.4-1       RcppHNSW_0.6.0            pkgload_1.4.0            
[148] reshape2_1.4.4            ScaledMatrix_1.14.0       XML_3.99-0.17             evaluate_1.0.3            scran_1.34.0              BiocManager_1.30.25       httpuv_1.6.15            
[155] batchelor_1.22.0          RANN_2.6.2                tidyr_1.3.1               purrr_1.0.4               polyclip_1.10-7           future_1.49.0             scattermore_1.2          
[162] rsvd_1.0.5                xtable_1.8-4              RSpectra_0.16-2           later_1.3.2               viridisLite_0.4.2         ggpointdensity_0.2.0      tibble_3.2.1             
[169] memoise_2.0.1             beeswarm_0.4.0            AnnotationDbi_1.68.0      cluster_2.1.6             globals_0.18.0            GSEABase_1.68.0     

[ncborcherding]

@cathalgking

Thanks for the excellent run down and issue formatting!

Your visualization issues have identified 2 big frustrations I have with Seurat:

Feature Names in Seurat can only have "-"

Seurat prevents any feature names in assays from having the "_" character and automatically converts them into "-". As far as I know there is no way around this, but I will include some additional notes in the documentation and vignette to explain this.

# seu
geyserEnrichment(song_seu, 
                 assay = "escape.ssGSEA",
                 gene.set = "BMP_geneset", 
                 group.by = "malignancy")
# sce
geyserEnrichment(song_sce, 
                 assay = "escape.ssGSEA",
                 gene.set = "BMP_geneset", 
                 group.by = "malignancy")

Duplicate feature names

In this instance, Seurat gets confused because there exist 2 features "BMPgeneset" in 2 assays. Seurat is agnostic when it comes to features that are not in your default assay and therefore errors out. To fix this, you need to declare the assay in your FeaturePlot(...,assay = "escape.GSEA") call.

FeaturePlot(song_seu, 
             "BMPgeneset", 
             reduction = "umap") + 
  scale_color_gradientn(colors = colorblind_vector) + 
  theme(plot.title = element_blank())

How can I plot the same feature plot for a SCE object?

You can used plotReducedDim() from the scater package to plot SCE objects

AUCell

This is strange - although there is a plotting function within AUCell::AUCell_buildRankings(), it has always been set to FALSE and no where should there be a return of a ggplot object. It looks like the authors of AUCell have not been maintaining there code base via github and right now the bioconductor code is not searchable. As of now, I cannot reproduce this issue and I can't troubleshoot unfortunately.

song_seu <- runEscape(song_seu, 
                           method = "AUC",
                           gene.sets = all_geneset, 
                           groups = 5000, 
                           min.size = 0,
                           new.assay.name = "escape.AUC")

[1] "Using sets of 5000 cells. Running 3 times."
Error in as.data.frame.default(a) :
cannot coerce class ‘c("gg", "ggplot")’ to a data.frame

[cathalgking]

Hi @ncborcherding

Thanks for your reply. Thats ok re Seurat feature names and the AUC issue.

However the argument assay = does not seem to be possible with Seurat::FeaturePlot(). Do you know why that might be?

> song_seu
An object of class Seurat 
21881 features across 13322 samples within 3 assays 
Active assay: RNA (21877 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 2 other assays present: escape.ssGSEA, escape.UCell
 3 dimensional reductions calculated: pca, tsne, umap
> song_seu@assays
$RNA
Assay data with 21877 features for 13322 cells
Top 10 variable features:
 SCGB1A1, OLFM4, LTF, KRT13, HBB, MMP7, TGM4, KRT4, CST1, PIGR 

$escape.ssGSEA
Assay (v5) data with 2 features for 13322 cells
First 2 features:
 PUFA-geneset, BMP-geneset 
Layers:
 data 

$escape.UCell
Assay (v5) data with 2 features for 13322 cells
First 2 features:
 PUFA-geneset, BMP-geneset 
Layers:
 data 

> # feature plot
> FeaturePlot(song_seu, 
+             "BMP-geneset", 
+             reduction = "umap", 
+             assay = "escape.ssGSEA") + 
+   scale_color_gradientn(colors = colorblind_vector) + 
+   theme(plot.title = element_blank())
Error in FeaturePlot(song_seu, "BMP-geneset", reduction = "umap", assay = "escape.ssGSEA") : 
  unused argument (assay = "escape.ssGSEA")

The plot scater::plotReducedDim does not seem to work. Where is the equivalent assays stored for a SCE object?

> assays(song_sce)
List of length 2
names(2): counts logcounts
> song_sce
class: SingleCellExperiment 
dim: 21877 13322 
metadata(0):
assays(2): counts logcounts
rownames(21877): A1BG-AS1 A1CF ... ZFHX4-AS1 C9orf171
rowData names(0):
colnames(13322): AUG_PB1 AUG_PB1.1 ... PR5199.1017 PR5199.1019
colData names(12): orig.ident nCount_RNA ... malignancy ident
reducedDimNames(3): PCA TSNE UMAP
mainExpName: RNA
altExpNames(1): escape.ssGSEA
> # feat plot
> plotReducedDim(song_sce, 
+                "UMAP", 
+                colour_by = "BMP_geneset", 
+                assay = "escape.ssGSEA")
Error in assay(current, assay.type, withDimnames = FALSE) : 
  'assay(<SummarizedExperiment>, i="character", ...)' invalid subscript 'i'
'logcounts' not in names(assays(<SummarizedExperiment>))