Seurat Weekly NO.06 || Scanpy2Seurat

天子呼来不上船，

自称臣是菜鸟团

在这里，和国际同行一起学习单细胞数据分析。

其实，我们在2019年的时候就介绍过单细胞转录组数据分析||Seurat3.1教程:Interoperability between single-cell object formats，讲了单细胞转录组数据对象的转化。对R语言境内的Seurat，CellDataSet，SingleCellExperiment，loom的格式转化起来还是比较方便的，但是对于异域的anndata转化一直不是很友好，所以我借此机会学会了python（在等短信验证码的那六十秒之内）。anndata的数据就在python中分析，完事。

但是这样的转化总有需求，于是，Seurat团队开发了SeuratDisk包，希望满足数据在Seurat和Scanpy之间快速搬家的需求。遇到规范的数据，当然可以一键搬家，但是如果有一点不同，就会带来各种Error。这里我们就以此为契机，看看遇到Error该如何处理，这属于进阶课程了，遇到问题Google解决不了了，我们怎么办？

我从网上下了一个不知道做了什么处理的单细胞数据，只知道是h5ad格式的，当我们用ReadH5AD读取时

library(Seurat)
cellxgene1 <- ReadH5AD(file = "some.processed.h5ad")

给出了这样的报错：

# 报错信息太长，我们只显示最有用的。In addition: Warning message:
Functionality for reading and writing H5AD files is being moved to SeuratDisk
For more details, please see https://github.com/mojaveazure/seurat-disk
and https://mojaveazure.github.io/seurat-disk/index.html

显然，作者这是在提示我们安装新的R包：seurat-disk，于是我们挺听话地去安装了。

