scRNA分析|使用AddModuleScore 和 AUcell进行基因集打分，可视化
scRNA分析|使用AddModuleScore 和 AUcell进行基因集打分，可视化
有了基因集文件除了做scRNA分析|单细胞GSVA + limma差异分析-celltype分组？样本分组？GSVA分析，还可以计算每个细胞的目标基因集评分 。
方式有很多种，本文简单的介绍seurat的AddModuleScore函数 以及 AUCell包 2种计算方法。
一 载入R包 数据
载入R包，加载单细胞数据
通过BiocManager::install的方式安装一下AUCell包 ，后面会用到。
library(Seurat)
library(tidyverse)
library(msigdbr) #提供基因集
#BiocManager::install("AUCell",force = TRUE)
library(AUCell)
#读取数据
load("sce.anno.RData")
#载入颜色
colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",
"#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0",
"#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE",
"#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")数据准备：使用seurat的AddModuleScore函数进行基因集评分分析比较简单，只需要准备（1）单细胞矩阵 和 （2）目标基因集（通路）的基因list即可。
1.1 单细胞矩阵
#设置Assay
DefaultAssay(sce2) <- "RNA"1.2 基因集准备
（1）基因比较少可以直接手动输入，转为list；
（2）基因较多，通过文件的形式读入，然后转为list ，这里使用msigdbr包选择的通路为例
# 手动输入基因向量，并转为list形式
WNT_features <- list(c(
"gene1","gene2","gene3","gene4","gene5","gene6","gene7" #示例，需要真是的基因symbol
))
#直接使用文件中基因向量，并转为list形式
human_KEGG = msigdbr(species = "Homo sapiens",
category = "C2",
subcategory = "KEGG") %>%
dplyr::select(gs_name,gene_symbol)#这里可以选择gene symbol，也可以选择ID
human_KEGG_Set = human_KEGG %>% split(x = .$gene_symbol, f = .$gs_name)#基因集是list
#随意选择其中一条通路，转为list
WNT_features <- list(human_KEGG_Set$KEGG_WNT_SIGNALING_PATHWAY)二 AddModuleScore 打分
使用AddModuleScore函数计算KEGG_WNT_SIGNALING_PATHWAY基因集的打分
sce2 <- AddModuleScore(sce2,
features = WNT_features,
ctrl = 100,
name = "WNT_features")
head(sce2@meta.data)
#这里就得到了基因集评分结果，但是注意列名为 WNT_features1
colnames(sce2@meta.data)[16] <- 'WNT_Score'得到的score 类似 在每个细胞中算出来的我们感兴趣的基因的表达均值。
三 AUCell 计算
AUCell使用曲线下面积来计算输入基因集的一个关键子集是否在每个细胞的表达基因中富集。
AUCell需要使用AUCell_buildRankings函数进行排序（Building the rankings），然后使用 AUCell_calcAUC 函数 计算AUC值（Calculate the Area Under the Curve ) 。
#cells_AUC <- AUCell_run(exprMatrix, geneSets)
cells_rankings <- AUCell_buildRankings(sce2@assays$RNA@data,splitByBlocks=TRUE)
cells_AUC <- AUCell_calcAUC(human_KEGG_Set, cells_rankings,
aucMaxRank=nrow(cells_rankings)*0.1)
#也可以直接使用 AUCell_run 函数得到上述2个步骤相同的结果。
#cells_AUC2 <- AUCell_run(exprMatrix, geneSets)
提取目标通路的结果，也就是提取目标基因集所在的行
#提取KEGG_WNT_SIGNALING_PATHWAY
geneSet <- "KEGG_WNT_SIGNALING_PATHWAY"
AUCell_auc <- as.numeric(getAUC(cells_AUC)[geneSet, ])
#添加至metadata中
sce2$AUCell <- AUCell_auc
head(sce2@meta.data)
四 可视化
得到基因集分数后，可以使用seurat内置的函数进行可视化，或者提取数据使用ggplot2 或者 扩展包进行可视化展示。
4.1 seurat的绘图函数
和其他信息一样，直接使用seurat的一些函数绘制 小提琴图，点图，feature 图等等。
VlnPlot(sce2,features = 'WNT_Score',
pt.size = 0,group.by = "sample")
4.2 ggpubr绘制
1）使用ggpubr包，自定义颜色绘制箱线图
ggboxplot(sce2@meta.data, x="sample", y="WNT_Score", width = 0.6,
color = "black",#轮廓颜色
fill="sample",#填充
palette = colour,
#palette =c("#E7B800", "#00AFBB"),#分组着色
xlab = F, #不显示x轴的标签
bxp.errorbar=T,#显示误差条
bxp.errorbar.width=0.5, #误差条大小
size=1, #箱型图边线的粗细
outlier.shape=NA, #不显示outlier
legend = "right") #图例
2）使用ggpubr，按照文献配色绘制小提琴图
#使用npg颜色
ggviolin(df, x="celltype", y="AUC", width = 0.6,
color = "black",#轮廓颜色
fill="celltype",#填充
palette = "npg",
add = 'mean_sd',
xlab = F, #不显示x轴的标签
bxp.errorbar=T,#显示误差条
bxp.errorbar.width=0.5, #误差条大小
size=1, #箱型图边线的粗细
outlier.shape=NA, #不显示outlier
legend = "right")
4.3 绘制umap图
提取基因集评分结果与umap的坐标 ，使用ggplot2 绘制umap点图，将基因集评分映射到umap图中 。
library(ggrepel)
#提取umap坐标数据
umap <- data.frame(sce2@meta.data, sce2@reductions$umap@cell.embeddings)
head(umap)
#1）计算每个celltype的median坐标位置
cell_type_med <- umap %>%
group_by(Anno) %>%
summarise(
UMAP_1 = median(UMAP_1),
UMAP_2 = median(UMAP_2)
)
#2）使用ggrepel包geom_label_repel 添加注释
ggplot(umap, aes(UMAP_1, UMAP_2)) +
geom_point(aes(colour = AUCell)) +
ggrepel::geom_label_repel(aes(label = Anno),fontface="bold",data = cell_type_med,
point.padding=unit(0.5, "lines")
) + theme_bw()
一些美化方式可以参考跟SCI学umap图| ggplot2 绘制umap图，坐标位置 ，颜色 ，大小还不是你说了算
4.4 绘制热图
如果展示多条通路的话，可以使用热图的方式
library(pheatmap)
aucMat <- getAUC(cells_AUC)
pheatmap(aucMat[1:20,], show_colnames = F, scale = "row")
◆ ◆ ◆ ◆ ◆
精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）
腾讯云开发者
