史上最全GSEA可视化教程，今天让你彻底搞懂GSEA！

探索生信之美，解构每段代码背后的故事

不知道你有没有过这样的体验，好不容易收集了心心念念已久的样品送去测序，从送出去的那一刻就开始想象自己即将拿到一堆宝藏数据，按照36策中讲到的“差异基因、正反回复、细胞动物”或者“多元分子交互”打出一套，就算不是飞黄腾达发表一篇nature、science，最起码3、5分保底吧。然而，现实总是残酷的，数据分析显示实验组和对照组之间竟然没有显著的差异基因！瞬间五雷轰顶，这该怎么办呢？别急，今天介绍的基因集富集分析（Gene Set Enrichment Analysis，GSEA）就可以作为一个替代方案，挖掘不同样品不同状态潜在的分子机制。

▌话说GSEA

通俗来讲，GSEA首先将基因在两种样品中的差异表达程度（如logFC）或者表型相关度进行排序（并不是计算差异基因），然后判断来自功能注释等预定义的基因集或自定义的基因集是否倾向于落在有序列表的顶部或底部，也就是说预定义基因集中大部分的基因是否在样品中高表达或者低表达，从而判断功能注释基因集在不同的样品中有没有发生显著变化。

▌GSEA结果解读

 第1部分是Enrichment Score的折线图 

横轴排序后的基因，纵轴为对应的Running ES, 在折线图中出现的峰值就是这个基因集的富集分数（Enrichment Score，ES）。ES是从排序后的表达基因集的第一个基因开始，如果排序表达基因集中的基因出现在功能基因数据集中则加分，反之则减分。正值说明在顶部富集，峰值左边的基因为核心基因，负值则相反。

 第2部分为基因位置图 

黑线代表排序后表达数据集中的基因存处于当前分析的功能注释基因集的位置，红蓝相间的热图是表达丰度排列，红色越深的表示该位置的基因logFC越大 ，蓝色越深表示logFC越小。如果研究的功能注释基因集的成员显著聚集在表达数据集的顶部或底部，则说明功能基因数据集中的基因在数据集中高表达或低表达，若随机分配，则说明表达数据集与该通路无关。

 第3部分表示每个基因对应的信噪比（Signal2noise）

以灰色面积图显展示。灰色阴影的面积比，可以从整体上反映组间的Signal2noise的大小。

除了会看图形，我们也要关注一些主要的统计学指标：

NES（校正后的富集分数）

FDR（错误发现率）、p值：FDR<0.25,p<0.05且|NES|>1则表示结果有意义。

▌GSEA的优点

GSEA是个非常强大的富集分析方法，可以针对多种数据库中的数据进行GSEA分析，包括常见的GO数据库，KEGG数据库，Reactome数据库以及MSigDb数据库或其他自定义数据集。并且与传统的GO/KEGG只对显著的差异基因进行功能注释相比，GSEA

1. 考虑了基因的表达水平，不需要指定明确的差异基因阈值，检验的是基因集合而非单个基因的表达变化，算法会根据实际数据的整体趋势进行分析。

2. GO和KEGG只能定位到功能，而GSEA通过检验预定义的基因集的基因在某种状态下表达水平而提示该通路是否被激活还是抑制。

3. 以待测功能基因集为对象来进行检验,使得检验结果针对性和灵敏性提高。

▌GSEA分析实战

讲了这么多的优点，那GSEA分析该怎么做呢？有些经验的小伙伴可能会说GSEA软件，但是GSEA软件需要输入特定的文件格式，图片也只能是png格式，远远达不到发表级图片的要求。然而，R语言的clusterprofile包几条代码轻松解决所有的难题，还有各种各样可视化方式，你不来看一看吗？

