scRNA分析|使用AddModuleScore 和 AUcell进行基因集打分，可视化-阿里云开发者社区

一载入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")

AGS无服务化分析基因数据 - mutect2 肿瘤样本分析

通过调用AGS的远程任务，可以完成一序列的基因数据的二级分析，不需要申请和持有云计算资源，就可以完成对海量数据的批量处理，目前可以支持人类全基因组，外显子，基因比对，宏基因组比对，Somatic胚系变异发现等业务场景的加速和低成本处理。通过AGS调用mutect2任务来检测体细胞短突变, 短突变包括单核苷酸（SNV）以及插入和缺失（Indel)的改变。本文介绍如何通过AGS分析肿瘤样本。场景实践 - 基于阿里云PAI机器学习平台使用时间序列分解模型预测商品销量【医学搜索Query相关性判断】赛题及baseline解读【算法实战】12. 利用 Apriori 算法进行关联分析机器学习基础与回归算法南瓜书《机器学习公式推导》【算法实战】7. 集成方法-随机森林

一 载入R包 数据