t/PopGen/MK.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 # This will outline many tests for the population genetics
   5 # objects in the Bio::PopGen namespace
   6
   7 use strict;
   8
   9 BEGIN {
  10     use lib '.';
  11     use Bio::Root::Test;
  12
  13     test_begin(-tests => 46);
  14     use_ok('Bio::AlignIO');
  15     use_ok('Bio::PopGen::Statistics');
  16     use_ok('Bio::PopGen::Utilities');
  17
  18 }
  19 my $verbose = test_debug();
  20
  21
  22 # - McDonald Kreitman tests -
  23
  24 my $stats = Bio::PopGen::Statistics->new(-verbose => $verbose);
  25 isa_ok($stats, 'Bio::PopGen::Statistics');
  26 my $alnio = Bio::AlignIO->new(-format => 'fasta',
  27                             -file   => test_input_file('CG2865.fasaln'));
  28 my $aln = $alnio->next_aln;
  29 isa_ok($aln,'Bio::SimpleAlign');
  30 my $population = Bio::PopGen::Utilities->aln_to_population(-alignment => $aln,
  31                                                         -site_model=> 'codon');
  32 isa_ok($population,'Bio::PopGen::Population');
  33
  34 my @marker_names = $population->get_marker_names;
  35 my @inds = $population->get_Individuals;
  36
  37 is(scalar @marker_names, 434, 'Marker Names');
  38 is(scalar @inds, 6,'Number of Inds');
  39
  40 my (@ingroup_seqs,@outgroup_seqs1, @outgroup_seqs2);
  41
  42 for my $ind ( $population->get_Individuals ) {
  43     my $id = $ind->unique_id;
  44     # do we allow ingroup to be a list as well?
  45     my $pushed = 0;
  46     if( $id =~ /sim/ ) {
  47         push @ingroup_seqs, $ind;
  48         $pushed++;
  49     }
  50     if( ! $pushed ) {
  51         if( $id =~ /mel/ ) {
  52             push @outgroup_seqs1, $ind;
  53             $pushed++;
  54         } elsif( $id =~ /yak/ ) {
  55             push @outgroup_seqs2, $ind;
  56             $pushed++;
  57         }
  58     }
  59     #    if( ! $pushed ) {
  60     #   warn("sequence $id was not grouped, ignoring...\n");
  61     #       push @ingroup_seqs, $ind;
  62     #}
  63 }
  64 is(scalar @ingroup_seqs, 4, 'number of ingroup sequences');
  65 is(scalar @outgroup_seqs1, 1, 'number of outgroup1 sequences');
  66 is(scalar @outgroup_seqs2, 1, 'number of outgroup2 sequences');
  67
  68 my $polarized = 0;
  69
  70 my @counts = $stats->mcdonald_kreitman(-ingroup   => \@ingroup_seqs,
  71                                        -outgroup  => \@outgroup_seqs1,
  72                                        -polarized => $polarized);
  73 is($counts[0], 0, 'NSpoly');
  74 is($counts[1], 1, 'NSfixed');
  75 is($counts[2], 3, 'Spoly');
  76 is($counts[3], 7, 'Sfixed');
  77 my $mk;
  78  SKIP: {
  79      test_skip(-tests => 1,
  80                -requires_module => 'Text::NSP::Measures::2D::Fisher2::twotailed');
  81      skip "Some problem with Bio::PopGen::Statistics::has_twotailed", 1 unless $Bio::PopGen::Statistics::has_twotailed;
  82
  83      $mk = $stats->mcdonald_kreitman_counts(@counts);
  84      is($mk, 1, 'McDonald Kreitman');
  85  }
  86 @counts = $stats->mcdonald_kreitman(-ingroup   => \@ingroup_seqs,
  87                                     -outgroup  => \@outgroup_seqs2,
  88                                     -polarized => $polarized);
  89 is($counts[0], 0, 'NSpoly');
  90 is($counts[1], 6, 'NSfixed');
  91 is($counts[2], 3, 'Spoly');
  92 is($counts[3], 16, 'Sfixed');
  93
  94 SKIP: {
  95         test_skip(-tests => 1,
  96                           -requires_module => 'Text::NSP::Measures::2D::Fisher2::twotailed');
  97         skip "Some problem with Bio::PopGen::Statistics::has_twotailed", 1 unless $Bio::PopGen::Statistics::has_twotailed;
  98
  99     $mk = $stats->mcdonald_kreitman_counts(@counts);
 100     is(sprintf("%.2f",$mk), 0.55, 'McDonald Kreitman');
 101 }
 102
 103 @counts = $stats->mcdonald_kreitman(-ingroup  => \@ingroup_seqs,
 104                                     -outgroup => [@outgroup_seqs1,
 105                                                   @outgroup_seqs2],
 106                                     -polarized=> 1);
 107 is($counts[0], 0, 'NSpoly');
 108 is($counts[1], 1, 'NSfixed');
 109 is($counts[2], 3, 'Spoly');
 110 is($counts[3], 1, 'Sfixed');
 111
 112
 113 # test 2nd aln file
 114 $alnio = Bio::AlignIO->new(-format => 'fasta',
 115                            -file   => test_input_file('CG11099.fasaln'));
 116 $aln = $alnio->next_aln;
 117 isa_ok($aln,'Bio::SimpleAlign');
 118 $population = Bio::PopGen::Utilities->aln_to_population(-alignment => $aln,
 119                                                         -site_model=> 'codon');
 120 isa_ok($population,'Bio::PopGen::Population');
 121
 122 @marker_names = $population->get_marker_names;
 123 @inds = $population->get_Individuals;
 124
 125 is(scalar @marker_names, 378, 'Marker Names');
 126 is(scalar @inds, 7,'Number of Inds');
 127
 128 @ingroup_seqs = ();
 129 @outgroup_seqs1 = ();
 130 @outgroup_seqs2 = ();
 131
 132 for my $ind ( $population->get_Individuals ) {
 133     my $id = $ind->unique_id;
 134     # do we allow ingroup to be a list as well?
 135     my $pushed = 0;
 136     if( $id =~ /sim/ ) {
 137         push @ingroup_seqs, $ind;
 138         $pushed++;
 139     }
 140     if( ! $pushed ) {
 141         if( $id =~ /mel/ ) {
 142             push @outgroup_seqs1, $ind;
 143             $pushed++;
 144         } elsif( $id =~ /yak/ ) {
 145             push @outgroup_seqs2, $ind;
 146             $pushed++;
 147         }
 148     }
 149     #    if( ! $pushed ) {
 150     #   warn("sequence $id was not grouped, ignoring...\n");
 151     #       push @ingroup_seqs, $ind;
 152     #}
 153 }
 154 is(scalar @ingroup_seqs, 5, 'number of ingroup sequences');
 155 is(scalar @outgroup_seqs1, 1, 'number of outgroup1 sequences');
 156 is(scalar @outgroup_seqs2, 1, 'number of outgroup2 sequences');
 157
 158 $polarized = 0;
 159
 160 @counts = $stats->mcdonald_kreitman(-ingroup   => \@ingroup_seqs,
 161                                     -outgroup  => \@outgroup_seqs1,
 162                                     -polarized => $polarized);
 163 is($counts[0], 9, 'NSpoly');
 164 is($counts[1], 1, 'NSfixed');
 165 is($counts[2], 26, 'Spoly');
 166 is($counts[3], 17, 'Sfixed');
 167
 168
 169 SKIP: {
 170         test_skip(-tests => 1,
 171                           -requires_module => 'Text::NSP::Measures::2D::Fisher2::twotailed');
 172         skip "Some problem with Bio::PopGen::Statistics::has_twotailed", 1 unless $Bio::PopGen::Statistics::has_twotailed;
 173
 174     $mk = $stats->mcdonald_kreitman_counts(@counts);
 175     is(sprintf("%.2f",$mk), 0.14, 'McDonald Kreitman');
 176 }
 177
 178 @counts = $stats->mcdonald_kreitman(-ingroup  => \@ingroup_seqs,
 179                                     -outgroup => \@outgroup_seqs2,
 180                                     -polarized=> $polarized);
 181 is($counts[0], 9, 'NSpoly');
 182 is($counts[1], 10, 'NSfixed');
 183 is($counts[2], 26, 'Spoly');
 184 is($counts[3], 42, 'Sfixed');
 185
 186 SKIP: {
 187         test_skip(-tests => 1,
 188                           -requires_module => 'Text::NSP::Measures::2D::Fisher2::twotailed');
 189         skip "Some problem with Bio::PopGen::Statistics::has_twotailed", 1 unless $Bio::PopGen::Statistics::has_twotailed;
 190
 191     $mk = $stats->mcdonald_kreitman_counts(@counts);
 192     is(sprintf("%.2f",$mk), '0.60', 'McDonald Kreitman');
 193 }
 194
 195 @counts = $stats->mcdonald_kreitman(-ingroup  => \@ingroup_seqs,
 196                                     -outgroup => [@outgroup_seqs1,
 197                                                   @outgroup_seqs2],
 198                                     -polarized=> 1);
 199 is($counts[0], 6, 'NSpoly');
 200 is($counts[1], 0, 'NSfixed');
 201 is($counts[2], 17, 'Spoly');
 202 is($counts[3], 1, 'Sfixed');