单细胞分析实录(14): 细胞类型注释的另一种思路 — CellID
之前写过一篇细胞类型注释的推文:单细胞分析实录(7): 差异表达分析/细胞类型注释,里面介绍了细胞类型注释的两种思路。在谈到自动注释时,我们基本都会用SingleR这个软件,类似的软件还有很多,原理上可能五花八门,但表现不一定有SingleR好。很多软件在做注释时会基于某个细胞和参考细胞表达谱的相似性,今天看到的这个软件使用的是另一种思路:将某个细胞的特征基因集与表征细胞类型的参考基因集做富集分析,当在某个细胞类型的基因集上显著富集时,就定义为这个细胞类型(和GO富集很像)。所以这种方法的关键点是:找单个细胞的特征基因集,以及找细胞类型的参考基因集。
CellID发表在NBT上面:Gene signature extraction and cell identity recognition at the single-cell level with Cell-ID。解决的主要是找单个细胞的特征基因集这个问题。
软件的官方文档:https://bioconductor.org/packages/devel/bioc/vignettes/CelliD/inst/doc/BioconductorVignette.html 有详细的演示,并且在选取细胞类型参考基因集时,使用的是网站PanglaoDB:https://panglaodb.se/index.html 这个网站有人和小鼠各种细胞类型的marker,这是下载链接:https://panglaodb.se/markers/PanglaoDB_markers_27_Mar_2020.tsv.gz 比我之前整理的详细得多。
题外话:在单细胞分析实录(7)发布之后,公众号后台断断续续收到几十位读者询问marker清单的信息,希望这个更全面的清单能帮助到你。
下面我简单演示一下CellID的用法:
library(CellID) library(Seurat) library(tidyverse) #获取原始表达矩阵 BaronMatrix <- readRDS(url("https://storage.googleapis.com/cellid-cbl/BaronMatrix.rds")) #仅考虑编码蛋白质的基因 data("HgProteinCodingGenes") BaronMatrixProt <- BaronMatrix[rownames(BaronMatrix) %in% HgProteinCodingGenes,] #这几步类似Seurat的标准流程 Baron <- CreateSeuratObject(counts = BaronMatrixProt, project = "Baron", min.cells = 5) Baron <- NormalizeData(Baron) Baron <- ScaleData(Baron, features = rownames(Baron)) Baron <- RunMCA(Baron) #该软件的主要分析函数,将细胞和基因同时降维到一个空间,离细胞近的基因被定义为细胞的特征基因——令人窒息的操作 #下载参考基因集 panglao <- read_tsv("https://panglaodb.se/markers/PanglaoDB_markers_27_Mar_2020.tsv.gz") #根据器官过滤,示例数据就是胰腺的 panglao_pancreas <- panglao %>% filter(organ == "Pancreas") #物种包含人 panglao_pancreas <- panglao_pancreas %>% filter(str_detect(species,"Hs")) #下面两步会得到一个列表,列表的每一个元素是由基因集构成的字符向量,且每个元素被命名为对应的细胞类型的名字 panglao_pancreas <- panglao_pancreas %>% group_by(`cell type`) %>% summarise(geneset = list(`official gene symbol`)) pancreas_gs <- setNames(panglao_pancreas$geneset, panglao_pancreas$`cell type`)
#富集分析,用的超几何检验 HGT_pancreas_gs <- RunCellHGT(Baron, pathways = pancreas_gs, dims = 1:50, n.features = 200) #每个细胞选200个特征基因
富集分析之后会返回一个矩阵,矩阵中的每一个值表示p值在校正之后,求-log10,即
-log10 corrected p-value,当然这个值越大越好。
> HGT_pancreas_gs[1:4,1:4] 4 x 4 sparse Matrix of class "dgCMatrix" human1_lib1.final_cell_0001 human1_lib1.final_cell_0002 Acinar cells 60.6897442 63.2674918 Alpha cells 0.1241303 0.1241303 Beta cells . . Delta cells . . human1_lib1.final_cell_0003 human1_lib1.final_cell_0004 Acinar cells 58.1720116 63.2674918 Alpha cells 0.1910771 0.1770659 Beta cells 0.4326425 0.1770659 Delta cells . . pancreas_gs_prediction <- rownames(HGT_pancreas_gs)[apply(HGT_pancreas_gs, 2, which.max)] #矩阵的每一列依次:判断是否是最大值,得到一列布尔值,结合矩阵的行名会返回行名中的一个元素(也就是最有可能的细胞类型)。 #所有列运行完了之后会得到所有细胞最可能的注释 pancreas_gs_prediction_signif <- ifelse(apply(HGT_pancreas_gs, 2, max)>2, yes = pancreas_gs_prediction, "unassigned") #如果`-log10 corrected p-value`的值小于等于2,则认为不显著,注释更正为"unassigned" Baron$pancreas_gs_prediction <- pancreas_gs_prediction_signif head(Baron@meta.data,2) orig.ident nCount_RNA nFeature_RNA human1_lib1.final_cell_0001 human1 21956 3362 human1_lib1.final_cell_0002 human1 27309 4009 pancreas_gs_prediction human1_lib1.final_cell_0001 Acinar cells human1_lib1.final_cell_0002 Acinar cells
到这儿,细胞类型预测算结束了。再来用真正的细胞类型(也可以是手动注释的结果),评估一下软件预测出来的准不准。走Seurat标准流程:
#加载预先知道的注释信息 BaronMetaData <- readRDS(url("https://storage.googleapis.com/cellid-cbl/BaronMetaData.rds")) Baron@meta.data=merge(Baron@meta.data,BaronMetaData,by="row.names") rownames(Baron@meta.data)=Baron@meta.data$Row.names Baron@meta.data$Row.names=NULL head(Baron@meta.data,2) orig.ident nCount_RNA nFeature_RNA human1_lib1.final_cell_0001 human1 21956 3362 human1_lib1.final_cell_0002 human1 27309 4009 pancreas_gs_prediction donor cell.type human1_lib1.final_cell_0001 Acinar cells GSM2230757 Acinar cells human1_lib1.final_cell_0002 Acinar cells GSM2230757 Acinar cells Baron <- FindVariableFeatures(Baron, selection.method = "vst", nfeatures = 2000) Baron <- ScaleData(Baron) Baron <- RunPCA(Baron, npcs = 50, verbose = FALSE) Baron <- FindNeighbors(Baron, dims = 1:30) Baron <- FindClusters(Baron, resolution = 0.5) Baron <- RunUMAP(Baron, dims = 1:30) Baron <- RunTSNE(Baron, dims = 1:30) library(cowplot) p1 <- DimPlot(Baron, reduction = "tsne", group.by = "cell.type", pt.size=0.5, label = TRUE,repel = TRUE) p2 <- DimPlot(Baron, reduction = "tsne", group.by = "pancreas_gs_prediction", pt.size=0.5, label = TRUE,repel = TRUE) fig_tsne <- plot_grid(p1, p2, labels = c('cell.type','pancreas_gs_prediction'),align = "v",ncol = 2) ggsave(filename = "tsne.pdf", plot = fig_tsne, device = 'pdf', width = 40, height = 15, units = 'cm')
这种细胞注释方法非常依赖于细胞类型的参考基因集的质量,我拿自己的数据集跑这个软件,并没有该软件的示例数据效果好,原因就是前面提到的PanglaoDB网站中我这个方向的参考数据质量不高。大家还知道类似的提供了细胞类型及对应的基因集网站吗?欢迎讨论,集思广益,我会在之后的推文中列出大家推荐的网站。
另外,这种用富集来注释的思路可以延伸一下,比如对于某个cluster不太清楚是什么类型,可以找完差异基因后,用细胞类型参考基因集做富集分析,来定义cluster类型。
因水平有限,有错误的地方,欢迎批评指正!
- 初识python(注释,数据类型,标识符和关键字,运算符,分支语句等)01/14
- 单细胞文章解读——用单细胞RNA测序技术分析与肿瘤转移相关的细胞间通讯
- 停车场模拟问题的分析思路,解决细节。加了关键注释。
- 专题二 符号的技巧---- 14.优先级和类型转换分析
- Java重修之路(二)第一个程序分析,注释,关键字,标识符,常量,进制转换,类型转换
- 单细胞分析实录(9): 展示marker基因的4种图形(二)
- 单细胞分析实录(2): 使用Cell Ranger得到表达矩阵
- 单细胞分析实录(8): 展示marker基因的4种图形(一)
- 单细胞分析实录(5): Seurat标准流程
- 链表的顺序插入 先分析有几种类型再编写 思路要清晰
- 韦东山嵌入式第一期学习笔记DAY_15——14_8改进上一节代码与用C实现面向对象编程思路分析
- 单细胞分析实录(3): Cell Hashing数据拆分
- 单细胞分析实录(6): 去除批次效应/整合数据
- 单细胞分析实录(15): 基于monocle2的拟时序分析
- 命名空间“Microsoft.Office.Interop”中不存在类型或命名空间名称“Excel”。是否缺少程序集引用 的另一种解决方案
- 判断cpu类型的另一种方案
- Behavior Designer中节点类型分析
- ASP.NET绘制折线图---(1)思路分析
- linux-0.11调试教程,fflush()函数实现分析和FILE类型定义
- 模拟电子技术14 22------——放大电路 静态分析续