bedapub/besca

The R code in function `valOutlier` needs modification

Closed this issue · 5 comments

I'm using Seurat v4.0.5 and it seems that the function valOutlier doesn't work with this version of Seurat:

rpy2.rinterface_lib.embedded.RRuntimeError: Error in normarg_assays(assays, as.null.if.no.assay = TRUE) : 
  trying to get slot "counts" from an object of a basic class ("NULL") with no slots

According to the bug report, I checked the seurat_obj in RStudio and found that the data was saved in originalexp assay, not RNA (in Seurat 3.0). And after I replaced the RNA with originalexp in valOutlier:

ro.r(
        "dat =  SingleCellExperiment(assays = list(counts=seurat_obj@assays$originalexp@counts) )"
    )

everything was alright. So maybe you can add if and else statements to judge the version of Seurat and determine which name of the assay you should use.

Dear @ToryDeng,
you say that the internal data format of the Seurat Object has changed from v3 to v4. But for me all worked fine. Maybe I miss something, but I tried to reproduce your issue using the following code:

# Install specific version of a package
require(devtools)
install_version("Seurat", version = "4.1.1", repos = "http://cran.us.r-project.org")
# install_version("Seurat", version = "3.2.3", repos = "http://cran.us.r-project.org")
install.packages("hdf5r", configure.args = c(hdf5r = "--with-hdf5=/opt/homebrew/bin/h5cc"))

# Download Seurat example data here:
# https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
# and stored it to ~/Downloads/

require(Seurat)
pbmc.data <- Read10X(data.dir = "~/Downloads/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)

#Works
pbmc@assays$RNA@counts

# Not working
pbmc@assays$originalexp@counts
# Error: trying to get slot "counts" from an object of a basic class ("NULL") with no slots

sessionInfo()

# Note: Based on this tutorial:
# https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

Can you try this little code snippet yourself and give some feedback?

Dear @kohleman,

I have tested your code. Yes, if you directly import data using R and convert it to the Seurat object, the count matrix is in the RNA assay. However, in Besaca you actually convert the AnnData to SingleCellExperiment by the Python package anndata2ri at first, and then you convert the SingleCellExperiment to Seurat object in R:

ro.globalenv["dat"] = adata
ro.globalenv["sym"] = adata.var["SYMBOL"]
ro.r('seurat_obj = as.Seurat(dat, counts="X", data = NULL)')
ro.r(
       "dat =  SingleCellExperiment(assays = list(counts=seurat_obj@assays$RNA@counts) )"
   )

So if you want to reproduce the error, you can try the following code:

# in Python
import anndata2ri
from rpy2.robjects import r, globalenv

anndata2ri.activate()
globalenv['dat'] = adata
r("save(dat, file = '/path/to/Rdata')")
anndata2ri.deactivate()
# in R
require(Seurat)

load("/path/to/Rdata")
seurat_obj = as.Seurat(dat, counts="X", data = NULL)
dat =  SingleCellExperiment(assays = list(counts=seurat_obj@assays$RNA@counts))

and you should see an error report:

Error in normarg_assays(assays, as.null.if.no.assay = TRUE) : 
  trying to get slot "counts" from an object of a basic class ("NULL") with no slots

and if you check the assay in seurat_obj, you should see similar output like this:

> seurat_obj@assays
$originalexp
Assay data with 772 features for 3331 cells
First 10 features:
 ISG15, AURKAIP1, SSU72, RPL22, PARK7, ENO1, CAPZB, MINOS1, DDOST,
CDC42

I guess the conversion between SingleCellExperiment and Seurat object using the function as.Seurat makes the error happen. Actually, I checked the source code of as.Seurat.SingleCellExperiment in Seurat v4.1.1 and found the following code:

if (packageVersion(pkg = "SingleCellExperiment") >= "1.14.0") {
    orig.exp <- SingleCellExperiment::mainExpName(x = x) %||% "originalexp"
  } else {
    orig.exp <- "originalexp"
  }
  if (!is.null(SingleCellExperiment::altExpNames(x = x))) {
    assayn <- assay %||% SingleCellExperiment::altExpNames(x = x)
    if (!all(assay %in% SingleCellExperiment::altExpNames(x = x))) {
      stop("One or more of the assays you are trying to convert is not in the SingleCellExperiment object")
    }
    assayn <- c(orig.exp, assayn)
  } else {
    assayn <- orig.exp
  }

I hope this will help you.

Thanks for your explanations and code samples. I will try it out.

Dear @ToryDeng,
I will change the code to this:

if (packageVersion("Seurat") >= "4.0" ) {
  print("Detected Seurat Version >= 4")
  dat =  SingleCellExperiment(assays = list(counts=seurat_obj@assays$originalexp@counts))
  
} else {
  print("Detected Seurat Version < 4")
  dat =  SingleCellExperiment(assays = list(counts=seurat_obj@assays$RNA@counts))
}

fixed in 9945d73