富集分析Enrich me again!

一般提到富集分析，首先想到的就是GO、KEGG这两把刷子，然后还需要知道两个重要概念：前景基因、背景基因

前景基因：你关注的要重点研究的基因集；
背景基因：所有的基因集
比如做转录组测序，一般都要设置处理组和对照组，前景基因是处理和对照的差异基因，背景基因就是两组样本的全部表达基因。

GO富集分析

Gene Ontology: 描述基因的层级关系【基于ORA算法】可以算得上是高通量数据分析的标配，转录组、甲基化、ChIP-seq、重测序等，都会用到对一个或多个集合的基因进行功能富集分析，来找这个基因集的功能偏好性

刚开始看到这两个词汇，我就看第一个单词比较熟，至于第二个嘛，“本体论”什么鬼？

本体论，既然是本体，那么就要求全，并且权威。因此，它就是某个科学领域中的专业词汇的集合，并且经过多方考证，能被公认的词汇表
GO发展
生物数据库有很多，但是不同的数据库之间可能会使用不同的术语，大量的研究时间都被浪费在查找信息上。于是基因本体联合会组织（Gene Ontology Consortium）建立了GO数据库，规范统一了对于不同物种的基因和蛋白描述。整个数据库前身就是1988年三大数据库的整合（果蝇数据库、酵母数据库、小鼠基因组数据库），现如今，已经支持几十种物种
发展30年，条目一定很多，有没有什么命名规律？

条目标准定义：

id：也就是GO编号，如：GO：0031985
name：全称，Golgi cisterna
ontology：命名分类，cellular_component
definition：定义，Any of the thin, flattened membrane-bounded compartments that form the central portion of the Golgi complex. Source: GOC:mah

除了这些标针对每个条目term的标准定义，还要定义条目之间的关系，让他们形成一个有结构的框架，而不是一盘散沙。于是GO采用了“有向无环图”的模式（有向指的是term之间的单向指向性关系，比如termA是内质网，termB是细胞器，规定A是B，却不能说B是A；无环指的是从任何一点开始沿着规定的指向都不能回到原点）表示这种关系的方式有三种

is_a/part_of/regulates
【详见官网http://www.geneontology.org/page/ontology-relations】is_a和part_of都具有传递性，比如A is a term of B / A is part of B；regulate 表示一个过程受另一个过程的影响，前者调控后者
能简单点吗？主要包括啥？

它将基因分门别类放入细胞组分CC、分子功能MF和生物过程BP三个功能类别中（分别对应基因产物在哪里发挥功能、怎样发挥功能、发挥什么样的功能）

基因产物被尽可能的富集到最低层的功能term上。如果一个基因产物被定为在Golgi cisterna（高尔基池），那么它只会定位到这个特定的条目，而不会定位到它的父节点或其他上级条目中
GO的具体作用？

帮助我们找到提交的基因集中各个基因是否有共同的GO条目，或者有没有共同的上级GO条目，可以发现具有某些共同特点的基因（比如说是某类细胞器的组成成分，或在某个催化反应起作用，又或者在某条通路中起作用）

GO气泡图可以表示显著性不同的GO条目对应的基因数量；

GO条形图可以表示三个分类的前多少条目对应基因的数量，可以每个分类单独做，也可以组合在一起；

GO网络图可以表示显著性较强的GO条目之间的相互关系（点大小代表条目内基因数目多少，颜色代表P-value，值越小越红）；

GOmap图可以表示显著性较强的GO条目的层级关系，以树状图的形式展现
GO数据库是不断更新的，因此GO富集结果和数据库使用关系密切，有时候公司给你返回来结果，你却找不到想要的基因，除了实验因素，有可能是注释环节采用的工具没有及时更新，采用了老版本的软件和数据库。比如Nature methods就在2016年发表了一篇Impact of outdated gene annotations on pathway enrichment analysis，吐槽DAVID在线工具的数据库老旧，更新慢的问题。再比如Y叔曾经吐槽过华大的WEGO数据库还是09年的。虽然做出来图是那么回事，但是解释不了当前的问题。其实数据库并没有错，GO注释每天都在更新，Pathway数据库例如Reactome和PathwayCommons每个季度都在更新。
因此，建议不要太依赖在线工具，因为它的维护不一定那么及时

KEGG

Kyoto Encyclopedia of Genes and Genomes: 系统分析基因产物和化合物在细胞中的代谢途径以及这些基因产物的功能的数据库【基于ORA算法】

