其他
读取loom格式的单细胞文件
万事开头难,考虑到很多小伙伴在做单细胞公共数据分析的时候往往是在第一个步骤读取作者上传的表达量矩阵去构建seurat对象就各种屏蔽,非常有必要把18种单细胞数据格式文件都给大家梳理一下 。
现在我们来演示一下如何读取loom格式的单细胞文件,首先需要安装并且加载一些包:
library(hdf5r)
library(loomR)
library(LoomExperiment)
# remotes::install_github("aertslab/SCopeLoomR")
# remotes::install_local('mojaveazure-loomR-0.2.0-beta-54-g1eca16a.tar.gz')
library(SCopeLoomR)
值得注意的是,有一些包其实是在GitHub上面哦,如果你网络比较差,需要自己想办法解决,如果连包读无法安装,不妨试试看我们的**马拉松授课(直播一个月互动教学) ,可以看完我们从2000多个提问互动交流里面精选的200个问答!2021第二期_生信入门班_微信群答疑整理,以及 2021第二期_数据挖掘班_微信群答疑笔记
与十万人一起学生信,你值得拥有下面的学习班:
接下来我们以 GSE160756 ,为例子,进入https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE160756 可以看到,其数据集的7个样品,都是以loom格式文件分享给大家的。
我们的示例代码如下所示 ;
###### step1:导入数据 ######
path='GSE160756_RAW/'
samples=list.files(path )
samples
sceList = lapply(samples,function(pro){
# pro=samples[1]
print(pro)
folder=file.path(path,pro)
print(pro)
print(folder)
# 导入loom文件
lfile <- connect(filename =folder, mode = "r+")
ct <- as.data.frame(lfile[["matrix"]][,])
colnames(ct) <- lfile[["row_attrs/gene_names"]][] # 添加gene_name
rownames(ct) <- lfile[["col_attrs/cell_names"]][] # 添加cell_name
ct[1:5,1:5]
ct <- t(ct) # 注意转置
sce=CreateSeuratObject(counts = ct,
project = pro )
return(sce)
})
names(sceList)
samples
读取全部的loom文件后,成为了一个列表,里面是每个样品独立的seurat对象,需要合并,代码如下所示:
sce.all=merge(x=sceList[[1]],
y=sceList[ -1 ],
add.cell.ids = samples)
as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all@meta.data$orig.ident)
library(stringr)
phe=str_split(rownames(sce.all@meta.data),'_',simplify = T)
head(phe)
table(phe[,4])
table(phe[,3])
sce.all@meta.data$rep=phe[,4]
sce.all@meta.data$group=phe[,3]
sce.all@meta.data$orig.ident = paste0(
sce.all@meta.data$group,
sce.all@meta.data$rep
)
table(sce.all@meta.data$orig.ident)
因为是多个样品,通常是建议走harmnoy流程进行整合哦。
大家在下面的文章里面可以搜索到10x单细胞转录组数据的文章公布在geo数据库的链接:
如果你对单细胞数据分析还没有基础认知,可以看基础10讲:
最基础的往往是降维聚类分群,参考前面的例子:人人都能学会的单细胞聚类分群注释
如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程
看完记得顺手点个“在看”哦!
长按扫码可关注