AggregateExpression() does not work with runEscape() output

[toobiwankenobi]

Hi there,

Thanks a lot for the amazing package!

I have a question regarding the interoperability of escape and Seurat. Somehow, I can't use the AggregateExpression function on a Seurat object created by runEscape. Ideally, I want to aggregate the enrichment scores using this function for downstream analysis.

GS <- list(Bcells = c("MS4A1", "CD79B", "CD79A", "IGH1", "IGH2"),
           Tcells = c("CD3E", "CD3D", "CD3G", "CD7","CD8A"))
pbmc_small <- SeuratObject::pbmc_small

pbmc_small <- runEscape(pbmc_small, 
                        gene.sets = GS, 
                        min.size = NULL)

pbmc_aggr <- AggregateExpression(pbmc_small, 
                           assays = "escape",
                           return.seurat = T, 
                           group.by = "groups")

Error in data.return[[1]] : subscript out of bounds

Am I doing something wrong? Thanks for your feedback.

Best wishes,
Tobi

[ncborcherding]

Hey @toobiwankenobi

Thanks for reaching out and the reproducible example - intuitively this should work, escape is using SeuratObject-based formatting.

Can you give me the output of sessionInfo() so I can ensure we get the right version for me to troubleshoot with?

Thanks,
Nick

[toobiwankenobi]

Hi @ncborcherding
Yea, sorry, forgot about that.

