读取loom格式的单细胞文件

2022-01-17 17:31:24 浏览数 (1)

万事开头难,考虑到很多小伙伴在做单细胞公共数据分析的时候往往是在第一个步骤读取作者上传的表达量矩阵去构建seurat对象就各种屏蔽,非常有必要把18种单细胞数据格式文件都给大家梳理一下 。

现在我们来演示一下如何读取loom格式的单细胞文件,首先需要安装并且加载一些包:

代码语言:javascript复制
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第二期_数据挖掘班_微信群答疑笔记

与十万人一起学生信,你值得拥有下面的学习班:

  • 数据挖掘(GEO,TCGA,单细胞)2021第11期(年度收官)
  • 生信入门课-2021第11期(年度收官)

接下来我们以 GSE160756 ,为例子,进入https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE160756 可以看到,其数据集的7个样品,都是以loom格式文件分享给大家的。

以loom格式文件分享给大家的

我们的示例代码如下所示 ;

代码语言:javascript复制

###### 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对象,需要合并,代码如下所示:

代码语言:javascript复制

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数据库的链接:

  • 我的课题只有一个10x样本肿么办?
  • 两个样品的10x单细胞转录组数据分析策略
  • 三个10X单细胞转录组样本CCA整合
  • 多个单细胞转录组样本的数据整合之CCA-Seurat包

如果你对单细胞数据分析还没有基础认知,可以看基础10讲:

  • 01. 上游分析流程
  • 02.课题多少个样品,测序数据量如何
  • 03. 过滤不合格细胞和基因(数据质控很重要)
  • 04. 过滤线粒体核糖体基因
  • 05. 去除细胞效应和基因效应
  • 06.单细胞转录组数据的降维聚类分群
  • 07.单细胞转录组数据处理之细胞亚群注释
  • 08.把拿到的亚群进行更细致的分群
  • 09.单细胞转录组数据处理之细胞亚群比例比较

最基础的往往是降维聚类分群,参考前面的例子:人人都能学会的单细胞聚类分群注释

0 人点赞