关于使用超几何分布做基因富集分析的一次讨论

biostar

转载请注明：来自生物统计家园：http://www.biostatistic.net

一次关于超几何分布的讨论，希望能和大家一起交流，共同讨论学习，可以持续跟帖~~·

A
22:14:27

10000基因，500个差异

其中一个PATHWAY有800个基因，那有50个基因来自这个PATHWAY的显著性，用超几何分布 怎么计算显著性

A  22:15:27
你知道超几何分布吗
B  22:15:50
你不就做富集分析吗
A  22:16:02
是的
B  22:16:50
背景10000
其中500个差异
一个通路相当于一次抽样，抽出800个
其中差异的50个
B  22:16:54
就这四个数
B  22:17:24
服从超几何分布，对应下公式计算就ok
A  22:17:24
不用减800
B  22:17:49
减掉800啥意思
A  22:18:35
就是把背景分成两部分了呀
抽出的800，和没有抽出的9200
A  22:19:13
这样才符合超几何分布呀
A  22:19:49
那你说
呢？
B  22:19:51
怎么可能，超几何分布，你再看看里了
B  22:19:55
理论
B  22:19:58
就一个背景
B  22:20:37
例如在有N个样本，其中m个是不及格的。超几何分布描述了在该N个样本中抽出n个，其中k个是无效的的机率：
B  22:21:00
在有N个样本，其中m个是不及格的，
注意，是其中，
B  22:21:14
认不是有N个好的，m个不好的
A  22:21:22
==
A  22:23:00
GO功能显著性富集分析给出与基因组背景相比，在差异表达基因中显著富集的GO功能条目，从而给出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向Gene Ontology数据库（http://www.geneontology.org/）的各个term映射，计算每个term的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的GO条目，其计算公式为
其中，N为所有Unigene中具有GO注释的基因数目；n为N中差异表达基因的数目；M为所有Unigene中注释为某特定GO term的基因数目；m为注释为某特定GO term的差异表达基因数目。计算得到的pvalue通过Bonferroni校正之后，以corrected-pvalue≤0.05为阈值，满足此条件的GO term定义为在差异表达基因中显著富集的GO term。通过GO功能显著性富集分析能确定差异表达基因行使的主要生物学功能。
A  22:23:14


A  22:23:24
你收不到图片
A  22:23:25
？
A  22:23:41
N为所有Unigene中具有GO注释的基因数目；n为N中差异表达基因的数目；M为所有Unigene中注释为某特定GO term的基因数目；m为注释为某特定GO term的差异表达基因数目
B  22:23:50
是啊
B  22:24:03
就一个背景N
A  22:24:42

那刚才的那个
就这样
1-phyper(50,800,10000,500)
A  22:24:47
我怎么感觉不对呀
B  22:24:51
咋
A  22:25:20
我在看看超几何分布
B  22:25:29
恩，再看看吧
A  22:25:57
超几何分布中，M,N 是独立的
B  22:26:13
本来就是独立的
A  22:26:16
这个 10000包括了800，所以要减去吧
B  22:26:24
独立又不是不相容
B  22:26:41
你把独立的概念和不相容搞混了
B  22:26:43
呵呵
B  22:26:52
独立是可以存在的
A  22:27:53
产品抽样检查中经常遇到一类实际问题，假定在N件产品中有M件不合格品，即不合格率p=M/N.在产品中随机抽n件做检查，发现X件是不合格品，可知X的概率函数为P(X=k)=C(k,M)*C(n-k,N-M)/C(n,N)
B  22:27:58
p(AB)=p(A)P(B)这叫独立
p(AB)=0这叫不相容
A  22:28:02
看到了吧
A  22:28:08
N-M
A  22:28:24
http://baike.baidu.com/view/984313.htm
B  22:28:59
一共你看这个公式的字母里面自己算的
B  22:29:17
这是公式的构架，
A  22:29:49
那和刚才的那个例子是一样的呀
B  22:30:20
N-M没错，
但这是公式自己算的
之前N就是所有的，M就是所有的里面的M个
不用你算，公式会自己算
你要算完就不对了

