点点点 | 真香！Simple GO GSEA 富集分析 ~ / 开普饭

写在前面

时间拨回去 2015年，那时我接触生信已有一年，TBtools 开发尚在萌芽阶段。那会，我写了几款小的软件，包括 “blast3go”，为的是应对即将收费的“blast2go”。当然，后来相关功能都整合到 TBtools 中。而其中有一个重点功能，即 GO 富集分析。那会在 Bioinformatics中国社群，我们开始了理论上是国内最早的公开社群学术Seminar（网络直播），我在其上也分享了相关学习经验和知识。主要介绍了最常用的基于超几何分布的富集分析方法，事实上，可以理解为简单的 Fisher's Extact Test。同时，也介绍了 GSEA，全称 Gene Set Enrichment Analysis。这个有趣。因为对于新手来说，概念似乎有点混乱。

常见的富集分析逻辑

目前，在科研工作中最为常见的三种富集分析策略，分别是：

ORA，全称，Over Represent Analysis，即前述基于超几何分布的统计检验。TBtools当前使用的也是，这个方法，最大的优点是，非常快！可以说，适用于几乎所有情况，毕竟生物学问题，我们常常希望的是抓主要矛盾，而我们能做的，也是抓主要矛盾。
卡方检验，这个比较简单，针对任何一个GO条目或者KEGG通路，更或者其他注释，很容易拿到背景属于和不属于该条目的两个数字，选择集属于和不属于该条目的两个数字，于是就是列联表检验。这个方法其实用得不过，尽管也可以看到。
GSEA，翻译，基因集合富集分析。Emmm，这个有点意思。一般来说，我们在现实生活中总是说 GO富集分析，KEGG富集分析，统称基因集合富集分析。而GSEA这个方法本身翻译过来，也是基因集合富集分析。但其中指代区别，还是要知悉。

既生ORA，何生GSEA？

事实上，关于富集分析方法的选择和优化，是一个持久不衰的热门话题，包括近期更新的国人佳作 KOBAS-i 也提出了新的更完善的方法，感兴趣的朋友，建议了解。那为啥如此？原因不用说，肯定是旧方法还是有一些值得改进的地方嘛。其中ORA最明显的问题，就是“硬阈值”。ORA要求的输入是，背景注释+选择集。而选择集要么就是差异表达基因，要么可能就是GWAS关联分析出来的基因。具体到差异表达基因上，往往要求 pvalue必须低于 0.05？或者0.01？同时差异表达倍数要达到 2倍以上。这是比较有趣的，因为我们都知道 1.99999.... 倍的差异其实也算差异。这是 GSEA 方法解决的主要问题之一。
对于 GSEA 方法，一般我们可以提供背景注释+基因排序信息表。后者或者物种所有基因的表达矩阵，由软件计算差异表达倍数或其他来排序以及作为权重；或者是直接输入物种所有基因在某个两两比较的差异表达倍数，Log2FoldChange值；其实还可以是其他信号强度数值。于是不存在“硬阈值”。解决了ORA的问题。同时其实还解决了另外一个问题，这里埋伏笔，下述给出。

GSEA~真香定律

事实上，15年我用过两三次，做完一个医学项目之后，就再也没碰过 GSEA，因为我几乎在所有项目中发现，ORA的表现已经足够好，根本不需要GSEA。尤其是 GSEA 计算起来，真的非常非常非常慢（当然，非常感兴趣的可以了解一下大湿兄的R语言包 ClusterProfiler，据说很快）。由于我确实不懂 R语言，所以一直也没去折腾，索性就一直没用。
但最近做到一组数据，做的是果皮软化。很明显，果皮肯定是软化的，表型非常明显，细胞壁代谢相关基因必然也差异表达，但是在常规富集分析中，没有任何相关显著富集结果。

这个不应该出现。我回头看了下差异表达基因列表，发现其实倍数差异排前的，几乎都与细胞壁代谢相关。直觉告诉我，果实软化过程，基因表达变化过大（对应差异表达基因数目偏多），导致没有合理的富集结果。于是，我想起来 GSEA。前述提及，他解决“硬阈值”的问题，但这个问题在本例其实没啥用，因为都差异表达的。而他其实同时还解决了一个问题，那么就是权重！只要我们使用每个基因的差异表达倍数作为权重，那么差异倍数大的，自然权重就高，可以被正确反映到富集分析结果上。