适用背景
昨天写的文章介绍了Anndata对象转成Seurat对象的流程，但是比较稀碎，所以这篇文章就补充了一个整理成函数的版本。
为什么要转格式
目前使用最为广泛的2个单细胞分析的软件包是基于R语言的Seurat和基于Python的Scanpy。Seurat的分析对象是一个S4 method的Seurat对象，S4 method可以简单理解为可以存储很多种格式的对象类型，不用过于深究，里面存了原始矩阵，各种处理过后的矩阵以及降维信息，样本信息等等，可以通过操作符@和$进行访问，如果还不熟悉其结构可以通过Rstudio使用Tab键查看其有关属性，也可以使用str函数快速查看其拓扑结构。而Scanpy的存储对象则为一个Anndata，其结构也类似Seurat对象，但其属性命名规则则与Seurat几乎完全不同，例如其meta.data信息存在obs里，矩阵存在属性X里。
之所以要转格式，是因为目前绝大多数单细胞分析工具都是基于这两种对象进行开发的，有些工具是R语言的，有些是Python语言的，所以常常需要转格式之后才能使用相关分析工具。
更新2023-06-29
偶然发现，SeuratDisk包能进行这种格式转换。教程网址
library(Seurat)
library(SeuratData)
library(SeuratDisk)
InstallData("pbmc3k")
data("pbmc3k.final")
#从Seurat转成h5ad
SaveH5Seurat(pbmc3k.final, filename = "pbmc3k.h5Seurat")
Convert("pbmc3k.h5Seurat", dest = "h5ad")
#从h5ad转成Seurat
Convert("pbmc3k_final.h5ad", dest = "h5seurat", overwrite = TRUE)
pbmc3k <- LoadH5Seurat("pbmc3k_final.h5seurat")
出现报错“Ambiguous assays”怎么办？
#添加assay参数
Convert("pbmc3k_final.h5ad", dest = "h5seurat", overwrite = TRUE,assay='Spatial')
ob1 <- LoadH5Seurat("pbmc3k_final.h5seurat",assays='Spatial')
Python 函数
我把函数写进一个function_anndata.py脚本里，这是脚本里的内容：
import os
import sys
import scanpy as sc
import anndata as ad
import numpy as np
import pandas as pd
import h5py
def extra_adata_info(infile,cord_name='spatial'):
ob1=sc.read(infile)
mat=pd.DataFrame(data=ob1.X.todense(),index=ob1.obs_names,columns=ob1.var_names)
mat.to_hdf("mat.h5","mat")
meta=pd.DataFrame(data=ob1.obs)
meta.to_csv('metadata.tsv',sep="\t")
cord=pd.DataFrame(data=ob1.obsm[cord_name],index=ob1.obs_names,columns=['x','y'])
cord.to_csv('position_'+cord_name+'.tsv',sep="\t")
然后在另一个脚本文件里加载这个脚本，网上收到加载脚本可以用import，但是需要先进入加载脚本所在路径，并且import的时候不用加.py后缀，然后再调用里面的函数extra_adata_info，需注意运行函数前需把目录转回输出目录：
import os
outdir=os.getcwd()
infile='tmp.h5ad'
cord_name='spatial'
os.chdir('/path/to/scripts/')
import function_anndata
from function_anndata import extra_adata_info
os.chdir(outdir)
extra_adata_info(infile,cord_name)
函数extra_adata_info只有两个参数，infile是Scanpy存储的Anndata对象，而cord_name是需保存坐标的slot名，默认是‘spatial’，因为我主要是想做空间组信息，如果要保留UMAP坐标，可以填‘X_umap’。
函数运行好后可以得到矩阵mat.h5，metadata信息metadata.tsv和坐标信息position_spatial.tsv共3个文件。
R函数
我把重建Seurat对象的函数rebulid_seob封装到一个脚本function_transfer_adata.R里：
library(rhdf5)
library(dplyr)
library(data.table)
library(Matrix)
library(rjson)
library(Seurat)
rebulid_seob <- function(fmat="mat.h5",cord_name='spatial',fmeta='metadata.tsv',fpos='position_spatial.tsv',assay='Spatial'){
mydata <- h5read(fmat,"mat")
mat <- mydata$block0_values
rownames(mat) <- mydata$axis0
colnames(mat) <- mydata$axis1
mat <- Matrix(mat, sparse = TRUE)
meta <- read.table(fmeta,sep="\t",header=T,row.names=1)
pos <- read.table(fpos,sep="\t",header=T,row.names=1)
obj <- CreateSeuratObject(mat, project = assay, assay = assay, min.cells=5, min.features=5,meta.data=meta)
##以下脚本将重建空间组对象
if (assay=='Spatial') {
tissue_lowres_image <- matrix(1, max(pos$y), max(pos$x))
tissue_positions_list <- data.frame(row.names = colnames(obj),
tissue = 1,
row = pos$y, col = pos$x,
imagerow = pos$y, imagecol = pos$x)
scalefactors_json <- toJSON(list(fiducial_diameter_fullres = 1,
tissue_hires_scalef = 1,
tissue_lowres_scalef = 1))
seurat_spatialObj <- obj
generate_spatialObj <- function(image, scale.factors, tissue.positions, filter.matrix = TRUE)
{
if (filter.matrix) {
tissue.positions <- tissue.positions[which(tissue.positions$tissue == 1), , drop = FALSE]
}
unnormalized.radius <- scale.factors$fiducial_diameter_fullres * scale.factors$tissue_lowres_scalef
spot.radius <- unnormalized.radius / max(dim(image))
return(new(Class = 'VisiumV1',
image = image,
scale.factors = scalefactors(spot = scale.factors$tissue_hires_scalef,
fiducial = scale.factors$fiducial_diameter_fullres,
hires = scale.factors$tissue_hires_scalef,
lowres = scale.factors$tissue_lowres_scalef),
coordinates = tissue.positions,
spot.radius = spot.radius))
}
spatialObj <- generate_spatialObj(image = tissue_lowres_image,
scale.factors = fromJSON(scalefactors_json),
tissue.positions = tissue_positions_list)
spatialObj <- spatialObj[Cells(seurat_spatialObj)]
DefaultAssay(spatialObj) <- 'Spatial'
seurat_spatialObj[['slice1']] <- spatialObj
return(seurat_spatialObj)
}else{
return(obj)
}
}
函数rebulid_seob需要传入上面Python函数生成的3个文件，对应三个参数：fmat，fmeta，fpos，如果是在同一个目录下，则可以不用传入这3个参数。另外还有2个参数：cord_name与之前一致是表示坐标的slot名，assay决定是返回什么类型的对象，默认是‘Spatial’，会返回空间组数据对象，如果不是‘Spatial’则返回普通的单细胞对象。
函数使用示例：
返回单细胞对象
que <- rebulid_seob(assay="RNA")
等价于
que <- rebulid_seob(fmat="mat.h5",cord_name='spatial',fmeta='metadata.tsv',fpos='position_spatial.tsv',assay='RNA')返回空间组对象
que <- rebulid_seob()
等价于
que <- rebulid_seob(fmat="mat.h5",cord_name='spatial',fmeta='metadata.tsv',fpos='position_spatial.tsv',assay='Spatial')
小结与补充
重构的Seurat对象是没有降维和聚类等信息的，需要根据具体需求进行后续分析。
