查看原文
其他

GSEA分析一文就够(单机版+R语言版)

jimmy 生信技能树 2022-06-06

通过前面的讲解,我们顺利的了解了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的算法,大概如下:


1

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

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

  • java版本GSEA软件的ES score图片的修改 http://www.bio-info-trainee.com/2105.html

  • GSEA的统计学原理试讲 http://www.bio-info-trainee.com/2102.html

  • 制作自己的gene set文件给gsea软件 http://www.bio-info-trainee.com/2144.html


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存