注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

云之南

风声,雨声,读书声,声声入耳;家事,国事,天下事,事事关心

 
 
 

日志

 
 
关于我

专业背景:计算机科学 研究方向与兴趣: JavaEE-Web软件开发, 生物信息学, 数据挖掘与机器学习, 智能信息系统 目前工作: 基因组, 转录组, NGS高通量数据分析, 生物数据挖掘, 植物系统发育和比较进化基因组学

网易考拉推荐

超几何分布和fisher精确检验---富集分析原理  

2016-10-21 14:39:08|  分类: 生信分析软件 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

超几何分布和fisher精确检验

https://en.wikipedia.org/wiki/Fisher%27s_exact_test

http://wangyinanchina.github.io/blog/2014/09/03/fishers-exact-test/

http://blog.sina.com.cn/s/blog_670445240101m4z3.html

http://www.bakerwm.com/r/2015/01/23/hypergeometric-analysis--enrichment-analysis/

Fisher精确检验(fisher’s exat test)是进行统计分析时经常碰到的一种检验方法,它基于超几何分布,作用于离散变量,用于检测两种分类方法的结果是否独立。

Fisher精确检验(fisher’s exat test)是进行统计分析时经常碰到的一种检验方法,它基于超几何分布,作用于离散变量,用于检测两种分类方法的结果是否独立。

首先,我们介绍超几何分布。超几何分布用来模拟这样的过程:将有限的总体分为两类A和B,从中不放回的抽样n次,结果中A的个数符合超几何分布。所以使用古典概型的方法,假设N个总体中有A和B两类,其中A有K个,从中不放回的抽样n次,我们可以推导出n中为A的数目x,即超几何分布的pmf:

P(X=x)=(Kx)(N?Kn?x)/(Nn)

回到fisher精确检验,fisher检验要回答的问题是,对数据进行两种分类,这两种分类是否独立?即在第一种分类条件下分为某一类的数据是否更倾向于在第二种分类中归于某类。举例说明(例子来源于wiki):

我们有24位测试对象,根据其性别和是否在节食,将其分为四类,分类结果如下:


男性 女性
节食 1 9
不节食 11 3

现在我们的问题是,是否女性更喜欢节食?从数据直观上来看,男性和女性都是12人,但是节食的女性是男性的9倍,似乎女性的确更容易节食,但是我们如何定量的去描述这件事呢?可以看出,我们要解决的是一个假设检验问题,我们将零假设设定为:是否节食和性别无关。那么,在零假设下,根据超几何分布我们观察到表格中数据的概率是:

(12,1)(12,9)/(24,10)=0.001346076
而如果我们做一个单尾的检验,那么我们观测到表格中数据或者比表格中数据更极端(即节食男性为0)的概率为:
(121)(129)/(2410)+(120)(1210)/(2410)=0.001379728

即fisher精确检验的p-value。我们用R去检验我们计算的结果

fisher.test(matrix(c(1, 11, 9, 3), 2, 2), alternative = 'less')

# 结果如下,可以看出与我们计算的结果相同
Fisher's Exact Test for Count Data

data:  matrix(c(1, 11, 9, 3), 2, 2)
p-value = 0.00138
alternative hypothesis: true odds ratio is less than 1
95 percent confidence interval:
 0.0000000 0.3260026
sample estimates:
odds ratio 
0.03723312 

对于双尾的fisher检验,目前还没有好的计算方法,因为对于两个极端并不好定义,目前最简单的方法就是对两个极端,分别计算观测到观测值或者更极端的观测值得概率(例如上例中节食男性为1或者0的情况(一个极端)以及节食男性为9,10,11,12(另一个极端)的情况),并将所得到的概率相加,即为最终的p-value。


超几何分析和GO富集分析

http://www.bakerwm.com/r/2015/01/23/hypergeometric-analysis--enrichment-analysis/

  1. 关于超几何分布(Hypergeometric distribution)

超几何分布是统计学上一种离散概率分布,它描述了由有限个物件 (N) 中抽取 (n) 个物件,成功抽出制定种类物件的个数 (k) (不归还)。
例如在有N个样品,其中m个是不及格的。超几何分布描述了在该N个样本中抽出n个,其中k个是不及格的几率。

f(k;n,m,N) = C(m,k)*C(N-m, n-k)/C(N,n)

