软件测评：百万级单细胞数据的Anndata和Seurat对象互转

生信菜鸟团

发布于 2024-03-25 17:02:24

2870

发布于 2024-03-25 17:02:24

最近终于搞完了博士预答辩，大修后送盲审，祝我好运~

随着单细胞相关研究成果的井喷式爆发，单细胞领域已进入百万级甚至千万级细胞量的时代。因此有不少R语言党（包括我）开始学习Python，使用Scanpy流程。但是，由于习惯了Seurat流程，有些时候需要把Anndata对象的单细胞数据转为Seurat对象，然后使用R语言进行一些分析。而最大的问题在于，如何丝滑的将Anndata对象的h5ad格式与Seurat对象相互转换。本文基于一个百万级的单细胞测试数据，对多种互转软件进行测评并总结。希望能够帮助到大家~

一. 测试数据

首先，测试数据来自这篇系统性红斑狼疮的PBMC单细胞测序文章《Single-cell RNA-seq reveals cell type–specific molecular and genetic associations to lupus》。数据储存于GEO数据库GSE174188，储存格式为h5ad，包含了约120w单细胞。点击下载之后，我们开始测评互转软件。

image-20240321174743755

下载后的h5ad数据解压后命名为GSE174188_CLUES1_adjusted.h5ad，约为10.6GB大小。

二. 互转测评

我收集了4种nndata对象的h5ad格式与Seurat对象互转方法，包括

R包SeuratDisk，这个是Seurat配套的算法（如果好用的话，我们就没必要折腾别的包了，所以可想而知，这个算法肯定不好使...）；
R包sceasy；
R包MuDataSeurat；
R包和Python包scDIOR，这个算法在我之前的帖子里经常用到，非常推荐，但是环境部署有点麻烦。

下面进行一一测评：

1. R包SeuratDisk

library(Seurat)
library(SeuratDisk)

# 1.1 h5ad转为Seurat
time.py2R = system.time({
  Convert('./Rawdata/GSE174188_CLUES1_adjusted.h5ad', "h5seurat",
          overwrite = TRUE,assay = "RNA")
})
print(time.py2R)

time.read = system.time({
  seurat_obj <- LoadH5Seurat("./Rawdata/GSE174188_CLUES1_adjusted.h5ad")
})
print(time.read)

image-20240321181135420

R语言下的转换运行奔溃了，没出结果...

这里再试试Seurat转为h5ad：

seurat_object = read_rds("./Outdata/Step1.seurat.data.Clean.rds")

# 1.2 Seurat转为h5ad
time.R2py = system.time({
  seu = DietSeurat(
    seurat_object,
    counts = TRUE, # so, raw counts save to adata.raw.X
    data = TRUE, # so, log1p counts save to adata.X
    scale.data = FALSE, # set to false, or else will save to adata.X
    features = rownames(seurat_object), # export all genes, not just top highly variable genes
    assays = "RNA",
    dimreducs = c("pca","umap"),
    graphs = c("RNA_nn", "RNA_snn"), # to RNA_nn -> distances, RNA_snn -> connectivities
    misc = TRUE
  )
  
  # step 2: factor to character, or else your factor will be number in adata
  i <- sapply(seu@meta.data, is.factor)
  seu@meta.data[i] <- lapply(seu@meta.data[i], as.character)
  
  # step 3: convert
  SaveH5Seurat(seu, filename = "./Outdata/SeuratDisk.h5seurat", overwrite = TRUE)
  Convert("./Outdata/SeuratDisk.h5seurat", "./Outdata/SeuratDisk.h5ad", assay="RNA", overwrite = TRUE)
})
print(time.R2py)

user system elapsed 577.649 32.573 611.104

使用SeuratDisk包进行Seurat转h5ad，百万级单细胞花了600多秒，速度比较慢。

2. R包sceasy

R包sceasy的安装有点麻烦，需要先在linux下安装一些软件：

conda install anndata==0.6.19 scipy==1.2.1 -c bioconda
conda install loompy -c biocond

然后在R语言中安装：

devtools::install_github("cellgeni/sceasy")
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("LoomExperiment", "SingleCellExperiment"))
install.packages('reticulate')

然后就可以运行sceasy了：

library(sceasy)
library(reticulate)
library(Seurat)

# 2.1 h5ad转为Seurat
time.py2R = system.time({
  sceasy::convertFormat(obj = "./Rawdata/GSE174188_CLUES1_adjusted.h5ad", from="anndata",
                        to="seurat",outFile = 'sceasy.rds')
})
print(time.py2R)

image-20240321194652273

似乎是单细胞量太大了，报错了。

试试Seurat转为h5ad：

# 2.2 Seurat转为h5ad
time.R2py = system.time({
  sceasy::convertFormat(seurat_object, from="seurat", to="anndata",
                        outFile='./Outdata/sceasy.h5ad')
})
print(time.R2py)

user system elapsed 243.205 22.359 267.202

使用sceasy包进行Seurat转h5ad，百万级单细胞花了267多秒，速度比较快。

3. R包MuDataSeurat

R包MuDataSeurat安装比较简单：

remotes::install_github("zqfang/MuDataSeurat", ref='dev', force = T)

然后开始转换：

library(MuDataSeurat)

# 3.1 h5ad转为Seurat
time.py2R = system.time({
  seurat.data = MuDataSeurat::ReadH5AD(file = "./Outdata/GSE174188_CLUES1_adjusted.h5ad")
})
print(time.py2R)

image-20240321202345088

报错了。

试试Seurat转为h5ad：

# 3.2 Seurat转为h5ad
time.R2py = system.time({
  MuDataSeurat::WriteH5AD(seurat_object, "MuDataSeurat.h5ad", assay="RNA")
})
print(time.R2py)

image-20240321202804306

还是报错。

4. R包dior和Python包scDIOR

这个算法在我之前的帖子里经常用到，安装比较麻烦，依赖包有版本的限制，否则就会报错。所以我的建议是设置一个scdiopy的Conda小环境：

conda create -n scdiopy python==3.8
conda activate scdiopy
# for python
pip install diopy
pip install scanpy==1.9.6 numpy==1.21.3 numba==0.57.1

然后在R语言中安装：

devtools::install_github('JiekaiLab/dior')

接下来，在R语言中测序转换：

library(Seurat)
library(dior)

# 4.1 h5ad转为Seurat
time.py2R = system.time({
  seurat.data <- read_h5ad(file = './Rawdata/GSE174188_CLUES1_adjusted.h5ad',
                           assay_name = 'RNA', 
                         target.object = 'seurat')
})
print(time.py2R)

user system elapsed 239.829 82.604 322.756

h5ad转Seurat花了322多秒。

# 4.2 Seurat转为h5ad
time.R2py = system.time({
  sceasy::convertFormat(seurat_object, from="seurat", to="anndata",
                        outFile='./Outdata/scDIOR.h5ad')
})
print(time.R2py)