#清空环境变量rm(list = ls())options(stringsAsFactors = F)#加载相关的包library(data.table)library(clusterProfiler)library(dplyr)library(org.Hs.eg.db)#加载所需物种基因组注释信息，可有20个物种选择，http:///packages/release/BiocViews.html#___OrgDblibrary(ggplot2)library(enrichpl)准备文件# 读入差异分析后文件，包含基因列（gene symbol或entrizid）和logFC列data <- read.table('D:/R study/data_input.txt',sep='\t',header=T,check.names=F)#查看文件内容head(data)# SYMBOL logFC#1 FAM182B -5.612055#2 MEP1B -5.134941#3 TMEM167B -3.383565#4 GRIK5 -4.903500#5 XKR9 -5.718563#6    NDST2 -6.843956#本次文件基因名为gene symbol，将gene symbol转换为ENTREZID#如果你的文件是entrizid则跳过该条语句gene.df <- bitr(data$SYMBOL, fromType = 'SYMBOL', toType = c('ENTREZID'), OrgDb = org.Hs.eg.db)#查看前6行数据head(gene.df) SYMBOL ENTREZID#1 FAM182B 728882#2 MEP1B 4225#3 TMEM167B 56900#4 GRIK5 2901#5 XKR9 389668#6    NDST2     8509data2=merge(data,gene.df,by.y='SYMBOL')#合并数据geneList<-data2$logFC #第二列可以是foldchange，也可以是logFCnames(geneList)=data2$ENTREZID #使用转换好的ENTREZIDgeneList=sort(geneList,decreasing = T) #从高到低排序接下来开始正式的GSEA分析啦#使用GO数据库进行GSEA分析ego <- gseGO( geneList = geneList, OrgDb = org.Hs.eg.db, ont = 'CC', nPerm = 1000, minGSSize = 100, maxGSSize = 500, pvalueCutoff = 0.05, verbose = FALSE, seed = FALSE, by = 'fgsea')##geneList是前面整理好的排序后的ENTREZID##OrgDb选择物种注释的数据库##ont同enrichGO分析，可选3大分类BP、CC、MF或ALL##nPerm置换检验的次数，默认为1000##minGSSize富集到某个条目的最小包含基因数，如果基因数小于该值则这个条目将被过滤掉，默认为10##maxGSSize富集到某个条目的最大包含基因数，如果基因数大于该值则这个条目将被过滤掉，默认为500##verbose是否输出提示信息##seed是否使结果具有可重复性##by选择使用的统计学方法，默认为fgseahead(ego)[1:6][1:6]## ID Description setSize enrichmentScore NES##GO:0098687 GO:0098687 chromosomal region 320 0.4187365 1.667091##GO:0000793 GO:0000793 condensed chromosome 209 0.4459184 1.713171##GO:0000775 GO:0000775 chromosome, centromeric region 191 0.4502879 1.722934##GO:0000776 GO:0000776 kinetochore 133 0.4537928 1.677579##GO:0000779 GO:0000779 condensed chromosome, centromeric region 118 0.5091699 1.850798##GO:0072562 GO:0072562 blood microparticle 114 0.4655046 1.680016## pvalue##GO:0098687 0.001052632##GO:0000793 0.001113586##GO:0000775 0.001123596##GO:0000776 0.001154734##GO:0000779 0.001187648##GO:0072562 0.001197605##保存结果write.table(ego,file='goGSEA.txt',sep='  ',quote=F,row.names = F)#使用KEGG数据库进行GSEA分析kk <- gseKEGG(geneList, organism = 'hsa',pvalueCutoff = 0.2, nPerm = 1000, minGSSize = 10, maxGSSize = 500, verbose = TRUE, seed = FALSE, by = 'fgsea')##参数同上##按照enrichment score从高到低排序，便于查看富集通路kk_order<-kk[order(kk$enrichmentScore,decreasing=T)]head(kk_order)[1:6][1:6]## ID Description setSize enrichmentScore NES pvalue##hsa03030 hsa03030 DNA replication 36 0.5600885 1.664196 0.008415147##hsa04740 hsa04740 Olfactory transduction 287 0.5005390 1.968686 0.001054852##hsa04110 hsa04110 Cell cycle 124 0.4983236 1.807226 0.001165501##hsa04610 hsa04610 Complement and coagulation cascades 84 0.4868227 1.685794 0.004944376##hsa03008 hsa03008 Ribosome biogenesis in eukaryotes 78 0.4607935 1.581720 0.009962640##hsa05206 hsa05206 MicroRNAs in cancer 160 0.3695695 1.371388 0.019252548##保存结果write.table(kk_order,file='keggGSEA.txt',sep='  ',quote=F,row.names = F)#使用MsigDb数据库进行GSEA分析##下载离线GMT文件，可以根据你的需要下载，也可以是全部，网址如下##https://www./gsea/downloads.jspall<-read.gmt('msigdb.v7.4.entrez.gmt') #读gmt文件gsea<-GSEA(geneList,TERM2GENE = all) ##参数大致同上，不同是多了TERM2GENE和TERM2NAME##TERM2GENE数据集中条目名称TERM与基因的对应关系,一般是两列的数据框形式##TERM2NAME条目名称与名称的对应关系，可用于个性化的GSEA分析，如系定义的基因集，可以为空 write.table(gsea,file='GSEA.txt',sep='  ',quote=F,row.names = F)#使用reactome数据库进行GSEA分析library(ReactomePA)pathway <- gsePathway(geneList)write.table(pathway,file='pathwayGSEA.txt',sep=' ',quote=F,row.names = F)

