GSEA分析之画图

作者: 疯子眼中的世界yia | 来源:发表于2021-03-09 11:51 被阅读0次

GSEA分析之画图
专题：富集分析
GASE相关知识（Linux）
GSEA的分析汇总-转载
单细胞之富集分析-5：一张图展示多组富集分析结果
单细胞之富集分析-6：PROGENy
单细胞之富集分析-7：Pseudobulk differenti
GSEA基因集富集分析
GSEA分析
第三步：简单的GEO数据的KEGG分析和GSEA分析

GSEA（Gene Set Enrichment Analysis）基因集合富集分析是一种计算方法，能够判定一个预先定义的基因集合（比如一个GO term或者一条通路）在两种生物学状态间是否呈现统计学上的显著的一致的差别。GSEA的主要有3个步骤：

基因集合（S）中的基因出现在排好的基因列表（L）的上端或者下端。

rank.png
计算感兴趣的基因集的富集得分(ES)，对排序后的基因列表，每遇到一个基因集S中有的基因，则增加其分值，如果遇到一个非基因集S的基因则降低其分值。
$ES(s) = max(\sum_{i=1}^jXi)$ $注：1<j<N$
$Xi = -\sqrt{\frac{G}{N-G}}$
随机扰动，评估ES的显著性（p值）
$pvalue = \frac{p(随机)>p(s)}{扰动次数(1000)}$

最近在做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