首先,让我们再简单回顾下GSEA的操作过程,(1)我们需要按顺序排列好的gene list用于分析,(2)需要参考基因集pre-defined gene set,那么这个从哪里来呢?这么跟大家说吧,在GSEA中富集出来的基因功能类或者基因集合都是提前定义好的,谁定义的呢?当然是GSEA官方或者一些权威数据库(比如KEGG通路数据库,Gene Ontology数据库等)。举个例子,有哪些基因隶属于p53 signaling pathway或者MAPK singaling pathway是不需要我们操心的,有专家已经帮我们定义好了。那么,如何去查看或者下载这些预定义好的基因集合呢?打开如下链接:
上面就是GSEA的官方网站的主页,如何找到我们感兴趣的基因集呢?如下图操作,找到MsigDB(Molecular Signatures Databases),这里面包含了GSEA定义的所有基因集,分为八个大类,分别由“H"和“C1”-“C7”开头,每个大类中都包含了哪些基因集在网页中有详细解释。比如C6是肿瘤相关,C7是免疫相关。
这个gmt文件里面是什么内容呢?我们可以把它用Excel表格打开看一眼。可以看到这个文件中每一行是一个基因集,第一列是基因集的名字,第二列是官网链接,后面的所有列就是该集合中包含的基因了(当然,这里是基因的Entrez ID)。
实例演练
现在给大家演示如何用特定的.gmt文件(基因集合)进行GSEA分析,我们还是用上次的数据集(没有测试数据的同学,可以在文末联系客服小姐姐领取)。我们先讲第一种方法:还是使用clusterProfiler包进行分析,代码如下,在计算出来的结果中,我们可以选择前10个富集出来结果绘制气泡图,同时选择第一个基因集合INTERFERON GAMMA RESPONSE绘制gsea图:
#clusterprofiler
pathway<-read.gmt("h.all.v7.0.entrez.gmt")
y <- GSEA(gene,TERM2GENE =pathway)
dotplot(y,showCategory=10)
gseaplot2(y,"HALLMARK_INTERFERON_GAMMA_RESPONSE",color = "red",pvalue_table = T)
如果你看倦了气泡图和GSEA图,在这里我们也给大家提供了另外一种方式,通过fgsea包选择特定基因集合进行分析,代码如下:
#fgsea
pathways<-gmtPathways("h.all.v7.0.entrez.gmt")
fgseaRes <- fgsea(pathways, gene, nperm=1000,minSize=15, maxSize=500)
topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n=10), pathway]
plot.new()
plotGseaTable(pathways[topPathwaysUp], gene, fgseaRes, gseaParam=0.5)
绘制出来的图形和结果如下所示,个人认为这种结果看上去更加简洁明了,我们最关心的NES、p value、adjusted p value以及Gene Ranks的缩略图都在图表中了。