GSEA（Gene Set Enrichment Analysis）基因集合富集分析是一种计算方法，能够判定一个预先定义的基因集合（比如一个GO term或者一条通路）在两种生物学状态间是否呈现统计学上的显著的一致的差别。GSEA的主要有3个步骤：
-
基因集合（S）中的基因出现在排好的基因列表（L）的上端或者下端。
rank.png 计算感兴趣的基因集的富集得分(ES)，对排序后的基因列表，每遇到一个基因集S中有的基因，则增加其分值，如果遇到一个非基因集S的基因则降低其分值。
随机扰动，评估ES的显著性（p值）
最近在做GSEA分析，但是R包直接得到的GSEA的图片实在是丑了点，于是小编决定自己重新绘制GSEA图片。
step1. 首先导入GSEA结果数据，选取 top15 的通路进行画图。
library(fgsea)
library(customLayout)
library(msigdbr)
library(cowplot)
gsea.result <- readRDS("D:/gsea_result.Rds")
gsea.result.top <- gsea.result %>% filter(pval < 0.05) %>% top_n(n = 15, wt = NES)
step2. 获取每条通路动态的ES得分结果
ToPlot <- sapply(gsea.result.top$pathway, function(Pathway){
pathway <- geneSet[[Pathway]]
stats <- useGene
rnk <- rank(-stats)
ord <- order(rnk)
statsAdj <- stats[ord]
statsAdj <- sign(statsAdj) * (abs(statsAdj)^gseaParam)
statsAdj <- statsAdj/max(abs(statsAdj))
pathway <- unname(as.vector(na.omit(match(pathway, names(statsAdj)))))
pathway <- sort(pathway)
gseaRes <- calcGseaStat(statsAdj, selectedStats = pathway, returnAllExtremes = TRUE)
bottoms <- gseaRes$bottoms
tops <- gseaRes$tops
n <- length(statsAdj)
xs <- as.vector(rbind(pathway - 1, pathway))
ys <- as.vector(rbind(bottoms, tops))
toPlot <- data.frame(x = c(0, xs, n + 1), y = c(0, ys, 0), pathway = Pathway)
return(list(toPlot))
})
ToPlot$GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
x y pathway
1 0 0.0000000 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
2 0 0.0000000 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
3 1 0.1866981 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
4 1 0.1866981 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
5 2 0.3477841 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
6 4 0.3425893 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
7 5 0.4525161 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
plot.data <- do.call(rbind, ToPlot)
plot.data <- plot.data[plot.data$y > 0,]
P1 <- ggplot(plot.data, aes(x = x, y = y, group = pathway, color = pathway)) + geom_point(aes(fill = pathway), size = 0) +
geom_line(size = 0.8, show.legend = FALSE) +
guides(fill = guide_legend(nrow = 15),
colour = guide_legend(override.aes = list(shape = 15, size = 2))) +
labs(x = "", y = "Enrichment score", title = "GSEA") +
theme_bw() +
theme(panel.grid =element_blank()) +
theme(axis.ticks.x = element_blank()) +
theme(panel.border = element_blank()) +
theme(axis.line = element_line(size = 0.1),
axis.text.x = element_blank(),
plot.title = element_text(hjust = 0.5, size = 15),
legend.title = element_blank(),
legend.text = element_text(size = 4),
legend.margin = margin(0,0,0,0),
legend.box.margin = margin(0,0,0,0))
GSEA.png
