答读者问专栏（答读者问43 怎么拓展思维的深度?吴军）

最近遇到很多朋友提到Scanpy和Seurat对象的转换问题，而大部分问题都在于如何把Scanpy的h5ad文件转为Seurat对象。比如：

问题当时在群里我提出了我自己的解决方案：在R中构建seurat对象，转为h5ad对象py中跑流程有需要时，把metadata和降维信息拿回来，加入之前的seurat对象收到群主邀请，为了帮助更多朋友解决这一问题，我把我自己所用的方法进一步总结成这一篇笔记，希望对大家有用。

1. Seurat对象转h5ad文件这一步通常问题不大，因为网上的解决方案已经有很多，包括MuDataSeurat、SeuratDisk等，就不列举了，大家有兴趣可以网上很容易搜索到这里我给出我自己所用的代码，。

同时我们分享一下如何在服务器的R中实现并行计算（并不清楚Windows是否能实现），以提高效率：首先我们构建Seurat对象# 清空环境变量和释放内存rm(list = ls())gc()# 这里是因为Rscripts需要，这行代码识别的是Rproj文件所在的路径，并将路径设置到此处

# 如果你已经通过Rproj打开，那么可以忽略这行setwd(here::here())# 在这里我们采用的数据集是包含了多个样本数据的，因此采用并行计算的方法加速# 需要注意的是，如果你的是共享服务器，我建议你把核数控制在10个以内，避免影响其他朋友

# 加载R包library(parallel)detectCores()# [1] 128# 这里初始化10个核cl <- makeCluster(10)# 构建样品名的向量# 在这里，把matrix/features/barcodes三个文件放在一个文件夹，并以样品名命名

# 而这些文件夹都统一放在data/human_EC_raw/下，因此用list.file()可以获取样品名path = data/human_EC_raw/sample = list.files(path) 

# 开启并行计算# 并行计算需要初始化，因此在此前所有的R包或变量都将不存在于新的核中# 而在运行结束后，所有的变量也将清空# 因此，我们需要在运行代码中重新加载R包和赋予变量，并及时返回结果# 与lappy一样，我们通过parLapply将seurat对象组成一个list

# 第一个参数是我们初始化的核# 第二个参数是我们输入的向量(样品名)# 并通过function定义我们的流程sc_list = parLapply(cl, sample, function(i){# 这里推荐用绝对路径

  path = /home/data/gz123/JupyterHub/PAEC-scRNA/data/human_EC_raw/  file_path <- paste0(path, i)# 读入矩阵

library(Matrix)library(readr)  counts = readMM(paste(file_path, "matrix.mtx.gz", sep = "/"))  barcodes = read_tsv(paste(file_path, 

"barcodes.tsv.gz", sep = "/"), col_names = F)  features = read_tsv(paste(file_path, "features.tsv.gz"

, sep = "/"), col_names = F)# 在这里我希望在barcodes上标记我的样品来源# 当你的数据量非常大的时候，或者数据来源非常多的时候，有可能会出现barcodes重复的现象，虽然seurat本身会处理这类问题，但是他的方式并不优雅

# 在这里我们可以预防这类情况出现，保证每个细胞的barocdes唯一性  colnames(counts) = paste(i, barcodes$X1, sep = "-")  rownames(counts) = features$X2

  counts = counts[!duplicated(rownames(counts)), ]# 上述步骤你也可以通过Read10X()来实现# 构建seurat对象library(Seurat)

library(SeuratObject)  scobj = CreateSeuratObject(counts = counts,                             project = i, 

# project写上样品名                             min.cells = 10,                             min.feature = 

500# 这里你可以根据实际情况调整                            )# 在这里我们把数据集来源也加上# 如果你有其他信息也可以一并添加  scobj@meta.data$source

 = GSE123456# 下面的代码是根据样品名来添加分组信息，如果你没有这个需求，可以忽略这里# 判断是否为Cont  iscontrol = grepl("^Cont", i)if (grepl(

"^Cont", i)) {    scobj@meta.data$group = Healthy  }else{    scobj@meta.data$group = Disease  }# 返回seurat结果至sc_list

return(scobj)})# 归还核和内存给系统stopCluster(cl)# 给seurat对象命名，防止错乱names(sc_list) = samplenames(sc_list)# 合并seu <- base::Reduce(f = merge, x = sc_list)

# 保存，这里推荐qs这个包，又快又好，占空间又小qs::qsave(seu.backup, "data/raw_backup.qs")# 读入代码是：# seu = qs::qread("data/raw_backup.qs")

接下来我们把Seurat对象保存成h5ad文件# remotes::install_github("zqfang/MuDataSeurat", ref=dev, force = T)# 单模态数据MuDataSeurat::WriteH5AD(seu, 

