细胞注释全流程解析
小伙伴们大家好,我是晨曦,上一期的scRNA-seq标准流程不知道大家学习的如何?其实自从Seurat包被开发出来以后,scRNA-seq的分析流程一直到细胞注释这里都有了很好的标准化流程,类似芯片分析或者bulk-seq一样,在理解完每一步的具体含义和做法后,在往后分析自己项目数据的时候,就是跑标准化的流程。
这其中唯一比较困难且具有一定挑战的就是细胞注释,怎么样可以快速帮助我们明确细胞类型一直是一个不小的难点,各种各样的细胞注释工具不断的被发表,据统计现在已经存在并且被证明可以实际应用的细胞注释工具就至少有30种,在这里晨曦将为大家展开scRNA-seq细胞注释全流程解析,完全手把手教学,相信大家一定可以快速掌握这门技术,为自己的scRNA-seq增光添彩~
在学习细胞注释之前,让我们再跑一下上次教程中的标准流程以获得细胞的marker基因
#标准流程(还有不会的快去看上一篇的推文教程哦~)library(dplyr)library(Seurat)library(patchwork)pbmc.data<- Read10X(data.dir = "hg19")pbmc<- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)pbmc[["percent.mt"]]<- PercentageFeatureSet(pbmc, pattern = "^MT-")pbmc<- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)pbmc<- NormalizeData(pbmc)pbmc<- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)pbmc<- ScaleData(pbmc)pbmc<- RunPCA(pbmc, features = VariableFeatures(object = pbmc))pbmc<- JackStraw(pbmc, num.replicate = 100)pbmc<- ScoreJackStraw(pbmc, dims = 1:20)ElbowPlot(pbmc)pbmc<- FindNeighbors(pbmc, dims = 1:10)pbmc<- FindClusters(pbmc, resolution = 0.5)pbmc<- RunTSNE(pbmc, dims = 1:10)DimPlot(pbmc,reduction = "tsne", label = TRUE)
晨曦解读
到这一步,我们已经完成了细胞亚群的分类,可视化如下
#细胞注释的方法在上一个教程中已经写明,但是这里我想做一个归类,其实无外乎两大类#第一类:根据每个细胞的marker基因进行注释#第二类:通过已知细胞类型的基因表达情况与未知的细胞类型的基因表达情况进行比对,相似即推测该未知细胞类型是其已知细胞类型#然后这两种思想下又会有很多种方法R包、文献、数据库等等#寻找细胞的marker基因pbmc.markers <- FindAllMarkers(pbmc)
晨曦解读
下面就是我们获得的all.marker基因的数据框,下面我们需要对这个数据框进行一些筛选
#筛选P值小于0.05(必须)all.markers = pbmc.markers %>% subset(p_val<0.05)
晨曦解读
流程是没有问题,但是晨曦看这个数据框会很难受,因为它和我们平时看的表达矩阵不一样,所以哦我们需要把这个数据框修饰一下~
all.markers=pbmc.markers %>% select(gene,everything()) %>% subset(p_val<0.05)#这样就好看不少,果然是这样看的习惯
晨曦解读
这里选择的是pvalue,为什么不选择adjp呢?
首先P值小于0.05是基本,adjP是校正后的P值,可以更好的控制假阳性率,所以也比P值更加严格,所以这里我选择P值是为了纳入更多的差异基因来增加下一步的输入数据~
#然后进行FC的筛选(必须)(选择每个亚群最大的前十个基因)top10 <- all.markers %>%group_by(cluster) %>%top_n(n = 10, wt = avg_log2FC)#如果想做的更详细一点,可以筛选pct.1>0.5的基因作为差异基因的候选,大家可以想想是为什么?#当然这一步骤是选做的,筛选FC更重要~
#然后导出文件,作为我们后续分析的输入文件write.csv(all.markers, "all.markers.csv",row.names = F)write.csv(top10, "top10.csv",row.names = F)
晨曦解读
这里我们就获得了每个细胞亚群的差异基因,且筛选出了Top10的DEGs
但是这里面仍然有我们可以探究的东西,那就让我们来探究一下FindAllMarkers这个函数吧
#FindAllMarkers函数?FindAllMarkers#我们可以调取这个函数的帮助文档,然后我们可以看到,这个函数其实提供了很多获取差异基因的方式,当然探索算法其实算是一个比较高的要求了,我们这里只是知道获得的结果有什么不同,不要过度探究算法,那样本身其实收益不是很大#pbmc.markers <- FindAllMarkers(pbmc,test.use = 'MAST')#pbmc.markers <- FindAllMarkers(pbmc,test.use = 'DESeq2',slot = 'counts')
晨曦解读
我们通过参数test.use可以选择不同的算法,默认是wilcox算法,这个其实就已经可以了,让我感兴趣的是里面还有几种我们以前见过的老朋友,比如"DESeq2"(如果选择这个算法, slot参数得选择“counts”),其中还有Mast算法(这个是专门给scRNA-seq定制的),运算后发现单纯取Top10来说结果其实差距不大(默认设置即可,这里只不过提供一下更多可能)
以上为止,我们就获得了每个cluster的marker,下面我们将通过在线数据库——Cellmarker数据库对细胞亚群进行注释并且把注释结果标注到我们的tSNE上
Cellmarker数据库
晨曦解读
上面这个就是Cellmarker数据库的主页面,下面我们将简单介绍一下这个数据库,并且会一步步完成细胞注释的过程(这个可能是网上比较全的保姆级教程了吧~骄傲)
首先需要了解一下这个数据库的基本知识,这个数据库来自哈尔滨医科大学,该团队通过梳理100,000+发表的文献,梳理出人的158个组织 (亚组织)的467个细胞类型的13,605个Marker基因,和鼠的81个组织 (亚组织)的389个细胞类型的9, 148个Marker基因,都是可搜索,可浏览,可下载(总结一句:数据很多,很靠谱)
首页(Home):提供人和小鼠的全局视图(右上角小图标可以切换)
晨曦解读
这其实就是从另一方面告诉我们,我们要想进行细胞注释,首先我们要明确我们的组织来源以及对我们组织内的细胞类型有一个大致的了解,有了这样的大致了解,我们就可以通过这个界面进行快速检索
浏览(Browse):该页面展示了细胞和组织的层次,我们可以通过点击左边的导航栏来选择我们想要的组织和细胞类型(对比上一个,像是餐厅点菜和餐厅自助之间的区别)
Search(搜索):这个界面分为两类,一个是Cell search,一个是Marker search
Cell search支持物种-组织-细胞类型来搜索目标细胞的细胞marker
Marker search支持Gene symbol、Gene ID或蛋白名称来检索以上是哪个组织那个细胞的marker基因
下载(Download):如果想用R语言操作,可以下载细胞类型注释(可以简单理解就像是探针注释表格那样)然后进入R语言进行匹配操作
晨曦解读
检索的方式就是以上这些,下面我们拿我们的数据来进行实际操作,看看我们是如何使用并且完成细胞注释的吧~
首先需要明确,我们的数据是外周血单核细胞(PBMC),所以说我们检索到的细胞类型最起码是要来源于外周血,然后我们进入实际操作
首先我们以一个亚群为例子,后续亚群只要重复进行相同的操作即可
晨曦解读
然后我们记录我的cluster0的marker基因为:CCR7、LDLRAP1、LDHB、LEF1、PRKCQ-AS1、PIK3IP1、MAL、NOSIP、CD3D、FHIT,然后我们把这些Gene marker键入到我们的搜索框,然后通过限定组织来源,我们得到了下面的数据
晨曦解读
然后这时候会有小伙伴犯难了,我这个只是检索了一个基因CCR7,就出现了这些细胞类型我该如何选择呢?
我们再检索一下下个marker基因,一般来说我们需要分别键入这10个marker基因,然后我们大概会得到很多细胞类型,然后我们需要结合我们这个领域的文献还有其它数据库,去进一步缩小这个范围,相信同学们也发现了,细胞注释是一个因人而异的事情,大方向对,有的人注释细,有的人注释粗,这个就是看每个人的习惯了
接下来我们检索LDLRAP1这个基因,然后限定组织类型(这个时候你就会知道预先知道组织类型是有多么重要?!)
晨曦解读
交集产生了,Naive CD4+ T cell和Naive CD8+ T cell,也就是说如果我做到这里,我可以直接把这个细胞注释成免疫细胞也是可以的(经验之谈,虽然我们选择了Top10,但是一般来说,我们注释Top10中的5个左右Gene即可得到细胞注释,当然这5个基因需要是Top10中也是处在前5的)
然后下面我们加快一下脚步,把Top5中的剩下三个注释出来,如下
晨曦解读
综上所述,表达量靠前的Gene基本上都是在Naive CD8+ T cell,所以我们可以认为Cluster0是Naive CD8+ T cell(这里其实和官网产生了一定的冲突,官网给Cluster0注释为Naive CD4+ T cell,这个可能是考虑到生物学背景,比如看文献或者联合多个数据库得到的结论,但是咱们这个以一个数据库进行的注释也不是错误答案,毕竟细胞亚群的注释本来就没有一个明确的答案)
但是毕竟产生了分歧,那么我们就要去解决,所以Cluster0我们可以注释为Naive T cell,这样的答案想必肯定是不会有错误的~
首先我这里个人认为细胞注释存在着太多主观的思路,参考一个数据库和多个数据库的结果可能会不一样,这里哦我认为我们需要以我们的Idea出发,比如说我们写作的出发点是免疫方向,那么我们只需要知道这个细胞亚群是免疫细胞,然后后续我们通过把这个亚群中的marker基因进行富集然后回过头来再定义这个细胞亚群的功能然后再次考虑细胞亚型的定义是否合理,个人感觉这是一个比较合理的思路~
然后有些文献中也会存在着作者整理好的细胞类型对应着marker基因,比如下面这个表格(这类文献如何检索,就是查找本领域内的单细胞文章,高分的文章往往作者会在图表部分附上这类表格)(参考文献:dropClust: efficient clustering of ultra-large scRNA-seq data)
晨曦解读
看起来我们的注释还是十分可靠的,然后这里说一下,上面这种表格有一些也会被我们的数据库收入,例如我们的Cellmarker数据库就收录的大量的文献,但是毕竟不是全部,所以这就是我们细胞注释需要多工具,多文献然后再结合自己的生物学知识进行综合注释~
#然后我们对tSNE重新进行注释,剩下细胞类群只需要重复上面的操作即可完成#注释这种方式属于人工注释,优点是准确性高,自由,缺点则是耗时费力new.cluster.ids <- c("Naive T cell", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono","NK", "DC", "Platelet")names(new.cluster.ids) <- levels(pbmc)pbmc <- RenameIdents(pbmc, new.cluster.ids)DimPlot(pbmc, reduction = "tsne", label = TRUE, pt.size = 0.5) + NoLegend()
晨曦解读
这样我们运用cellmarker数据库进行细胞注释的流程到这里就告一段落了,当然今天这次的侧重主要是在线数据库的讲解与使用,下次给大家介绍基于R语言进行细胞注释的方法,其代表就是Singler包,那么我们就下次再见啦~
内心OS:本来今天想增加Singler包的讲解,但是奈何展开篇幅实在是太长,可能会影响观感(绝对不是参考数据集没有下载成功QAQ,仍在不断调试中......)
我是晨曦,我们下次再见~
回复“晨曦04,即可获得今天的范例数据和代码哦~
END

撰文丨晨   曦
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
继续阅读
阅读原文