富集分析中的超几何分布

2022-06-18

当我们从测序数据中得到了一系列的基因后,我们需要研究这些基因都属于那些通路,那些通路是跟我们研究课题相关的。为了解决这一问题,最常用的做法就是做通路注释然后做通路富集分析,看看我们得到的基因都分布在哪些通路。 超几何分布是富集分析的常用方法,常用的GO富集分析都是用超几何分布计算的。下面将浅显的探讨一下超几何分布的原理。

1、超几何分布

超几何分布是一种非常常见的分布,常用来表示在N个物品中有指定商品M个,不放回抽取n个,抽中指定商品的个数,即X~H(N,n,M),则抽中k件M商品的概率为:

img

在这里我们做一个简单的概念转换即可知道软件是如何做GO富集分析的:

  1. N为GO注释数据库中的总基因数;
  2. M为数据库中属于某个GO子类的基因数;
  3. n为我们得到的需要进行GO富集分析的基因的总数目4;
  4. k为n中属于M的数目。

因此我们就可以计算基因集n是否在M类中富集的概率。 但是知道这个概率后并不能直接用来作为富集分析的结果,必须要对其进行一个评估,因为我们必须要考虑到随机情况,如果随机从N中抽取n个基因,其中k个在M中的概率很高的话,那我们富集得到的通路意义就是极小的。这时候我们引入p值对富集分析的概率结果进行分析。

2、p-value检验

P值就是当原假设为真时所得到的样本观察结果或更极端结果出现的概率。如果P值很小,说明这种情况的发生的概率很小,而如果出现了,根据小概率原理,我们就有理由拒绝原假设,P值越小,我们拒绝原假设的理由越充分。通俗的讲,p值就是指随机出现的概率,p值越小说明越不可能随机出现,也就是说我们得到的结果越具有显著性。 总之,P值越小,表明结果越显著。但是检验的结果究竟是“显著的”、“中度显著的”还是“高度显著的”需要我们自己根据P值的大小和实际问题来解决。 在我们的富集分析中,p值是由下面这个式子计算得到的: img

上面式子的意思是: 从总N个基因抽n个基因, 作为分母,分子是M个基因有i个落在通路里,有n-i个不落在通路里。 p-value是指你观察到m个基因落在通路里,比这还要更极端的概率之和,所以i是从m到M。 就是说看到更多的基因落在这个通路里的所有可能。所以超几何检验很方便地 可以给你算一个p-value,最后得到p-value<0.01或者0.05,你的结果如果定义p-value<0.05 那就有5%的概率看到是一个假阳性,这里我们只是在谈拿一个通路来做检测, KEGG现在大概有360多个通路,每一个通路都做一个超几何检验,每一次有5%的概率出错,一共进行360次, 那出错的概率就很会大很多, 所以怎么评估最终看到的结果是真的而不是被误导的呢? 你就要算一个叫多假设检验的矫正,只要做了多次的statistical test, 就要做多假设矫正。矫正有多种方法,现在大家用的最多的是FDR校正。

3、FDR校正

FDR矫正的是false discovery rate, 也就是FP/(TP+FP)的期望值,看这个期望值是多少。 如果这个期望值小于0.05,大家就认为有可能是 一个真实的有生物学意义的结果。

img

参考

补充1. 组合数计算公式

从 n 个不同元素中每次取出 m 个不同元素 (0≤m≤n),不管其顺序合成一组,称为从 n 个元素中不重复地选取 m 个元素的一个组合。所有这样的组合的种数称为组合数。 组合数的计算公式为

组合数计算公式

n 元集合 A 中不重复地抽取 m 个元素作成的一个组合实质上是 A 的一个 m 元子集和。如果给集 A 编序 成为一个序集,那么 A 中抽取 m 个元素的一个组合对应于数段 到序集 A 的一个确定的严格保序映射。

在线性写法中被写作C(n,m),如下:

组合数符号线性写法

补充2. 超几何分布

(不放回抽样) 设一批产品共有N个,其中有M个次品。每次从这批产品中随机地抽出一件来检查,检查后不放回,共取n次(相当于一次同时取n件产品),试求在n次检查中有k次是次品的概率 img 。从N件产品中抽取n件共有 img

种不同的取法,现要求在抽取的一组n件产品中,有k件次品和n-k件合格品。因为这k件次品有 img 种不同的取法,n-k件合格品有 img 种不同的取法,故得

上式确定的分布律称做超几何分布,这是因为 img 可看成是所谓超几何级数的通项的系数。

超几何分布是统计学上一种离散概率分布。它描述了由有限个物件中抽出n个物件,成功抽出指定种类的物件的个数(不归还)。 超几何分布和Fisher’s Exact Test是完全一模一样的原理,只是两种不同的称谓。

补充3. R包中实现超几何分布

R包中实现 R中自带超几何分布的检验(stats包)

4.1 方法1 phyper(q-1, m, n, k, lower.tail=F)

4.2 方法2 1 - phyper(q-1, m, n, k)

Note:两种方法的参数如下: q = the number of white balls drawn from the urn (without replacement) q对应到抽样问题,为k

m = the number of white balls in the urn m对应到抽样问题,为M

n = the number of black balls in the urn n对应到抽样问题,为N-M

k = the number of balls drawn from the urn (sample size) k对应到抽样问题,为n

Note: 加上lower.tail=F参数时,算的是X > q的概率 (不包括等于的情况),此时要用q-1代替q计算,以便概率包括X=q的情形;不加lower.tail=F参数时,算的是X <= q的概率(包括等于),此时也要用q-1来替代q,以便概率包括X=q的情形。

这里有一个对基因集用R语言做超几何分布检验的代码例子,解释得很清楚,请点击这里here