B  22:31:18
你自己吧M的值变成N-M了，带到公式里面公式还要再减一下那不就变成N-M-M了吗
A  22:31:47
哦
B  22:32:04
假定在N件产品中有M件不合格品，机抽n件做检查，发现X件是不合格品
这几个数就不要再动了
A  22:32:07
1-phyper(50,800,10000,500)
A  22:32:30
计算时，它会自动减的？
B  22:32:51
不信你看看原代码
B  22:32:52
哈哈
A  22:33:06
R的，你会吧
B  22:33:22
或者你再看看参数说明
A  22:34:01
Hypergeometric {stats} R Documentation
The Hypergeometric Distribution
Description

Density, distribution function, quantile function and random generation for the hypergeometric distribution.
Usage

dhyper(x, m, n, k, log = FALSE)
phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)
qhyper(p, m, n, k, lower.tail = TRUE, log.p = FALSE)
rhyper(nn, m, n, k)

Arguments
x, q  vector of quantiles representing the number of white balls drawn without replacement from an urn which contains both black and white balls.
m  the number of white balls in the urn.
n  the number of black balls in the urn.
k  the number of balls drawn from the urn.
p  probability, it must be between 0 and 1.
nn  number of observations. If length(nn) > 1, the length is taken to be the number required.
log, log.p  logical; if TRUE, probabilities p are given as log(p).
lower.tail  logical; if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].
A  22:34:22
m  the number of white balls in the urn.
n  the number of black balls in the urn.
A  22:34:34
m ,n相互独立的
A  22:37:26
怎么解释呢
B  22:37:38
稍等，我看看
A  22:38:09
B  22:32:04
假定在N件产品中有M件不合格品，机抽n件做检查，发现X件是不合格品
这几个数就不要再动了
是不是要 N-M
A  22:39:06
背景10000
其中500个差异
一个通路相当于一次抽样，抽出800个
其中差异的50个
应该是 10000-500
1-phyper(50,800,10000-500,500)

A  22:39:55
？
B  22:40:02
。。。。
A  22:40:25
哦，我以为你走了
是不是
背景10000
其中500个差异
一个通路相当于一次抽样，抽出800个
其中差异的50个
应该是 10000-500
1-phyper(50,800,10000-500,500)

A  22:41:06
1-phyper(50,500,10000-500,800)？
A  22:41:56
1-phyper(50,500,10000-500,800) 应该是这样才对
A  22:41:59
哈哈
B  22:43:53
m  the number of white balls in the urn.
n  the number of black balls in the urn.
k  the number of balls drawn from the urn.
B  22:44:23
k  the number of balls drawn from the urn.
这个肯定是抽样
A  22:44:28
是的
A  22:45:08
背景10000
其中500个差异
一个通路相当于一次抽样，抽出800个
其中差异的50个
1-phyper(50,500,10000-500,800) 应该是这样才对

A  22:45:09
对吧
A  22:45:25
这样算出P-VALUE
B  22:46:33
别着急，我理顺下
x, q  vector of quantiles representing the number of white balls drawn without replacement from an urn which contains both black and white balls.
这个是抽样中白的

B  22:46:37
对吧
A  22:46:43
好的
B  22:47:03
那现在就差一个背景和一个背景中白球
B  22:47:57
m  the number of white balls in the urn.
n  the number of black balls in the urn.

这个按字面理解应该是你说的意思
A  22:47:58
http://stat.ethz.ch/R-manual/R-patched/library/stats/html/Hypergeometric.html
A  22:48:02
你看这
B  22:49:02
我觉得你还是用超几何原始公式算下，
然后再用这个命令算下，对下结果
A  22:49:20
好的
B  22:49:26
R包本来就是开源的，谁传上的也不知道
B  22:50:41
而且，我觉得对于一个自然的程序来说，方便用的话不应该用你自己计算一些东西
B  22:50:43
呵呵
B  22:51:09
保险起见，自己用公式算下，再用他命令对下，免的错了

wanp_cnu

“计算得到的pvalue通过Bonferroni校正之后，以corrected-pvalue≤0.05为阈值，满足此条件的GO term定义为在差异表达基因中显著富集的GO term。”
Bonferroni校正如何做？
谢谢！

xiaoliu

bonferroni校正 ，来源于如下错误叠加，  1-（1-a)^n=a'
本来单次检验的时候，a为0.05，但次数多了导致总体的错误率多于0.05，所以要想让a'=0.05，反推a就是你校正后的显著水平

还有一种更简单的，就是   校正后的显著水平=0.05/检验次数