t/PopGen.t

   1 # -*-Perl-*- mode for emacs
   2 # $Id$
   3
   4 # This will outline many tests for the population genetics
   5 # objects in the Bio::PopGen namespace
   6
   7 my $error;
   8
   9 use vars qw($SKIPXML $LASTXMLTEST);
  10 use strict;
  11 use lib '.';
  12
  13 BEGIN {
  14     # to handle systems with no installed Test module
  15     # we include the t dir (where a copy of Test.pm is located)
  16     # as a fallback
  17     eval { require Test; };
  18     if( $@ ) {
  19         use lib 't';
  20     }
  21     use vars qw($NTESTS);
  22     $NTESTS = 89;
  23     $error = 0;
  24
  25     use Test;
  26     plan tests => $NTESTS;
  27
  28 }
  29
  30 if( $error == 1 ) {
  31     exit(0);
  32 }
  33
  34
  35 use Bio::PopGen::Individual;
  36 use Bio::PopGen::Genotype;
  37 use Bio::PopGen::Population;
  38 use Bio::PopGen::IO;
  39 use Bio::PopGen::PopStats;
  40 use Bio::AlignIO;
  41 use Bio::PopGen::Statistics;
  42 use Bio::PopGen::Utilities;
  43
  44
  45 # test Fu and Li's D using data from the paper
  46
  47 ok(sprintf("%.3f",
  48            Bio::PopGen::Statistics->fu_and_li_D_counts(24, 18, 9)),
  49    -1.529);
  50
  51 ok(sprintf("%.3f",
  52            Bio::PopGen::Statistics->fu_and_li_D_star_counts(24, 18, 10)),
  53    -1.558);
  54
  55 ok(sprintf("%.3f",
  56            Bio::PopGen::Statistics->fu_and_li_F_counts(24, 3.16, 18, 9)),
  57    -1.735);
  58
  59 ok(sprintf("%.2f",
  60            Bio::PopGen::Statistics->fu_and_li_F_star_counts(24, 3.16, 18, 10)),
  61    -1.71);
  62
  63 my ($FILE1) = qw(popgentst1.out);
  64
  65 END {
  66     # unlink($FILE1);
  67 }
  68 my @individuals = ( new Bio::PopGen::Individual(-unique_id => '10a'));
  69 ok($individuals[0]);
  70
  71 my @genotypes = ( new Bio::PopGen::Genotype(-marker_name    => 'Mkr1',
  72                                             -individual_id  => '10a',
  73                                             -alleles => [ qw(A a)]),
  74                   new Bio::PopGen::Genotype(-marker_name    => 'Mkr2',
  75                                             -individual_id  => '10a',
  76                                             -alleles => [ qw(B B)]),
  77                   new Bio::PopGen::Genotype(-marker_name    => 'Mkr3',
  78                                             -individual_id  => '10a',
  79                                             -alleles => [ qw(A a)]));
  80 ok(($genotypes[1]->get_Alleles)[0], 'B');
  81
  82 $individuals[0]->add_Genotype(@genotypes);
  83 ok($individuals[0]->get_Genotypes,3);
  84 ok($individuals[0]->get_Genotypes(-marker => 'Mkr3')->get_Alleles(),2);
  85 my @alleles = $individuals[0]->get_Genotypes(-marker => 'Mkr2')->get_Alleles();
  86 ok($alleles[0], 'B');
  87
  88
  89 my $population = new Bio::PopGen::Population(-name        => 'TestPop1',
  90                                              -source      => 'testjasondata',
  91                                              -description => 'throw away example',
  92                                              -individuals => \@individuals);
  93
  94 ok(scalar ($population->get_Individuals()), 1);
  95 ok($population->name, 'TestPop1');
  96 ok($population->source, 'testjasondata');
  97 ok($population->description, 'throw away example');
  98
  99 my @genotypes2 = ( new Bio::PopGen::Genotype(-marker_name   => 'Mkr1',
 100                                              -individual_id => '11',
 101                                              -alleles       => [ qw(A A)]),
 102                    new Bio::PopGen::Genotype(-marker_name   => 'Mkr2',
 103                                              -individual_id => '11',
 104                                              -alleles       => [ qw(B B)]),
 105                    new Bio::PopGen::Genotype(-marker_name   => 'Mkr3',
 106                                              -individual_id => '11',
 107                                              -alleles       => [ qw(a a)]),
 108                    new Bio::PopGen::Genotype(-marker_name   => 'Mkr4',
 109                                              -individual_id => '11',
 110                                              -alleles       => [ qw(C C)])
 111                    );
 112 push @individuals, new Bio::PopGen::Individual(-genotypes   => \@genotypes2,
 113                                                -unique_id   => '11');
 114 $population->add_Individual($individuals[1]);
 115
 116 ok(scalar ($population->get_Individuals()), 2);
 117 my ($found_ind) = $population->get_Individuals(-unique_id => '10a');
 118 ok($found_ind->unique_id, '10a');
 119 ok(scalar($population->get_Individuals(-marker => 'Mkr4')) , 1);
 120 ok(scalar($population->get_Individuals(-marker => 'Mkr3')) , 2);
 121
 122 my @g = $population->get_Genotypes(-marker => 'Mkr4');
 123
 124 ok($g[0]->individual_id, '11');
 125 ok(($g[0]->get_Alleles())[0], 'C');
 126
 127 my $marker = $population->get_Marker('Mkr3');
 128 ok($marker);
 129
 130 @alleles = $marker->get_Alleles;
 131 ok(@alleles,2);
 132 my %af = $marker->get_Allele_Frequencies();
 133 ok($af{'a'}, 0.75);
 134 ok($af{'A'}, 0.25);
 135
 136
 137 # Read in data from a file
 138 my $io = new Bio::PopGen::IO(-format => 'csv',
 139                              -file   => Bio::Root::IO->catfile(qw(t data
 140                                                                   popgen_saureus.dat)));
 141
 142 my @inds;
 143 while( my $ind = $io->next_individual ) {
 144     push @inds, $ind;
 145 }
 146
 147 my @mrsainds = grep { $_->unique_id =~ /^MRSA/ } @inds;
 148 my @mssainds = grep { $_->unique_id =~ /^MSSA/ } @inds;
 149 my @envinds = grep { $_->unique_id =~ /^NC/ } @inds;
 150
 151 ok(scalar @mrsainds, 9);
 152 ok(scalar @mssainds, 10);
 153 ok(scalar @envinds, 5);
 154
 155 my $mrsapop = new Bio::PopGen::Population(-name        => 'MRSA',
 156                                           -description => 'Resistant S.aureus',
 157                                           -individuals => \@mrsainds);
 158
 159 my $mssapop = new Bio::PopGen::Population(-name        => 'MSSA',
 160                                           -description =>'Suceptible S.aureus',
 161                                           -individuals => \@mssainds);
 162
 163 my $envpop = new Bio::PopGen::Population(-name        => 'NC',
 164                                          -description => 'WT isolates',
 165                                           -individuals => \@envinds);
 166
 167 my $stats = new Bio::PopGen::PopStats(-haploid => 1);
 168 my $fst = $stats->Fst([$mrsapop,$mssapop],[qw(AFLP1)]);
 169 # We're going to check the values against other programs first
 170 ok(sprintf("%.3f",$fst),0.077,'mrsa,mssa aflp1');
 171
 172 $fst = $stats->Fst([$envpop,$mssapop,$mrsapop],[qw(AFLP1 )]);
 173 ok(sprintf("%.3f",$fst),0.035,'all pops, aflp1');
 174
 175 $fst = $stats->Fst([$mrsapop,$envpop],[qw(AFLP1 AFLP2)]);
 176 ok(sprintf("%.3f",$fst),0.046,'mrsa,envpop aflp1,aflp2');
 177
 178 # Read in data from a file
 179 $io = new Bio::PopGen::IO(-format => 'csv',
 180                           -file   => Bio::Root::IO->catfile
 181                           (qw(t data popgen_saureus.multidat)));
 182
 183 @inds = ();
 184 while( my $ind = $io->next_individual ) {
 185     push @inds, $ind;
 186 }
 187
 188 @mrsainds = grep { $_->unique_id =~ /^MRSA/ } @inds;
 189 @mssainds = grep { $_->unique_id =~ /^MSSA/ } @inds;
 190 @envinds = grep { $_->unique_id =~ /^NC/ } @inds;
 191
 192 ok(scalar @mrsainds, 7);
 193 ok(scalar @mssainds, 10);
 194 ok(scalar @envinds, 5);
 195
 196 $mrsapop = new Bio::PopGen::Population(-name        => 'MRSA',
 197                                        -description => 'Resistant S.aureus',
 198                                        -individuals => \@mrsainds);
 199
 200 $mssapop = new Bio::PopGen::Population(-name        => 'MSSA',
 201                                        -description =>'Suceptible S.aureus',
 202                                        -individuals => \@mssainds);
 203
 204 $envpop = new Bio::PopGen::Population(-name        => 'NC',
 205                                       -description => 'WT isolates',
 206                                       -individuals => \@envinds);
 207
 208 $stats = new Bio::PopGen::PopStats(-haploid => 1);
 209 my @all_bands = map { 'B' . $_ } 1..20;
 210 my @mkr1     = map { 'B' . $_ } 1..13;
 211 my @mkr2     = map { 'B' . $_ } 14..20;
 212
 213 # still wrong ?
 214 $fst = $stats->Fst([$mrsapop,$mssapop],[@all_bands ]);
 215 skip(sprintf("%.3f",$fst),'-0.001','mssa,mrsa all_bands'); # We're going to check the values against other programs first
 216 $fst = $stats->Fst([$envpop,$mssapop],[ @mkr1 ]);
 217 ok(sprintf("%.3f",$fst),0.023,'env,mssa mkr1'); # We're going to check the values against other programs first
 218
 219 $fst = $stats->Fst([$envpop,$mssapop,$mrsapop],[ @all_bands ]);
 220 ok(sprintf("%.3f",$fst),0.071,'env,mssa,mrsa all bands'); # We're going to check the values against other programs first
 221
 222 $fst = $stats->Fst([$envpop,$mssapop,$mrsapop],[ @mkr2 ]);
 223 ok(sprintf("%.3f",$fst),0.076, 'env,mssa,mrsa mkr2'); # We're going to check the values against other programs first
 224
 225 $fst = $stats->Fst([$mrsapop,$envpop],[@all_bands ]);
 226 ok(sprintf("%.3f",$fst),0.241,'mrsa,nc all_bands'); # We're going to check the values against other programs first
 227
 228 # test overall allele freq setting for a population
 229
 230 my $poptst1 = new Bio::PopGen::Population(-name => 'tst1');
 231 my $poptst2 = new Bio::PopGen::Population(-name => 'tst2');
 232
 233 $poptst1->set_Allele_Frequency(-frequencies =>
 234                                { 'marker1' => { 'a' => '0.20',
 235                                                 'A' => '0.80' },
 236                                  'marker2' => { 'A' => '0.10',
 237                                                 'B' => '0.20',
 238                                                 'C' => '0.70' }
 239                              });
 240
 241 my $mk1 = $poptst1->get_Marker('marker1');
 242 my %f1 = $mk1->get_Allele_Frequencies;
 243 ok($f1{'a'}, '0.20');
 244 ok($f1{'A'}, '0.80');
 245 my $mk2 = $poptst1->get_Marker('marker2');
 246 my %f2 = $mk2->get_Allele_Frequencies;
 247 ok($f2{'C'}, '0.70');
 248
 249 $poptst2->set_Allele_Frequency(-name      => 'marker1',
 250                                -allele    => 'A',
 251                                -frequency => '0.60');
 252 $poptst2->set_Allele_Frequency(-name      => 'marker1',
 253                                -allele    => 'a',
 254                                -frequency => '0.40');
 255
 256 #$fst = $stats->Fst([$poptst1,$poptst2],[qw(marker1 marker2) ]);
 257 skip('Fst not calculated yet for just allele freqs',1,'marker1 test'); #
 258
 259 $io = new Bio::PopGen::IO(-format => 'csv',
 260                           -file   => ">$FILE1");
 261
 262 $io->write_individual(@inds);
 263 $io->close();
 264 ok( -e $FILE1);
 265 unlink($FILE1);
 266 $io = new Bio::PopGen::IO(-format => 'csv',
 267                           -file   => ">$FILE1");
 268
 269 $io->write_population(($mssapop,$mrsapop));
 270 $io->close();
 271 ok( -e $FILE1);
 272 unlink($FILE1);
 273
 274 $io = new Bio::PopGen::IO(-format => 'prettybase',
 275                           -file   => ">$FILE1");
 276
 277 $io->write_individual(@inds);
 278 $io->close();
 279 ok( -e $FILE1);
 280 unlink($FILE1);
 281
 282 $io = new Bio::PopGen::IO(-format => 'prettybase',
 283                           -file   => ">$FILE1");
 284
 285 $io->write_population(($mssapop,$mrsapop));
 286 $io->close();
 287 ok( -e $FILE1);
 288 unlink($FILE1);
 289
 290
 291 # Let's do PopGen::Statistics tests here
 292
 293 $io = new Bio::PopGen::IO(-format          => 'prettybase',
 294                           -no_header       => 1,
 295                           -file            => Bio::Root::IO->catfile
 296                           (qw(t data popstats.prettybase )));
 297 my (@ingroup,@outgroup);
 298 my $sitecount;
 299 while( my $ind = $io->next_individual ) {
 300     if($ind->unique_id =~ /out/) {
 301         push @outgroup, $ind;
 302     } else {
 303         push @ingroup, $ind;
 304         $sitecount = scalar $ind->get_marker_names() unless defined $sitecount;
 305     }
 306 }
 307 $stats = new Bio::PopGen::Statistics();
 308
 309 # Real data and values courtesy M.Hahn and DNASP
 310
 311 ok($stats->pi(\@ingroup),2);
 312 ok(Bio::PopGen::Statistics->pi(\@ingroup,$sitecount),0.4);
 313
 314 ok(Bio::PopGen::Statistics->theta(\@ingroup),1.92);
 315 ok(Bio::PopGen::Statistics->theta(\@ingroup,$sitecount),0.384);
 316
 317 # Test with a population object
 318 my $ingroup  = new Bio::PopGen::Population(-individuals => \@ingroup);
 319 my $outgroup = new Bio::PopGen::Population(-individuals => \@outgroup);
 320
 321 ok($stats->pi($ingroup),2);
 322 ok(Bio::PopGen::Statistics->pi($ingroup,$sitecount),0.4);
 323
 324 ok(Bio::PopGen::Statistics->theta($ingroup),1.92);
 325 ok(Bio::PopGen::Statistics->theta($ingroup,$sitecount),0.384);
 326
 327 my $haploidpop = $ingroup->haploid_population;
 328 ok(sprintf("%.5f",Bio::PopGen::Statistics->tajima_D($haploidpop)), 0.27345);
 329
 330 # to fix
 331 ok(sprintf("%.5f",Bio::PopGen::Statistics->tajima_D(\@ingroup)),0.27345);
 332 ok(sprintf("%.5f",Bio::PopGen::Statistics->tajima_D($ingroup)),0.27345);
 333
 334 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_D_star(\@ingroup)),
 335    0.27345);
 336 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_D_star($ingroup)),
 337    0.27345);
 338
 339 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_F_star(\@ingroup)),
 340      0.27834);
 341 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_F_star($ingroup)),
 342    0.27834);
 343
 344 ok((Bio::PopGen::Statistics->derived_mutations(\@ingroup,\@outgroup))[0], 1);
 345 ok((Bio::PopGen::Statistics->derived_mutations($ingroup,\@outgroup))[0], 1);
 346 ok((Bio::PopGen::Statistics->derived_mutations(\@ingroup,$outgroup))[0], 1);
 347 ok((Bio::PopGen::Statistics->derived_mutations($ingroup,$outgroup))[0], 1);
 348
 349 # expect to have 1 external mutation
 350 @ingroup = $haploidpop->get_Individuals;
 351
 352 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_D(\@ingroup,1)),0.75653);
 353 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_D($ingroup,1)),0.75653);
 354
 355 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_D(\@ingroup,
 356                                                        \@outgroup)),0.75653);
 357 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_D($ingroup,
 358                                                        \@outgroup)),0.75653);
 359 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_D($ingroup,
 360                                                        $outgroup)),0.75653);
 361 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_D(\@ingroup,
 362                                                        $outgroup)),0.75653);
 363
 364 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_F(\@ingroup,1)),
 365      0.77499);
 366 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_F($haploidpop,1)),0.77499);
 367 ok(sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_F($ingroup,
 368                                                        \@outgroup)),0.77499);
 369 ok( sprintf("%.5f",Bio::PopGen::Statistics->fu_and_li_F($ingroup,
 370                                                          $outgroup)),0.77499);
 371
 372
 373 # Test composite LD
 374
 375 $io = new Bio::PopGen::IO(-format => 'prettybase',
 376                           -file   => Bio::Root::IO->catfile
 377                           (qw(t data compLD_test.prettybase)));
 378
 379 my $pop = $io->next_population;
 380
 381 my %LD = $stats->composite_LD($pop);
 382
 383 ok($LD{'01'}->{'02'}->[1], 10);
 384 ok($LD{'01'}->{'03'}->[1], 0);
 385 ok($LD{'02'}->{'03'}->[1], 0);
 386
 387 # Test composite LD
 388
 389 $io = new Bio::PopGen::IO(-format => 'prettybase',
 390                           -file   => Bio::Root::IO->catfile
 391                           (qw(t data compLD_missingtest.prettybase)));
 392
 393 $pop = $io->next_population;
 394
 395 %LD = $stats->composite_LD($pop);
 396
 397 ok(sprintf("%.4f",$LD{'ProC9198EA'}->{'ProcR2973EA'}->[0]), -0.0375);
 398 ok(sprintf("%.2f",$LD{'ProC9198EA'}->{'ProcR2973EA'}->[1]), 2.56);
 399
 400
 401
 402 # build a population from an alignment
 403
 404 my $alnin = Bio::AlignIO->new(-format => 'clustalw',
 405                               -file   => Bio::Root::IO->catfile(qw(t data T7.aln)));
 406 my $aln = $alnin->next_aln;
 407 $population = Bio::PopGen::Utilities->aln_to_population(-alignment => $aln);
 408 ok($population->get_number_individuals,9);
 409 #warn($aln->match_line,"\n");
 410 my $matchline = $aln->match_line;
 411 ok( $population->get_marker_names, $matchline =~ tr/ //);
 412 for my $name ( $population->get_marker_names ) {
 413     my $marker = $population->get_Marker($name);
 414 #    warn("$name ",join(" ",$marker->get_Alleles()),"\n");
 415 }
 416
 417
 418 # test Rich's phase and hap parsers
 419
 420 $io = new Bio::PopGen::IO(-format   => 'hapmap',
 421                           -verbose  => 1,
 422                           -no_header=> 1,
 423                           -starting_column => 10,
 424                           -file     => Bio::Root::IO->catfile(qw(t data
 425                                                                 example.hap)));
 426
 427 # Some IO might support reading in a population at a time
 428
 429 my @population;
 430 while( my $ind = $io->next_individual ) {
 431     push @population, $ind;
 432 }
 433 ok(@population, 90);
 434 ok($population[3]->unique_id, 'NA06994');
 435 ok($population[3]->get_Genotypes, 34);
 436 $population = Bio::PopGen::Population->new(-individuals => \@population);
 437
 438 ok(sprintf("%.3f",$stats->pi($population)),12.335);
 439 # if forced haploid population is called within pi
 440 # need to decide about that...
 441 # ok(sprintf("%.3f",$stats->pi($population)),12.266);
 442
 443 ok(sprintf("%.3f",$stats->theta($population)),5.548);
 444 skip(1,'tjd inconsistency, need to recalculate');
 445 skip(1,'tjd inconsistency, need to recalculate');
 446 #ok(sprintf("%.3f",$stats->tajima_D($population)),2.926);
 447 #ok(sprintf("%.3f",$stats->tajima_D($population->haploid_population)),3.468);
 448
 449 $io = new Bio::PopGen::IO(-format => 'phase',
 450                           -file   => Bio::Root::IO->catfile(qw(t data
 451                                                                example.phase)));
 452
 453 # Some IO might support reading in a population at a time
 454
 455 @population = ();
 456 while( my $ind = $io->next_individual ) {
 457     push @population, $ind;
 458 }
 459 ok(@population, 4);
 460
 461
 462 # test diploid data
 463