若n=1,超几何分布还原为伯努利分布。 若n接近∞,超几何分布可视为二项式分布。
From: Wikipedia

转换到生物学分析(GO)中,某样品具有GO注释的基因数 (N), 其中属于TermA的基因数目为 (M), (举例)差异表达分析后得到的基因数量为 (n),其中有 (k) 个基因属于TermA,现在需要使用超几何分析判断该TermA是否富集?!

 而对于别的注释分析,如:所有物种总的基因数是N个,其中注释到这个TermA的基因数是M个,对于某一个物种的基因数是n个,其中有k个基因属于TermA, 现在问这个物种在这个TermA是否富集呢?

f(k,N,M,n) = C(k,M)*C(n-k,N-M)/C(n,N)
N = size of population
M = # of items in population with property "E"
N-M = # of items in population without property "E"
n = number of items sampled
k = number of items in sample with property "E"

解读为:先计算出k-1的累积可能性,再用1减去该数,即可得到获得 >=k 的可能性。(个人理解)

R中的函数

1-phyper(k-1,M,N-M,n)
  1. 关于富集性分析(Enrichment analysis)

需要知道有几种方法执行富集分析:
1. 超几何分布(常用软件都使用该方法)。
2. fisher's exact test (也是使用超几何分布).
3. chi-sauared test.


1.《超几何分析》
Posted by ygc on April 28, 2012
原文链接:http://ygc.name/2008/08/20/hypergeometric-distribution/

简单点说,超几何分布就是有限样本的无放回抽样。不同于有放回抽样的二项分布(每次贝努里试验成功概率是一样的),每次的概率不相等。 随机变量X的超几何概率分布:
f(k,N,M,n) = C(k,M)*C(n-k,N-M)/C(n,N)
N = size of population
M = # of items in population with property "E"
N-M = # of items in population without property "E"
n = number of items sampled
k = number of items in sample with property "E"

这个公式可以理解为有C(n,N)种可能的样本,有C(k,M)种方法得到k个属于M的抽样、有C(n-k,N-M)种方法得到n-k个不属于M的抽样。
X服从参数n,N,M的超几何分布记为 X~H(n,N,M).

参考:http://en.wikipedia.org/wiki/Hypergeometric_distribution

对于基因进行GO注释,看基因集在某个GO子类中是否富集,富集的概率服从超几何分布。 * N为GO注释的总基因数。 * M为属于某个GO子类的基因数。 * n为进行GO富集分析的基因集的数目。 * k为n中属于M的数目。 * 基因集n是否在M类中富集的概率

1-phyper(k-1,M,N-M,n)  ##R代码

或者是

phyper(k-1,M,N-M,n, lower.tail=FALSE)

在已知总体分布下,抽样n个中出现M类的个数是k以及k以上个数的概率。

lower.tail: logical; if TRUE (default), probabilities are P[X <= x],
          otherwise, P[X > x]. 

2.《富集性分析》
Posted by ygc on April 28, 2012

原文链接: http://ygc.name/2012/04/28/enrichment-analysis/

当我们用组学测定了一大堆分子之后,我们希望站在更高的角度去看这些分子和那些生物学过程相关。那么通常各种注释,对这些基因/蛋白进行分类,那么从分类的比例上,是不能草率下结论,正如上面有钱人学历分布的例子一样。我们需要把总体的分布考虑进去。
和某个注释/分类是否有相关性,把基因分成属于这一类,和不属于这一类两种,这就好比经典统计学中的白球和黑球的抽样问题。也可以列一个2x2的表,进行独立性分析。
以文章Gene Expression in Ovarian Cancer Reflects Both Morphology and Biological Behavior, Distinguishing Clear Cell from Other Poor-Prognosis Ovarian Carcinomas: URL 所鉴定的差异基因为例。

73个差异基因的Symbol,我把它转为 entrezgene ID得到57个(漏掉的不管它,只是做为一个例子):

eg [1] 

#下面测试一下这些基因和化学刺激响应的相关性。

goid <- “GO:0042221# response to chemical stimulus

#那么做为背景,总体基因为N,属于“化学刺激响应”这个分类的基因有M个。

allgeneInCategory <- unique(get(goid, org.Hs.egGO2ALLEGS))
M <- length(allgeneInCategory)
N <- length(mappedkeys(org.Hs.egGO))

