/usr/bin/bp_heterogeneity_test is in bioperl 1.6.924-1.
This file is owned by root:root, with mode 0o755.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 | #!/usr/bin/perl
eval 'exec /usr/bin/perl -S $0 ${1+"$@"}'
if 0; # not running under some shell
# -*-Perl-*- (for my emacs)
use strict;
use warnings;
=head1 NAME
bp_heterogeneity_test - a test for distinguishing between selection and population expansion.
=head1 SYNOPSIS
heterogenetity_test -mut_1/--mutsyn synonymous_mut_count -mut_2/--mutnon nonsyn_mut_count -s/--smaplesize sample_size [-i/--iterations iterations] [-o/--observed observed_D] [-v/--verbose] [--silent ] [-m/--method tajimaD or fuD] [--precision]
=head2 DESCRIPTION
This is an implementation of the Heterogenetity test as described in
Hahn MW, Rausher MD, and Cunningham CW. 2002. Genetics 161(1):11-20.
=head2 OPTIONS
Options in brackets above are optional
-s or --samplesize samplesize
-mut_1 or --mutsyn synonymous mutation count
-mut_2 or --mutnon nonsynonmous mutation count
-i or --iterations number of iterations
-o or --observed observed D
-m or --method tajimaD or fuD for Tajima's D or Fu and Li's D
-v or --verbose print out extra verbose messages
--silent Be extra quiet
--precision Level of precision - specify the number of digits
(default 4)
=head2 AUTHOR Matthew Hahn E<lt>matthew.hahn-at-duke.eduE<gt>
For more information contact:
Matthew Hahn, E<lt>matthew.hahn-at-duke.eduE<gt>
Jason Stajich E<lt>jason-at-bioperl-dot-orgE<gt>
=cut
use Getopt::Long;
use Bio::PopGen::Simulation::Coalescent;
use Bio::PopGen::Statistics;
use Bio::PopGen::Individual;
use Bio::PopGen::Genotype;
my $sample_size = 4;
my $mut_count_1 = 10; # synonymous
my $mut_count_2 = 20; # non-synonymous
my $iterations = 1;
my $verbose = 0;
my $observedD = undef;
my $method = 'fuD';
my $help = 0;
my $precision = '4'; # Let's make the random precision between
# 0->1 to 1000th digits
GetOptions(
's|samplesize|samp_size:i' => \$sample_size,
'mut_1|mutsyn:i' => \$mut_count_1,
'mut_2|mutnon:i' => \$mut_count_2,
'i|iterations:i' => \$iterations,
'o|obsered|observedD:f' => \$observedD,
'v|verbose' => \$verbose,
'm|method:s' => \$method,
'h|help' => \$help,
'silent' => sub { $verbose = -1; },
'p|precision:i' => \$precision,
);
if( $help ) {
system("perldoc",$0);
exit(0);
}
if( $method ne 'fuD' and $method ne 'tajimaD' ) {
die("available methods are [fu and li's D] (fuD) and Tajima's D (tajimaD)");
}
my @D_distribution;
printf("sample size is %d iteration count = %d\n", $sample_size,
$iterations);
my $sim = new Bio::PopGen::Simulation::Coalescent
(-sample_size => $sample_size);
for(my $iter = 0; $iter < $iterations; $iter++ ) {
my $tree = $sim->next_tree($sample_size);
my $f1 = 0;
if( $mut_count_1 > 0 ) {
$sim->add_Mutations($tree,$mut_count_1,$precision);
my @leaves = $tree->get_leaf_nodes;
# the outgroup is just an individual with the ancestral state
# (no mutations)
my $outgroup = new Bio::PopGen::Individual();
foreach my $m ( $leaves[0]->get_marker_names ) {
$outgroup->add_Genotype(Bio::PopGen::Genotype->new
(-marker_name=> $m,
-alleles => [ 0 ]));
}
if( $method eq 'fuD' ) {
$f1 = Bio::PopGen::Statistics->fu_and_li_D(\@leaves,[$outgroup]);
} elsif( $method eq 'tajimaD' ) {
$f1 = Bio::PopGen::Statistics->tajima_D(\@leaves);
}
print "(mutation count = $mut_count_1) D=$f1\n"
if( $verbose >= 0);
}
my $f2 = 0;
if( $mut_count_2 > 0 ) {
$sim->add_Mutations($tree,$mut_count_2,$precision);
my @leaves = $tree->get_leaf_nodes;
# the outgroup is just an individual with the ancestral state
# (no mutations)
my $outgroup = new Bio::PopGen::Individual();
foreach my $m ( $leaves[0]->get_marker_names ) {
$outgroup->add_Genotype(Bio::PopGen::Genotype->new
(-marker_name=> $m,
-alleles => [ 0 ]));
}
if( $method eq 'fuD' ) {
$f2 = Bio::PopGen::Statistics->fu_and_li_D(\@leaves,[$outgroup]);
} elsif( $method eq 'tajimaD' ) {
$f2 = Bio::PopGen::Statistics->tajima_D(\@leaves);
}
print "(mutation count = $mut_count_2) D=$f2\n" if( $verbose >= 0);
}
my $deltaD = ( $f1 - $f2 );
push @D_distribution, $deltaD;
if( $iter % 10 == 0 && $iter > 0 ) {
print STDERR "iter = $iter\n";
}
}
if( defined $observedD && $iterations > 1 ) {
my @sortedD = sort { $a <=> $b } @D_distribution;
my $i;
for($i = 0; $i < scalar @sortedD; $i++ ) {
if( $sortedD[$i] > $observedD ) {
last;
}
}
printf( "index %d value=%.4f out of %d total (obs=%.4f)\n",
$i, $sortedD[$i], scalar @sortedD, $observedD);
}
|