bioperl 计算 pi 和 Theta
2012-10-03 06:22:06| 分类:
生物信息编程
| 标签:
|举报
|字号大中小 订阅
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
评论这张
转发至微博
转发至微博
评论