本文学习单细胞转录组分析过程中找marker基因、根据marker基因确定细胞、只标定特定细胞群的聚类分群图
六、找marker基因
marker基因和差异基因里面的上调基因有点类似：某个基因在某一簇细胞里表达量都很高（或低），在其他簇表达量很低（或高），那么这个基因就是这簇细胞的象征
1.找全部cluster的marker基因
下面这行代码运行也很耗时，尤其是细胞数量多的时候，最好设置“存在即跳过”这个选项，得到了2889个，两千多个已经需要等十秒多了
pbmc.markers only.pos = TRUE, # 只返回positive基因
min.pct = 0.25)
图29
我们点开来看一下，如图30，会发现一个问题：例如ABI3、ABI31、ABI32这三个行名对应的gene列的内容都是ABI3，这里可能踩坑
首先为什么会这样呢？因为在数据框的行名里面不允许重复，所以会被加上1、2；那后续用行名来取基因会有什么后果呢？虽然不会报错，但如果正好筛选到了ABI31，那么它就没法和表达矩阵里面的ABI3匹配；所以按照行名来筛选就会踩坑，没被匹配上这个基因就丢失了，后续得到基因的数量就会比别人少
图30
说完坑，我们再来看一下图30得到的结果是什么意思，以第一行为例
在7这一簇里面，第一组中有45.9%的细胞表达了AAMP基因，第二组中有11.3%的细胞表达了AAMP基因
图31
接着，只计算至少在(两簇细胞总数的)25%的细胞中有表达的基因
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)#将每个cluster里面logFC排名前2的基因筛选出来
有些p值很小，甚至是0。这不是出错，而是它们小到逼近于0了
图32
2.比较某个基因在几个cluster之间的表达量
我们拿图32结果中最后两个基因画个图看看
# 小提琴图 默认拿ScaleData来做的
VlnPlot(pbmc, features = c("PF4", "PPBP"))
从图33能一眼看出这两个基因确实是第8簇的marker gene，它们在第8簇中的表达量很高
图33
# 也可以拿count数据画，标准流程里面有所以才在此展示，一般不用这个
# VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
3.FeaturePlot
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
FeaturePlot是一种经典图，实现了在umap图上进行标记。从深紫→银白，数值越来越小；哪里颜色深，就表明基因在哪簇表达量高，
这张图经常出现在数据挖掘文章，比如预后分析得到了几个基因，想在单细胞的数据上也进行投影，期望得到的现象是：单细胞水平上这些基因也是在tumor细胞中表达量比较高，在normal里面表达量比较低
图34
4.marker基因的热图
# 从pbmc.markers的每簇里面选出top 10，共有9个簇，就会得到90个基因
top10 % group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
黄色高表达，紫色低表达。可以看到在对角线位置，都是黄色块
图35
七、根据marker基因确定细胞
1.DotPlot
markers.txt这个文件是手动创建的，如果换一个数据了就要重新找了，比如在数据库找、在文献里查，后续介绍
a = read.delim("supp/markers.txt",header = F)#得到一个数据框，V1是细胞名，V2是对应的marker gene名
36
下面这一步格式转换，把数据框转为列表，是为了接下来在图38中，我们能清楚地看到细胞的名字
gt = split(a[,2],a[,1])#把a里面的V2列按照V1列来拆分
