gsea富集分析结果怎么看_“便携式”GSEA分析 Do GSEA without "GSEA software"

接触生物信息有段日子了，自己也发表了几篇数据挖掘的文章，感觉数据挖掘很大程度上来说是在做两件事：1.比较(异同) 2.富集(特征)。

举个例子来说，如果我们对control-treatment做差异表达分析，算法会给出的差异表达基因list，按照某个统计量，比如fold change，也就是control相较于treatment的变化倍数，从小到大排序，得到一个rank list,怎么从这个list中解读出有值得研究的信息呢？也就是如何富集呢？比较常见的办法就是选取fold change(变化倍数)最大值和最小值的基因按顺序一个一个查文献，如果恰好有课题需要的蛋白，就算是中奖了，然而更多的时候我们对于差异基因两眼一抹黑，不知道它位于什么通路，是否重要。那第二种方法就是GO/KEGG了，我们可以选择foldchange>2或者重要的不同就被掩盖掉了。GSEA(Gene Set Enrichment Analysis， 基因集富集分析)就是为了解决这样的问题而产生的。

基本概念

GSEA的想法其实很简单，就是看这些差异表达的基因在一些既往的通路中的富集情况。原假设是，某个通路的所有基因，在list中是随机的分布的，假如我们能观测到某个通路的所有基因富集到list中的一端，计算其富集程度，计算其统计显著性，如果小于某个cutoff，那么我们就可以拒绝原假设，认为该通路在list中富集，并且通过富集程度的打分，如果为正，则该通路倾向于在上调的基因中富集，如果为负，则该通路倾向于在下调的基因中富集。与差异表达不同，GSEA不需要指定明确的差异基因阈值，算法会根据实际整体趋势进行分析。

对于GSEA的结果，我们一般关注ES值(enrichment score)，峰出现在前端还是后端(ES值大于0在前端，小于0在后端)以及Leading-edge subset(即对富集贡献最大的部分，领头亚集)；在ES图中出现领头亚集的形状，表明这个功能基因集在某处理条件下具有更显著的生物学意义。具体的算法原理也不是很复杂，这里给大家简单介绍下。

一大堆公式是不是？哈哈，看似复杂，其实很简单啦，假如某个基因属于这一通路，那么它的fold change值除以fold change总和就是它的running ES，假如它不属于这个通路，那么它的running ES就等于-1除以(基因总数-这一通路基因的个数)，然后按照顺序不断累加，当达到最大值是就是这条通路的ES。

值得注意的是，官方GSEA软件默认使用signal2noise，还可以选择其他值，而今天我们的方法则采用fold change。

我们在之前的一篇文章中曾经介绍过GSEA软件的用法，可以查阅GSEA分析全攻略(视频教程)。在这篇文章中，我们通过视频详细介绍了GSEA官方软件的使用方法，但是在实际操作的过程中，很多同学从软件安装都开始报错，还有就是文件格式的整理，由于GSEA需要读入特定格式的表达谱文件和表型文件，很多同学在实际操作中可能要花费很多时间来制作这两个文件。按照上次的教程，我们通过GSEA软件可以分析得到如下图：

有没有什么快捷的方法，可以不要在制作各种不同格式的文件，也不要按照Java运行环境等，直接在R语言中运行GSEA呢？当然有！同时，在这个颜值取胜的时代，上面这种灰头土脸的图怎么能配得上你的CNS呢！有没有更好看的图呢？当然也有！

同样是GSEA，这张图有没有让人眼前一亮呢？讲完了大致原理，今天就来教给大家如何用快速做出上面的GSEA图。

实例操作

首先是导入数据，正如前面所说，我们的测试数据其实就是一个基因列表的Fold change值。不知道怎么产生这个列表的同学，可以参考我们的GEO数据挖掘系列【干货】数据分析教程-往期汇总。好了，我们的测试数据长这样：

然后是ID转换，将Gene symbol转换成通用的entrze ID，最后把输入数据转换成一个包含基因fold change的数值型向量。R代码如下：

 1# read in data 2df"DEG.txt",header=T) 3head(df) 4 5# annotation 6df.id"SYMBOL", toType = "ENTREZID",OrgDb = "org.Hs.eg.db") 7name.df"SYMBOL",all=F) 8head(name.df) 910# sort 11df.sortT),]12gene = df.sort$foldChange13names(gene) 14head(gene)

然后就是GSEA的过程，在这里我们使用clusterProfiler包进行分析。这个R包整合了常规KEGG富集和GO富集以及KEGG、GO进行GSEA富集的函数。如发表文章，请尊重作者工作并引用文献：

“Guangchuang Yu, Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.”

一个看似简单的工作，通过这个R包都可以快速解决，代码如下：

1# gsea2kegg "hsa",pvalueCutoff = 1)
3sortkeggT),]

这个计算过程要耗费我们几分钟的时间，大家耐心等待一下，计算的结果我们可以把结果输出到一个csv表格当中去。最后就是画图了，如何才能绘制出高级的GSEA矢量图呢？用gseaplot()函数就可以绘制一条或者多条GSEA的曲线。如下图所示：

                 

← 左右滑动图片 →

是不是感觉可以给文章增色不少呀？当然有的小伙伴可能会想，我是搞免疫方向的，我只想盯着免疫相关通路去研究，我该怎么去富集呢？这个问题其实也很好解决，有需要的小伙伴可以给客服留言，下一次我带你3分钟解决这个问题。

本期教程，关于通过R语言进行GSEA分析，我们讲到这里，需要完整代码的同学，关注“科研猫”，联系文末客服领取。

本期干货

·

- GSEA分析R代码 -

关注“科研猫”公众号，联系客服

胖雨小姐姐

or

折耳猫小姐姐

领取

更多科研新鲜资讯、文献精读和生物信息技能

请关注科研猫公众号

未经许可请勿随意转载，

版权事宜由上海辰明律师事务所提供法务支持