文章预览

简介

首先想跟各位小伙伴说的是漂亮的图形是次要的，正确的分析才是第一位，不要代码出了结果就忘记了回头思考。本篇文章将为大家分享使用R语言实现GSEA的方法以及若干GSEA（Gene Set Enrichment Analysis，基因集富集分析）结果的可视化方法，文章的结构如下：

包含图形

GSEA介绍

1. GSEA原理
GSEA文献来源https://www.pnas.org/doi/10.1073/pnas.0506580102, 关于GSEA的原理解释可以移步 生信技能树详解GSEA
简单来说，GSEA的开发是为了更好的从生物学功能角度来解释基因组数据。之前所分享的GO、KEGG富集分析也是旨在解释不同实验分组之间的功能差异，但相信很多小伙伴在实践的时候也深有体会GO、KEGG方法得出的结果存在着诸多限制，比如技术噪声、可重复性差、依赖于生物背景、忽略有意义通路等等，GSEA在一定程度上打破了这些限制，并且可以更好的进行个性化分析。

2. GSEA准备
清楚基本原理之后，就可以上手分析了。进行GSEA需要准备如下：

1. 基因排序信息
   其实GSEA可以根据自己的需求定制基因排序信息，这里为了方便使用的是之前差异分析结果的logFC排序。需要注意的是这个基因排序在R语言中是一个带名字的向量。
   
2. 功能基因集
   一些已经预定义的功能基因集可以在MSigDB:http://www.gsea-msigdb.org/gsea/msigdb/index.jsp中下载，自定义基因集需要有较强的背景知识。以HALLMARK基因集为例，一般基因集都有两个常用版本，一种是SYMBOL ID版本，另一种是ENtrez ID版本，基因中所包含的基因名需要与第一部分文件中的基因名属于同一ID类型才可以进行后续分析，在本篇文章中对logFC基因排序命名时使用的是SYMBOL ID，这里也下载SYMBOl ID版本的hallmark gene sets。
3. GSEA实现方法
   这里小编常用的方法有两种，一种是clusterProfiler包中的GSEA函数，另一种则是fgsea包中的fgsea方法。
   

获得GSEA分析结果

1. 基因排序
清楚了上述知识之后我们就可以开始代码实操了，首先我们加载之前的差异分析结果，然后得到根据logFC的基因排序，以SYMBOL ID给排序向量命名。

rm(list = ls())

library(gggsea)
library(ggplot2)
library(GSVA)
library(enrichplot)
library('GSEABase')
library(fgsea)
library(clusterProfiler)

###定义排序
load('deg.RData')
alldiff <- deg[order(deg$logFC,decreasing = T),]
id <- alldiff$logFC
names(id) <- rownames(alldiff)

2. 加载预定义基因集
使用clusterProfiler中的read.gmt函数读取下载的gmt基因集文件。

gmtfile <- "./h.all.v7.4.symbols.gmt"
hallmark <- read.gmt(gmtfile)

3. 使用不同方法得到结果
使用两种方法得到各自结果并存储，具体函数的参数设置可参阅帮助文档，另外小编这里因为是连续的流程所以所有RData存放在同一个目录，实际项目的操作中最好分模块存储，以防出错。

gsea.re1<-clusterProfiler::GSEA(id,TERM2GENE=hallmark,verbose = T)
g1<-as.data.frame(gsea.re1)
g1<-subset(g1,p.adjust<0.05)
g1<-g1[order(g1$NES,decreasing = T),]


## 这里去掉了基因集前缀
hallmark$term <- gsub('HALLMARK_','',hallmark$term)
hallmark.list <- hallmark %>% split(.$term) %>% lapply( "[[", 2)

gsea.re2 <- fgsea(pathways = hallmark.list,
                  stats = id,
                  minSize=1,
                  maxSize=10000,
                  nperm=1000)
colnames(gsea.re2)
#
g2 <- gsea.re2[gsea.re2$padj<0.05,]
g2 <- g2[order(g2$NES,decreasing = T),]

save(gsea.re1,g1,gsea.re2,g2,file = 'gsea.RData')

可视化