> sessionInfo()
R version 4.4.2 (2024-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.3.2

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: Europe/Zurich
tzcode source: internal

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

other attached packages:
 [1] BiocParallel_1.40.2 gridExtra_2.3       msigdbr_24.1.0      escape_2.2.3        tidyr_1.3.1         rstatix_0.7.2       ggpubr_0.6.0        ggbeeswarm_0.7.2   
 [9] ggplot2_3.5.2       stringr_1.5.1       dplyr_1.1.4         reticulate_1.42.0   Seurat_5.3.0        SeuratObject_5.1.0  sp_2.2-0           

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22            splines_4.4.2               later_1.4.2                 R.oo_1.27.1                 tibble_3.2.1               
  [6] polyclip_1.10-7             graph_1.84.1                XML_3.99-0.18               fastDummies_1.7.5           lifecycle_1.0.4            
 [11] globals_0.18.0              lattice_0.22-7              MASS_7.3-65                 ggdist_3.3.3                backports_1.5.0            
 [16] magrittr_2.0.3              limma_3.62.2                plotly_4.10.4               rmarkdown_2.29              yaml_2.3.10                
 [21] httpuv_1.6.16               sctransform_0.4.2           spam_2.11-1                 spatstat.sparse_3.1-0       cowplot_1.1.3              
 [26] pbapply_1.7-2               DBI_1.2.3                   RColorBrewer_1.1-3          pkgload_1.4.0               abind_1.4-8                
 [31] zlibbioc_1.52.0             GenomicRanges_1.58.0        Rtsne_0.17                  R.utils_2.13.0              purrr_1.0.4                
 [36] BiocGenerics_0.52.0         GenomeInfoDbData_1.2.13     IRanges_2.40.1              S4Vectors_0.44.0            ggrepel_0.9.6              
 [41] irlba_2.3.5.1               listenv_0.9.1               spatstat.utils_3.1-4        GSVA_2.0.7                  goftest_1.2-3              
 [46] RSpectra_0.16-2             spatstat.random_3.3-3       annotate_1.84.0             fitdistrplus_1.2-2          parallelly_1.44.0          
 [51] DelayedMatrixStats_1.28.1   codetools_0.2-20            DelayedArray_0.32.0         tidyselect_1.2.1            UCell_2.10.1               
 [56] UCSC.utils_1.2.0            farver_2.1.2                ScaledMatrix_1.14.0         matrixStats_1.5.0           stats4_4.4.2               
 [61] spatstat.explore_3.4-2      jsonlite_2.0.0              BiocNeighbors_2.0.1         progressr_0.15.1            Formula_1.2-5              
 [66] ggridges_0.5.6              survival_3.8-3              tools_4.4.2                 ica_1.0-3                   Rcpp_1.0.14                
 [71] glue_1.8.0                  SparseArray_1.6.2           xfun_0.52                   distributional_0.5.0        MatrixGenerics_1.18.1      
 [76] GenomeInfoDb_1.42.3         AUCell_1.28.0               HDF5Array_1.34.0            withr_3.0.2                 fastmap_1.2.0              
 [81] ggpointdensity_0.2.0        rhdf5filters_1.18.1         rsvd_1.0.5                  digest_0.6.37               R6_2.6.1                   
 [86] mime_0.13                   colorspace_2.1-1            scattermore_1.2             tensor_1.5                  spatstat.data_3.1-6        
 [91] RSQLite_2.3.11              R.methodsS3_1.8.2           generics_0.1.4              data.table_1.17.2           S4Arrays_1.6.0             
 [96] httr_1.4.7                  htmlwidgets_1.6.4           uwot_0.2.3                  pkgconfig_2.0.3             gtable_0.3.6               
[101] blob_1.2.4                  lmtest_0.9-40               SingleCellExperiment_1.28.1 XVector_0.46.0              htmltools_0.5.8.1          
[106] carData_3.0-5               dotCall64_1.2               GSEABase_1.68.0             scales_1.4.0                Biobase_2.66.0             
[111] png_0.1-8                   SpatialExperiment_1.16.0    spatstat.univar_3.1-3       knitr_1.50                  rstudioapi_0.17.1          
[116] rjson_0.2.23                reshape2_1.4.4              curl_6.2.2                  nlme_3.1-168                rhdf5_2.50.2               
[121] zoo_1.8-14                  cachem_1.1.0                KernSmooth_2.23-26          parallel_4.4.2              miniUI_0.1.2               
[126] vipor_0.4.7                 AnnotationDbi_1.68.0        pillar_1.10.2               grid_4.4.2                  vctrs_0.6.5                
[131] RANN_2.6.2                  promises_1.3.2              BiocSingular_1.22.0         car_3.1-3                   beachmat_2.22.0            
[136] xtable_1.8-4                cluster_2.1.8.1             beeswarm_0.4.0              evaluate_1.0.3              magick_2.8.6               
[141] cli_3.6.5                   compiler_4.4.2              rlang_1.1.6                 crayon_1.5.3                future.apply_1.11.3        
[146] ggsignif_0.6.4              labeling_0.4.3              plyr_1.8.9                  stringi_1.8.7               viridisLite_0.4.2          
[151] deldir_2.0-4                babelgene_22.9              assertthat_0.2.1            Biostrings_2.74.1           lazyeval_0.2.2             
[156] spatstat.geom_3.3-6         Matrix_1.7-3                RcppHNSW_0.6.0              patchwork_1.3.0             sparseMatrixStats_1.18.0   
[161] bit64_4.6.0-1               future_1.49.0               Rhdf5lib_1.28.0             statmod_1.5.0               KEGGREST_1.46.0            
[166] shiny_1.10.0                SummarizedExperiment_1.36.0 ROCR_1.0-11                 igraph_2.1.4                broom_1.0.8                
[171] memoise_2.0.1               bit_4.6.0    

[ncborcherding]

Hey @toobiwankenobi,

It looks like AggregateExpression() does not work on Seurat Objects with v3 assay information (this is because there is no slot argument). If you update the Seurat object and then run it, it should work:

GS <- list(Bcells = c("MS4A1", "CD79B", "CD79A", "IGH1", "IGH2"),
           Tcells = c("CD3E", "CD3D", "CD3G", "CD7","CD8A"))
pbmc_small <- UpdateSeuratObject(SeuratObject::pbmc_small) #Here is my change

pbmc_small <- runEscape(pbmc_small, 
                        gene.sets = GS, 
                        min.size = NULL)

pbmc_aggr <- AggregateExpression(pbmc_small, 
                                 assays = "escape",
                                 return.seurat = T, 
                                 group.by = "groups")

Hope that helps and let me know if you have any other questions,

Nick

[toobiwankenobi]

Hi @ncborcherding

Thank you for the feedback. Now it works, but the returned layer does not contain any values. I had the same problem with my data. Furthermore, the row/column names for the output are not indicated - I don't see a way to deconvolute what the elements of the output correspond to. AggregateExpression() would normally name the rows and columns correctly. Any ideas why this is not happening?

GS <- list(Bcells = c("MS4A1", "CD79B", "CD79A", "IGH1", "IGH2"),
           Tcells = c("CD3E", "CD3D", "CD3G", "CD7","CD8A"))
pbmc_small <- UpdateSeuratObject(SeuratObject::pbmc_small) #Here is my change
pbmc_small <- runEscape(pbmc_small, 
                        gene.sets = GS, 
                        min.size = NULL)

pbmc_aggr <- AggregateExpression(pbmc_small, 
                                 assays = "escape",
                                 return.seurat = T, 
                                 group.by = "groups")

pbmc_aggr@assays$escape$data

>
         Cell_1 Cell_2
Feature1    NaN    NaN
Feature2    NaN    NaN

Example with correct naming of aggregated data

GS <- list(Bcells = c("MS4A1", "CD79B", "CD79A", "IGH1", "IGH2"),
           Tcells = c("CD3E", "CD3D", "CD3G", "CD7","CD8A"))
pbmc_small <- UpdateSeuratObject(SeuratObject::pbmc_small) #Here is my change
pbmc_aggr <- AggregateExpression(pbmc_small,  
                                 assays = "RNA",
                                 return.seurat = T, 
                                 group.by = "groups")
head(pbmc_aggr@assays$RNA$data)
> 
              g1       g2
MS4A1    2.424756 3.283687
CD79B    3.072645 3.477518
CD79A    3.271565 3.561953
HLA-DRA  6.443298 5.999495
TCL1A    2.188596 3.043027
HLA-DQB1 4.374933 4.010973

Thank you for your help!

Tobi

[ncborcherding]

Apologies, you need to declare the slot:

> pbmc_aggr <- AggregateExpression(pbmc_small, 
                                 assays = "escape",
                                 return.seurat = T, 
                                 slot = "data",
                                 group.by = "groups")

> pbmc_aggr@assays$escape$data

pbmc_aggr@assays$escape$data
g1 g2
Bcells 9.210264 9.189247
Tcells 1.015096 5.350531