#样本的大小是n,属于“化学刺激响应”这个分类的基因有k个。

n <- length(eg)
k <- sum(eg %in% allgeneInCategory)

#白球黑球问题,最简单的莫过于用二项式分布,从总体上看,要拿到一个基因属于“化学刺激响应”这个分类的概率是M/N。那么现在抽了n个基因,里面有k个基于这个分类,p值为:

> 1-sum(sapply(0:k-1, function(i) choose(n,i) * (M/N)^i * (1-M/N)^(n-i)))
[1] 8.561432e-10

#二项式分布,是有放回的抽样,你可以多次抽到同一基因,这是不符合的。所以这个计算只能说是做为近似的估计值,无放回的抽样,符合超几何分布,通过超几何分布的计算,p值为:

> phyper(k-1,M, N-M, n, lower.tail=FALSE)
[1] 7.879194e-10

#如果用2x2表做独立性分析,表如下:

> d <- data.frame(gene.not.interest=c(M-k, N-M-n+k), gene.in.interest=c(k, n-k))
> row.names(d) <- c("In_category", "not_in_category")
> d
                gene.not.interest gene.in.interest
In_category                  2613               28
not_in_category             15310               29

#这个也有很多方法可以做检验,经典的有卡方检验和fisher's exact test。
如果用卡方检验来做。  

> chisq.test(d, )

    Pearson's Chi-squared test with Yates' continuity correction

data:  d 
X-squared = 51.3849, df = 1, p-value = 7.592e-13

对于2x2表来说,卡方检验通常也只能做为近似估计值,特别是当sample size或expected ell count比较小的时候,计算并不准确。 fisher's exact test,名副其实,真的就比较exact,因为它使用的是超几何分布来计算p值。这也是为什么fisher's exact test和超几何模式计算的p-值是一样的,

> fisher.test(d)
    Fisher's Exact Test for Count Data
data:  d 
p-value = 7.879e-10
alternative hypothesis: true odds ratio is not equal to 1 
95 percent confidence interval:
 0.1013210 0.3089718 
sample estimates:
odds ratio 
 0.1767937 

通常各种软件做GO富集性分析,都是使用超几何分布进行计算。 IPA软件则是使用fisher's exact test来检验基因在某个网络中是否富集。
从这个例子上看,chi-squared test的表现最差,binomial做为p值的近似估计还是不错的,因为计算较为简单。 富集性分析应用范围非常广,从Disease Ontology, Gene Ontology, KEGG, 到Reactome Pathway等等。





利用超几何概率分布做富集分析

http://www.chenlianfu.com/?p=1122

1. 了解超几何概率分布

参考:The Hypergeometric Distribution
超几何分布的公式为:

p(x) = choose(m, x) choose(n, k-x) / choose(m+n, k)
for x = 0, …, k.
其中, m 是桶里面白球的个数, n 是黑球的个数, k 是从桶中随机取出的球数, x 是取出球
中白球的个数.

Fisher‘s Exact Test就是依据超几何概率分布得到的。

2. 超几何分布运算方法

使用R的运算方法:

dhyper(x, m, n, k, log = FALSE)
    计算某一个点的p值
phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)
    计算一个分布的p值,默认下计算P(X<=x)
qhyper(p, m, n, k, lower.tail = TRUE, log.p = FALSE)
    计算某一个p值对应的分位数
rhyper(nn, m, n, k)
    按超几何分布给出nn的可能的模拟结果值

3. 对一组p值进行校正

参考:Adjust P-values for Multiple Comparisons
使用R的运算方法:

p.adjust(p, method = p.adjust.methods, n = length(p))

p.adjust.methods
# c(“holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”,
# “fdr”, “none”)

4. 利用上述原理来编写perl程序,来进行GO或PATHWAY富集分析

输入两个文件:第1个文件是各个基因对应的GO分类或KEGG分类;第2个文件是差异表达基因的名称。以下为KEGG富集分析的perl程序。

 #!/usr/bin/perl

=head1 Name

  kegg_enrich_analysis.pl

=head1 Description

  This program is designed to do enrichment analysis by kegg results.

=head1 Version

  Author: Chen Lianfu, chenllianfu@foxmail.com
  Version: 1.0,  Date: Sat Apr 27 02:06:02 2013