走完上述流程，就可以来到最喜欢的可视化环节了，这里介绍三种不同的方法，应该足够简单的发表使用了。
1. enrichplot
enrichplot中的gseaplot2函数可以展示单个以及多个基因集的富集结果。我们先来看看单个基因集的结果，设置geneSetID参数为感兴趣基因集名称即可，这里以第一个作为示例。

load('gsea.RData')

library(ggsci)
col_gsea1<-pal_simpsons()(16)

num1=1
gseaplot2(gsea.re1,geneSetID = rownames(g1)[1:num1],
              title = "",#标题
              color = col_gsea1[1:num1],#颜色
              base_size = 14,#基础大小
              rel_heights = c(1, 0.2, 0.4),#小图相对高度
              subplots = 1:3,#展示小图
              pvalue_table = FALSE,#p值表格
              ES_geom = "line"#line or dot
          )

出图如下：

展示多个基因集，设置多个名称即可，比如这里展示4个，如果有特定基因集，定义完成后传入geneSetID参数即可。

num2=4
gseaplot2(gsea.re1,geneSetID = rownames(g1)[1:num2],
          title = "",#标题
          color = col_gsea1[1:num2],#颜色
          base_size = 14,#基础大小
          rel_heights = c(1, 0.2, 0.4),#小图相对高度
          subplots = 1:3,#展示小图
          pvalue_table = FALSE,#p值表格
          ES_geom = "line"#line or dot
)

使用illustrator再稍作编辑即可得到下图：2. fgsea
这篇cell也用到了fgsea中的结果，使用plotGseaTable函数即可得到。

plotGseaTable(hallmark.list[g2$pathway],
              id, 
              gsea.re2,gseaParam = 0.5,
              colwidths = c(0.5,0.2,0.1,0.1,0.1)
)

稍作编辑，出图如下：3. gggsea
最后，总会一直想找ggplot2的实现方法，发现其实是有包的，它就是gggsea，其github链接为：https://github.com/NicolasH2/gggsea。具体原理就不过多介绍了，主要是利用三个图层组合得到的，可以用来展示单基因集也可以用来展示多个基因集，只是多个基因集的结果是分面展示的形式，这里选择显著的4个基因集作为示例进行展示。

names(hallmark.list)
se_hall<-c(head(g2$pathway,2),tail(g2$pathway,2))

#3 所选中的基因集
sig1<-g2[g2$pathway%in%se_hall,]

hallmark.se<-hallmark.list[sig1$pathway]
df.new <- gseaCurve(id, hallmark.se, sig1)

#主要由三部分组成
# ggplot() + geom_gseaLine(df.new) + theme_gsea() #曲线
# ggplot() + geom_gseaTicks (df.new) + theme_gsea()#竖线
# ggplot() + geom_gseaGradient(df.new) + theme_gsea()#色块

pal_line<-pal_lancet()(9)#各曲线颜色

ggplot() + 
  geom_gsea(df.new,
            prettyGSEA=T,
            tickcolor='grey30',#竖线颜色
            ticksize=0.1,#线条粗细
            #colour='grey',#色块边框
            linecolor=pal_line,
            linesize=2,lty=1,
            #多个图形时使用，设定每行列图形个数
            nrow = 2,
            ncol = 2
  ) + 
  theme_bw()+
  theme(strip.text = element_text(size = 7,face = 'italic'),#增大分面标题字体
        strip.background = element_rect(fill = 'white'))+
  xlab(bquote(italic('Rank')))+ylab(bquote(italic('Enrichment Score')))+
  theme(axis.text.x = element_text(size=8,angle = -30,face = 'plain',hjust = 0.5),
        axis.text.y = element_text(size=8,angle = 0,face = 'plain',vjust = 0.5))

具体参数的作用就不多说了，有一部分注释代码，另外帮助文档也很好用，可以加深理解。出图如下：到这里本篇的分享就结束了，感谢大家的耐心！

总结

这篇文章和大家分享了用R语言进行GSEA的方法，以及GSEA几种可视化图形的制作，内容还是比较多，各位觉得有用的话可以点个小赞，非常感谢！公众号回复GSEA分析即可获得本篇示例数据和代码，希望可以有所收获~