GEOdatabase下载了一个.loom文件,想要用一下这个数据,折磨了一周,终于成功的将它转换成了seurat对象,现记录如下:
数据为GSE162183,地址:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE162183
关于loomR的简介:
随着单细胞数据量的增长,计算要求成指数增长,当数据量大于10万个细胞的时候,seurat包分析就显得非常有压力了,因为在实时内存中储存数据就变得非常困难,HDF5数据格式提供了高效的磁盘存储,而不是在内存中存储数据,这就将分析扩展到大规模数据集,甚至可以达到大于100万细胞的级别 ,Linnarson实验室开发了一种基于hdf5的数据结构,loom,可以方便地存储单细胞基因组数据集和元数据。他们还发布了一个名为loompy的Python API来与loom文件交互,而loomR为基于R的与loom交互的R包,详细的信息可以参考: https://satijalab.org/loomr/loomr_tutorial
安装并加载loomR包
devtools::install_github(repo = "hhoeflin/hdf5r")
devtools::install_github(repo = "mojaveazure/loomR", ref = "develop")
library(loomR)
loomR包基于R6对象与seurat基于R4对象不同,R6对象用$而不是@来取值
与loom建立连接
与将包含在其内部的所有数据加载到内存中的标准R对象不同,loom对象仅仅是到磁盘上的一个文件的连接,这使得可以在低内存消耗的情况下扩展到大量数据集。可以通过 connect 连接到一个现有的loom文件,使用loom::create从表达式矩阵创建您自己的文件,或者使用Convert从现有的Seurat对象创建一个loom文件。
>hc.2=connect('D:/Google/GSE162183_Raw_gene_counts_matrix_LoomFile.loom/GSE162183_Raw_gene_counts_matrix_LoomFile.loom',mode = 'r+',skip.validate = TRUE)
Warning message:
In initialize(...) :
Skipping validation step, some fields are not populated
Class: loom
Filename: D:\Google\GSE162183_Raw_gene_counts_matrix_LoomFile.loom\GSE162183_Raw_gene_counts_matrix_LoomFile.loom
Access type: H5F_ACC_RDWR
Attributes: last_modified, version
Listing:
name obj_type dataset.dims
attrs H5I_GROUP <NA>
col_attrs H5I_GROUP <NA>
col_graphs H5I_GROUP <NA>
layers H5I_GROUP <NA>
matrix H5I_DATASET 24234 x 19968
row_attrs H5I_GROUP <NA>
row_graphs H5I_GROUP <NA>
dataset.type_class
H5T_INTEGER
在没加参数skip.validate之前,我遇到了一个error(Error in validateLoom(object = self) :
There can only be 5 groups in the loom file: 'row_attrs', 'col_attrs', 'layers', 'row_graphs', 'col_graphs'),github上找到了解决方法,貌似是版本的问题,加上那个参数设置便可以了
但是我直接as.seurat之后,显示error如下:
hc.2=as.Seurat(hc.2)
Error in UseMethod(generic = "as.Seurat", object = x) :
"as.Seurat"没有适用于"c('loom', 'H5File', 'H5RefClass', 'R6')"目标对象的方法
仍然不行,观察.loom文件中的matrix即为创建seurat对象中的单细胞表达矩阵,因此考虑将matrix提取出来,再将barcode和feature提出来,即可创建seurat对象,如下所示:
psoriasis=hc.2[["matrix"]][,] #提取.loom文件的matrix
psoriasis=t(psoriasis) #发现那个matrix的gene和barcode颠倒了,换过来才符合seurat对象中的矩阵行为基因,列为barcode
dim(psoriasis)
[1] 19968 24234
gene=hc.2$row.attrs$Gene[] #提取基因名
barcode=hc.2$col.attrs$CellID[]# 提取barcode
length(gene)
[1] 19968
length(barcode)
[1] 24234
colnames(psoriasis)= barcode
row.names(psoriasis)= gene
dim(psoriasis)
[1] 19968 24234
####创建seurat对象####
psoriasis=CreateSeuratObject(counts = psoriasis,project = 'psoriasis',min.cells = 3, min.features = 200)
psoriasis
An object of class Seurat