生信学习之通路富集二(通路打分):
生信学习之通路富集一(GO分析)
对富集通路打分:
AUC打分代码
。
// 用来看AUC对某一个通路的评分
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())
###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(harmony)
library(cowplot)
library(ggplot2)
load("scRNA1.rdata")
# BiocManager::install("AUCell")
library(AUCell)
# BiocManager::install("clusterProfiler")
library(clusterProfiler)
sc.id=sample(colnames(scRNA1),1500)
sc2=scRNA1[,sc.id]
# install.packages("doParallel")
# install.packages("doRNG")
cells_rankings <- AUCell_buildRankings(sc2@assays$RNA@data, splitByBlocks=TRUE, nCores=6, plotStats=TRUE)
cells_rankings
c2 <- read.gmt("c2.cp.kegg.v2023.1.Hs.symbols.gmt")
geneSets <- lapply(unique(c2$term), function(x){print(x);c2$gene[c2$term == x]})
names(geneSets) <- unique(c2$term)
?AUCell_calcAUC
cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings,nCores =1, aucMaxRank=nrow(cells_rankings)*0.1)
# install.packages("msigdbr")
library(msigdbr)
# msigdbr_species()
#
# m_df<- msigdbr(species = "Mus musculus", category = "C2", subcategory = "KEGG")
# fgsea_sets<- m_df %>% split(x = .$gene_symbol, f = .$gs_name)
#
#
# cells_AUC <- AUCell_calcAUC(fgsea_sets, cells_rankings,nCores =1, aucMaxRank=nrow(cells_rankings)*0.1)
# grep("OX",rownames(cells_AUC@assays@data$AUC),value = T)
geneSet <- "KEGG_OXIDATIVE_PHOSPHORYLATION"
aucs <- as.numeric(getAUC(cells_AUC)[geneSet, ])
sc2$AUC <- aucs
df<- data.frame(sc2@meta.data, sc2@reductions$umap@cell.embeddings)
colnames(df)
class_avg <- df %>%
group_by(seurat_clusters) %>%
summarise(
UMAP_1 = median(UMAP_1),
UMAP_2 = median(UMAP_2)
)
ggplot(df, aes(UMAP_1, UMAP_2)) +
geom_point(aes(colour = AUC)) + viridis::scale_color_viridis(option="A") +
ggrepel::geom_label_repel(aes(label = seurat_clusters),
data = class_avg,
size = 3,
label.size = 1,
segment.color = NA
)+ theme(legend.position = "none") + theme_bw() + facet_grid(.~orig.ident)
