Error in enrichIt()

[ncborcherding]

Error in enrichIt()

[ncborcherding]

Hey Khushbu,

What version of escape are you using?

I recently implemented a filtering step that checks the count matrix to see what genes are expressed and automatically remove gene sets with less than a specific number of genes.

enrichIt <- function(obj, gene.sets = NULL, 
                     method = "ssGSEA", groups = 1000, cores = 2,
                     min.size = 5)

Here min.size = 5 means there needs to be at least 5 genes in the count matrix. This could be giving the error above.

If not, please let me know and I'd be happy to troubleshoot with you further.

Nick

[kpatel427]

Hi Nick,

Thank you for responding.

I am using escape_1.3.1.

I tried giving a min.size = 5 parameter, but it throws me this error - Error in enrichIt(obj = seurat_integrated, gene.sets = gene.sets, groups = 1000, : unused argument (min.size = 5)

I tried re-running with even smaller min.sizes but got the same error.

Thanks,
Khushbu

[ncborcherding]

Hey Khushbu,

So the version you are using does not have the min.size parameter. However the error could also be occurring at the time of the enrichment scoring - have you check your count matrix to see if the MES and ADRN genes are expressed? If there is 0 variation in the genes or 0 count values, the above error could be outputted.

Nick

[kpatel427]

Hi Nick,

I checked my count matrices and I do see counts for these genes in the cells.
Here's a subset of count matrix for these genes -

MES genes:

PRDM6  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
LIPA   . 1 . . . . . . . . . . . . . . . . . 2 . . 1 . . . . 1 . . 1 . . . . . . . 2 . . 1 1 . 1 . . . 1 . ......
ANTXR1 1 . 1 . 5 1 3 . . 1 . . 1 1 . . . . 2 1 2 . 1 . 1 . . . . 1 . . . 1 3 . 2 . 1 2 . 2 2 1 . 1 . . . 1 ......
SMAD3  . . . 1 . 1 1 . . . . . . . 1 . 2 . . 1 . . 1 . . . . . . 1 1 1 . . . . . . . . . . . 1 1 . . . . . ......```


ADRN genes:

GAP43   . . 1 . . . . . . 1 . . . . 1 . . 1 1 . . . 1 1 1 . . . . . . 3 . . . . 2 1 . . . . . . . 1 1 . 2 2 ......
DCX     1 . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . ......
NOL4    . 1 2 1 1 1 2 . 2 . . 2 . 3 3 2 1 2 . 4 . 1 . 3 1 . . . . . . . . 4 . 1 2 1 . . . . 2 . . . 1 2 3 . ......
BIRC5   1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . ......

Could you please tell me how do I install the newest version of escape? I installed using devtools::install_github("ncborcherding/escape")

Thanks,
Khushbu

[ncborcherding]

Hey Khushbu,

You can install the most-up-to-date version using the code below - thinking about it, this is probably why you got the error for UCell as it is only implemented right now in the dev version.

devtools::install_github("ncborcherding/escape@dev")

Going over your initial error - Do you have GSEAbase installed? What is the output of your sessionInfo()? Maybe the error could be arising from the pre-enrichment checks in the function?

Nick

[kpatel427]

Hi Nick,

Unfortunately, my R version does not support dev version. It requires R >=4.1.

And yes, I have GSEAbase installed. Here's the output from my sessionInfo():

R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

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

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

other attached packages:
 [1] GSEABase_1.52.1      graph_1.68.0         annotate_1.68.0      XML_4.0-0            escape_1.3.1         fgsea_1.16.0         org.Hs.eg.db_3.12.0  reactome.db_1.74.0  
 [9] AnnotationDbi_1.52.0 Biobase_2.50.0       GenomicRanges_1.42.0 GenomeInfoDb_1.26.7  IRanges_2.24.1       png_0.1-7            apeglm_1.12.0        S4Vectors_0.28.1    
[17] BiocGenerics_0.36.1  reshape2_1.4.4       Matrix_1.3-4         magrittr_2.0.1       edgeR_3.32.1         limma_3.46.0         cowplot_1.1.1        singscore_1.13.1    
[25] dittoSeq_1.2.6       viridis_0.6.1        viridisLite_0.4.0    RColorBrewer_1.1-2   pheatmap_1.0.12      forcats_0.5.1        stringr_1.4.0        dplyr_1.0.7         
[33] purrr_0.3.4          readr_2.0.2          tidyr_1.1.4          tibble_3.1.5         ggplot2_3.3.5        tidyverse_1.3.1      Seurat_4.0.4         SeuratObject_4.0.2  

loaded via a namespace (and not attached):
  [1] utf8_1.2.2                  reticulate_1.22             tidyselect_1.1.1            RSQLite_2.2.8               htmlwidgets_1.5.4           BiocParallel_1.24.1        
  [7] grid_4.0.3                  Rtsne_0.15                  devtools_2.4.2              munsell_0.5.0               codetools_0.2-18            ragg_1.1.3                 
 [13] ica_1.0-2                   future_1.22.1               miniUI_0.1.1.1              withr_2.4.2                 colorspace_2.0-2            rstudioapi_0.13            
 [19] SingleCellExperiment_1.12.0 ROCR_1.0-11                 tensor_1.5                  listenv_0.8.0               MatrixGenerics_1.2.1        labeling_0.4.2             
 [25] bbmle_1.0.24                GenomeInfoDbData_1.2.4      polyclip_1.10-0             bit64_4.0.5                 farver_2.1.0                rprojroot_2.0.2            
 [31] coda_0.19-4                 parallelly_1.28.1           vctrs_0.3.8                 generics_0.1.0              R6_2.5.1                    msigdbr_7.4.1              
 [37] locfit_1.5-9.4              bitops_1.0-7                spatstat.utils_2.2-0        cachem_1.0.6                DelayedArray_0.16.3         assertthat_0.2.1           
 [43] promises_1.2.0.1            scales_1.1.1                gtable_0.3.0                globals_0.14.0              processx_3.5.2              goftest_1.2-2              
 [49] rlang_0.4.11                genefilter_1.72.1           systemfonts_1.0.2           splines_4.0.3               lazyeval_0.2.2              spatstat.geom_2.2-2        
 [55] broom_0.7.9                 BiocManager_1.30.16         yaml_2.2.1                  abind_1.4-5                 modelr_0.1.8                backports_1.2.1            
 [61] httpuv_1.6.3                tools_4.0.3                 usethis_2.0.1               ellipsis_0.3.2              spatstat.core_2.3-0         sessioninfo_1.1.1          
 [67] ggridges_0.5.3              Rcpp_1.0.7                  plyr_1.8.6                  zlibbioc_1.36.0             RCurl_1.98-1.5              ps_1.6.0                   
 [73] prettyunits_1.1.1           rpart_4.1-15                deldir_0.2-10               pbapply_1.5-0               zoo_1.8-9                   SummarizedExperiment_1.20.0
 [79] haven_2.4.3                 ggrepel_0.9.1               cluster_2.1.2               fs_1.5.0                    data.table_1.14.2           scattermore_0.7            
 [85] lmtest_0.9-38               reprex_2.0.1                RANN_2.6.1                  mvtnorm_1.1-2               fitdistrplus_1.1-6          matrixStats_0.61.0         
 [91] pkgload_1.2.2               GSVA_1.38.2                 hms_1.1.1                   patchwork_1.1.1             mime_0.12                   xtable_1.8-4               
 [97] emdbook_1.3.12              readxl_1.3.1                gridExtra_2.3               bdsmatrix_1.3-4             testthat_3.1.0              compiler_4.0.3             
[103] KernSmooth_2.23-20          crayon_1.4.1                htmltools_0.5.2             mgcv_1.8-37                 later_1.3.0                 tzdb_0.1.2                 
[109] geneplotter_1.68.0          lubridate_1.7.10            DBI_1.1.1                   dbplyr_2.1.1                MASS_7.3-54                 babelgene_21.4             
[115] cli_3.0.1                   igraph_1.2.6                pkgconfig_2.0.3             numDeriv_2016.8-1.1         plotly_4.9.4.1              spatstat.sparse_2.0-0      
[121] xml2_1.3.2                  XVector_0.30.0              rvest_1.0.1                 callr_3.7.0                 digest_0.6.28               sctransform_0.3.2          
[127] RcppAnnoy_0.0.19            spatstat.data_2.1-0         fastmatch_1.1-3             cellranger_1.1.0            leiden_0.3.9                uwot_0.1.10                
[133] curl_4.3.2                  shiny_1.7.1                 lifecycle_1.0.1             nlme_3.1-153                jsonlite_1.7.2              desc_1.4.0                 
[139] fansi_0.5.0                 pillar_1.6.3                lattice_0.20-45             fastmap_1.1.0               httr_1.4.2                  pkgbuild_1.2.0             
[145] survival_3.2-13             glue_1.4.2                  remotes_2.4.1               bit_4.0.4                   stringi_1.7.5               blob_1.2.2                 
[151] textshaping_0.3.5           DESeq2_1.30.1               memoise_2.0.0               irlba_2.3.3                 future.apply_1.8.1  

[ncborcherding]

Hey Khushbu,

Sorry for the delay - everything looks right to me from you session info. I have not been able to recreate the error unfortunately from my end.

I do have a follow up question - what is the exact code you are using to create the gene set list ?

gene.sets <- list(MES=c(genes1),
                  ADRN=c(genes2))

In version 1.3.2, there is no filtering of the gene sets before the enrich function and in fact, the only evaluation of the gene set object is where you are getting an error (below) which transforms a gene set collection object into a list (which is what you are entering). But the error implies you gene set (egc in the function call) has a length of 0?

#Line60
if(inherits(egc, what = "GeneSetCollection")){
        egc <- GSEABase::geneIds(egc) # will return a simple list, 
        #which will work if a matrix is supplied to GSVA
    }

Nick

[kpatel427]

Hi Nick,

Following is the code I am using to create the gene.sets:

genes1 <- c("BNC2", "ROBO1",  "COL1A1", "COL3A1", "GPC6", "PTPRK","APOE", "PDE3A", "DMD", "PTPRG")
genes2 <- c("CRH", "SYNPO2",  "TMEM108", "BMPR1B",  "RRM2", "ADGRB3",  "ESRRG", "PRSS12",  "EYA1", "HS6ST2")

gene.sets <- list(MES=c(genes1),
                  ADRN=c(genes2))

And here's a peak into gene.sets list:

gene.sets[['MES']]
  [1] "BNC2"      "ROBO1"     "COL1A1"    "COL3A1"    "GPC6"      "PTPRK"     "APOE"      "PDE3A"     "DMD"    "PTPRG"
  
gene.sets[['ADRN']]
 [1] "CRH"     "SYNPO2"  "TMEM108" "BMPR1B"  "RRM2"    "ADGRB3"  "ESRRG"   "PRSS12"  "EYA1"    "HS6ST2" 

Thanks,
Khushbu

[ncborcherding]

Hey Khushbu,

Your gene set list looks formed right, not sure what is going wrong. I've installed the 1.3.2 version and matching packages based on your session info and the gene sets you used above and cannot reproduce the error.

And the output of class(gene.sets) is a list? Do you have anything else in your global environment named egc? This should not effect it, but maybe?

Nick