|
有了基因集文件除了做scRNA分析|单细胞GSVA + limma差异分析-celltype分组？样本分组？GSVA分析，还可以计算每个细胞的目标基因集评分 。 方式有很多种，本文简单的介绍seurat的AddModuleScore函数 以及 AUCell包 2种计算方法。
载入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 单细胞矩阵#设置AssayDefaultAssay(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，也可以选择IDhuman_KEGG_Set = human_KEGG %>% split(x = .$gene_symbol, f = .$gs_name)#基因集是list
#随意选择其中一条通路，转为listWNT_features <- list(human_KEGG_Set$KEGG_WNT_SIGNALING_PATHWAY)
使用AddModuleScore函数计算KEGG_WNT_SIGNALING_PATHWAY基因集的打分 sce2 <- AddModuleScore(sce2, features = WNT_features, ctrl = 100, name = "WNT_features")head(sce2@meta.data)#这里就得到了基因集评分结果，但是注意列名为 WNT_features1colnames(sce2@meta.data)[16] <- 'WNT_Score'
得到的score 类似 在每个细胞中算出来的我们感兴趣的基因的表达均值。 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_PATHWAYgeneSet <- "KEGG_WNT_SIGNALING_PATHWAY"AUCell_auc <- as.numeric(getAUC(cells_AUC)[geneSet, ])#添加至metadata中sce2$AUCell <- AUCell_auchead(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绘图，生信图形可视化汇总）
|