它的产生就是因为大量的数据被获得，却埋没在浩如烟海的文章之中，没有人去归档整理，但是使用这些数据的人还不少，他们每次都在检索上花费了大量时间。而日本的商人看到了商机，他们的收纳功力世界文明，于是把很多数据分类整理，做成了一个付费的平台，这开始就是一个“卖方市场”，十分抢手。但是到如今，其实大量的开源软件对它构成了威胁
包括什么？整合了基因组、化学分子和生化系统等方面的数据,包括代谢通路（KEGG PATHWAY）、药物（KEGG DRUG）、疾病（KEGG DISEASE）、功能模型（KEGG MODULE）、基因序列（KEGG GENES）及基因组（KEGG GENOME）等等。KEGG 有一套完整KO注释的系统，可完成新测序物种的基因组或转录组的功能注释（KO是蛋白质或酶的一个分类体系，将同一条通路上功能相似、序列相似的蛋白质归为一类）。因此它可以将基因一个个归置到代谢网络指定位置上。
命名规则：
K（大写） +num基因ID号，表示所有同源物种中具有相似结构或功能的一类同源蛋白，如：K04456表示丝氨酸蛋白激酶；
ko+num代谢通路，表示特定的生物路径，如：ko04151表示PI3K-Akt信号通路【也是我们常用的代谢通路】；

M+num表示模块，如：M00676表示PI3K-Akt信号模块
C+num表示化合物，如：C00533表示一氧化碳
ECx.x.x.x表示酶，如：EC2.7.11.1表示丝氨酸
R+num表示反应名称

关于物种缩写：https://www.genome.jp/kegg/catalog/org_list.html

举个例子：三个字母表示物种，hsa表示Homo sapiens；具体的KO号，如K12407表示和葡糖激酶glucokinase序列和功能相近的蛋白质/酶类，当然一个KO号有可能会对应好几个数字（基因登录号），表示细胞中存在几种不同的葡糖激酶，分别由以上几种数字表示的基因编码
怎么使用该图？

image

这个图可以用pathview函数获得。先看看有没有关注的基因能注释到通路上，主要看颜色：红色表示上调的差异基因，绿色表示下调，白色的表示没有差异基因成功注释。
如果自己期望的一些基因没有成功注释，就代表没有吗？其实也不是，可以再试试KEGG官网数据库，使用BlastKOALA，上传蛋白序列-选择物种-填写邮箱
当然，KEGG还可以做条形图、气泡图等，但得到数据是最重要的，详细内容可以参考clusterProfiler

GSEA

Gene Set Enrichment Analysis 基因集富集分析，用于评估一个基因集的基因在表型相关度排序中的分布趋势，进而判断它们对表型的贡献

特点：

第一，它与GO和KEGG在算法层面上就是不同的，它采用的是FCS算法；

第二，与GO不同的地方在于，GO是先筛选基因（需要人为设定阈值），再判断差异基因在哪些通路有富集；GSEA可以考虑那些表达差异不大却功能重要的基因对通路影响，相比GO和KEGG能保留更多信息；

比如，你现在要研究某种昆虫的抗性，用传统的富集分析发现10个基因在细胞凋亡通路里有显著差异，同时还有10个基因在细胞扩增通路中有显著差异，那么到底是由于阻遏抗性的细胞凋亡引起抗性还是由于本身突变的细胞扩增产生了抗性？
用GSEA分析就可能会发现在细胞凋亡通路里，除了10个差异显著的基因表达量上升外，还有一些其他的相关基因也有整体上升的趋势；而扩增通路就不存在，因此可以进一步说明抗性和细胞凋亡的相关性更强

第三，研究表型的差异就要找到相应的基因，但是“孤木不成林”，即便是有一两个差异较大的基因，也不会直接导致表型变化。相对应的还要有一大帮相关的上游或者下游基因的变化，但它们的强度可能没有主效的那几个基因大。用GSEA就可以全面地进行分析，分析基因集中整体表达量得到的结果更可靠；

第四，你可能会想，为什么这个东西这么好，但是日常研究见不到呢？因为GSEA的数据库还并不完善，目前主要支持人类相关的研究，动植物还欠缺http://software.broadinstitute.org/gsea/msigdb/collections.jsp#C3
原理： 使用一个排序的基因列表genelist（使用gseKEGG函数需要降序排列）、一个基因集geneset（比如可以自己选择某个代谢通路的基因、同一个GO注释的基因），判断geneset中的基因在genelist中是随机分布还是聚集在genelist的顶部或底部【落在顶部或底部表示geneset中的基因对表型差异用贡献】

image

ES是Enrichment score 富集得分：表示基因集的基因在基因列表两端的富集程度。开始从左到右基因集的第一个基因开始，每次计算统计值并且逐次累加，有一个在基因列表中的基因，就增加计算的统计值，反之就减小。并且每次的增加和减小都被记录在上图红线的位置，它代表了基因和表型的相关性。结果得到的峰值计作富集得分ES，正值表示基因集的基因主要在基因列表顶部富集，负值表示在底部富集
绿色圆圈表示的Leading edge subset：表示对ES富集得分最大的基因集中的基因，ES为正值表示是峰左侧的基因，负值表示右侧的基因。这个图中的每一条线表示基因集中一个基因