關(guān)鍵字
Anndata對象轉(zhuǎn)成Seurat對象 h5文件讀寫 空間組格式轉(zhuǎn)換 已補(bǔ)充快速使用的函數(shù)整理版本，如果不想看細(xì)節(jié)可以直接看已整理好的版本。
適用背景
眾所周知，單細(xì)胞數(shù)據(jù)分析有兩大軟件：基于R語言的Seurat和基于Python的Scanpy，在平時的分析中常常需要把Seurat對象轉(zhuǎn)成Scanpy的Anndata對象，這已經(jīng)有比較成熟的流程了。但是，如果反過來把Anndata對象轉(zhuǎn)成Seurat對象 ，網(wǎng)上搜到的方案寥寥無幾，而且在本人親測之下均報錯無法成功實(shí)現(xiàn)。再加上我需要轉(zhuǎn)的是空間組對象，結(jié)構(gòu)比單細(xì)胞的更為復(fù)雜，只好自己想法從Anndata對象提取信息重新構(gòu)建出一個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格式的矩陣，可以運(yùn)行下面代碼：
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")
3 提取坐標(biāo)信息（UMAP坐標(biāo)或空間組坐標(biāo)）
#保存UMAP坐標(biāo)
cord=pd.DataFrame(data=ob1.obsm['X_umap'],index=ob1.obs_names,columns=['x','y'])
cord.to_csv('position_X_umap.tsv',sep="\t")
#保存空間坐標(biāo)
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】方式則直接運(yùn)行上面一行代碼即可得到矩陣的值，但之后要手動添加行名和列名，下面的代碼不適用，但如果用的【to_hdf】則運(yùn)行下面代碼即可。
#獲取矩陣值
mat <- mydata$block0_values
#添加行名
rownames(mat) <- mydata$axis0
#添加列名
colnames(mat) <- mydata$axis1
#轉(zhuǎn)成稀疏矩陣
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對象進(jìn)行簡單聚類后才能添加UMAP坐標(biāo)信息，可略過直接進(jìn)行常規(guī)分析
obj <- seob_cluster(obj)
obj@reductions$umap@cell.embeddings[,1]<-pos[,1]
obj@reductions$umap@cell.embeddings[,2]<-pos[,2]
如果只是單細(xì)胞數(shù)據(jù)，以上步驟已經(jīng)足夠了，但如果是重建一個Seurat空間組對象，其實(shí)也就是填補(bǔ)image的slice1槽，就需要運(yùn)行以下腳本。
library(dplyr)
library(data.table)
library(Matrix)
library(rjson)
#以下腳本參考了部分華大生命科學(xué)研究院的余浩師兄和鄒軒軒師姐寫的腳本，在此表示感謝
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
看一下新建的對象結(jié)構(gòu)，說明已經(jīng)是一個標(biāo)準(zhǔn)的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_"
可以往后面進(jìn)行分析了。
小結(jié)與補(bǔ)充
希望Scanpy或Seurat官方能出一下相關(guān)函數(shù)吧，涉及到空間組數(shù)據(jù)時，常規(guī)的轉(zhuǎn)換流程也容易報錯，有很多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
總結(jié)
以上是生活随笔 為你收集整理的从Scanpy的Anndata对象提取信息并转成Seurat对象（适用于空间组且涉及h5文件读写） 的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔 網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔 推薦給好友。
