關鍵字
Anndata對象轉成Seurat對象 h5文件讀寫 空間組格式轉換 已補充快速使用的函數(shù)整理版本，如果不想看細節(jié)可以直接看已整理好的版本。
適用背景
眾所周知，單細胞數(shù)據(jù)分析有兩大軟件：基于R語言的Seurat和基于Python的Scanpy，在平時的分析中常常需要把Seurat對象轉成Scanpy的Anndata對象，這已經有比較成熟的流程了。但是，如果反過來把Anndata對象轉成Seurat對象 ，網上搜到的方案寥寥無幾，而且在本人親測之下均報錯無法成功實現(xiàn)。再加上我需要轉的是空間組對象，結構比單細胞的更為復雜，只好自己想法從Anndata對象提取信息重新構建出一個Seurat對象了。
步驟一 從Scanpy的Anndata對象中提取信息
import os
import sys
import scanpy as sc
import anndata as ad
import numpy as np
import pandas as pd
import h5pyob1=sc.read('tmp.h5ad')
mat=pd.DataFrame(data=ob1.X.todense(),index=ob1.obs_names,columns=ob1.var_names)
mat.to_hdf("mat.h5","mat")
如果要在Python里重新讀取h5格式的矩陣，可以運行下面代碼：
mat=ob1.X.todense().T
with h5py.File('mat.h5','w') as f1:f1.create_dataset("mat",data=mat)##在Python里讀取h5格式矩陣的方法
with h5py.File('mat.h5','r') as f1:mat=f1['mat'][:]
為什么用h5保存矩陣？
pd.DataFrame(data=ob1.X.todense().T, index=ob1.var_names,columns=ob1.obs_names).to_csv('raw_mat.csv',sep="\t",float_format='%.0f')
meta=pd.DataFrame(data=ob1.obs)
meta.to_csv('metadata.tsv',sep="\t")
#保存UMAP坐標
cord=pd.DataFrame(data=ob1.obsm['X_umap'],index=ob1.obs_names,columns=['x','y'])
cord.to_csv('position_X_umap.tsv',sep="\t")
#保存空間坐標
cord_name='spatial'
cord=pd.DataFrame(data=ob1.obsm[cord_name],index=ob1.obs_names,columns=['x','y'])
cord.to_csv('position_'+cord_name+'.tsv',sep="\t")
提取完以上信息之后，就可以在R里面重建Seurat對象了。
步驟二 重建Seurat對象
library(rhdf5)
library(Matrix)
mydata <- h5read("mat.h5","mat")
#如果在Python保存h5文件是用【with h5py.File】方式則直接運行上面一行代碼即可得到矩陣的值，但之后要手動添加行名和列名，下面的代碼不適用，但如果用的【to_hdf】則運行下面代碼即可。
#獲取矩陣值
mat <- mydata$block0_values
#添加行名
rownames(mat) <- mydata$axis0
#添加列名
colnames(mat) <- mydata$axis1
#轉成稀疏矩陣
mat <- Matrix(mat, sparse = TRUE)
#讀取metadata信息
cord_name='X_umap'
meta <- read.table('metadata.tsv',sep="\t",header=T,row.names=1)
pos <- read.table(paste0('position_',cord_name,'.tsv'),sep="\t",header=T,row.names=1)
obj <- CreateSeuratObject(mat,assay='Spatial',meta.data=meta)
需要對seurat對象進行簡單聚類后才能添加UMAP坐標信息，可略過直接進行常規(guī)分析
obj <- seob_cluster(obj)
obj@reductions$umap@cell.embeddings[,1]<-pos[,1]
obj@reductions$umap@cell.embeddings[,2]<-pos[,2]
如果只是單細胞數(shù)據(jù)，以上步驟已經足夠了，但如果是重建一個Seurat空間組對象，其實也就是填補image的slice1槽，就需要運行以下腳本。
library(dplyr)
library(data.table)
library(Matrix)
library(rjson)
#以下腳本參考了部分華大生命科學研究院的余浩師兄和鄒軒軒師姐寫的腳本，在此表示感謝
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))
mat <- obj@assays$Spatial@countsseurat_spatialObj <- CreateSeuratObject(mat, project = 'Spatial', assay = 'Spatial',min.cells=5, min.features=5)
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_scalefspot.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
看一下新建的對象結構，說明已經是一個標準的Seurat空間組對象了：
> str(seurat_spatialObj@images$slice1)
Formal class 'VisiumV1' [package "Seurat"] with 6 slots..@ image : num [1:6, 1:9] 1 1 1 1 1 1 1 1 1 1 .....@ scale.factors:List of 4.. ..$ spot : num 1.. ..$ fiducial: num 1.. ..$ hires : num 1.. ..$ lowres : num 1.. ..- attr(*, "class")= chr "scalefactors"..@ coordinates :'data.frame': 71668 obs. of 5 variables:.. ..$ tissue : num [1:71668] 1 1 1 1 1 1 1 1 1 1 ..... ..$ row : num [1:71668] -2.904 -2.379 -0.798 -2.162 -1.436 ..... ..$ col : num [1:71668] 5.84 5.58 8.63 7.08 7.99 ..... ..$ imagerow: num [1:71668] -2.904 -2.379 -0.798 -2.162 -1.436 ..... ..$ imagecol: num [1:71668] 5.84 5.58 8.63
7.08 7.99 ... ..@ spot.radius : num 0.111 ..@ assay : chr "Spatial" ..@ key : chr "slice1_"
可以往后面進行分析了。
小結與補充
希望Scanpy或Seurat官方能出一下相關函數(shù)吧，涉及到空間組數(shù)據(jù)時，常規(guī)的轉換流程也容易報錯，有很多bugs。
以下為艱辛的報錯心路歷程，可借鑒可忽略……
error1
> library(Seurat)
library(SeuratDisk)Attaching SeuratObject
> library(SeuratDisk)
Registered S3 method overwritten by 'cli':method fromprint.boxx spatstat.geom
Registered S3 method overwritten by 'SeuratDisk':method fromas.sparse.H5Group Seurat
> Convert("cell_adata.h5ad",'h5SeuratWarning: Unknown file type: h5ad
Error in H5File.open(filename, mode, file_create_pl, file_access_pl) :HDF5-API Errors:error #000: H5F.c in H5Fopen(): line 793: unable to open fileclass: HDF5major: File accessibilityminor: Unable to open fileerror #001: H5VLcallback.c in H5VL_file_open(): line 3500: open failedclass: HDF5major: Virtual Object Layerminor: Can't open objecterror #002: H5VLcallback.c in H5VL__file_open(): line 3465: open failedclass: HDF5major: Virtual Object Layerminor: Can't open objecterror #003: H5VLnative_file.c in H5VL__native_file_open(): line 100: unable to open fileclass: HDF5major: File accessibilityminor: Unable to open fileerror #004: H5Fint.c in H5F_open(): line 1622: unable to lock the fileclass: HDF5major: File accessibilityminor: Unable to open fileerror #005: H5FD.c in H5FD_lock(): line 1675: driver lock request failedclass: HDF5major: Virtual File Layer
error2
> Convert("cell_adata.h5ad",'h5SeuraWarning: Unknown file type: h5ad
Creating h5Seurat file for version 3.1.5.9900
Adding X as data
Adding X as counts
Adding meta.features from var
Adding bbox as cell embeddings for bbox
Adding contour as cell embeddings for contour
Error: Not a matrix dataset
error3
> library(reticulate)
> scanpy <- import("scanpy")
*** Error in `/share/app/R/4.0.2/lib64/R/bin/exec/R': free(): invalid pointer: 0x000000000aa5f3b8 ***
======= Backtrace: =========
/usr/lib64/libc.so.6(+0x81329)[0x7fb79efc4329]
error4
> fmat <- npyLoad("fmat.npy")
> fmat <- npyLoad("fmat.npy")*** caught segfault ***
address 0x7f3bf5770000, cause 'invalid permissions'Traceback:1: .External(list(name = "InternalFunction_invoke", address = <pointer: 0x421bb10>, dll = list(name = "Rcpp", path = "/jdfssz1/software/R4.0/lib/Rcpp/libs/Rcpp.so", dynamicLookup = TRUE, handle = <pointer: 0x957cc10>, info = <pointer: 0x1d5ba90>), numParameters = -1L), <pointer: 0x140526f0>, filename, type, dotranspose)2: npyLoad("/jdfssz3/fmat.npy")Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection: 3
rm: cannot remove ‘/hwfssz5/tmp/RtmpRSe1IF’: Directory not empty
總結
以上是生活随笔 為你收集整理的从Scanpy的Anndata对象提取信息并转成Seurat对象（适用于空间组且涉及h5文件读写） 的全部內容，希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔 網站內容還不錯，歡迎將生活随笔 推薦給好友。
