8f37c46e064b71c7fc55ec17806e6e2.png
1.R包的下载
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
# To support paralell execution:
BiocManager::install(c("doMC", "doRNG","doSNOW"))
# For the main example:
BiocManager::install(c("mixtools", "GEOquery", "SummarizedExperiment"))
# For the examples in the follow-up section of the tutorial:
BiocManager::install(c("DT", "plotly", "NMF", "d3heatmap", "shiny", "rbokeh",
"dynamicTreeCut","R2HTML","Rtsne", "zoo"))
browseVignettes("AUCell")
2.需要两个文件，一个是seruat的表达矩阵，一个是基因集（MsigDB数据库)
#加载R包
library(AUCell)
library(clusterProfiler)
library(ggplot2)
library(Seurat)
library(SeuratData)
library(msigdbr)
library(patchwork)
rm(list=ls())
#导入数据
#导入seurat对象
scRNA <- readRDS("scRNAsub.rds") #这个是我自己的单细胞数据
dim(scRNA)
cells_rankings <- AUCell_buildRankings(scRNA@assays$RNA@data) # 关键一步 细胞的表达矩阵
#制作基因集
#选择GO-BP
a=msigdbr(species = "Homo sapiens",category = "C5",subcategory = c("BP"))
a = subset(a,select =c("gs_name","gene_symbol")) %>% as.data.frame()
#选择GO-MF
b=msigdbr(species = "Homo sapiens",category = "C5",subcategory = c("MF"))
b = subset(b,select =c("gs_name","gene_symbol")) %>% as.data.frame()
#将二者合并
c <- rbind(a,b)
genesets = split(c$gene_symbol,c$gs_name)
#这样就制作好了关于GO BP与MF的基因集（唉，从这里就可以看出我的水平了，一言难尽......）
3.计算AUC： 耗时长
cells_AUC <- AUCell_calcAUC(genesets, cells_rankings, aucMaxRank=nrow(cells_rankings)*0.1)
4.找一些通路（该找哪些通路呢？） 用grep函数找
length(rownames(cells_AUC@assays@data$AUC)) #查看一下一共有多少条
data.frame(grep("METABOLIC",rownames(cells_AUC@assays@data$AUC),value = T)) #找与代谢相关的term
#选择其中一条，进行plot ##set gene set of interest here for plotting
geneSet <- "GOBP_VITAMIN_B6_METABOLIC_PROCESS"
aucs <- as.numeric(getAUC(cells_AUC)[geneSet, ]) #提取这个通路在每一个细胞的得分
scRNA$AUC <- aucs #将得分添加入scRNA（seruat）对象
saveRDS(scRNA,"scRNAsub.rds") #保存数据
#选择细胞展示的维度
df<- data.frame(scRNA@meta.data, scRNA@reductions$umap@cell.embeddings) #选择用UMAP维度看 也可以选择TSNE
head(df)
#我们看到每个细胞现在都加上AUC值了，下面做一下可视化。
5.plot
#做一个注释文件 就好比要在什么层次对细胞进行注释
class_avg <- df %>%
group_by(celltype) %>% #这里可以改成cluster seurat_clusters/或者其他的annotation
summarise(
UMAP_1 = median(UMAP_1),
UMAP_2 = median(UMAP_2)
)
#通过ggplot画图
ggplot(df, aes(UMAP_1, UMAP_2)) +
geom_point(aes(colour = AUC)) + viridis::scale_color_viridis(option="A") +
ggrepel::geom_label_repel(aes(label = celltype),
data = class_avg,
size = 6,
label.size = 0,
segment.color = NA
)+ theme(legend.position = "none") + theme_bw()
6.references：
https://www.jianshu.com/p/784a04c49873
https://www.jianshu.com/p/2fb20f44da67
