Cannot use expressed_genes_receiver with assay_oi = 'SCT'

[jaquol]

I am trying to run nichenetr on a Seurat object resulting from the Seurat integration workflow. When I first try to determine which genes are expressed with:

expressed_genes_receiver = get_expressed_genes(receiver, seuratObj, pct = 0.10)

I get this warning:

Seurat object is result from the Seurat integration workflow. Make sure that this way of defining expressed genes is appropriate for your integrated data. The expressed genes are defined based on the integrated slot

My reasoning was then to use the SCT assay instead with:

expressed_genes_receiver = get_expressed_genes(receiver, seuratObj, pct = 0.10, assay_oi = 'SCT')

But then I get this error:

Error in match(x, table, nomatch = 0L) : 
  'match' requires vector arguments

I looked into the get_expressed_genes code and I narrowed the error to this part:

    if (!is.null(assay_oi)) {
        if (!assay_oi %in% Assays(seurat_obj)) {
            stop("assay_oi should be an assay of your Seurat object")
        }
    }

!assay_oi %in% Assays(seurat_obj) cannot be assessed because, I understand, a vector is expected but the output of Assays(seurat_obj) is a list.

My first question is whether my reasoning of forcing the use of the SCT assay is correct.

My second question is if you could please help in overcoming the issue here reported.

Thank you in advance!

Javier

[browaeysrobin]

Hi @jaquol,

I would suggest setting assay_oi = 'RNA' to define expressed genes. Can you check whether this works?

[jaquol]

Hi @browaeysrobin,

Thank you for your quick reply! I am afraid that assay_oi = 'RNA' does not work either on my side; I get the same error as in my previous message.

[jaquol]

Hi @browaeysrobin,

I thought that I might be able to overcome the issue above by deleting the "SCT" and/or "integrated" assays. Inspired by this, I deleted both the "SCT" and "integrated" assays with:

seuratObj@assays$integrated <- NULL
seuratObj@assays$SCT <- NULL

Unlike previously, this way I do not get any warning when running the get_expressed_genes function:

expressed_genes_receiver = get_expressed_genes(receiver, seuratObj, pct = 0.10)

My questions are then:

1. Does deleting the "SCT" and "integrated" assays makes sense?
2. Is "RNA" the preferred assay to run nichenetr on?

Thank you!

[browaeysrobin]

Hi @jaquol

Deleting the SCT and integrated assays is only a temporary solution to this error of course.
I tried to reproduce the error, but again, everything runs fine with me and I can select different assays (RNA, SCT) via the get_expressed_genes function. In my case Assays(seurat_obj) returns a vector and not a list, in contrast to what you mentioned before. Maybe there is again an issue with a conflict with another package (as for your other issue)?

The selection of which slot (RNA/SCT/integrated) in this get_expressed_genes() function will only effect which genes will be considered to be expressed (which affects the background of genes for the ligand activity analysis and the selection of ligands/receptors used as input for this). This warning that this function throws when having a SCT and/or integrated assay in addition to the RNA assay is just to let users think about whether the way get_expressed_genes() defines expressed genes is a correct way to do it for transformed data. get_expressed_genes() will look at the fraction of cells having a non-zero expression value for a specific gene. Because data transformation could change zeroes to non-zeroes and vice versa, I would recommend using the unaffected RNA slot. But, out of experience I can say that choosing SCT / RNA / integrated does not make much difference in which genes you get out as expressed in the end. Typically, choosing RNA will result in a slightly higher number of expressed genes, but differences in number are minor (eg 5950 vs 6000 expressed genes).

[tkapello]

I started implementing the NicheNet code to my sc dataset and I wanted to report the same issue. Switching the argument to "RNA" does not work:

>expressed_genes_receiver = get_expressed_genes(receiver, baldataset, pct = 0.10, assay_oi = 'RNA')
Error in .normarg_assays(assays) : 
  'assays' must be a SimpleList, list or array

In my case, the integrated gene list is a lot smaller so I d rather use the RNA assay to have a bigger expressed gene list. If deleting the integrated slot works for some people, I would do it, but of course, an alternative would be more preferable.

Thanks in advance,
Theo

[browaeysrobin]

