本文是参考学习单细胞转录组基础分析四：细胞类型鉴定的学习笔记。可能根据学习情况有所改动。
细胞完成聚类和降维可视化之后，细胞群的异质性已经体现出来了，那么各个cluster的细胞如何鉴定细胞类型呢？通用的方法有两种：一种是通过细胞类型特异性表达的marker基因识别，就好像流式细胞仪用特定的抗体筛选细胞一样；第二种是建立已知细胞类型的转录谱数据库，将未知细胞类型的转录谱数据与之比较相似性，就知道它最有可能是哪种细胞了。第一种方法需要人工收集marker基因比对各个cluster的显著高表达基因综合分析，第二种方法可以使用SingleR包自动识别细胞类型。建议两种方法结合起来进行细胞鉴定。Cluster差异****基因通过marker基因鉴定细胞类型目前仍是最常用的方法，应用此方法的基础是找到各个cluster的显著高表达的基因，这就需要差异分析提供候选基因列表。Seura提供多种差异分析的方法，默认wilcox方法，MASK方法和DESeq2方法也可以留意一下。MASK是专门针对单细胞数据差异分析设计的，DESeq2是传统bulkRNA数据差异分析的经典方法。Seurat函数通过test.use参数确定分析方法，如果选用 "negbinom", "poisson", or "DESeq2", 需要将slot参数设为"counts"。
##以下方法三选一，建议第一种
以上三种方法我都运行了一遍，top10显著高表达基因差别很小，cluster0的结果如下图所示：
图片
鉴定细胞类型****人工鉴定细胞类型人工鉴定细胞类型，首先要明确自己研究组织可能有哪些细胞类型，各种细胞的marker基因是什么，然后对比差异分析得到的各个cluster的显著高表达基因，综合分析就可以判断细胞类型了。细胞的marker基因主要通过相关领域的文献收集，也可以通过专门的数据库查找，推荐两个比较常用的数据库：
CellMarker：http://biocc.hrbmu.edu.cn/CellMarker/index.jsp
PanglaoDB：https://panglaodb.se/index.html
使用top10显著高表达基因，结合CellMarker数据库，粗略鉴定结果如下：
图片
选择基因作图展示
#挑选部分基因
图片
SingleR鉴定细胞类型
library(SingleR)
鉴定结果与手工鉴定的基本一致：
图片
鉴定结果展示
p1 = DimPlot(scRNA, group.by="celltype", label=T, label.size=5, reduction='tsne')
图片
# 细胞完成聚类和降维可视化之后，细胞群的异质性已经体现出来了，那么各个cluster的细胞如何鉴定细胞类型呢？通用的方法有两种：一种是通过细胞类型特异性表达的marker基因识别，就好像流式细胞仪用特定的抗体筛选细胞一样；第二种是建立已知细胞类型的转录谱数据库，将未知细胞类型的转录谱数据与之比较相似性，就知道它最有可能是哪种细胞了。第一种方法需要人工收集marker基因比对各个cluster的显著高表达基因综合分析，第二种方法可以使用SingleR包自动识别细胞类型。建议两种方法结合起来进行细胞鉴定。
#
# Cluster差异基因
# 通过marker基因鉴定细胞类型目前仍是最常用的方法，应用此方法的基础是找到各个cluster的显著高表达的基因，这就需要差异分析提供候选基因列表。Seura提供多种差异分析的方法，默认wilcox方法，MASK方法和DESeq2方法也可以留意一下。MASK是专门针对单细胞数据差异分析设计的，DESeq2是传统bulkRNA数据差异分析的经典方法。Seurat函数通过test.use参数确定分析方法，如果选用 "negbinom", "poisson", or "DESeq2", 需要将slot参数设为"counts"。
##以下方法三选一，建议第一种
rm(list=ls())
options(stringsAsFactors = F)
dir.create("cell_identify")#
#默认wilcox方法
scRNA <- readRDS(file = "scRNA.rds")
library(Seurat)
diff.wilcox = FindAllMarkers(scRNA)
library(tidyverse)
all.markers = diff.wilcox %>% dplyr::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)
plot1#先查看再保存
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)
p1 <- VlnPlot(scRNA, features = select_genes, pt.size=0, ncol=2)
p1
ggsave("cell_identify/selectgenes_VlnPlot.png", p1, width=6 ,height=8)
#featureplot展示
p2 <- FeaturePlot(scRNA, features = select_genes, reduction = "tsne", label=T, ncol=2)
p2
ggsave("cell_identify/selectgenes_FeaturePlot.png", p2, width=8 ,height=12)
p3=p1|p2
p3
ggsave("cell_identify/selectgenes.png", p3, width=10 ,height=8)
# 图片
#
#
#
# SingleR鉴定细胞类型
library(SingleR)
refdata <- HumanPrimaryCellAtlasData()
testdata <- GetAssayData(scRNA, slot="data")
clusters <- scRNA@meta.data$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)
scRNA@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
scRNA@meta.data[which(scRNA@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
# 鉴定结果与手工鉴定的基本一致：
# 图片
#
# 鉴定结果展示
library(patchwork)
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')
p1
p2
p3
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)
#图片
