|
有了基因集文件除了做scRNA分析|單細胞GSVA + limma差異分析-celltype分組？樣本分組？GSVA分析，還可以計算每個細胞的目標基因集評分 。 方式有很多種，本文簡單的介紹seurat的AddModuleScore函數(shù) 以及 AUCell包 2種計算方法。
載入R包，加載單細胞數(shù)據(jù)通過BiocManager::install的方式安裝一下AUCell包 ，后面會用到。library(Seurat) library(tidyverse)library(msigdbr) #提供基因集#BiocManager::install("AUCell",force = TRUE)library(AUCell)#讀取數(shù)據(jù)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")
數(shù)據(jù)準備：使用seurat的AddModuleScore函數(shù)進行基因集評分分析比較簡單，只需要準備（1）單細胞矩陣 和 （2）目標基因集（通路）的基因list即可。
1.1 單細胞矩陣#設(shè)置AssayDefaultAssay(sce2) <- "RNA"
1.2 基因集準備（1）基因比較少可以直接手動輸入，轉(zhuǎn)為list； （2）基因較多，通過文件的形式讀入，然后轉(zhuǎn)為list ，這里使用msigdbr包選擇的通路為例 # 手動輸入基因向量，并轉(zhuǎn)為list形式 WNT_features <- list(c( "gene1","gene2","gene3","gene4","gene5","gene6","gene7" #示例，需要真是的基因symbol))#直接使用文件中基因向量，并轉(zhuǎn)為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
#隨意選擇其中一條通路，轉(zhuǎn)為listWNT_features <- list(human_KEGG_Set$KEGG_WNT_SIGNALING_PATHWAY)
使用AddModuleScore函數(shù)計算KEGG_WNT_SIGNALING_PATHWAY基因集的打分 sce2 <- AddModuleScore(sce2, features = WNT_features, ctrl = 100, name = "WNT_features")head(sce2@meta.data)#這里就得到了基因集評分結(jié)果，但是注意列名為 WNT_features1colnames(sce2@meta.data)[16] <- 'WNT_Score'
得到的score 類似 在每個細胞中算出來的我們感興趣的基因的表達均值。 AUCell使用曲線下面積來計算輸入基因集的一個關(guān)鍵子集是否在每個細胞的表達基因中富集。 AUCell需要使用AUCell_buildRankings函數(shù)進行排序（Building the rankings），然后使用 AUCell_calcAUC 函數(shù) 計算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 函數(shù)得到上述2個步驟相同的結(jié)果。#cells_AUC2 <- AUCell_run(exprMatrix, geneSets)
提取目標通路的結(jié)果，也就是提取目標基因集所在的行 #提取KEGG_WNT_SIGNALING_PATHWAYgeneSet <- "KEGG_WNT_SIGNALING_PATHWAY"AUCell_auc <- as.numeric(getAUC(cells_AUC)[geneSet, ])#添加至metadata中sce2$AUCell <- AUCell_auchead(sce2@meta.data)
得到基因集分數(shù)后，可以使用seurat內(nèi)置的函數(shù)進行可視化，或者提取數(shù)據(jù)使用ggplot2 或者 擴展包進行可視化展示。
4.1 seurat的繪圖函數(shù)和其他信息一樣，直接使用seurat的一些函數(shù)繪制 小提琴圖，點圖，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圖提取基因集評分結(jié)果與umap的坐標 ，使用ggplot2 繪制umap點圖，將基因集評分映射到umap圖中 。 library(ggrepel)#提取umap坐標數(shù)據(jù)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學(xué)umap圖| ggplot2 繪制umap圖，坐標位置 ，顏色 ，大小還不是你說了算
4.4 繪制熱圖如果展示多條通路的話，可以使用熱圖的方式 library(pheatmap)aucMat <- getAUC(cells_AUC)pheatmap(aucMat[1:20,], show_colnames = F, scale = "row")
精心整理（含圖PLUS版）|R語言生信分析，可視化（R統(tǒng)計，ggplot2繪圖，生信圖形可視化匯總）
|