Hi @tkapello

Is it possible you loaded the SummarizedExperiment package? If yes, could you try to run this code again without loading SummarizedExperiment?

[tkapello]

@browaeysrobin

SummarizedExperiment seems to be needed, otherwise:

Error in Assays(seurat_obj) : could not find function "Assays"

[browaeysrobin]

Hi @tkapello
Normally, SummarizedExperiment is not needed, because Assays() is also function of the Seurat package. So: did you load the Seurat package (should be yes)?

I also was able to replicate the error you both mentioned here above if SummarizedExperiment was loaded after loading nichenetr/Seurat/tidyverse (because it overwrites some functions of Seurat/base R/tidyverse and loads other packages that also overwrite functions of Seurat/base R/tidyverse). However, no error was thrown when SummarizedExperiment was loaded before.

So: I don't think there is an issue if you do the following

library(nichenetr)
library(Seurat)
library(tidyverse)
seurat_obj = readRDS("data/seuratObj.rds")
expressed_genes_receiver = get_expressed_genes("celltype X", seurat_obj, pct = 0.10)

Here, Assays should be used from Seurat.
If you do want to load other packages, the following should work:

library(SummarizedExperiments) # and other packages if really necessary
library(nichenetr)
library(Seurat)
library(tidyverse)
seurat_obj = readRDS("data/seuratObj.rds")
expressed_genes_receiver = get_expressed_genes("celltype X", seurat_obj, pct = 0.10)

In the near future, I will change the underlying code to avoid these conflicts such that the order in which you load the packages should not matter.

[tkapello]

Hi @browaeysrobin

I tried the order of the packages you suggested, but the error remains:

Error in .normarg_assays(assays) : 
  'assays' must be a SimpleList, list or array

Does the Seurat object need to be loaded as a .rds? I load a class Seurat object from my environment.

[browaeysrobin]

Hi @tkapello

Could you once try this to see whether this resolves the issue?

[tkapello]

@browaeysrobin

I tried everything you suggested and still I get the same error :-(

Error in .normarg_assays(assays) : 
  'assays' must be a SimpleList, list or array

[browaeysrobin]

Hi @tkapello

Could you try to

1. re-install nichenetr and run your script again ?
2. let me know whether it works now and if not: provide me with your sessionInfo() ?

[tkapello]

Hi @browaeysrobin

I uninstalled and re-installed nichenetr and now I have the following error:

> # Define a receiver cell type
> receiver = "Neutrophil_0"
> expressed_genes_receiver = get_expressed_genes(receiver, baldataset, pct = 0.10, assay_oi = 'RNA')
Error: 'Assays' is not an exported object from 'namespace:Seurat'

As you requested, my SessionInfo():

> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS: /usr/lib/openblas-base/libblas.so.3
LAPACK: /usr/lib/libopenblasp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] circlize_0.4.6              forcats_0.4.0               stringr_1.4.0              
 [4] dplyr_0.8.0.1               purrr_0.3.2                 readr_1.3.1                
 [7] tidyr_0.8.3                 tibble_2.1.1                ggplot2_3.1.0              
[10] tidyverse_1.2.1             Seurat_3.0.0                nichenetr_0.1.0            
[13] SummarizedExperiment_1.12.0 DelayedArray_0.8.0          BiocParallel_1.16.6        
[16] matrixStats_0.54.0          Biobase_2.42.0              GenomicRanges_1.34.0       
[19] GenomeInfoDb_1.18.2         IRanges_2.16.0              S4Vectors_0.20.1           
[22] BiocGenerics_0.28.0         BiocManager_1.30.4         

