对每个单细胞亚群进行独立分析，以确定两个亚群的差异

最编程 2024-03-05 21:20:29

...

前些天我们的单细胞学徒培养有小伙伴分享了文章；在 JCI Insight 2022 https://doi.org/10.1172/jci.insight.152616 ，里面对第一次降维聚类分群后的各个单细胞亚群独立在两个分组做差异分析，如下所示：

各个单细胞亚群独立在两个分组做差异分析

可以看到，每个单细胞亚群都有自己的差异分析火山图，会议上有人提问这个分析如何做。其实主要是大家可能是初次接触生物信息学就是单细胞数据处理，所以基础知识有点欠缺。它就是普通的表达量矩阵分析而已，我七八年前就写过系列笔记，公众号推文在：

解读GEO数据存放规律及下载，一文就够
解读SRA数据库规律一文就够
从GEO数据库下载得到表达矩阵一文就够
GSEA分析一文就够（单机版+R语言版）
根据分组信息做差异分析- 这个一文不够的
差异分析得到的结果注释一文就够

我们这里以大家熟知的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" 在这里基本上没有统计学显著性。

上一篇：简单了解 k8s 资源限制参数中的 CPU 计量单位

下一篇：压力测试指标（QPS、TPS、PV、RT）