關(guān)鍵字

• Anndata對(duì)象轉(zhuǎn)成Seurat對(duì)象
• h5文件讀寫
• 空間組格式轉(zhuǎn)換

已補(bǔ)充快速使用的函數(shù)整理版本，如果不想看細(xì)節(jié)可以直接看已整理好的版本。

適用背景

眾所周知鸡捐，單細(xì)胞數(shù)據(jù)分析有兩大軟件：基于R語(yǔ)言的Seurat和基于Python的Scanpy嗅绸，在平時(shí)的分析中常常需要把Seurat對(duì)象轉(zhuǎn)成Scanpy的Anndata對(duì)象，這已經(jīng)有比較成熟的流程了断国。但是贤姆，如果反過(guò)來(lái)把Anndata對(duì)象轉(zhuǎn)成Seurat對(duì)象，網(wǎng)上搜到的方案寥寥無(wú)幾稳衬，而且在本人親測(cè)之下均報(bào)錯(cuò)無(wú)法成功實(shí)現(xiàn)霞捡。再加上我需要轉(zhuǎn)的是空間組對(duì)象，結(jié)構(gòu)比單細(xì)胞的更為復(fù)雜薄疚，只好自己想法從Anndata對(duì)象提取信息重新構(gòu)建出一個(gè)Seurat對(duì)象了碧信。
這個(gè)步驟主要分為2步：

步驟一 從Scanpy的Anndata對(duì)象中提取信息

• 1提取矩陣

import os
import sys
import scanpy as sc
import anndata as ad
import numpy as np
import pandas as pd
import h5py

ob1=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=pd.read_hdf("mat.h5",key="mat")
上面的腳本是我測(cè)試過(guò)比較好的保存矩陣的方案街夭，下面代碼塊則是最初的版本砰碴，但是在R里面讀入之后會(huì)缺失行名和列名，所以還要額外保存行名與列名板丽，之后加上呈枉，特別麻煩，所以采用上面的代碼塊更為簡(jiǎ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保存矩陣猖辫？
理論上是可以轉(zhuǎn)成數(shù)據(jù)框之后保存為csv或tsv文件，但這樣保存很慢乃正，讀取數(shù)據(jù)也很慢住册，因此存為h5文件更加方便。但對(duì)于小文件瓮具，h5文件反而占的空間更大荧飞，但h5文件應(yīng)該是壓縮過(guò)的凡人，就很奇怪，不太懂其中原理叹阔，但如果是大文件則能有效進(jìn)行壓縮使得所占空間更小挠轴。

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') 

• 2 提取metadata信息

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對(duì)象了耳幢。

步驟二 重建Seurat對(duì)象

• 1 讀取上面提取的信息

library(rhdf5)
library(Matrix)
mydata <- h5read("mat.h5","mat")
#如果在Python保存h5文件是用【with h5py.File】方式則直接運(yùn)行上面一行代碼即可得到矩陣的值岸晦，但之后要手動(dòng)添加行名和列名，下面的代碼不適用睛藻，但如果用的【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)

• 2 重建單細(xì)胞Seurat對(duì)象

obj <- CreateSeuratObject(mat,assay='Spatial',meta.data=meta)

需要對(duì)seurat對(duì)象進(jìn)行簡(jiǎn)單聚類后才能添加UMAP坐標(biāo)信息，可略過(guò)直接進(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)足夠了冈在，但如果是重建一個(gè)Seurat空間組對(duì)象，其實(shí)也就是填補(bǔ)image的slice1槽按摘，就需要運(yùn)行以下腳本包券。

• 3 重建空間組Seurat對(duì)象

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@counts

seurat_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_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

看一下新建的對(duì)象結(jié)構(gòu)炫贤，說(shuō)明已經(jīng)是一個(gè)標(biāo)準(zhǔn)的Seurat空間組對(duì)象了：

> 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ù)時(shí)兰珍，常規(guī)的轉(zhuǎn)換流程也容易報(bào)錯(cuò)侍郭，有很多bugs。

以下為艱辛的報(bào)錯(cuò)心路歷程俩垃，可借鑒可忽略……

error1
網(wǎng)上一查励幼，基本都是以下教程，很簡(jiǎn)單嘛口柳，然而……

> library(Seurat)
library(SeuratDisk)

Attaching SeuratObject
> library(SeuratDisk)
Registered S3 method overwritten by 'cli':
  method     from
  print.boxx spatstat.geom
Registered S3 method overwritten by 'SeuratDisk':
  method            from
  as.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 file
        class: HDF5
        major: File accessibility
        minor: Unable to open file

    error #001: H5VLcallback.c in H5VL_file_open(): line 3500: open failed
        class: HDF5
        major: Virtual Object Layer
        minor: Can't open object

    error #002: H5VLcallback.c in H5VL__file_open(): line 3465: open failed
        class: HDF5
        major: Virtual Object Layer
        minor: Can't open object

    error #003: H5VLnative_file.c in H5VL__native_file_open(): line 100: unable to open file
        class: HDF5
        major: File accessibility
        minor: Unable to open file

    error #004: H5Fint.c in H5F_open(): line 1622: unable to lock the file
        class: HDF5
        major: File accessibility
        minor: Unable to open file

    error #005: H5FD.c in H5FD_lock(): line 1675: driver lock request failed
        class: HDF5
        major: Virtual File Layer

error2
修復(fù)上面的bug后（export HDF5_USE_FILE_LOCKING=FALSE）苹粟，還是報(bào)錯(cuò)

>  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
根據(jù)某個(gè)教程，可以直接在Rstudio里使用Python跃闹，加載半天后直接報(bào)錯(cuò)嵌削，可能是因?yàn)槲矣玫氖羌骸?/p>

> 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
嘗試在Python導(dǎo)出矩陣后再用R讀入，R直接崩潰退出……

> 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

erro n++ ……