loaded via a namespace (and not attached):
  [1] reticulate_1.12        R.utils_2.8.0          tidyselect_0.2.5       RSQLite_2.1.1         
  [5] AnnotationDbi_1.44.0   htmlwidgets_1.3        grid_3.5.3             Rtsne_0.15            
  [9] munsell_0.5.0          codetools_0.2-16       ica_1.0-2              future_1.12.0         
 [13] withr_2.1.2            GOSemSim_2.8.0         colorspace_1.4-1       rgexf_0.15.3          
 [17] knitr_1.22             rstudioapi_0.10        ROCR_1.0-7             DOSE_3.8.2            
 [21] gbRd_0.4-11            listenv_0.7.0          Rdpack_0.11-0          GenomeInfoDbData_1.2.0
 [25] bit64_0.9-7            downloader_0.4         generics_0.0.2         ipred_0.9-8           
 [29] xfun_0.6               randomForest_4.6-14    R6_2.4.0               rsvd_1.0.0            
 [33] fgsea_1.8.0            bitops_1.0-6           assertthat_0.2.1       SDMTools_1.1-221      
 [37] scales_1.0.0           nnet_7.3-12            gtable_0.3.0           npsurv_0.4-0          
 [41] globals_0.12.4         timeDate_3043.102      rlang_0.3.3            GlobalOptions_0.1.0   
 [45] splines_3.5.3          lazyeval_0.2.2         ModelMetrics_1.2.2     acepack_1.4.1         
 [49] broom_0.5.1            brew_1.0-6             checkmate_1.9.1        yaml_2.2.0            
 [53] reshape2_1.4.3         modelr_0.1.4           backports_1.1.3        qvalue_2.14.1         
 [57] Hmisc_4.2-0            caret_6.0-82           DiagrammeR_1.0.0       tools_3.5.3           
 [61] lava_1.6.5             influenceR_0.1.0       gplots_3.0.1.1         RColorBrewer_1.1-2    
 [65] ggridges_0.5.1         Rcpp_1.0.1             plyr_1.8.4             base64enc_0.1-3       
 [69] visNetwork_2.0.6       zlibbioc_1.28.0        RCurl_1.95-4.12        rpart_4.1-13          
 [73] pbapply_1.4-0          viridis_0.5.1          cowplot_0.9.4          zoo_1.8-5             
 [77] haven_2.1.0            ggrepel_0.8.0          cluster_2.0.7-1        magrittr_1.5          
 [81] data.table_1.12.0      DO.db_2.9              lmtest_0.9-36          RANN_2.6.1            
 [85] fitdistrplus_1.0-14    hms_0.4.2              lsei_1.2-0             XML_3.98-1.19         
 [89] readxl_1.3.1           gridExtra_2.3          shape_1.4.4            compiler_3.5.3        
 [93] KernSmooth_2.23-15     crayon_1.3.4           R.oo_1.22.0            htmltools_0.3.6       
 [97] Formula_1.2-3          lubridate_1.7.4        DBI_1.0.0              MASS_7.3-51.3         
[101] Matrix_1.2-17          cli_1.1.0              R.methodsS3_1.7.1      gdata_2.18.0          
[105] metap_1.1              gower_0.2.0            igraph_1.2.4           pkgconfig_2.0.2       
[109] foreign_0.8-71         plotly_4.9.0           recipes_0.1.5          xml2_1.2.0            
[113] foreach_1.4.4          XVector_0.22.0         prodlim_2018.04.18     bibtex_0.4.2          
[117] rvest_0.3.2            digest_0.6.18          sctransform_0.2.0      tsne_0.1-3            
[121] fastmatch_1.1-0        cellranger_1.1.0       htmlTable_1.13.1       Rook_1.1-1            
[125] gtools_3.8.1           nlme_3.1-137           jsonlite_1.6           viridisLite_0.3.0     
[129] limma_3.38.3           pillar_1.3.1           lattice_0.20-38        GO.db_3.6.0           
[133] httr_1.4.0             survival_2.44-1.1      glue_1.3.1             fdrtool_1.2.15        
[137] png_0.1-7              iterators_1.0.10       bit_1.1-14             class_7.3-15          
[141] stringi_1.4.3          blob_1.1.1             memoise_1.1.0          latticeExtra_0.6-28   
[145] caTools_1.17.1.2       irlba_2.3.3            future.apply_1.2.0     ape_5.3 

[browaeysrobin]

Hi @tkapello

It seems that your Seurat version is not the most recent. Could you try to reinstall the newest version?
Can you then first try to run Seurat::Assays() on your seurat object, and then run get_expressed_genes() from nichenetr again? And let me know which of these tests did work (and which not).

[tkapello]

@browaeysrobin

you were right, it was the Seurat version. There was no Assays function