單細(xì)胞測(cè)序數(shù)據(jù)也可以做gsea，步驟跟用RNAseq的數(shù)據(jù)差不多藏澳，主要是要用到差異基因并且根據(jù)Fold change來(lái)排序。
library(msigdbr)
library(fgsea)
library(dplyr)
library(ggplot2)
選擇自己數(shù)據(jù)的物種以及要做的GSEA的數(shù)據(jù)庫(kù)類型
##查看物種的數(shù)據(jù) msigdbr_show_species()
m_df<- msigdbr(species = "Homo sapiens", category = "H")
#將m_df的基因與通路取出并改成一個(gè)通路對(duì)應(yīng)相應(yīng)基因的格式
fgsea_sets<- m_df %>% split(x = .$gene_symbol, f = .$gs_name)
這里選用seurat求差異基因耀找，并將差異基因按顯著性排序
#每一個(gè)細(xì)胞類型的GSEA按顯著性進(jìn)行降序排序
gesa_TvsC_allgenes<-FindMarkers(seurat.combined, ident.1 = "TCF7KO", ident.2 = "NC", verbose = FALSE,test.use ="roc",logfc.threshold = 0.01,only.pos =F)
gesa_TvsC_allgenes$gene<-rownames(gesa_TvsC_allgenes)
gsea_genes<-gesa_TvsC_allgenes %>%
arrange(desc(myAUC), desc(avg_diff)) %>%
dplyr::select(gene,avg_diff)
library(tibble)
ranks <- deframe(gsea_genes)
fgseaRes <- fgsea(pathways = fgsea_sets,
stats = ranks ,
minSize=5,
maxSize=500,
nperm=10000)
對(duì)結(jié)果進(jìn)行整理
library(dplyr)
fgseaResTidy <- fgseaRes %>%
as_tibble() %>%
arrange(desc(NES))
fgseaResTidy %>%
dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>%
arrange(padj) %>%
head()
ggplot(fgseaResTidy %>% filter(padj < 0.5) %>% head(n= 15), aes(reorder(pathway, NES), NES)) +
geom_col(aes(fill= NES < 0)) +
coord_flip() +
labs(x="Pathway", y="Normalized Enrichment Score",
title="Hallmark pathways NES from GSEA") +
theme_minimal()
plotEnrichment(fgsea_sets[["HALLMARK_APICAL_JUNCTION"]],
ranks) + labs(title="HALLMARK_APICAL_JUNCTION")