▌GSEA分析可视化

除了上方介绍的经典的GSEA结果展示图以外，GSEA和GO、KEGG一样，都有很多让人眼花缭乱的可视化方法，让我们来看一看把。

#新建文件夹储存图片dir.create('GSEA图片集')#最简单的点图dotplot(kk,x='GeneRatio',showCategory = 10)##x为横坐标的变量，除了GeneRatio' or 'Count'以外，还有~GeneRatio/BgRatio，~-log(qvalue)##showCategory为展示前10个条目

#分面点图激活和抑制dotplot(kk,split='.sign') facet_grid(~.sign)

# 通过改变scales改变Y轴展示dotplot(kk,split='.sign') facet_wrap(~.sign,scales = 'free')

# 条形图##加载需要的R包library(forcats)library(ggstance)##arrange去给NES的绝对值从大到小排序y <- arrange(x, abs(NES)) %>% group_by(sign(NES))##绘图ggplot(y, aes(NES, fct_reorder(Description, NES), fill=qvalues), showCategory=10) geom_barh(stat='identity') scale_fill_continuous(low='red', high='blue') theme_minimal() ylab(NULL)

#网络图cnetplot(kk, showCategory = 3,foldChange = geneList,node_label='category',          colorEdge = TRUE)##colorEdge不同的term展示不同的颜色##node_label 可以选择label的内容，category', 'gene', 'all' and 'none', 默认 'all'.##也可以把想要展示的通路名称赋值给y#，然后showCategory = y##categorySize='pvalue'##circular=true以环形展示

#网络图library(enrichplot)ego2 <- pairwise_termsim(ego)##直接运行emapplot会报错，要加上这条语句emapplot(ego2,showCategory=10,color='pvalue',cex_category=1,layout='kk')##cex_category调节节点的大小##showCategory展示显示的条目个数

#高级韦恩图upsetplot(kk)

# 经典的gsea结果图library(enrichplot)gseaplot2(gsea,1,color='red',pvalue_table = T,title='',base_size=10,ES_geom='line') # 按第一个条目做二维码图，并显示p值#color是enrichment score线的颜色#base_sizexy轴标题字的大小#ES_geom是enrichment score线用线或点显示

##当然也可以在一张图上展示多个条目gseaplot2(gsea,1:3,color='red',pvalue_table = T,title=gsea$Description[1],base_size=10,ES_geom='line')

#山脊图ridgeplot(kk,fill = 'pvalue',5) scale_fill_continuous(type = 'viridis')

好了，今天的分享就到这结束啦，换上自己的数据趁热赶紧试一下吧！发表文章时千万别忘记引用R包文献哦。

—END—

撰文丨念   念

排版丨四金兄

编辑丨小雪球

欢迎大家关注解螺旋生信频道-挑圈联靠公号~