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"; } }
写的好!
不知道如何区分dhyper和phyper两个函数?
dhyper是计算点的概率,phyper是计算的面积
这个超几何分布函数用phyper计算的是累计概率,如果我没理解错的话就是比如说 差异基因总数为100 其中某个GO下有6个基因 背景是10000个基因 在这个GO下有100个基因,此时phyper计算的p值就是 少于6个基因的所有概率之和吧?那p矫正之后很小 也就是说明少于6个基因的情况存在的可能性很小,说明大于6个基因富集到该GO才正常 这和富集怎么联系起来呢?
你的理解很对。如果“少于6个基因的所有概率之和”的概率很大,这个概率之和超过了95%。那么结果 6 就属于小概率事件了;这个小概率事件发生了,就表示这个事件可能性太低了,原因是存在着富集。
你好,我刚刚开始做RNA-seq的数据分析,由于自学,所以经常参考你网站的一些的资料,非常感谢你的分享。
现在我刚刚好做到富集分析,使用的数学模型中大家现在好像比较接受goseq的方法,因为它把基因的长度差异进行了考虑,采用了类似超几何分布的模型,所以打算开始用。不知到你是否有使用goseq进行过分析,使用你提到的超几何分布和goseq的分析区别大吗?(一些测序公司似乎也是采用超几何分布的方法。)
不知能否上传两个样本文件~因为不知道阁下设定的 第一个文件 (即GO 分类 )格式到底是怎样的。我用blast2go做的annotion,不知道格式符不符合。