单细胞转录组基础分析四：细胞类型鉴定

单细胞测序技术是近年最大的生命科学突破之一，相关文章频繁发表于各大顶级期刊，然而单细胞数据的分析依然是大家普遍面临的障碍。本专题将针对10X Genomics单细胞转录组数据演示各种主流分析，包括基于Seurat的基础分析、以及基于clusterProfiler、Monocle、SingleR等R包的延伸分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨！（扫描文末二维码）

细胞完成聚类和降维可视化之后，细胞群的异质性已经体现出来了，那么各个cluster的细胞如何鉴定细胞类型呢？通用的方法有两种：一种是通过细胞类型特异性表达的marker基因识别，就好像流式细胞仪用特定的抗体筛选细胞一样；第二种是建立已知细胞类型的转录谱数据库，将未知细胞类型的转录谱数据与之比较相似性，就知道它最有可能是哪种细胞了。第一种方法需要人工收集marker基因比对各个cluster的显著高表达基因综合分析，第二种方法可以使用SingleR包自动识别细胞类型。建议两种方法结合起来进行细胞鉴定。

Cluster差异基因

通过marker基因鉴定细胞类型目前仍是最常用的方法，应用此方法的基础是找到各个cluster的显著高表达的基因，这就需要差异分析提供候选基因列表。Seura提供多种差异分析的方法，默认wilcox方法，MASK方法和DESeq2方法也可以留意一下。MASK是专门针对单细胞数据差异分析设计的，DESeq2是传统bulkRNA数据差异分析的经典方法。Seurat函数通过test.use参数确定分析方法，如果选用 "negbinom", "poisson", or "DESeq2",  需要将slot参数设为"counts"。

##以下方法三选一，建议第一种
#默认wilcox方法
diff.wilcox = FindAllMarkers(scRNA)
all.markers = diff.wilcox %>% select(gene, everything()) %>% subset(p_val<0.05)
top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
write.csv(all.markers, "cell_identify/diff_genes_wilcox.csv", row.names = F)
write.csv(top10, "cell_identify/top10_diff_genes_wilcox.csv", row.names = F)
#专为单细胞设计的MAST
diff.mast = FindAllMarkers(scRNA, test.use = 'MAST')
all.markers = diff.mast %>% select(gene, everything()) %>% subset(p_val<0.05)
top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
write.csv(all.markers, "cell_identify/diff_genes_mast.csv", row.names = F)
write.csv(top10, "cell_identify/top10_diff_genes_mast.csv", row.names = F)
#bulkRNA经典方法DESeq2
diff.deseq2 = FindAllMarkers(scRNA, test.use = 'DESeq2', slot = 'counts')
all.markers = diff.deseq2 %>% select(gene, everything()) %>% subset(p_val<0.05)
top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
write.csv(all.markers, "cell_identify/diff_genes_deseq2.csv", row.names = F)
write.csv(top10, "cell_identify/top10_diff_genes_deseq2.csv", row.names = F)
##top10基因绘制热图
top10_genes <- read.csv("cell_identify/top10_diff_genes_wilcox.csv")
top10_genes = CaseMatch(search = as.vector(top10_genes$gene), match = rownames(scRNA)) 
plot1 = DoHeatmap(scRNA, features = top10_genes, group.by = "seurat_clusters", group.bar = T, size = 4)
ggsave("cell_identify/top10_markers.pdf", plot=plot1, width=8, height=6) 
ggsave("cell_identify/top10_markers.png", plot=plot1, width=8, height=6)

以上三种方法我都运行了一遍，top10显著高表达基因差别很小，cluster0的结果如下图所示：

鉴定细胞类型

人工鉴定细胞类型

人工鉴定细胞类型，首先要明确自己研究组织可能有哪些细胞类型，各种细胞的marker基因是什么，然后对比差异分析得到的各个cluster的显著高表达基因，综合分析就可以判断细胞类型了。细胞的marker基因主要通过相关领域的文献收集，也可以通过专门的数据库查找，推荐两个比较常用的数据库：

CellMarker：http://biocc.hrbmu.edu.cn/CellMarker/index.jsp

PanglaoDB：https://panglaodb.se/index.html

使用top10显著高表达基因，结合CellMarker数据库，粗略鉴定结果如下：

选择基因作图展示

#挑选部分基因
select_genes <- c('LYZ','CD79A','CD8A','CD8B','GZMB','FCGR3A')
#vlnplot展示
p1 <- VlnPlot(scRNA, features = select_genes, pt.size=0, group.by="celltype", ncol=2)
ggsave("cell_identify/selectgenes_VlnPlot.png", p1, width=6 ,height=8)
#featureplot展示
p2 <- FeaturePlot(scRNA, features = select_genes, reduction = "tsne", label=T, ncol=2)
ggsave("cell_identify/selectgenes_FeaturePlot.png", p2, width=8 ,height=12)
p3=p1|p2
ggsave("cell_identify/selectgenes.png", p3, width=10 ,height=8)

SingleR鉴定细胞类型

library(SingleR)
refdata <- HumanPrimaryCellAtlasData()
testdata <- GetAssayData(scRNA, slot="data")
clusters <- [e-mail]$seurat_clusters
cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.main, 
                     method = "cluster", clusters = clusters, 
                     assay.type.test = "logcounts", assay.type.ref = "logcounts")
celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels, stringsAsFactors = F)
write.csv(celltype,"cell_identify/celltype_singleR.csv",row.names = F)
[e-mail]$celltype = "NA"
for(i in 1:nrow(celltype)){
[e-mail][which([e-mail]$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}

鉴定结果与手工鉴定的基本一致：

鉴定结果展示

p1 = DimPlot(scRNA, group.by="celltype", label=T, label.size=5, reduction='tsne')
p2 = DimPlot(scRNA, group.by="celltype", label=T, label.size=5, reduction='umap')
p3 = plotc <- p1+p2+ plot_layout(guides = 'collect')
ggsave("cell_identify/tSNE_celltype.pdf", p1, width=7 ,height=6)
ggsave("cell_identify/UMAP_celltype.pdf", p2, width=7 ,height=6)
ggsave("cell_identify/celltype.pdf", p3, width=10 ,height=5)
ggsave("cell_identify/celltype.png", p3, width=10 ,height=5)

获取帮助

本教程的目的在于把常用的单细胞分析流程串起来，适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释，官网有详细的教程和权威的解释，翻译好的中文教程也有多个版本，有兴趣的可以搜索一下。如果您学习本教程有一定困难，可以分享此篇文章到朋友圈，截图微信发给Kinesin（文末二维码），我会抽时间组群直播讲解单细胞数据分析的全过程。本专题所用的软件、数据、代码脚本和中间结果，我打包放在了百度云上，需要的朋友添加Kinesin微信索取。