单细胞GSEA进阶分析之画图改造

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

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

• 计算感兴趣的基因集的富集得分(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))

（向左查看更多）

step3.获取GSEA图下方的排秩后的基因信息，用于画图。

rankPathway  <- 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)
  rank.data <- data.frame(rank = pathway, pathway = Pathway)
  return(list(rank.data))
})
rankPathway$GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
   rank                                         pathway
1     1 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
2     2 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
3     5 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
4    15 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
5    17 GO_POSITIVE_REGULATION_OF_INFLAMMATORY_RESPONSE
plot.rank <- do.call(rbind, rankPathway)
ticksSize <- 0.8
P2 <- ggplot() + geom_segment(data = plot.rank, mapping = aes(x = rank, y = -0.2, xend = rank, yend = -0.1, color = pathway), size = ticksSize) + 
  facet_wrap(~ pathway, ncol = 1) + 
  theme(
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    panel.spacing.y = unit(0, "mm"),
    panel.border = element_blank(), 
    axis.text = element_blank(), 
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank(), 
    panel.grid.major = element_blank(),
    panel.background = element_rect(fill = "transparent",colour = NA)) + labs(x = "", y = "") &NoLegend()

（向左查看更多）

将通路的动态ES得分图与基因排秩后的图拼接在一起。R语言里图片的拼接有很多方式，plot_grid()、par()、layout()、grid.arrange()，arrangeGrob()、marrangeGrob()等，感兴趣的伙伴儿可以都试一试。

p.list <- list()
p.list[[1]] <- P1
p.list[[2]] <- P2
patchwork::wrap_plots(plotlist = p.list, ncol = 1)

更多好玩的单细胞数据分析挖掘尽在百迈客云，让我们一起期待单细胞百迈客云的上线，当然强大的百迈客云还提供了52款APP实现高度自由个性化分析，118款工具、交互的个性化分析提供无上限分析实操，云平台课堂高效助您方案设计和信息分析（http://www.biocloud.net/），百迈客云成为您的私人信息分析平台，让基因分析更简单。

文：peach

排版：市场部

往期阅读

百迈客生物基于高通量测序技术、生物信息分析技术和生物云计算技术，为广大科研工作者提供以综合技术服务、生物云分析、三代高通量测序以及试剂、仪器等科研周边业务。

公司拥有Nanopore、PacBio、Illumina、Waters、10XGenomics等主流服务平台，以及基于云架构的生物云计算平台—百迈客云，提供涵盖人重外显子、三维基因组、单细胞与空间转录组、基因组组装、转录调控、微生物、群体遗传、质谱及表观遗传等研究方向的技术服务。目前百迈客云平台拥有200多款基因分析工具，分析结果可直接用于文章发表，更有近百部科研相关视频和8大基因数据库助力科研工作者深度数据挖掘。

自公司成立起先后在《Cell》、《Nature》、《Nature Genetics》、《Nature Communications》、《Plant Cell》等学术刊物发表论文数千篇，拥有国家发明专利技术40余项，软件著作权近200余项。

我们一直秉承”生物科技创新，服务社会，造福人民”的企业使命，致力于打造“生物科技创新中心”的发展愿景，让生物科技更快，更好的提高人类生活质量。