富集分析
基因富集分析(gene set enrichment analysis)是在一组基因或蛋白中找到一类过表达的基因或蛋白。研究方法可分为三种:Over-Repressentation Analysis(ORA),Functional Class Scoring(FCS)和Pathway Topology。ORA是目前应用最多的方法,GO富集分析和KEGG富集分析就是使用的这种方法;FCS这种方法应用于GSEA分析。
功能分析(functional analysis)/ 通路分析(pathway analysis)是将一堆基因按照基因的功能/通路来进行分类。换句话说,就是把一个基因列表中,具有相似功能的基因放到一起,并和生物学表型关联起来。GO分析是将基因分门别类放入一个个功能类群,而pathway则是将基因一个个具体放到代谢网络中的指定位置。
为了解决将基因按照功能进行分类的问题,科学家们开发了很多基因功能注释数据库。这其中比较有名的就是Gene Ontology(基因本体论,GO)和Kyoto Encyclopedia of Genes and Genomes(京都基因与基因组百科全书,KEGG)。
GO
GO数据库是基因本体论联合会(Gene Ontology Consortium)建立的一个数据库(官网http://geneontology.org/),旨在建立一个适用于各种物种的、对基因和蛋白功能进行限定和描述的,并能随着研究不断深入而更新的语义词汇标准。分别从分子功能、参与的生物途径及细胞中的定位对基因产物进行了标准化描述,一个基因对应有一个或多个GO term(GO功能),一个term对应多个gene。
GO注释分为三大类,分别是:分子生物学功能(Molecular Function,MF)、生物学过程(Biological Process,BP)和细胞学组分(Cellular Components,CC),通过这三个功能大类,对一个基因的功能进行多方面的限定和描述。
Cellular component解释的是基因产物在哪里,在细胞质还是在细胞核,如果存在细胞质那在哪个细胞器上,如果是在线粒体中那是存在线粒体膜上还是在线粒体的基质中。
Biological process该基因参与了哪些生物学过程,比如参与了rRNA的加工或参与了DNA的复制。
Molecular function该基因在分子层面的功能是什么,它是催化什么反应的。
通常在得到差异表达基因后,可通过功能富集分析进一步筛选对生物体具有重要意义的基因。将筛选得到的基因分门别类放入细胞组分CC、分子功能MF和生物过程BP三个功能类别中,基因产物被尽可能的富集到最低层的功能term上。寻找各个基因是否有共同的GO条目,或者有没有共同的上级GO条目,可以发现具有某些共同特点的基因。根据超几何分布关系,GO分析会对涉及的GO返回一个p-value,小的p值表示差异基因在该GO 中出现了富集。GO 分析对实验结果有提示的作用,通过差异基因的GO 分析,可以找到富集差异基因的GO分类条目,寻找不同样品的差异基因可能和哪些基因功能的改变有关。
补充:
GO是Gene Ontology的缩写。本体论是哲学概念,它是研究存在的本质的哲学问题。后来这个词被应用到计算机界,定义为概念化的详细说明。在实现上,一个ontology往往就是一个正式的词汇表,其核心作用就在于定义某一领域或领域内的专业词汇以及他们之间的关系,是领域内部不同主体之间进行交流的一种语义基础。
使用GO的时候一般需要GO定义文件和GO关联文件。GO定义文件存放GO词条的定义,而GO关联文件则是不同命名体系与GO词条的映射关系。条目标准定义:
id:GO编号,如:GO:0031985
name:全称,Golgi cisterna
ontology:命名空间namespace,cellular_component
definition:定义,Any of the thin, flattened membrane-bounded compartments that form the central portion of the Golgi complex. Source: GOC:mah
条目之间的关系,采用有向无环图(Directed Acyclic Graphs,DAG)的形式。注释系统中每一个节点就代表了一个基本描述单元(term),有向指的是term之间的单向指向性关系,比如termA是内质网,termB是细胞器,规定A是B,却不能说B是A;无环指的是从任何一点开始沿着规定的指向都不能回到原点。
KEGG
KEGG是一个整合了基因组、化学和系统功能信息的综合数据库。KEGG下属4个大类和17和子数据库,而其中有一个数据库叫做KEGG Pathway,专门存储不同物种中基因通路的信息,也是用的最多的一个,所以,久而久之,KEGG就被大家当做是一个通路数据库了。
GO分析好比是将基因分门别类放入一个个功能类群,而pathway则是将基因一个个具体放到代谢网络中的指定位置。根据挑选出的差异基因,计算这些差异基因同Pathway 的超几何分布关系,Pathway 分析会对每个有差异基因存在的pathway 返回一个p-value,小的p 值表示差异基因在该pathway 中出现了富集。pathway 分析对实验结果有提示的作用,通过差异基因的Pathway 分析,可以找到富集差异基因的Pathway 条目,寻找不同样品的差异基因可能和哪些细胞通路的改变有关。pathway 是蛋白质之间的相互作用,pathway 的变化可以由参与这条pathway 途径的蛋白的表达量或者蛋白的活性改变而引起,因此pathway 分析的结果更显得间接。
Over-Repressentation Analysis(ORA)
过表征分析,其实就是想看看某类功能或分类和随机事件相比是否有更明显的趋势。统计方法包括Fisher精确检验、卡方检验等。Fisher精确检验是基于超几何分布计算的,它分为两种,分别是单边检验(等同于超几何检验)和双边检验。超几何分布检验常用来对venn图两个圈overlap的显著性进行检验,Fisher精确检验常用来对2×2的列联表进行检验。
分析列联表中两个变量的关联,可以采取卡方检验(Chi-square test)。先假设两个变量之间没有关系(是否在这个GO term和是否在目标基因集没有关系,即目标基因集在特定GO term没有富集),计算统计量:Σ(实际值-理论值)^2/理论值,然后根据自由度(等于1)查表得到p值。如果p值小于0.05,说明原假设不成立,即目标基因集在特定GO term出现了富集。
卡方检验最大的优势在于计算比较简便,可以徒手计算,Fisher精确检验的计算要相对复杂些,但是现在实现起来也很容易了。对于2×2列联表来说,卡方检验通常只能做为近似估计值,特别是当总样本量或理论频数比较小的时候,计算并不准确。一般情况下,如果总样本量大于40,最小理论频数大于5,可以使用卡方检验。但是,如果采用卡方检验得到的P值在0.05附近时,应该用Fisher确切概率法。如果差异很大,采用卡方检验和Fisher确切概率法得到的结果相差不大。现在GO富集分析一般都是使用超几何分布进行计算的。
富集分析的超几何分布检验的p值计算如下。N为所有基因中具有pathway/GO term注释的基因数目;n为N中差异表达基因的数目;M为所有基因中注释为某特定pathway/GO term的基因数目;m为注释为某特定pathway/GO term的差异表达基因数目。通过计算得到的P value会进一步经过多重检验校正,通常应用的是BH方法,得到FDR值。然后以FDR≤0.05为阈值,满足此条件的pathway/GO term定义为在差异表达基因中显著富集的pathway/GO term。此外还有很多其他的算法来试图解决一个基因对应多个GO term、一个term对应多个gene的问题,但是本质上也是基于Fisher's exact test。
常见的富集结果描述包括:
RichFactor,富集因子,是指感兴趣基因列表中属于这个term的基因的数量/背景基因集中富集在这个term中所有基因的数量。
p值或q值:代表富集显著程度,可以映射到图形颜色。
GeneNumber:感兴趣基因列表中属于这个term的基因数量。
Gene Percent(%):感兴趣基因列表属于这个term的基因的数量占感兴趣基因列表所有基因数量的百分比
富集分析工具
推荐clusterProfiler,它支持ORA和FCS两类算法。函数为:enrichGO, gseGO: GO富集分析;enrichKEGG, gseKEGG: KEGG富集分析;enrichDAVID: DAVID富集分析。
DAVID(https://david.ncifcrf.gov/)是由美国Leidos 生物医学研究公司的LHRI团队开发的一个在线基因注释及功能富集网站,最为常用且权威,引用超高(>21000)。但是它的数据库版本比较老,目前最新版的DAVID 6.8还是在2016年更新的,而且基本只更新了GO和ID转换的数据,KEGG也没有更新。听说2016年Nature Methods 专门写了Impact of outdated gene annotations on pathway enrichment analysis 吐槽大家还在用老旧的DAVID。
GSEA
Gene Set Enrichment Analysis(基因集富集分析)用来评估一个预先定义的基因集S(已知功能的基因集)的基因在与表型相关度排序的基因列表L(按照logFC、Signal to Noise Ratio等排序的基因列表)中的分布趋势(是随机分布,还是主要分布在顶部或底部),从而判断其对表型的贡献。GSEA确定一个预先定义的基因集是否能在两个生物学状态中显示出显著的一致性的差异,通俗一点就是某个通路/GO条目中的基因集在实验组和对照组中呈现出一 致的上调或者下调趋势。
富集分数enrichment score(ES)代表集合S在排序列表L的顶部或底部被过表达的程度。这个分数是通过遍历列表L来计算的,当我们遇到一个在S中的基因时增加一个running-sum statistic(类Kolmogorov-Smirnovlike统计量),当遇到的基因不在S中时减少统计量。增量的大小取决于基因统计(例如基因与表型的相关性)。ES为random walk中遇到的与零的最大偏差(maximum deviation from zero)。GSEA的那条曲曲折折的线就是通过不断的加分减分做出来的,图中的每一条垂直线表示基因集S中一个基因。
利用置换检验(permutation test)计算ES的p值。具体地说,我们对基因列表L的gene labels进行重新排列(permute),并为排列后的数据重新计算基因集的ES(重复1000次),从而为ES生成一个null distribution。然后相对于这个零分布计算观察到的ES的p值。并使用FDR调整计算q值。
各种方法的特点
ORA方法存在一些问题:仅使用了基因数目信息,而没有利用基因表达水平或表达差异值,为了获得感兴趣或者差异表达基因,需要人为的设置阈值;ORA法通常仅使用最显著的基因,而忽略差异不显著的基因。在获得感兴趣的基因时, 往往需要选取合适的阈值, 有可能会丢失显著性较低但比较关键的基因, 导致检测灵敏性的降低;假设每个基因都是独立的,忽视了基因在通路内部生物学意义的不同(如调控和被调控基因的不同)及基因间复杂的相互作用;ORA假设通路与通路间是独立的,但这个前提假设是错误的。
FCS方法相较于ORA 法在理论上有明显突破,考虑到了基因表达值的属性信息, 以待测基因功能集为对象来进行检验, 也使得检验结果更加灵敏。认为虽然个体基因表达改变之后会更多在通路中体现,但是一些功能相关基因中较弱但协调的变化(small but consistent changes)也有明显的影响。但是仍独立分析每一条通路,但同一个基因可能涉及多条通路,所以不同通路间的基因出现重叠,所以别的通路可能由于重叠的基因,也出现显著富集;仍然把待测基因功能集中的每个基因作为独立的个体, 忽略了基因的生物学属性和基因间的复杂相互作用关系。