使用火山图呈现GSEA分析的结果

火山图本身很简单，最少需要两列数据就可以画出来，一列是横坐标，就是变化值logFC, 一列是纵坐标pvaue值。如果需要最终打出想要的标签，还需要一列gene symble，就是基因的名称。

还有一个重要的点是，火山图需要输入所有基因的信息，不仅仅是差异基因的。

当我一开始学习ggplot2作图的时候，就靠着网络，一句句检索，
画了人生第一张火山图。

library(ggplot2)
library(ggrepel)
load(file = "allDiff.Rda")
data <- allDiff
data$significant <- as.factor(data$adj.P.Val<0.05 & abs(data$logFC) > 1)
data$gene <- rownames(data)

ggplot(data=data, aes(x=logFC, y =-log10(adj.P.Val),color=significant)) +
  geom_point(alpha=0.8, size=1.2)+
  scale_color_manual(values =c("black","red"))+
  labs(x="log2 (fold change)",y="-log10 (adj.P.Val)")+
  theme(plot.title = element_text(hjust = 0.4))+
  geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
  geom_vline(xintercept = c(1,-1),lty=4,lwd=0.6,alpha=0.8)+
  #theme(legend.position='none')
  theme_bw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),   
        axis.line = element_line(colour = "black")) +
   geom_text_repel(data=subset(data, abs(logFC) > 3), 
                   aes(label=gene),col="black",alpha = 0.8)

自制火山图

当时觉得自己画的可以了，有一天我看到小Y在众筹火山图，是海涛实现的，没看代码的情况下，我觉得可以通过修改这个已有的代码实现。
改进的点在这里：

• 1.高表达的基因和低表达的基因可以用两个颜色表示。我考虑直接用geom_point分别画三个部分。
• 2.点的大小可以用另外一个变量表示，鉴于这里没有额外的变量，可以就用变化倍数来表示，最终越靠边，点越大。

改进后的代码如下：

library(ggplot2)
library(ggrepel)
library(dplyr)
load(file = "allDiff.Rda")
data <- allDiff
data$gene <- rownames(data)
## 仔细观察data数据，如果是你自己的数据，至少有三列，logFC，P.Value，gene
ggplot(data=data, aes(x=logFC, y =-log10(P.Value))) +
  ## 画出变化倍数少的点
  geom_point(data=subset(data,abs(data$logFC) <= 1),aes(size=abs(logFC)),color="black",alpha=0.1) +
  ## 画出上调的点
  geom_point(data=subset(data,data$P.Value<0.05 & data$logFC > 1),aes(size=abs(logFC)),color="red",alpha=0.1) +
  ## 画出下调的点
  geom_point(data=subset(data,data$P.Value<0.05 & data$logFC < -1),aes(size=abs(logFC)),color="darkgreen",alpha=0.1) +
  ## 画横线
  geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
  ## 画竖线
  geom_vline(xintercept = c(1,-1),lty=4,lwd=0.6,alpha=0.8)+
  ## 主题
  theme_bw()+
  ## 去边框
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),   
        axis.line = element_line(colour = "black"))+
  ## 修改坐标轴
  labs(x="log2 (fold change)",y="-log10 (q-value)")+
  ## 去掉注释
  theme(legend.position='none')+
  ## 打上标签
  geom_text_repel(data=subset(data, abs(logFC) > 3), aes(label=gene),col="black",alpha = 0.8)

修改后

这里面的点的大小还是用logFC这个已经显示在x轴的变量，我们可以考虑给点添加一个变量。
比如Y叔的神包clusterProfiler来做GSEA分析。
需要两个准备文件，

• 一个是geneList，这个可以用差异分析的结果来制作，Y叔的帖子讲过。
• 二是一个gene set集合，这个可以从board官网获取。

假设我们有了这两个文件，就可以十分便捷地做GSEA分析

## 读入hallmarks gene set
hallmarks <- read.gmt("h.all.v6.2.entrez.gmt")
# GSEA
y <- GSEA(geneList,TERM2GENE =hallmarks,pvalueCutoff = 1)

一般情况下，我们直接就会选一个富集分析图来展示，放到文章里面，但是这明显是不科学的，为什么值选取某一个gene set数据？为什么不给出其他所有集合的富集情况呢？就跟开题时专家经常问的问题那样，你研究的这个通路如果重要，究竟有多重要，跟你师兄那个比，谁更重要？
下面这行代码可以用分面的形式看激活的集合以及抑制的集合

dotplot(y,showCategory=12,split=".sign")+facet_grid(~.sign)

这时候你如果选择箭头所指的那个集合进行下一步研究，没有人会说你有问题，因为从图上看，这个是抑制通路里面，富集比率最高的。

但是，我自己觉得，这个图还有改进的空间，判断GSEA富集的结果，可以从富集分数NES，p值，以及GeneRatio三个维度去考虑。这时候我想到了用ggplot2画火山图。

首先我们要获取到刚才说的那几列,这时候我还增加了一列，就是gene set的名称。

### 选择需要呈现的来作图
yd <- data.frame(y)
### 增加GeneRatio这一列
enrich_num <- stringr::str_count(yd$core_enrichment,"/")+1
yd$GeneRatio <- as.numeric(enrich_num/yd$setSize)
yd$ID <- substring(yd$ID,10)
data <- yd

得到了数据，就可以依据火山图的原理作图了。

library(ggplot2)
library(ggrepel)
ggplot(data=data, aes(x=NES, y =-log10(pvalue))) +
  geom_point(data=subset(data,data$pvalue>=0.05),aes(size=GeneRatio),alpha=0.6)+
  ## 画出上调的点
  geom_point(data=subset(data,data$pvalue<0.05 & data$NES > 1),aes(size=abs(GeneRatio)),color="red",alpha=0.6) +
  ## 画出下调的点
  geom_point(data=subset(data,data$pvalue<0.05 & data$NES < -1),aes(size=abs(GeneRatio)),color="darkgreen",alpha=0.6) +
  ## 画横线
  geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
  ## 画竖线
  geom_vline(xintercept = c(1,-1),lty=4,lwd=0.6,alpha=0.8)+
  ## 主题
  theme_bw()+
  ## 去边框
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),   
        axis.line = element_line(colour = "black"))+
  ## 打上标签
  geom_text_repel(data=subset(data, abs(NES) > 2), 
                   aes(label=ID),col="black",alpha = 0.8,
                   ylim = c(0.5,1.2),
                  direction    = "x",
                  angle        = 90,
                  vjust        = 0)

GSEA富集分布图

因为我们选取的是board公司总结的50条hallmarks gene set，所以看起来你不怎么像火山图。不过如果尝试GO 这样包含很多数据集的gene set集合，就会很像。比如下面这个：

GO的GSEA富集分布图

回到我们本次的富集分布图，x坐标代表富集分数，y是pvalue的负log10，点的大小代表富集的基因数目占gene set总数的比例，红色是激活的，绿色是抑制的。激活的明显的是，

HALLMARK_INTERFERON_ALPHA_RESPONSE(有上角)

这时候可以画出对应gene set 富集分析图

library(enrichplot)
gseaplot2(y,"HALLMARK_INTERFERON_ALPHA_RESPONSE")

HALLMARK_INTERFERON_ALPHA_RESPONSE

对应的，我们可以画出抑制最明显的通路

HALLMARK_E2F_TARGETS

library(enrichplot)
gseaplot2(y,"HALLMARK_E2F_TARGETS")

HALLMARK_E2F_TARGETS

更多帖子，关注“果子学生信”微信公众号。