Ⅰ GO富集分析
GO富集是组学数据分析常用的手段,通常用来挖掘差异基因中GO term的富集程度。Fisher's exact test是常用的统计检验方法,但这种方法存在明显的缺点。很多公司提供的测序分析结果都普遍使用这样的方法,导致很多后续的分析与实验结果不一致的情况。对于这种情况,目前还有其他算法来弥补这些缺点。(文中例子来源于《the Gene Ontology handbook》)
fisher's exact test是基于超几何分布来计算的,单边检验就是超几何检验。通常用来检验两组分类是否有显著差异。
m:研究物种的基因数;
n:研究的样本中基因数;
mt:总体中被注释到term t(GO 词条t)的基因数;
nt:样本中被注释到term t的基因数。
随机变量Xt表示样本中被观察到的term t 的数目,所有根据超几何分别,观察到k个term t 的概率P(Xt)是:
零假设H0:样本中出现的term t 的数目与总体中总的term t数目没有正关联。也就是说样本中的term t数目的比例与总体中term t的接近。
为了拒绝H0,使用单尾检验:
一个简单的例子:假设总体中有18个基因,其中有5个被注释到binding这个term,转录组分析发现有5个差异表达的基因,其中有3个被注释到binding这个term,为了说明binding这个term是否是overrepresentation,用上面的Fisher's exact test计算p值:
在现实中,我们不可能只对某一个term进行检验,而是对很多term进行检测,即多重检验,但这样就会导致假阳性的term数目非常高。所以我们需要对p值进行校正。通常是使用Benjamini-Hochberg校正方法来控制预期的错误发现率(false discoveries rate-FDR)进行校正。( 如何通俗地解释错误发现率(FDR) )。尽管多重检验的校正可以减少假阳性,但并不能从根本上解决GO(或KEGG)富集的问题。
GO富集的根本问题在于一个基因对应的GO term有多个,一个term对应多个gene,同时还有层级关系。这样导致如果一个term显著富集,那和它共享很多基因的term也会显著富集。
有很多其他的算法来试图解决这个问题,其中包括parent-child approach、topology-based algorithms、model-based approaches和gene set enrichment analysis。下面是对这些算法的简单介绍:
该算法还是基于Fisher's exact test,只不过考虑了term的父节点。在计算概率时,不再是在总体m中取样,而是从term的父节点中取样,所以计算公式变成了:
当一个term有多个父节点时计算就变得复杂了,具体方法还得参考原始文献(improved detection of overrepresentation of Gene-Ontology annotatins with parent child analysis)
这两种方法原理反正我没看懂,有兴趣的可以看原始文献:
1、Improved scoring of functional groups from gene expression data by decorrelating GO graph structure.
2、GOing Bayesian: model-based gene set analysis of genome-scale data
该算法首先根据感兴趣的特征(比如差异基因的表达量)对基因进行排序,形成一个列表。零假设是某个基因集(genes encoding procts in a metabolic pathway, located in the same cytogenetic band, or sharing the sam GO category)里基因顺序与这个列表没有关联,即排序是随机的。对应的备择假设是它们之间有关联。如果基因集里的基因都聚集在基因列表的前端或底端或者非随机分布,我们就倾向于相信它们之间有关联。
S:想研究的基因集;
L:整个排序的基因列表;
统计量:Kolmogorov-Smirnov(KS)
step 1:计算富集得分(Enrichment Score)。按顺序从头到尾逐个比较L中的基因与S中的基因,加和统计量,如果两者相同就增加KS统计量,反之就减少KS统计量。增加的多少与这个基因和表型的相关性有关。最后ES就是KS的最大方差值。
step 2、检验ES的显著性。重复k次随机选择的大小为nt的基因集(Nt1,...,Ntk),p值计算公式为:
step3、使用FDR进行多重检验的校正。
相关软件:GSEA-P software
1、Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practercal and powerful approch to multiple testing.
2、 https://www.hu.com/question/3560619
3、Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles
4、 《the Gene Ontology handbook》
Ⅱ 「GO富集分析」从原理到实践 ~ 零基础掌握
原本,我并无写这一稿件的想法。主要原因有二:
如果要找合理解释,那么针对第一点,就是每天仍然有大量新接触生信数据分析的朋友;针对第二点,......在前两天我推的文稿《零基础快速完成基因功能注释 / GO / KEGG / PFAM...》中,评论区答应了下,阅读过5000,那就写一写富集分析。于是,如果不写,总是不对。如果要写,只能现在写。毕竟有些事情,现在不做,以后真的不会做。
对于这一块,完全陌生的朋友,尤其是不少生物学背景朋友,有必要温习一下数理统计基础。这一稿件只做原理最简单的但使用最广泛其速度最快的Over-Represence Analysis模式的富集分析讲演。其他模式,不涉及。
回到主题,先举个经典的抽球例子:
小红小绿小蓝三个人自称有超能力,可以用手摸摸球就分辨出黑球白球,于是我们找来黑袋子,放100个球,其中20个白球80个黑球,让三人分别无放回地抽取。
小红随机抽出来10个球,其中2个白球8个黑球,情况即,
抽球中白球比例与背景白球比例完全一致,说明小红抽球结果随机。
球放回去,小绿来抽球,抽出来的10个球,其中3个白球7个黑球,情况即,
这是经典的抽球案例,抽取到的白球个数的概率分布为超几何分布。基于此,我们可以简单计算抽取到比小绿抽取到球个数(或更多即更极端)的概率如何,在 R语言中计算,即
而对于小蓝的情况,那么概率如何?
在 TBtools 中也可以计算,只是写法有点区别
可以看到,尽管这只是一次抽球,小绿抽球中白球比例(或更极端情况)出现的概率是31.88%+,还是挺高的,于是我们有较高的把握说,小绿嘛,只是走了狗屎运。相反,小蓝抽球中白球比例或更极端情况出现的概率几乎为 0 ,我们几乎没啥把握说,小蓝走狗屎运....换句话说,我们有理由相信,或许小蓝真有抽白球的超能力.....
说了这么多,那么跟基因集合富集分析有啥关系?....基因集合功能富集分析。那么我们就需要有一个基因集合(如差异表达基因集合或ChIP-seq的Peaks或GWAS定位的系列区间),还有一个功能标签(如 生长素信号转导相关 )。于是黑白球案例可以简单调整一下。假定现在这个物种一共有100个基因,其中20个基因与生长素信号转导相关,80个没有注释到与生长素信号转导相关(换句话说,约等于无关),我们做了对植株做了处理,和CK分别测定转录表达谱,通过差异表达分析,鉴定到10个差异表达基因,其中2个与生长素信号转导相关,而另外8个则没注释到生长素信号转导相关,简单画一下,即
好,剩下的两个就不替换了。整体上,ORA模式的富集分析,本身就是经典的抽球案例,感兴趣的自行替换就可以了。
基本原理,相信都搞清楚了。不过还是有两三点需要注意:
具体如何做物种所有基因的背景注释,请参考前述推文《零基础快速完成基因功能注释 / GO / KEGG / PFAM...》。
首先,打开 TBtools GO 富集分析界面
整体如上,一共三个文件:
具体示例如下
点击 Start ,随后等待即可。完成时会有弹窗提示。查看输出文件
(写到这里,突然觉得这些都没啥意思,不知为何....就不详细写了,大伙自己看看列名,猜猜吧)
很多时候,我们会选择,筛选第一列,只看 Biological Process。一般这些与我们的生物学认知会贴近一些。
基因集合功能富集分析,是一个常常被谈起的话题,甚至近期都有不少新方法或算法被提出。感兴趣的朋友可以去了解。这份教程,只与大伙说最简单,但也是使用最为广泛的一种富集分析模式。无论是不是 TBtools 用户,理论上来说,都可以轻松理解并掌握,从原理到实践。
写到一半,其实我已经不想写了。原因非常简单,这也是为什么在我之前,并没有一个人写出来 TBtools 类似的工具。不是写不了,而是不想写。有时候,随着能力增长和知识积累,往往不再愿意做一些简单的事情。或许这还涉及到年龄的增长,角色的转变,责任的变化....云云。
小时候,我以为写 TBtools 玩玩;
后来,我以为我会一直写下去;
现在,,,,,,
Ⅲ 什么是GO富集分析,常说的GO功能分析、功能分析、Pathway分析是什么意思
Gene Ontology可分为分子功能(Molecular Function),生物过程(biological process)和细胞组成(cellular component)三个部分。蛋白质或者基因可以通过ID对应或者序列注释的方法找到与之对应的GO号,而GO号可对于到Term,即功能类别或者细胞定位。
功能富集分析: 功能富集需要有一个参考数据集,通过该项分析可以找出在统计上显著富集的GO Term。该功能或者定位有可能与研究的目前有关。
GO功能分类是在某一功能层次上统计蛋白或者基因的数目或组成,往往是在GO的第二层次。此外也有研究都挑选一些Term,而后统计直接对应到该Term的基因或蛋白数。结果一般以柱状图或者饼图表示。
1.GO分析
根据挑选出的差异基因,计算这些差异基因同GO 分类中某(几)个特定的分支的超几何分布关系,GO 分析会对每个有差异基因存在的GO 返回一个p-value,小的p 值表示差异基因在该GO 中出现了富集。
GO 分析对实验结果有提示的作用,通过差异基因的GO 分析,可以找到富集差异基因的GO分类条目,寻找不同样品的差异基因可能和哪些基因功能的改变有关。
2.Pathway分析
根据挑选出的差异基因,计算这些差异基因同Pathway 的超几何分布关系,Pathway 分析会对每个有差异基因存在的pathway 返回一个p-value,小的p 值表示差异基因在该pathway 中出现了富集。
Pathway 分析对实验结果有提示的作用,通过差异基因的Pathway 分析,可以找到富集差异基因的Pathway 条目,寻找不同样品的差异基因可能和哪些细胞通路的改变有关。与GO 分析不同,pathway 分析的结果更显得间接,这是因为,pathway 是蛋白质之间的相互作用,pathway 的变化可以由参与这条pathway 途径的蛋白的表达量或者蛋白的活性改变而引起。而通过芯片结果得到的是编码这些蛋白质的mRNA 表达量的变化。从mRNA 到蛋白表达还要经过microRNA 调控,翻译调控,翻译后修饰(如糖基化,磷酸化),蛋白运输等一系列的调控过程,mRNA 表达量和蛋白表达量之间往往不具有线性关系,因此mRNA 的改变不一定意味着蛋白表达量的改变。同时也应注意到,在某些pathway 中,如EGF/EGFR 通路,细胞可以在维持蛋白量不变的情况下,通过蛋白磷酸化程度的改变(调节蛋白的活性)来调节这条通路。所以芯片数据pathway 分析的结果需要有后期蛋白质功能实验的支持,如Western blot/ELISA,IHC(免疫组化),over expression(过表达),RNAi(RNA 干扰),knockout(基因敲除),trans gene(转基因)等。
3.基因网络分析
目的:根据文献,数据库和已知的pathway 寻找基因编码的蛋白之间的相互关系(不超过1000 个基因)。
Ⅳ RNA-Seq(9):使用GSEA做GO/KEGG富集分析
最广为人知的富集分析做法是把上调、下调基因分别或者合并,拿来做GO和KEGG富集分析。经常有一些数据集,拿差异基因做得不到结果,那是因为确实富集不到任何通路,是正常的。不妨试试GSEA,不是拿差异基因,而是拿全部基因作为输入。
GSEA与GO,KEGG分析区别:GO,KEGG分析更加依赖差异基因,实则是对一部分基因的分析 (忽略差异不显著的基因),而GSEA是从全体基因的表达矩阵中找出具有协同差异 (concordant differences)的基因集,故能兼顾差异较小的基因
GO,KEGG富集是定性的分析,GSEA考虑到了表达或其它度量水平的值的影响。GSEA分析不需要指定阈值(p值或FDR)来筛选差异基因,在没有经验存在的情况下分析我们感兴趣的基因集,而这个基因集不一定是显著差异表达的基因。GSEA分析可以将那些GO/KEGG富集分信息中容易遗漏掉的差异表达不显著却有着重要生物学意义的基因包含在内。
另外,对于时间序列数据或样品有定量属性时,GSEA的优势会更明显,不需要每个分组分别进行富集,直接对整体进行处理。
数据准备,制作geneList
我们现在知道Cytokine-cytokine receptor interaction setSize enrichmentScore是被抑制的,如果还想看一下这个通路里面的基因是如何变化的,应该怎么办呢,pathview 可以帮到我们。