ncborcherding / escape

Easy single cell analysis platform for enrichment
https://www.borch.dev/uploads/screpertoire/articles/running_escape
MIT License
139 stars 20 forks source link

Error in enrichIt() #31

Closed kpatel427 closed 2 years ago

kpatel427 commented 2 years ago

Hello,

Thank you for creating such a useful package.

I am following your vignette and trying to perform ssGSEA on my single cell RNA-Seq data. However, I am running into this error - Error in if (attr(class(egc), "package") == "GSEABase") { : argument is of length zero

Here's my code:

gene.sets <- list(MES=c(genes1),
                  ADRN=c(genes2))
ES <- enrichIt(obj = seurat_integrated, 
               gene.sets = gene.sets, 
               groups = 1000, cores = 4)

Where genes1 and genes2 are vectors containing list of genes.

I also tried running ES <- enrichIt(obj = seurat_integrated, method = "UCell", gene.sets = gene.sets) but that throws an error too - Error in enrichIt(obj = seurat_integrated, method = "UCell", gene.sets = gene.sets) : unused argument (method = "UCell")

My R version is 4.0.3. I am not sure what I am missing here. Any help is greatly appreciated.

Thanks, Khushbu

ncborcherding commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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