缺氧跟其它生物学特征一样，本质上就是一个功能基因列表，比如缺氧基因列表就可以根据 GO:0001666 ‘‘Response to hypoxia’’ 获得，是 Hif1a, Slc2a1, Vegfa, Hmox1, Bnip3, Nos2, Mmp2, Sod3, Cited2, Ldha 这些基因。
评分的算法很多，gsea，gsva等等，单细胞领域比较出名的是Seurat包的AddModuleScore函数，以及UCell和AUCell等包，从代码的角度来看，都是一个函数即可。
有了评分这个武器，我们很容易对各种各样的基因列表进行自由自在的探索，我们分享过的就有：
- • 借鉴escape包的一些可视化GSVA或者ssGSEA结果矩阵的方法
- • 对单细胞表达矩阵做gsea分析
而且，这样的分析确实出现在几乎全部的单细胞数据分析文章里面，比如 2022年5月发表在CELL杂志的文章：《Quiescent cancer cells resist T cell attack by forming an immunosuppressive niche》，就有 we quantified a hypoxic signature comparing cells inside and outside QCC regions. （这个时候这个单细胞项目的数据分成了两个组，背后有一个差异分析的哦）
All sub- sets of conventional dendritic cells (cDCs) showed an enrichment in hypoxia-induced genes，主要是关心树突细胞，包括：
- • DC1 (identified by Clec9a, Xcr1, and Wdfy4),
- • DC2 (Itgax, Sirpa, and Cd209a),
- • mregDCs (Ccr7, Cd80, Cd200, and Cd247)
如下所示：
不同区域的各种树突细胞都缺氧富集
这个文章使用的方法就是Seurat包的AddModuleScore函数，如下所示：
Seurat包的AddModuleScore函数
那么，我们来复现一下，仍然是以大家熟知的pbmc3k数据集为例：
大家先安装这个数据集对应的包，并且对它进行降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释 ，而且每个亚群找高表达量基因，都存储为Rdata文件。
代码语言：javascript复制library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
table(Idents(sce))
DimPlot(sce,label = T)因为这个数据集已经进行了不同单细胞亚群的标记，就是完成了：人人都能学会的单细胞聚类分群注释 。我们就直接开始缺氧评分吧。
代码语言：javascript复制gs=list(
DC1 = c( 'Clec9a', 'Xcr1', 'Wdfy4'),
DC2 = c('Itgax', 'Sirpa', 'Cd209a'),
mregDCs= c('Ccr7', 'Cd80', 'Cd200', 'Cd247') ,
hypoxia=c('Hif1a', 'Slc2a1', 'Vegfa', 'Hmox1',
'Bnip3', 'Nos2', 'Mmp2', 'Sod3',
'Cited2', 'Ldha')
)
gs = lapply(gs, toupper)
sce = AddModuleScore(object = sce,gs)
colnames(sce@meta.data)
FeaturePlot(sce,'Cluster1')
VlnPlot(sce,'Cluster4')
ncol(sce@meta.data)
ac=sce@meta.data[,4,drop=F]
pheatmap::pheatmap(sce@meta.data[,8:11],
show_rownames = F,
annotation_row = ac)
因为我们这个示例数据是pbmc，里面的dc细胞非常少，跟原文不一样，所以看不出来很清晰的规律。需要把里面的三十多个DC细胞单独提取后继续降维聚类分群。
打分矩阵
要做出前面的条形图需要一个简单的转换，我随意写了一个代码，发现居然是 FCGR3A Mono 的缺氧打分略高于所有的其它免疫细胞，蛮让我费劲的，不知道是不是PBMC单细胞数据集的特异性呢，还是说，其它单细胞也是如此，代码是：
代码语言：javascript复制p=VlnPlot(sce,'Cluster4')
library(ggpubr)
df = aggregate(p$data$Cluster4,list(p$data$ident),median)
ggbarplot(df,'Group.1','x') coord_flip()出图如下所示；
FCGR3A Mono 的缺氧打分略高于所有的其它免疫细胞
其实关键是基因集合
虽然我演示的是缺氧评分，但是大家很容易根据我的示例代码进行复现，但是也很容易发散你的思维，理论上任何类似的基因列表，基因集合，都是可以同样的进行打分哦！
MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 包括至少几万 个基因列表，但是我们通常是不会每个都肉眼看或者说理解和记忆，只能是看看 hallmark里面的50个核心基因列表。
