ixxmu
commented
1 year ago 各个单细胞亚群的特异性基因集合的打分能准确划分其亚群吗？ by 生信技能树
绝大部分初学者在拿到了单细胞表达量矩阵后就直接开始降维聚类分群，都是有标准代码，遇到的第一个拦路虎就是亚群的生物学命名，因为这个步骤并不是代码技巧反而是生物学背景。往往是有一个错误的预期，就是各个单细胞亚群泾渭分明，非此即彼，区分的清清楚楚。
而实际情况下，不同层次的细胞亚群的界限容忍度就不一样。比如肿瘤相关单细胞数据集常规分析都是拿到表达量矩阵后的第一层次降维聚类分群通常是：
immune (CD45+,PTPRC), epithelial/cancer (EpCAM+,EPCAM), stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)
参考我前面介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这3大单细胞亚群构成了肿瘤免疫微环境的复杂。绝大部分文章都是抓住免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群。但是也有不少文章是抓住stromal 里面的fibo 和endo进行细分，并且编造生物学故事的。
在第一层次降维聚类分群，当然是免疫细胞和上皮细胞界限非常明晰，除非是背景污染没有去除或者少量的实验环节引入误差，比如：不知道为什么偶尔能看到T淋巴细胞跟肿瘤细胞在umap临近 就是这样的情况。而往往是在第3层次，很多亚群其实就本来是不应该是有清晰的界限，比如我们熟知的pbmc3k数据集里面的部分免疫细胞就比较类似：
# devtools::install_github('satijalab/seurat-data')
library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
# InstallData("pbmc3k")
data("pbmc3k")
library(Seurat)
sce <- pbmc3k.final
p1=DimPlot(sce, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
p1
# 需要合并一些单细胞亚群：
levels(Idents(sce))
## Assigning cell type identity to clusters
new.cluster.ids <- c("CD4 T","CD4 T", "Mono", "B", "CD8 T",
"Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(sce)
sce <- RenameIdents(sce, new.cluster.ids)
p2 = DimPlot(sce, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
library(patchwork)
p1+p2
如下所示，两个CD4的T细胞比较类似，两个单核细胞也是如此，CD8和NK也是如此：
因为比较类似，如果你没有十足的把握在这个umap把它们区分开来， 就可以不用那么细致的命名。这个时候我们来做一个有意思的探索， 就是各个单细胞亚群的特异性基因集合的打分能准确划分其亚群吗？
future::plan("multiprocess", workers = 4)
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE,
min.pct = 0.25,
thresh.use = 0.25)
save(sce.markers,file = 'sce.markers.FindAllMarkers.pbmc.Rdata')
library(dplyr)
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_log2FC)
DoHeatmap(sce ,unique(top10$gene))
ggplot2::ggsave('top10_DoHeatmap.pdf',height = 10)
可以看到，两个CD4的T细胞其实大量共享高表达量基因，两个单核细胞也是如此，而CD8和NK也是如此：
这本来就是符合生物学背景的，当然了如果你换一个算法，或者修改FindAllMarkers的参数，也可以达到每个单细胞亚群可以寻找到更加特异一点的高表达量基因列表。
如果我们替换成为了top基因的打分，然后可视化，代码如下所示：
top10 = split(top10$gene,top10$cluster)
names(top10)
# ?AddModuleScore()
sce = AddModuleScore(sce,top10)
names(sce@meta.data)
colnames(sce@meta.data)[8:16] = levels(Idents(sce))
VlnPlot(sce,features = levels(Idents(sce)) ,
stack = T)
ggplot2::ggsave('top10_AddModuleScore_VlnPlot.pdf' )
可以看到区分度好了一点：
如果把top10的基因替换成为top100，你会发现肉眼看上去，区分情况会有一点点好转。
如果我们把这个pbmc3k数据集区分的粗糙一点，效果会更好一点，说明这个时候的cd8和nk其实并不是很好的界限，本来就是可以提取子集后继续细分，这个时候的pbmc3k的示范命名反而并不是很合理。
写在文末
我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。
https://mp.weixin.qq.com/s/3DCdxI2rwdylYwub7Vqpaw