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

Bioinformatics home

 
 
 

日志

 
 

bioperl 计算 pi 和 Theta  

2012-10-03 06:22:06|  分类: 生物信息编程 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
use Bio::PopGen::IO;
use Bio::PopGen::Statistics;
use Bio::AlignIO;
use Bio::PopGen::Simulation::Coalescent;
 
 
my $stats = Bio::PopGen::Statistics->new();
my $io = Bio::PopGen::IO->new(-format => 'prettybase',
                              -file     =>  'data.txt');          
                              
                              
if( my $pop = $io->next_population ) {
    my $pi = $stats->pi($pop);
    my $theta  = $stats->theta($pop);
    print "pi is $pi   ; theta is $theta\n";

    # to generate pi just for 3 of the individuals;
    my @inds;
    for my $ind ( $pop->get_Individuals ) {
        if( $ind->unique_id =~ /A0[1-3]/ ) {
            push @inds, $ind;
        }
    }
    print "pi for inds 1,2,3 is ", $stats->pi(\@inds),"\n";

}


result:
pi is 2.26666666666667   ; theta is 2.18978102189781
pi for inds 1,2,3 is 2

##############################################
‘data.txt’ (tab-delimited)内容为:
1    A01    A
1    A02    A
1    A03    A
1    A04    A
1    A05    A
2    A01    A
2    A02    T
2    A03    T
2    A04    T
2    A05    T
4    A01    G
4    A02    G
4    A03    C
4    A04    C
4    A05    G
5    A01    T
5    A02    C
5    A03    T
5    A04    T
5    A05    T
11    A01    G
11    A02    G
11    A03    G
11    A04    A
11    A05    A

1    out    G
2    out    A
4    out    G
5    out    T
11    out    G
  评论这张
 
阅读(937)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

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

页脚

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