GSEA分析一文就够（单机版+R语言版）

生信技能树

发布于 2018-03-29 16:04:18

5.9K0

发布于 2018-03-29 16:04:18

文章被收录于专栏：生信技能树生信技能树

通过前面的讲解，我们顺利的了解了GEO数据库以及如何下载其数据，得到我们想要的表达矩阵，但，这只是分析的开始，最经典的分析就是GSEA了，看看基因全局表达量的变化是否有某些特定的基因集合的倾向性。

历史目录： 解读GEO数据存放规律及下载，一文就够解读SRA数据库规律一文就够从GEO数据库下载得到表达矩阵一文就够

GSEA软件的用法

这个是java软件，所以各个电脑操作系统都可以很容易安装及使用。我在生信菜鸟团博客也手把手讲解了详细操作过程，这里就不再赘述咯：

用GSEA来做基因集富集分析 http://www.bio-info-trainee.com/1282.html
批量运行GSEA，命令行版本 http://www.bio-info-trainee.com/1334.html

GSEA的原理

首先对每个样本里面的基因的表达值在样本内部进行排序，本质是是根据该基因在两个group之间的差异来排序！但是差异如何量化，就有多种方法了，可以是Signal2Noise 值，或者是Ttest值，或者是fold change，logFC等等。默认的，GSEA会根据signal-to-noise metric 来对基因进行排序。但是也可以选择其它metric。

如果是自己已经排序好了的基因，也可以直接拿来做GSEA分析了，见： GSEAPreranked Page in the GSEA User Guide.
如果是affymetrix的表达矩阵，不需要提前进行Present/Marginal/Absent Calls. 来过滤掉一些表达探针，GSEA需要各种情况的表达数据。
如果是gct and pcl 的表达矩阵，缺失值空着就好了。但是如果缺失值太多了，这样在计算signal-to-noise的时候，不同group的样本数就不一致了，mean和sd都会变好，最好是避免这样的情况，可以考虑进行插值，或者过滤掉这样的探针。

同时不需要提前过滤掉低表达量的探针或者低variance的探针。它们都会在我们算好的 ranked gene list 的中间部分，增强我们的统计效应。完全不用担心数据量计算时间的问题。

值得注意的是如果要想计算Signal2Noise ，每个group必须要有3个及以上的samples

值得一提的是除了两个group之间的比较可以做gsea之外，还可以针对连续性的phenotypes和time-course数据。参考：http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html

假设芯片或者其它测量方法测到了2万个基因，那么这两万个基因在case和control组的差异度量(六种差异度量，默认是signal 2 noise，GSEA官网有提供公式，也可以选择大家熟悉的foldchange)肯定不一样,那么根据它们的差异度量，就可以对它们进行排序，并且Z-score标准化，在下图的最底端展示的就是

img

那么图中间，就是我们每个gene set里面的基因在所有的2万个排序好基因的位置，如果gene set里面的基因集中在2万个基因的前面部分，就是在case里面富集，如果集中在后面部分，就是在control里面富集着。

而最上面的那个ES score的算法，大概如下：

仔细看，其实还是能看明白的，每个基因在每个gene set里面的ES score取决于这个基因是否属于该gene set，还有就是它的差异度量，上图的差异度量就是FC（foldchange）,对每个gene set来说，所有的基因的ES score都要一个个加起来，叫做running ES score，在加的过程中，什么时候ES score达到了最大值，就是这个gene set最终的ES score！

所谓的GSEA分析，就是一个个遍历探索已知的基因集合，在我们的表达矩阵里面是否出现了某种统计学显著的扰动，如上图所示，要深入理解，请看我在生信菜鸟团写的另外3个教程：