前些天我们的单细胞学徒培养有小伙伴分享了文章;在 JCI Insight 2022 https://doi.org/10.1172/jci.insight.152616 ,里面对第一次降维聚类分群后的各个单细胞亚群独立在两个分组做差异分析 ,如下所示:
各个单细胞亚群独立在两个分组做差异分析
可以看到,每个单细胞亚群都有自己的差异分析火山图,会议上有人提问这个分析如何做。其实主要是大家可能是初次接触生物信息学就是单细胞数据处理,所以基础知识有点欠缺。它就是普通的表达量矩阵分析而已,我七八年前就写过系列笔记,公众号推文在:
我们这里以大家熟知的pbmc3k数据集为例。大家先安装这个数据集对应的包,并且对它进行降维聚类分群,参考前面的例子:人人都能学会的单细胞聚类分群注释
# 0.安装R包 ----
# InstallData("pbmc3k")
library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
table(Idents(sce))
DimPlot(sce,label = T)
这个时候,因为它pbmc3k数据集并没有分组,所以我们没办法做差异分析。不过我们可以简单的模拟一个分组。如下所示:
sce$celltype = Idents(sce)
sce$group = sample(1:2,ncol(sce),replace = T)
table(sce$celltype,sce$group )
# 如下所示,两个随机赋予的分组,每个分组里面的都是有这些不同单细胞亚群
Naive CD4 T 349 348
Memory CD4 T 231 252
CD14+ Mono 224 256
B 186 158
CD8 T 135 136
FCGR3A+ Mono 81 81
NK 66 89
DC 12 20
Platelet 9 5
下面对它们进行批量差异分析:
Idents(sce) = paste0('c',sce$group )
table(Idents(sce))
degs = lapply(unique(sce$celltype), function(x){
FindMarkers(sce[,sce$celltype==x],ident.1 = 'c1',
ident.2 = 'c2')
})
x=degs[[1]]
do.call(rbind,lapply(degs, function(x){
table(x$avg_log2FC > 0 )
}))
值得注意的是这个FindMarkers函数并不是最好的单细胞转录组表达量矩阵的差异分析方法,我这里仅仅是举例哦!
可以看到,如果是以为 avg_log2FC 标准,这个时候很容易得到假阳性的差异基因:
FindMarkers函数的结果
其实我们更应该关心的是 "pct.1" 和 "pct.2" 的差异,这个也是各个单细胞亚群特异性高表达量基因的金标准,不过跟我们这个时候的差异分析不太一样的需求哦,需要自己多思考。
另外,这个时候因为我们是随机赋值的两个分组,所以差异分析理论上就没有意义,所以"p_val_adj" 在这里基本上没有统计学显著性。