remotes::install_github("mojaveazure/seurat-disk")
library(SeuratDisk)
Convert("some.processed.h5ad", dest = "h5seurat", overwrite = TRUE)`
Warning: Unknown file type: h5ad
Warning: 'assay' not set, setting to 'RNA'
Creating h5Seurat file for version 3.1.2
Adding X as data
Adding X as counts
Error: Cannot find feature names in this H5AD file

同样，转角处遇到Error，于是直接google Error信息，我们看到Seurat团队给出的答案：

For this specific H5AD file, it’s compressed using the LZF filter. This  filter is Python-specific, and cannot easily be used in R. To use this  file with Seurat and SeuratDisk, you’ll need to read it in Python and  save it out using the gzip compression

import anndata
adata = anndata.read("some.processed.h5ad")
adata.write("some.processed.gzip.h5ad", compression="gzip")

这显然是python语法，在R里面该如何操作呢？自然是在R里调用Python了，所以别问人家用的是R还是python了，这两个可以相互运行的语言，其实是一种语言。

library(reticulate)
reticulate::py_config()
Sys.which('python')  # 该python 下要安装了anndata 
# usethis::edit_r_environ()

filesmy2 ='some.processed.gzip.h5ad'
ad <- import("anndata", convert = FALSE)
some_ad <- ad$read_h5ad(filesmy)
adata = anndata.read("some.processed.h5ad")
some_ad$write("some.processed.gzip.h5ad", compression="gzip")

之后，我们终于可以转化了。

Convert("some.processed.gzip.h5ad", dest = "h5seurat", overwrite = TRUE)
Warning: Unknown file type: h5ad
Warning: 'assay' not set, setting to 'RNA' # 这里其实要格外小心，看看数据里面是不是 RNA啊
Creating h5Seurat file for version 3.1.5.9900
Adding X as data
Adding raw/X as counts
Adding meta.features from raw/var
Adding X_umap as cell embeddings for umap

Convert 成功后目录下多了一个文件：some.processed.gzip.h5seurat，就差一个Seurat对象了。

cellxgene <- LoadH5Seurat('some.processed.gzip.h5seurat',
                       assays ="RNA")

Validating h5Seurat file
Initializing RNA with data
Error in dimnamesGets(x, value) : 
  invalid dimnames given for “dgCMatrix” object

同样，在转角处遇到了Error ，于是我们再次Google 这个Error 。一番浏览，我们发现自己遇到了Google解决不了的问题。

决心debug这个函数。

debug(LoadH5Seurat)
cellxgene <- LoadH5Seurat('some.processed.gzip.h5seurat',
                       assays ="RNA")
# 一波回车
debug( as.Seurat) #因为LoadH5Seurat里面用了这个函数，所以在LoadH5Seurat的debug环境种再debug  as.Seurat
# 一波回车
debug(AssembleAssay) #因为as.Seurat里面用了这个函数，所以在as.Seurat的debug环境种再debug  AssembleAssay
# 一波回车

我们找到问题所在了：

slots.assay <- names(x = Filter(f = isTRUE, x = index[[assay]]$slots))
  slots <- slots %||% slots.assay
  slots <- match.arg(arg = slots, choices = slots.assay, several.ok = TRUE)
  if (!any(c("counts", "data") %in% slots)) {
    stop("At least one of 'counts' or 'data' must be loaded", 
      call. = FALSE)
  }
  assay.group <- file[["assays"]][[assay]]
  features <- assay.group[["features"]][]
  if ("counts" %in% slots && !"data" %in% slots) {
    if (verbose) {
      message("Initializing ", assay, " with counts")
    }
    counts <- as.matrix(x = assay.group[["counts"]])
    rownames(x = counts) <- features # 还是第一个features <- assay.group[["features"]][] 
    colnames(x = counts) <- Cells(x = file)
    obj <- CreateAssayObject(counts = counts, min.cells = -1, 
      min.features = -1)
  }
  else {
    if (verbose) {
      message("Initializing ", assay, " with data")
    }
    data <- as.matrix(x = assay.group[["data"]])
    rownames(x = data) <- features  # 还是第一个features <- assay.group[["features"]][]
    colnames(x = data) <- Cells(x = file)
    obj <- CreateAssayObject(data = data)
  }

也就是在这部分代码中作者是认为，slots 之counts和data 的行名是一样的，其实我们知道Seurat的每部分存的数据其实都是可以独立操作的，所以可能并不是都一样。

怀疑就要检查

Browse[8]> counts <- as.matrix(x = assay.group[["counts"]])
Browse[8]> dim(counts)
[1] 33567 10550
Browse[8]> data <- as.matrix(x = assay.group[["data"]])
Browse[8]> dim(data)
[1] 33421 10550
Browse[8]> length(features )
[1] 33567

果然不一样。

既然我们找到了invalid dimnames given for “dgCMatrix” objectError 的报错原因我们就可以针对counts来转化，data的部分我们在Seurat里面做。为什么不找到data的行名，赋值给data呢？这里留作思考题吧。坏。

我们开始改造原函数使之能够接受slots=’counts’ 这样的限定，于是我们找到 as.Seurat源代码中调用AssembleAssay的部分：

for (assay in names(x = assays)) {
    assay.objects[[assay]] <- AssembleAssay(
      assay = assay,
      file = x,
      slots =  assays[[assay]],  #这里强制对所有的slots进行转化，我们要让他接受传参
      verbose = verbose
    )
  }

除了在函数参数中加入slots = NULL,之外，这部分改为：

for (assay in names(x = assays)) {
    assay.objects[[assay]] <- AssembleAssay(
      assay = assay,
      file = x,
      slots =slots %||%  assays[[assay]],
      verbose = verbose
    )
  }

考虑到基本上要改动https://github.com/mojaveazure/seurat-disk/blob/master/R/LoadH5Seurat.R这个脚本的大部分函数，所以我们Fork一份seurat-disk 对应的我们改https://github.com/tuqiang2014/seurat-disk/blob/master/R/LoadH5Seurat.R

至于怎么改的，这里就略过了。如果遇到一模一样的问题，你可以安装这个改过了的。改完之后，我们卸载掉这个官方的seurat-disk ，安装自己改过的。

detach("package:SeuratDisk", unload = TRUE)
remove.packages('SeuratDisk')
remotes::install_github("tuqiang2014/seurat-disk")  # tuqiang2014就是我啦
library(SeuratDisk)
undebug(LoadH5Seurat) 
undebug( as.Seurat)
undebug(AssembleAssay) # 这里的提示忽略，不是同一个环境。

cellxgene <- LoadH5Seurat('some.processed.gzip.h5seurat',
                       assays ="RNA",slots='counts')

# 提示信息：
Validating h5Seurat file
Initializing RNA with counts
Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')
Adding feature-level metadata for RNA
Adding reduction umap
Adding cell embeddings for umap
Adding miscellaneous information for umap
Adding command information
Adding cell-level metadata

于是，一个标致的Seurat对象就呈现在眼前了，现在可以在Seurat里面尽情的分析啦。

cellxgene
An object of class Seurat 
XXXX  features across XXXXX samples within 1 assay 
Active assay: RNA (XXXX features, 0 variable features)
 1 dimensional reduction calculated: umap



library(ggplot2)
ggplot(DimPlot(cellxgene)$data,aes(umap_1        , umap_2   ,fill=ident)) + 
  geom_point(shape=21,colour="black",stroke=0.25,alpha=0.8) +
  DimPlot(cellxgene,label = T)$theme +
  theme_bw()+ NoLegend()

参考： 

[单细胞转录组数据分析||Seurat3.1教程:Interoperability between single-cell object formats](https://www.jianshu.com/p/396345566479)
[Convert AnnData to Seurat fails with raw h5ad](https://gitmemory.com/issue/satijalab/seurat/1021/479993549)
[conversion_vignette](https://satijalab.org/seurat/v3.1/conversion_vignette.html)
[convert-anndata RMD](https://github.com/mojaveazure/seurat-disk/blob/master/vignettes/convert-anndata.Rmd)
[Error: Cannot find feature names in this H5AD file #7](https://github.com/mojaveazure/seurat-disk/issues/7)