=head1 Usage

  --out         the output result name, default: enrichment_output
  --tmp         to keep the tmp files, default: false
  --fdr         set the false discovery rate, default: 0.05
  --kdn         When a kegg term was perpared for Fisher's exact 
test, the i(the num of this kegg terms have so many differential 
expression genes) should >= this value, default: 3
  --help        output help information to screen

=head1 Example

  perl ../kegg_enrich_analysis.pl annot_keggs genelist.txt

=cut

use strict;
use Getopt::Long;

my ($out,$tmp,$fdr_threshold,$kegg2degs_num_threshold,$Help);

GetOptions(
    "out:s"=>\$out,
    "fdr:s"=>\$fdr_threshold,
    "kdn:s"=>\$kegg2degs_num_threshold,
    "tmp"=>\$tmp,
    "help"=>\$Help
);

$out ||= 'enrichment_output';
$fdr_threshold ||= 0.05;
$kegg2degs_num_threshold ||= 3;

die `pod2text $0` if (@ARGV == 0 || $Help);

my $kegg2gene_file = shift @ARGV;
my $deglist = shift @ARGV;

open kegg2GENE, '<', $kegg2gene_file or die $!;
open DEGLIST, '<', $deglist or die $!;

my %kegg2genes;
my %gene2keggs;
while ( <kegg2GENE> ) {
    next unless /^(\S+).*?(K\d+)/;
    $gene2keggs{$1} .= "$2,";
    $kegg2genes{$2}{"genes"} .= "$1,";
    $kegg2genes{$2}{"num"} += 1;
}

my ( $total_deg_num, $deg_own_kegg_num );
my %deg_kegg2genes;
while ( <DEGLIST> ) {
    chomp;
    my $genename = $_;
    $total_deg_num ++;
    if ( $gene2keggs{$genename} ) {
        $deg_own_kegg_num ++;
        my @keggs = split /,/, $gene2keggs{$genename};

        foreach ( @keggs ) {
            $deg_kegg2genes{$_}{"genes"} .= "$genename,";
            $deg_kegg2genes{$_}{"num"} ++;
        }
    }
}

print "$total_deg_num differential expression genes in total, and 
$deg_own_kegg_num genes own kegg terms.\n";

my @for_enrich_keggs = keys %deg_kegg2genes;
open FOR_R, '>', "tmp_4value" or die $!;

my $kegg_num_for_fisher_test = @for_enrich_keggs;
print "$kegg_num_for_fisher_test kegg terms is perpared for Fisher
's exact test.\n";

foreach my $kegg ( @for_enrich_keggs ) {
    my $n_and_m = keys %gene2keggs;
    my $N = $deg_own_kegg_num;
    my $n = $kegg2genes{$kegg}{"num"};
    my $i = $deg_kegg2genes{$kegg}{"num"};
    my $m = $n_and_m -$n;

    print FOR_R "$kegg $i $n $m $N\n" if $i >= $kegg2degs_num_threshold;
}

open R_SCRIPT, '>', "phyper_adjust.R" or die $!;

print R_SCRIPT 
'phy <- read.table(file="tmp_4value")
pvalue <- phyper(phy$V2,phy$V3,phy$V4,phy$V5,lower.tail=FALSE)
qvalue <- p.adjust(pvalue,method=\'fdr\')
value <- cbind ( pvalue, qvalue )
rownames(value)=phy$V1
value
';

my @fdr = `cat phyper_adjust.R | R --vanilla --slave`;

#print @fdr;

`rm phyper_adjust.R tmp_4value` unless $tmp;

open OUTPUT, '>', $out or die $!;

foreach ( @fdr ) {
    next unless /(\S+)\s+(\S+)\s+(\S+)/;
    my $kegg = $1; my $pvalue = $2; my $fdr = $3;

    if ( $fdr <= $fdr_threshold ) {
        my $n_and_m = keys %gene2keggs;
        my $N = $deg_own_kegg_num;
        my $n = $kegg2genes{$kegg}{"num"};
        my $i = $deg_kegg2genes{$kegg}{"num"};
        my $m = $n_and_m -$n;

        my $genes = $deg_kegg2genes{$kegg}{"genes"};

        print OUTPUT "$kegg\t$i\t$N\t$n\t$n_and_m\t$pvalue\t$fdr\t$genes\n";
    }
}
  评论这张
 
阅读(146)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016