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