"data/seurat_to_scanpy.h5ad", assay="RNA")这里需要补充一下，通常而言笔者本人更喜欢在Python中跑流程，主要是因为超过10万细胞后，R就会显得特别慢通常来说我构建好了。

Seurat对象就会马上转到Scanpy中另外，当你的Seurat对象元素越少，出错的概率越低！（个人经验）那有人就要问了，为什么不直接从Scanpy开始，主要有两点原因，一是我更熟悉R，很多情况处理起来更顺手，希望在R这边保留一个备份；另一个原因是我可能会回到R中画图。

2. 从Scanpy回到Seurat这通常是令人头疼的问题，笔者浏览了很多帖子，尝试了各种方式，网上所谓的一键转换都是条件限制很多的直到我的一位生信专业的老师点醒了我：我认为你并不需要把这么多的东西全部转换过来，你的anndata中最重要的结果，无非就是矩阵、metadata以及降维数据。

至于差异分析那些结果，你把它保存成csv不就得了，反正你在R里面也是这样做的确实！如果我要画图，我只需要把坐标拿出来就行了！费这么大劲干嘛呢？于是我转换了思路，来看看我是如何实现的假设我们现在已经在Python中完成了分析，并保存成了文件

04_annotation.h5ad首先我们使用reticulate包在R中使用Python环境，这是为了完整地读入anndata对象的信息如果你不熟悉reticulate包的使用，可以自行搜索相关文档

# 使用reticulatelibrary(reticulate)# 导入scanpysc <- import(scanpy)# 读入h5ad文件scvi_h5ad <- sc$read_h5ad(data/GSE136831/04_annotation.h5ad

)scvi_h5ad## 在这里我只取需要的信息## 1.样本信息meta = scvi_h5ad$obs # obs储存细胞信息，对应seuratobj@meta.datavar = scvi_h5ad$var  

# var储存基因信息，这是由于anndata提取的矩阵没有行列名，所以这里是为了提取基因名## 2.scvi坐标X_umap_scVI <- scvi_h5ad$obsm[X_umap_scVI]rownames(X_umap_scVI) = rownames(meta)

colnames(X_umap_scVI) = c(scVI_UMAP-1, scVI_UMAP-2)head(X_umap_scVI)X_umap_Harmony <- scvi_h5ad$obsm[

X_umap_Harmony]rownames(X_umap_Harmony) = rownames(meta)colnames(X_umap_Harmony) = c(Harmony_UMAP-1, 

Harmony_UMAP-2)head(X_umap_Harmony)# 提取矩阵，这里需要行列转置，因为anndata矩阵是 cells x genes，而seurat是genes x cellsfilter_counts = t(scvi_h5ad$layers[

counts])dim(filter_counts)# 提取的矩阵缺少行列名colnames(filter_counts) = rownames(meta)# 这里的矩阵是只保留了高可变基因rownames(filter_counts) = rownames(var)[var$highly_variable_features]

filter_counts[1:5,1:3]# 如果你希望拿到原始矩阵，你可以从这里提取: adata$raw$X (前提是你在Scanpy中保存了原始矩阵)# 对应的，你应该从adata$raw$var提取基因名

# 如果你对anndata对象不熟悉，可以忽略这几行注释# 重新构建seurat对象scobj <- CreateSeuratObject(counts = filter_counts,                            meta.data = meta)

scobj然后我们简单走一个流程scobj <- NormalizeData(scobj)scobj <- FindVariableFeatures(scobj, selection.method = 

"vst", nfeatures = 1000)scobj <- ScaleData(scobj, features = VariableFeatures(object = scobj))现在我们把降维坐标信息添加

# 添加 reductionsscobj[[X_umap_scVI]] <- CreateDimReducObject(embeddings = X_umap_scVI)scobj[[X_umap_Harmony

]] <- CreateDimReducObject(embeddings = X_umap_Harmony)# 使用qs包保存qs::qsave(scobj, file = data/05_rebuilt_seurat.qs

)我们可以简单画个图来验证一下是否成功library(SCP)CellDimPlot(  srt = scobj,   group.by = c("celltype"),   reduction = "X_umap_scVI"

,  theme_use = "theme_blank",  xlab = "UMAP-1",  ylab = "UMAP-2",  label = F,  label.size = 3.5)# 这里只是给了代码，实际上图并不是这个图，我没有保存 T^T

scVI_umap往期答疑回顾答读者问(9)：Seurat报错object ‘CsparseMatrix_validate‘ not found答读者问(8)：如何批量查询marker基因(对应的蛋白)会不会在膜上表达？

答读者问(7)：关于doublet答读者问(6)：单细胞TPM矩阵如何分析？答读者问(5)：提取monocle2的拟时序/坐标重新画图答读者问(4)：inferCNV的几个问题如何加入讨论群最后，如果你想加入我们的讨论群，可以公众号后台回复：

加群