t/Tools/Phylo/PAML.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use lib '.';
   8     use Bio::Root::Test;
   9
  10     test_begin(-tests => 256,
  11                            -requires_module => 'IO::String');
  12
  13         use_ok('Bio::Tools::Phylo::PAML');
  14 }
  15
  16 my $inpaml = Bio::Tools::Phylo::PAML->new(-file => test_input_file('codeml.mlc'));
  17 ok($inpaml);
  18 my $result = $inpaml->next_result;
  19 ok($result);
  20 is($result->model, 'several dN/dS ratios for branches');
  21 like($result->version, qr'3\.12');
  22 my $MLmat = $result->get_MLmatrix;
  23 my $NGmat = $result->get_NGmatrix;
  24
  25 is($NGmat->[0]->[1]->{'omega'}, 0.2507);
  26 is($NGmat->[0]->[1]->{'dN'}, 0.0863);
  27 is($NGmat->[0]->[1]->{'dS'}, 0.3443);
  28
  29 is($NGmat->[2]->[3]->{'omega'}, 0.2178);
  30 is($NGmat->[2]->[3]->{'dN'}, 0.1348);
  31 is($NGmat->[2]->[3]->{'dS'}, 0.6187);
  32
  33 is($MLmat->[0]->[1]->{'omega'}, 0.19479);
  34 is($MLmat->[0]->[1]->{'dN'}, 0.0839);
  35 is($MLmat->[0]->[1]->{'dS'}, 0.4309);
  36 is($MLmat->[0]->[1]->{'lnL'}, -1508.607268);
  37 is($MLmat->[0]->[1]->{'t'}, 0.47825);
  38 is($MLmat->[0]->[1]->{'kappa'}, 2.29137);
  39
  40 is($MLmat->[2]->[3]->{'omega'}, 0.16114);
  41 is($MLmat->[2]->[3]->{'dN'}, 0.1306);
  42 is($MLmat->[2]->[3]->{'dS'}, 0.8105);
  43 is($MLmat->[2]->[3]->{'lnL'},-1666.440696);
  44 is($MLmat->[2]->[3]->{'t'}, 0.85281);
  45 is($MLmat->[2]->[3]->{'kappa'}, 2.21652);
  46
  47 my @codonposfreq = $result->get_codon_pos_basefreq();
  48 is($codonposfreq[0]->{'A'}, 0.23579);
  49 is($codonposfreq[0]->{'T'}, 0.14737);
  50 is($codonposfreq[1]->{'C'}, 0.25123);
  51 is($codonposfreq[2]->{'G'}, 0.32842);
  52
  53 # AAML parsing - Empirical model
  54 $inpaml = Bio::Tools::Phylo::PAML->new(-file => test_input_file('aaml.mlc'));
  55
  56 ok($inpaml);
  57 $result = $inpaml->next_result;
  58 ok($result);
  59 is($result->model, 'Empirical (wag.dat)');
  60 my @trees = $result->get_trees;
  61 is(@trees, 1);
  62 is($trees[0]->score, -1042.768973);
  63
  64 is((scalar grep { $_->is_Leaf } $trees[0]->get_nodes), $result->get_seqs);
  65
  66 my $aadistmat = $result->get_AADistMatrix();
  67 ok($aadistmat);
  68 is($aadistmat->get_entry('Cow', 'Horse'), 0.5462);
  69 is($aadistmat->get_entry('Baboon', 'Langur'), 0.1077);
  70
  71 my %aafreq = %{$result->get_AAFreqs()};
  72 ok(%aafreq);
  73 is($aafreq{'Human'}->{'N'}, 0.0769);
  74 is($aafreq{'Human'}->{'R'}, 0.1077);
  75
  76 my %ratfreqs = %{$result->get_AAFreqs('Rat')};
  77 is($ratfreqs{'R'},0.0923);
  78 is($ratfreqs{'F'},0.0154);
  79 my %avgfreqs = %{$result->get_AAFreqs('Average')};
  80 is($avgfreqs{'Q'},0.0411);
  81
  82 is($result->get_AAFreqs('Average')->{'I'},0.0424);
  83
  84 my $patterns = $result->patterns;
  85 my @pat = @{$patterns->{'-patterns'}};
  86 is(scalar @pat, 98);
  87 is($patterns->{'-ns'}, 6);
  88 is($patterns->{'-ls'}, 130);
  89
  90 is((sort $result->get_stat_names)[0], 'constant_sites');
  91 is($result->get_stat('constant_sites'), 46);
  92 is($result->get_stat('constant_sites_percentage'), 35.38);
  93
  94 # AAML parsing - pairwise model
  95 $inpaml = Bio::Tools::Phylo::PAML->new(-file => test_input_file('aaml_pairwise.mlc'));
  96
  97 ok($inpaml);
  98 $result = $inpaml->next_result;
  99 ok($result);
 100 is($result->model, 'Empirical_F (wag.dat)');
 101 is($result->get_stat('loglikelihood'),-1189.106658);
 102 is($result->get_stat('constant_sites'), 170);
 103 is($result->get_stat('constant_sites_percentage'), 59.65);
 104
 105 is($result->get_AAFreqs('Average')->{'R'},0.0211);
 106 is($result->get_AAFreqs('rabbit')->{'L'},0.1228);
 107
 108 $aadistmat = $result->get_AADistMatrix();
 109 ok($aadistmat);
 110 is($aadistmat->get_entry('rabbit', 'marsupial'), 0.2877);
 111 is($aadistmat->get_entry('human', 'goat-cow'), 0.1439);
 112
 113 $aadistmat = $result->get_AAMLDistMatrix();
 114 ok($aadistmat);
 115 is($aadistmat->get_entry('rabbit', 'marsupial'), 0.3392);
 116 is($aadistmat->get_entry('human', 'goat-cow'), 0.1551);
 117
 118 my @seqs = $result->get_seqs;
 119 is($seqs[0]->display_id, 'human');
 120
 121 # YN00 parsing, pairwise Ka/Ks from Yang & Nielsen 2000
 122 $inpaml = Bio::Tools::Phylo::PAML->new(-file => test_input_file('yn00.mlc'));
 123
 124 ok($inpaml);
 125 $result = $inpaml->next_result;
 126
 127 ok($result);
 128 $MLmat = $result->get_MLmatrix;
 129 $NGmat = $result->get_NGmatrix;
 130
 131 is($NGmat->[0]->[1]->{'omega'}, 0.251);
 132 is($NGmat->[0]->[1]->{'dN'}, 0.0863);
 133 is($NGmat->[0]->[1]->{'dS'}, 0.3443);
 134 is($NGmat->[2]->[3]->{'omega'}, 0.218);
 135 is($NGmat->[2]->[3]->{'dN'}, 0.1348);
 136 is($NGmat->[2]->[3]->{'dS'}, 0.6187);
 137
 138 is($MLmat->[0]->[1]->{'omega'}, 0.1625);
 139 is($MLmat->[0]->[1]->{'dN'}, 0.0818);
 140 is($MLmat->[0]->[1]->{'dS'}, 0.5031);
 141 is($MLmat->[2]->[3]->{'omega'}, 0.1262);
 142 is($MLmat->[2]->[3]->{'dN'}, 0.1298);
 143 is($MLmat->[2]->[3]->{'dN_SE'}, 0.0149);
 144 is($MLmat->[2]->[3]->{'dS'}, 1.0286);
 145 is($MLmat->[2]->[3]->{'dS_SE'}, 0.2614);
 146
 147 # codeml NSSites parsing
 148
 149 $inpaml = Bio::Tools::Phylo::PAML->new
 150     (-file => test_input_file('codeml_nssites.mlc'));
 151
 152 ok($inpaml);
 153 $result = $inpaml->next_result;
 154
 155 ok($result);
 156 is($result->model, 'One dN/dS ratio dGamma (ncatG=11)');
 157 is($result->version, 'paml 3.13, August 2002');
 158 $NGmat = $result->get_NGmatrix;
 159 ok($NGmat);
 160
 161 is($NGmat->[0]->[1]->{'omega'}, 0.2782);
 162 is($NGmat->[0]->[1]->{'dN'}, 0.0133);
 163 is($NGmat->[0]->[1]->{'dS'}, 0.0478);
 164 is($NGmat->[1]->[2]->{'omega'}, 1.1055);
 165 is($NGmat->[1]->[2]->{'dN'}, 0.0742);
 166 is($NGmat->[1]->[2]->{'dS'}, 0.0671);
 167           # this is
 168           #   model num  description
 169           #   kappa   log-likelihood tree length time used
 170           #   shape   alpha/gamma r          f
 171 my @tstr = ([qw(0 one-ratio 0
 172                 4.54006 -906.017440    0.55764
 173                 )],
 174             [qw(1 neutral 2
 175                 4.29790 -902.503869    0.56529
 176                 )],
 177             [qw(2 selection 3
 178                 5.12250 -900.076500    0.6032
 179                 )],
 180              );
 181 my $iter = 0;
 182 my $lastmodel;
 183 foreach my $model ( $result->get_NSSite_results ) {
 184     my $i = 0;
 185     my $r = shift @tstr;
 186     is($model->model_num, $r->[$i++]);
 187     like($model->model_description, qr/$r->[$i++]/);
 188     is($model->num_site_classes,$r->[$i++]);
 189     my $tree = $model->next_tree;
 190     is($model->kappa, $r->[$i++]);
 191     is($model->likelihood,$r->[$i]);
 192     is($tree->score, $r->[$i++]);
 193     is($tree->total_branch_length, $r->[$i++]);
 194     if( $iter == 0 ) {
 195         my $params = $model->shape_params;
 196         is($params->{'shape'}, 'alpha');
 197         is($params->{'gamma'},   '0.50000');
 198         is($params->{'r'}->[0], '1.00000');
 199         is($params->{'f'}->[0], '1.00000');
 200     } elsif( $iter == 2 ) {
 201         my $class = $model->dnds_site_classes;
 202         is($class->{'p'}->[0], '0.38160');
 203         is($class->{'w'}->[1], '1.00000');
 204     }
 205     $iter++;
 206     $lastmodel = $model;
 207 }
 208
 209 my ($firstsite) = $lastmodel->get_pos_selected_sites;
 210 is($firstsite->[0], 15);
 211 is($firstsite->[1], 'L');
 212 is($firstsite->[2], 0.6588);
 213
 214 # codeml NSSites parsing
 215 # for M0 model
 216
 217 my $codeml_m0 = Bio::Tools::Phylo::PAML->new
 218     (-file => test_input_file('M0.mlc'));
 219 ok($codeml_m0);
 220 my $result_m0 = $codeml_m0->next_result;
 221 my ($nssite_m0,$nssite_m1) = $result_m0->get_NSSite_results;
 222 is($nssite_m0->num_site_classes,1);
 223 my $class_m0 = $nssite_m0->dnds_site_classes;
 224 is($class_m0->{q/p/}->[0],q/1.00000/);
 225 is($class_m0->{q/w/}->[0],0.09213);
 226
 227 is($nssite_m0->model_num, "0");
 228 @trees= $nssite_m0->get_trees;
 229 is (@trees , 1 );
 230 # model 0
 231 is($trees[0]->score, -30.819156);
 232 is($nssite_m1->model_num, "1");
 233 @trees= $nssite_m1->get_trees;
 234 is($trees[0]->score, -30.819157);
 235
 236 # test BASEML
 237 # pairwise first
 238
 239 my $baseml_p = Bio::Tools::Phylo::PAML->new
 240     (-file => test_input_file('baseml.pairwise'));
 241 ok($baseml_p);
 242 my $baseml = $baseml_p->next_result;
 243 my @b_seqs =  $baseml->get_seqs;
 244 is($b_seqs[0]->seq, 'GTAGAGTACTTT');
 245 is($b_seqs[1]->seq, 'GTAAGAGACGAT');
 246
 247 my @otus = map { $_->display_id } @b_seqs;
 248 is(scalar @otus, 3);
 249 my $ntfreq = $baseml->get_NTFreqs;
 250 ok($ntfreq);
 251 is($ntfreq->{$otus[0]}->{'A'}, '0.3333');
 252 is($ntfreq->{$otus[1]}->{'G'}, '0.2105');
 253 my $kappaM = $baseml->get_KappaMatrix;
 254 ok($kappaM);
 255 is($kappaM->get_entry($otus[1],$otus[0]), '0.3240');
 256 is($kappaM->get_entry($otus[0],$otus[1]),
 257    $kappaM->get_entry($otus[1],$otus[0]));
 258 is($kappaM->get_entry($otus[1],$otus[2]), '0.1343');
 259 my $alphaM = $baseml->get_AlphaMatrix;
 260 ok($alphaM);
 261 is($alphaM->get_entry($otus[1],$otus[0]), '9.3595');
 262 is($alphaM->get_entry($otus[0],$otus[1]),
 263    $alphaM->get_entry($otus[1],$otus[0]));
 264 is($alphaM->get_entry($otus[1],$otus[2]), '1.1101');
 265 is($alphaM->get_entry($otus[0],$otus[2]), '33.1197');
 266
 267 # codeml NSSites parsing
 268 # for only 1 model
 269
 270 my $codeml_single = Bio::Tools::Phylo::PAML->new
 271     (-file => test_input_file('singleNSsite.mlc'));
 272 ok($codeml_single);
 273 my $result_single = $codeml_single->next_result;
 274 my ($nssite_single) = $result_single->get_NSSite_results;
 275 is($nssite_single->num_site_classes,q/3/);
 276 is($nssite_single->kappa, q/5.28487/);
 277 is($nssite_single->likelihood,q/-30.819156/);
 278
 279 is($baseml->get_stat('loglikelihood'),-110.532715);
 280 is($baseml->get_stat('constant_sites'),46);
 281 is($baseml->get_stat('constant_sites_percentage'),'80.70');
 282 is($baseml->model,'HKY85 dGamma (ncatG=5)');
 283
 284 # user trees
 285 $baseml_p = Bio::Tools::Phylo::PAML->new
 286     (-file => test_input_file('baseml.usertree'));
 287 $baseml = $baseml_p->next_result;
 288
 289 @trees = $baseml->get_trees;
 290 is(@trees, 1);
 291 is($trees[0]->score, -129.328757);
 292
 293 # codeml NSSites parsing
 294 # for branch site model/clade model
 295
 296 my $codeml_bs = Bio::Tools::Phylo::PAML->new
 297     (-file => test_input_file('branchSite.mlc'));
 298 ok($codeml_bs);
 299 my $result_bs = $codeml_bs->next_result;
 300 my ($nssite_bs) = $result_bs->get_NSSite_results;
 301 is($nssite_bs->num_site_classes,q/4/);
 302 my $class_bs = $nssite_bs->dnds_site_classes;
 303 is($class_bs->{q/p/}->[1],q/0.65968/);
 304 is($class_bs->{q/w/}->[1]->{q/background/},q/0.00000/);
 305 is($class_bs->{q/w/}->[2]->{q/foreground/},q/999.00000/);
 306
 307 # Let's parse the RST file
 308
 309 my $paml = Bio::Tools::Phylo::PAML->new
 310     (-file => test_input_file('codeml_lysozyme', 'mlc'),
 311      -dir  => test_input_file('codeml_lysozyme'));
 312
 313 $result = $paml->next_result;
 314
 315 my ($rst) = grep {$_->id eq 'node#8'} $result->get_rst_seqs;
 316 ok($rst);
 317 is($rst->seq, join('',qw(
 318 AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACTCTGAAAAGATTGGGACTGGATGGCTAC
 319 AGGGGAATCAGCCTAGCAAACTGGATGTGTTTGGCCAAATGGGAGAGTGATTATAACACA
 320 CGAGCTACAAACTACAATCCTGGAGACCAAAGCACTGATTATGGGATATTTCAGATCAAT
 321 AGCCACTACTGGTGTAATAATGGCAAAACCCCAGGAGCAGTTAATGCCTGTCATATATCC
 322 TGCAATGCTTTGCTGCAAGATAACATCGCTGATGCTGTAGCTTGTGCAAAGAGGGTTGTC
 323 CGTGATCCACAAGGCATTAGAGCATGGGTGGCATGGAGAAATCATTGTCAAAACAGAGAT
 324 GTCAGTCAGTATGTTCAAGGTTGTGGAGTG)),
 325    'node#8 reconstructed seq');
 326
 327 my ($first_tree) = $result->get_rst_trees;
 328 my ($node) = $first_tree->find_node(-id => '5_Mmu_rhesus');
 329 my @changes = $node->get_tag_values('changes');
 330 my ($site) = grep { $_->{'site'} == 94 } @changes;
 331 is($site->{'anc_aa'}, 'A','ancestral AA');
 332 is($site->{'anc_prob'}, '0.947','ancestral AA');
 333 is($site->{'derived_aa'}, 'T','derived AA');
 334
 335 ($node) = $first_tree->find_node(-id => '12');
 336 @changes = $node->get_tag_values('changes');
 337 ($site) = grep { $_->{'site'} == 88 } @changes;
 338 is($site->{'anc_aa'}, 'N','ancestral AA');
 339 is($site->{'anc_prob'}, '0.993','ancestral AA');
 340 is($site->{'derived_aa'}, 'D','derived AA');
 341 is($site->{'derived_prob'}, '0.998','derived AA');
 342
 343 my $persite = $result->get_rst_persite;
 344 # minus 1 because we have shifted so that array index matches site number
 345 # there are 130 sites in this seq file
 346 is(scalar @$persite -1, $result->patterns->{'-ls'});
 347 # let's score site 1
 348 $site = $persite->[2];
 349 # so site 2, node 2 (extant)
 350 is($site->[2]->{'codon'}, 'GTC');
 351 is($site->[2]->{'aa'}, 'V');
 352 # site 2, node 3
 353 is($site->[3]->{'codon'}, 'ATC');
 354 is($site->[3]->{'aa'}, 'I');
 355
 356 # ancestral node 9
 357 is($site->[9]->{'codon'}, 'GTC');
 358 is($site->[9]->{'aa'},    'V');
 359 is($site->[9]->{'prob'},  '1.000');
 360 is($site->[9]->{'Yang95_aa'},'V');
 361 is($site->[9]->{'Yang95_aa_prob'},'1.000');
 362
 363 # ancestral node 10
 364 is($site->[10]->{'codon'}, 'ATC');
 365 is($site->[10]->{'aa'},    'I');
 366 is($site->[10]->{'prob'},  '0.992');
 367 is($site->[10]->{'Yang95_aa'},'I');
 368 is($site->[10]->{'Yang95_aa_prob'},'0.992');
 369
 370
 371 ## PAML 3.15
 372 my $paml315 = Bio::Tools::Phylo::PAML->new(-file => test_input_file('codeml315.mlc'));
 373 $result = $paml315->next_result;
 374
 375 is($result->model, 'One dN/dS ratio');
 376 like($result->version, qr'3\.15');
 377 $MLmat = $result->get_MLmatrix;
 378 $NGmat = $result->get_NGmatrix;
 379
 380 is($NGmat->[0]->[1]->{'omega'}, 0.2264);
 381 is($NGmat->[0]->[1]->{'dN'}, 0.0186);
 382 is($NGmat->[0]->[1]->{'dS'}, 0.0821);
 383
 384 is($MLmat->[0]->[1]->{'omega'}, 0.32693);
 385 is($MLmat->[0]->[1]->{'dN'}, '0.0210');
 386 is($MLmat->[0]->[1]->{'dS'}, 0.0644);
 387
 388 ## PAML 4
 389 my $codeml4 = Bio::Tools::Phylo::PAML->new(-file => test_input_file('codeml4.mlc'));
 390 $result = $codeml4->next_result;
 391
 392 is($result->model, 'One dN/dS ratio');
 393 like($result->version, qr'4');
 394 $MLmat = $result->get_MLmatrix;
 395 $NGmat = $result->get_NGmatrix;
 396
 397 is($NGmat->[0]->[1]->{'omega'}, 0.2507);
 398 is($NGmat->[0]->[1]->{'dN'}, 0.0863);
 399 is($NGmat->[0]->[1]->{'dS'}, 0.3443);
 400
 401 is($MLmat->[0]->[1]->{'omega'}, 0.29075);
 402 is($MLmat->[0]->[1]->{'dN'}, '0.0874');
 403 is($MLmat->[0]->[1]->{'dS'}, 0.3006);
 404 is($MLmat->[0]->[1]->{'lnL'}, -1596.739984);
 405
 406 ## PAML 4.3a
 407 # codeml pairwise ML comparison (runmode=-2)
 408 my $codeml43 = Bio::Tools::Phylo::PAML->new(-file => test_input_file('codeml43.mlc'));
 409 $result = $codeml43->next_result;
 410
 411 is($result->model, 'One dN/dS ratio for branches');
 412 like($result->version, qr'4\.3', 'codeml 4.3 runmode=-2');
 413 $MLmat = $result->get_MLmatrix;
 414 $NGmat = $result->get_NGmatrix;
 415
 416 is($NGmat->[0]->[2]->{'omega'}, 0.2627);
 417 is($NGmat->[0]->[2]->{'dN'}, 0.0867);
 418 is($NGmat->[0]->[2]->{'dS'}, 0.3301);
 419
 420 is($MLmat->[0]->[2]->{'omega'}, 0.19819);
 421 is($MLmat->[0]->[2]->{'dN'}, '0.0842');
 422 is($MLmat->[0]->[2]->{'dS'}, 0.4247);
 423 is($MLmat->[0]->[2]->{'lnL'}, -1512.583367);
 424
 425 ## PAML 4.3a
 426 # codeml NSSites parsing (two NSSites models, 1 and 2)
 427 {
 428     my $codeml43_nssites = Bio::Tools::Phylo::PAML->new
 429         (-file => test_input_file('codeml43_nssites.mlc'));
 430     ok($codeml43_nssites);
 431
 432     my $result = $codeml43_nssites->next_result;
 433     ok($result);
 434
 435     is($result->model, 'One dN/dS ratio for branches');
 436     like($result->version, qr'4\.3', 'codeml 4.3 two NSSites models');
 437     my $NGmat = $result->get_NGmatrix;
 438     ok($NGmat);
 439
 440     is($NGmat->[0]->[1]->{'omega'}, 0.2507);
 441     is($NGmat->[0]->[1]->{'dN'}, 0.0863);
 442     is($NGmat->[0]->[1]->{'dS'}, 0.3443);
 443     is($NGmat->[1]->[2]->{'omega'}, 0.2943);
 444     is($NGmat->[1]->[2]->{'dN'}, 0.1054);
 445     is($NGmat->[1]->[2]->{'dS'}, 0.3581);
 446
 447     # these are
 448     # "model num" description "number of site classes" kappa log-likelihood "tree length" "time used"
 449     my @tstr = ([qw(1 NearlyNeutral     2 2.06684 -2970.527521 2.898 0:08)],
 450                 [qw(2 PositiveSelection 3 2.18136 -2965.809712 3.589 0:26)],);
 451     my $iter = 0;
 452     my $lastmodel;
 453     foreach my $model ( $result->get_NSSite_results ) {
 454         my $i = 0;
 455         my $r = shift @tstr;
 456         is($model->model_num, $r->[$i++]);
 457         like($model->model_description, qr/$r->[$i++]/);
 458         is($model->num_site_classes,$r->[$i++]);
 459         my $tree = $model->next_tree;
 460         is($model->kappa, $r->[$i++]);
 461         is($model->likelihood,$r->[$i]);
 462         is($tree->score, $r->[$i++]);
 463         like($tree->total_branch_length, qr/$r->[$i++]/);
 464         if( $iter == 1 ) {
 465             my $class = $model->dnds_site_classes;
 466             is($class->{'p'}->[0], '0.83347');
 467             is($class->{'w'}->[1], '1.00000');
 468         }
 469         $iter++;
 470         $lastmodel = $model;
 471     }
 472
 473     my @sites = $lastmodel->get_NEB_pos_selected_sites;
 474     my $firstsite = $sites[0];
 475     my $lastsite  = $sites[-1];
 476     is($firstsite->[0], 35, 'NEB positively selected sites');
 477     is($firstsite->[1], 'S');
 478     is($firstsite->[2], 0.643);
 479     is($firstsite->[4], '4.400');
 480     is( $lastsite->[0], 264);
 481     is( $lastsite->[1], 'P');
 482     is( $lastsite->[2], 0.971);
 483     is( $lastsite->[3], '*');
 484     is( $lastsite->[4], 6.134);
 485 }
 486
 487 # bug #3040
 488 {
 489     my $parser = Bio::Tools::Phylo::PAML->new
 490         (-file => test_input_file('codeml_nan.mlc'));
 491     ok($parser);
 492
 493     my $result = $parser->next_result;
 494     ok($result);
 495
 496     my $MLmatrix = $result->get_MLmatrix();
 497     ok($MLmatrix);
 498
 499     is($MLmatrix->[1]->[2]->{'dS'}, 'nan', 'bug 3040');
 500 }
 501
 502 # bugs 3365, 3366
 503 {
 504     my $parser = Bio::Tools::Phylo::PAML->new
 505         (-file => test_input_file('codeml45.mlc'));
 506
 507     my $result = $parser->next_result;
 508
 509         my @otus = $result->get_seqs();
 510         is(scalar @otus, 9, 'bug 3365');
 511
 512         my $MLmatrix = $result->get_MLmatrix();
 513         is($MLmatrix->[1]->[2]->{dN},0.0103,'bug 3366');
 514 }
 515
 516 # bug 3367
 517 {
 518     my $parser = Bio::Tools::Phylo::PAML->new
 519         (-file => test_input_file('yn00_45.mlc'));
 520
 521     my $result = $parser->next_result;
 522
 523         my @otus = $result->get_seqs();
 524         is(scalar @otus, 9, 'bug 3367');
 525 }
 526
 527 # bug 3332
 528 {
 529     my $parser = Bio::Tools::Phylo::PAML->new
 530         (-file => test_input_file('codeml45b.mlc'));
 531
 532         my $result = $parser->next_result;
 533         my $omega2 = $result->get_NGmatrix()->[0]->[1]->{'omega'};
 534         is($result->get_NGmatrix()->[0]->[1]->{'omega'}, '-1.0300', 'bug 3332');
 535 }
 536
 537 # bug 3331
 538 {
 539     my $parser = Bio::Tools::Phylo::PAML->new
 540         (-file => test_input_file('bug3331.mlc'));
 541         my $result = $parser->next_result;
 542         my $MLmatrix = $result->get_MLmatrix();
 543         my $kappa = $MLmatrix->[0]->[1]->{'kappa'};
 544         is ($kappa, '2.000', 'bug 3331');
 545 }