t/PAML.t

   1 # This is -*-Perl-*- code
   2 ## Bioperl Test Harness Script for Modules
   3 ##
   4 # $Id$
   5
   6 # Before `make install' is performed this script should be runnable with
   7 # `make test'. After `make install' it should work as `perl test.t'
   8
   9 use strict;
  10 BEGIN {
  11     use Bio::Root::Test;
  12     test_begin(-tests => 28,
  13                            -requires_module => 'IO::String');
  14         use_ok('Bio::Tools::Phylo::PAML'); # PAML parser
  15         use_ok('Bio::Tools::Run::Phylo::PAML::Codeml');
  16         use_ok('Bio::Tools::Run::Phylo::PAML::Yn00');
  17         use_ok('Bio::AlignIO');
  18 }
  19
  20 my $testnum;
  21 my $verbose = 0;
  22
  23 ## End of black magic.
  24 ##
  25 ## Insert additional test code below but remember to change
  26 ## the print "1..x\n" in the BEGIN block to reflect the
  27 ## total number of tests that will be run.
  28
  29 my $inpaml = Bio::Tools::Phylo::PAML->new(-file =>
  30                                          test_input_file('codeml.mlc'));
  31
  32 ok($inpaml);
  33
  34 my $codeml = Bio::Tools::Run::Phylo::PAML::Codeml->new
  35     (-params => {'runmode' => -2,
  36                  'seqtype' => 1,
  37                  'model'   => 0,
  38                  'alpha'   => '0',
  39                  'omega'   => 0.4,
  40                  'kappa'    => 2,
  41                  'CodonFreq'=> 2,
  42                  'NSsites'   => 0,
  43                  'model'    => 0,
  44              },
  45      -verbose => $verbose);
  46
  47 SKIP: {
  48         test_skip(-requires_executable => $codeml,
  49                   -tests => 23);
  50
  51         my $in = Bio::AlignIO->new(-format => 'phylip',
  52                                   -file   => test_input_file('gf-s85.phylip'));
  53         my $aln = $in->next_aln;
  54         $codeml->alignment($aln);
  55         my ($rc,$results) = $codeml->run();
  56
  57         is($rc,1);
  58
  59         if( ! defined $results ) {
  60                 skip('No results', 22);
  61         }
  62
  63         my $result = $results->next_result;
  64         if( ! defined $result ) {
  65                 skip('No result', 22);
  66         }
  67
  68         my $MLmatrix = $result->get_MLmatrix;
  69
  70         my ($vnum) = ($result->version =~ /(\d+(\.\d+)?)/);
  71         # PAML 2.12 results
  72         if( $vnum == 3.12 ) {
  73                 is($MLmatrix->[0]->[1]->{'dN'}, 0.0693);
  74                 is($MLmatrix->[0]->[1]->{'dS'},1.1459);
  75                 is($MLmatrix->[0]->[1]->{'omega'}, 0.0605);
  76                 is($MLmatrix->[0]->[1]->{'S'}, 273.5);
  77                 is($MLmatrix->[0]->[1]->{'N'}, 728.5);
  78                 is($MLmatrix->[0]->[1]->{'t'}, 1.0895);
  79
  80                 skip($MLmatrix->[0]->[1]->{'lnL'}, "I don't know what this should be, if you run this part, email the list so we can update the value");
  81
  82         } elsif( $vnum >= 3.13  && $vnum < 4) {
  83         # PAML 2.13 results
  84                 is($MLmatrix->[0]->[1]->{'dN'}, 0.0713);
  85                 is($MLmatrix->[0]->[1]->{'dS'},1.2462);
  86                 is(sprintf("%.4f",$MLmatrix->[0]->[1]->{'omega'}), 0.0572);
  87                 is($MLmatrix->[0]->[1]->{'S'}, 278.8);
  88                 is($MLmatrix->[0]->[1]->{'N'}, 723.2);
  89                 is(sprintf("%.4f",$MLmatrix->[0]->[1]->{'t'}), 1.1946);
  90                 skip($MLmatrix->[0]->[1]->{'lnL'}, "I don't know what this should be, if you run this part, email the list so we can update the value");
  91
  92         } elsif( $vnum == 4 ) {
  93                 is($MLmatrix->[0]->[1]->{'dN'}, 0.0713);
  94                 is($MLmatrix->[0]->[1]->{'dS'},1.2462);
  95                 is(sprintf("%.4f",$MLmatrix->[0]->[1]->{'omega'}), 0.0572);
  96                 is($MLmatrix->[0]->[1]->{'S'}, 278.8);
  97                 is($MLmatrix->[0]->[1]->{'N'}, 723.2);
  98                 is(sprintf("%.4f",$MLmatrix->[0]->[1]->{'t'}), 1.1946);
  99                 is($MLmatrix->[0]->[1]->{'lnL'}, -1929.935243);
 100
 101         } else {
 102                 for( 1..7) {
 103                 skip("Can't test the result output, don't know about PAML version ".$result->version,1);
 104                 }
 105         }
 106
 107         unlike($codeml->error_string, qr/Error/); # we don't expect any errors;
 108
 109         my $yn00 = Bio::Tools::Run::Phylo::PAML::Yn00->new();
 110         $yn00->alignment($aln);
 111         ($rc,$results) = $yn00->run();
 112         is($rc,1);
 113         if( ! defined $results ) {
 114                 exit(0);
 115         }
 116         $result = $results->next_result;
 117         $MLmatrix = $result->get_MLmatrix;
 118
 119         is($MLmatrix->[0]->[1]->{'dN'}, 0.0846);
 120         is($MLmatrix->[0]->[1]->{'dS'}, 1.0926);
 121         is($MLmatrix->[0]->[1]->{'omega'}, 0.0774);
 122         is($MLmatrix->[0]->[1]->{'S'}, 278.4);
 123         is($MLmatrix->[0]->[1]->{'N'}, 723.6);
 124         is($MLmatrix->[0]->[1]->{'t'}, 1.0941);
 125
 126         unlike($yn00->error_string, qr/Error/); # we don't expect any errors;
 127
 128         $codeml = Bio::Tools::Run::Phylo::PAML::Codeml->new
 129                 (-params => { 'alpha' => 1.53 },
 130                  -verbose => $verbose);
 131
 132         ok($codeml);
 133
 134
 135         # AAML
 136         my $cysaln = Bio::AlignIO->new(-format => 'msf',
 137                                            -file => test_input_file('cysprot.msf'))->next_aln;
 138
 139         my $cystre = Bio::TreeIO->new(-format => 'newick',
 140                                           -file  => test_input_file('cysprot.raxml.tre'))->next_tree;
 141         ok($cysaln);
 142         ok($cystre);
 143
 144         $codeml = Bio::Tools::Run::Phylo::PAML::Codeml->new
 145                 (
 146                  -verbose => 0,
 147                  -tree   => $cystre,
 148                  -params => { 'runmode' => 0, # provide a usertree
 149                           'seqtype' => 2, # AMINO ACIDS,
 150                           'model'   => 0, # one dN/dS rate
 151                           'NSsites' => 0, # one -- swap this with 1, 2, 3 etc
 152                           'clock'   => 0, # 0 = no clock
 153                           'getSE'   => 1, # get Standard Error
 154                           'fix_blength' => 0, # use initial BLengths
 155                           'ncatG' => 1, #increase approrpriately for NSsites,
 156                  },
 157                  -alignment => $cysaln,
 158                  -save_tempfiles => 1,
 159                 );
 160     test_skip(-requires_executable => $codeml,
 161               -tests => 14);
 162         ok($codeml);
 163
 164         ($rc,$results) = $codeml->run();
 165         is($rc,1);
 166
 167
 168         unless( defined $results ) {
 169                 warn($codeml->error_string, "\n");
 170                 skip('No results',1);
 171         }
 172
 173         $result = $results->next_result;
 174         unless( defined $result ) {
 175                 skip('No result',1);
 176         }
 177
 178         ($vnum) = ($result->version =~ /(\d+(\.\d+)?)/);
 179         for my $tree ( $result->get_trees ) {
 180                 my $node = $tree->find_node(-id => 'CATL_HUMAN');
 181                 if( $vnum == 4 ) {
 182                 is($node->branch_length, '0.216223');
 183                 } else {
 184                 is($node->branch_length, '0.216223');
 185                 }
 186         }
 187 }