ixxmu
commented
6 months ago 软件测评：百万级单细胞数据的Anndata和Seurat对象互转 by 生信菜鸟团
这里做一个勘误，dior包并不提供h5ad的转换方法，dior包可以将seurat转为h5文件然后用python读取，但是dior包的seurat转为h5不是特别快：dior::seurat_write_h5(seurat.data,file='./Outdata/scdata.h5',object.type = 'seurat')
因此，旧文文末给的代码实际是我总结的最佳读取和转存策略，也就是组合不同的R包，达到最快的读和存的目的，总结来说：用dior包快速读h5ad为seurat，用sceasy将seurat快速转为h5ad文件。
虽然旧文的小错误不影响代码的运行，但是仍给读者带来了一些误解，在此非常抱歉！
勘误版本如下：
随着单细胞相关研究成果的井喷式爆发，单细胞领域已进入百万级甚至千万级细胞量的时代。因此有不少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单细胞。点击下载之后，我们开始测评互转软件。
下载后的h5ad数据解压后命名为GSE174188_CLUES1_adjusted.h5ad，约为10.6GB大小。
二. 互转测评
我收集了XX种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)
R语言下的转换运行崩溃了，没出结果...
这里再试试Seurat转为h5ad：
seurat_object = qread("./Outdata/Step1.seurat.data.Clean.qs")
# 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，百万级单细胞花了611.104多秒，速度比较慢。
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)
似乎是单细胞量太大了，报错了。
试试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)
报错了。
试试Seurat转为h5ad：
# 3.2 Seurat转为h5ad
time.R2py = system.time({
MuDataSeurat::WriteH5AD(seurat_object, "MuDataSeurat.h5ad", assay="RNA")
})
print(time.R2py)
还是报错。
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转为h5
time.R2py = system.time({
dior::seurat_write_h5(seurat.data,file='./Outdata/scdata.h5',object.type = 'seurat')
})
print(time.R2py)
dior包的Seurat转h5运行不是特别快。
因此，如果需要将Seurat转h5ad数据用于python环境读取的话，非常推荐搭配sceasy包：
# 4.2 Seurat转为h5ad
time.R2py = system.time({
sceasy::convertFormat(seurat_object, from="seurat", to="anndata",
outFile='./Outdata/sceasy.h5ad')
})
print(time.R2py)
三. 总结
如果需要将百万级细胞数量的数据从Anndata/h5ad格式转换为Seurat对象，我非常推荐使用R包dior。它的优点是运行速度快，数据兼容性强；缺点是依赖包有版本限制，容易报错。
而对于将Seurat对象转换为h5ad格式，R包sceasy是一个不错的选择，转换速度非常快；相对而言，R包SeuratDisk的速度较慢，不推荐使用。
https://mp.weixin.qq.com/s/8fwJSc9Dnp8h_Suv76oXVA