Bio/Tools/Phylo/PAML/ModelResult.pm

   1 #
   2 # BioPerl module for Bio::Tools::Phylo::PAML::ModelResult
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Jason Stajich <jason@open-bio.org>
   7 #
   8 # Copyright Jason Stajich
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::Tools::Phylo::PAML::ModelResult - A container for NSSite Model Result from PAML
  17
  18 =head1 SYNOPSIS
  19
  20   # get a ModelResult from a PAML::Result object
  21   use Bio::Tools::Phylo::PAML;
  22   my $paml = Bio::Tools::Phylo::PAML->new(-file => 'mlc');
  23   my $result = $paml->next_result;
  24   foreach my $model ( $result->get_NSSite_results ) {
  25     print $model->model_num, " ", $model->model_description, "\n";
  26     print $model->kappa, "\n";
  27     print $model->run_time, "\n";
  28 # if you are using PAML < 3.15 then only one place for POS sites
  29    for my $sites ( $model->get_pos_selected_sites ) {
  30     print join("\t",@$sites),"\n";
  31    }
  32 # otherwise query NEB and BEB slots
  33    for my $sites ( $model->get_NEB_pos_selected_sites ) {
  34      print join("\t",@$sites),"\n";
  35    }
  36
  37    for my $sites ( $model->get_BEB_pos_selected_sites ) {
  38     print join("\t",@$sites),"\n";
  39    }
  40
  41   }
  42
  43 =head1 DESCRIPTION
  44
  45 Describe the object here
  46
  47 =head1 FEEDBACK
  48
  49 =head2 Mailing Lists
  50
  51 User feedback is an integral part of the evolution of this and other
  52 Bioperl modules. Send your comments and suggestions preferably to
  53 the Bioperl mailing list.  Your participation is much appreciated.
  54
  55   bioperl-l@bioperl.org                  - General discussion
  56   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  57
  58 =head2 Support
  59
  60 Please direct usage questions or support issues to the mailing list:
  61
  62 I<bioperl-l@bioperl.org>
  63
  64 rather than to the module maintainer directly. Many experienced and
  65 reponsive experts will be able look at the problem and quickly
  66 address it. Please include a thorough description of the problem
  67 with code and data examples if at all possible.
  68
  69 =head2 Reporting Bugs
  70
  71 Report bugs to the Bioperl bug tracking system to help us keep track
  72 of the bugs and their resolution. Bug reports can be submitted via
  73 email or the web:
  74
  75   https://github.com/bioperl/bioperl-live/issues
  76
  77 =head1 AUTHOR - Jason Stajich
  78
  79 Email jason@open-bio.org
  80
  81 =head1 APPENDIX
  82
  83 The rest of the documentation details each of the object methods.
  84 Internal methods are usually preceded with a _
  85
  86 =cut
  87
  88
  89 # Let the code begin...
  90
  91
  92 package Bio::Tools::Phylo::PAML::ModelResult;
  93 use strict;
  94
  95 # Object preamble - inherits from Bio::Root::Root
  96
  97
  98 use base qw(Bio::Root::Root);
  99
 100 =head2 new
 101
 102  Title   : new
 103  Usage   : my $obj = Bio::Tools::Phylo::PAML::ModelResult->new();
 104  Function: Builds  a new Bio::Tools::Phylo::PAML::ModelResult object
 105  Returns : an instance of Bio::Tools::Phylo::PAML::ModelResult
 106  Args    : -model_num           => model number
 107            -model_description   => model description
 108            -kappa               => value of kappa
 109            -time_used           => amount of time
 110            -pos_sites           => arrayref of sites under positive selection
 111            -neb_sites           => arrayref of sites under positive selection (by NEB analysis)
 112            -beb_sites           => arrayref of sites under positive selection (by BEB analysis)
 113            -trees               => arrayref of tree(s) data for this model
 114            -shape_params        => hashref of parameters
 115                                    ('shape' => 'alpha',
 116                                     'gamma' => $g,
 117                                     'r' => $r,
 118                                     'f' => $f
 119                                     )
 120                                     OR
 121                                     ( 'shape' => 'beta',
 122                                       'p' => $p,
 123                                       'q' => $q
 124                                      )
 125            -likelihood          => likelihood
 126            -num_site_classes    => number of site classes
 127            -dnds_site_classes   => hashref with two keys, 'p' and 'w'
 128                                    which each point to an array, each
 129                                    slot is for a different site class.
 130                                    'w' is for dN/dS and 'p' is probability
 131
 132 =cut
 133
 134 sub new {
 135   my($class,@args) = @_;
 136
 137   my $self = $class->SUPER::new(@args);
 138   my ($modelnum,$modeldesc,$kappa,
 139       $timeused,$trees,
 140       $pos_sites,$neb_sites,$beb_sites,
 141       $num_site_classes, $shape_params,
 142       $dnds_classes,
 143       $likelihood) =          $self->_rearrange([qw(MODEL_NUM
 144                                                       MODEL_DESCRIPTION
 145                                                       KAPPA
 146                                                       TIME_USED
 147                                                       TREES
 148                                                       POS_SITES
 149                                                       NEB_SITES BEB_SITES
 150                                                       NUM_SITE_CLASSES
 151                                                       SHAPE_PARAMS
 152                                                       DNDS_SITE_CLASSES
 153                                                       LIKELIHOOD)],
 154                                                  @args);
 155   if( $trees ) {
 156       if(ref($trees) !~ /ARRAY/i ) {
 157           $self->warn("Must provide a valid array reference to initialize trees");
 158       } else {
 159           foreach my $t ( @$trees ) {
 160               $self->add_tree($t);
 161           }
 162       }
 163   }
 164   $self->{'_treeiterator'} = 0;
 165   if( $pos_sites ) {
 166       if(ref($pos_sites) !~ /ARRAY/i ) {
 167           $self->warn("Must provide a valid array reference to initialize pos_sites");
 168       } else {
 169           foreach my $s ( @$pos_sites ) {
 170               if( ref($s) !~ /ARRAY/i ) {
 171                   $self->warn("Need an array reference for each entry in the pos_sites object");
 172                   next;
 173               }
 174               $self->add_pos_selected_site(@$s);
 175           }
 176       }
 177   }
 178   if( $beb_sites ) {
 179     if(ref($beb_sites) !~ /ARRAY/i ) {
 180           $self->warn("Must provide a valid array reference to initialize beb_sites");
 181       } else {
 182           foreach my $s ( @$beb_sites ) {
 183               if( ref($s) !~ /ARRAY/i ) {
 184                   $self->warn("need an array ref for each entry in the beb_sites object");
 185                   next;
 186               }
 187               $self->add_BEB_pos_selected_site(@$s);
 188           }
 189       }
 190   }
 191   if( $neb_sites ) {
 192     if(ref($neb_sites) !~ /ARRAY/i ) {
 193           $self->warn("Must provide a valid array reference to initialize neb_sites");
 194       } else {
 195           foreach my $s ( @$neb_sites ) {
 196               if( ref($s) !~ /ARRAY/i ) {
 197                   $self->warn("need an array ref for each entry in the neb_sites object");
 198                   next;
 199               }
 200               $self->add_NEB_pos_selected_site(@$s);
 201           }
 202       }
 203   }
 204
 205   defined $modelnum  && $self->model_num($modelnum);
 206   defined $modeldesc && $self->model_description($modeldesc);
 207   defined $kappa     && $self->kappa($kappa);
 208   defined $timeused  && $self->time_used($timeused);
 209   defined $likelihood  && $self->likelihood($likelihood);
 210
 211   $self->num_site_classes($num_site_classes || 0);
 212   if( defined $dnds_classes ) {
 213       if( ref($dnds_classes) !~ /HASH/i ||
 214           ! defined $dnds_classes->{'p'} ||
 215           ! defined $dnds_classes->{'w'} ) {
 216           $self->warn("-dnds_site_classes expects a hashref with keys p and w");
 217       } else {
 218           $self->dnds_site_classes($dnds_classes);
 219       }
 220   }
 221   if( defined $shape_params ) {
 222       if( ref($shape_params) !~ /HASH/i ) {
 223           $self->warn("-shape_params expects a hashref not $shape_params\n");
 224       } else {
 225           $self->shape_params($shape_params);
 226       }
 227   }
 228   return $self;
 229 }
 230
 231
 232 =head2 model_num
 233
 234  Title   : model_num
 235  Usage   : $obj->model_num($newval)
 236  Function: Get/Set the Model number (0,1,2,3...)
 237  Returns : value of model_num (a scalar)
 238  Args    : on set, new value (a scalar or undef, optional)
 239
 240
 241 =cut
 242
 243 sub model_num {
 244     my $self = shift;
 245     return $self->{'_num'} = shift if @_;
 246     return $self->{'_num'};
 247 }
 248
 249 =head2 model_description
 250
 251  Title   : model_description
 252  Usage   : $obj->model_description($newval)
 253  Function: Get/Set the model description
 254            This is something like 'one-ratio', 'neutral', 'selection'
 255  Returns : value of description (a scalar)
 256  Args    : on set, new value (a scalar or undef, optional)
 257
 258
 259 =cut
 260
 261 sub model_description{
 262     my $self = shift;
 263     return $self->{'_model_description'} = shift if @_;
 264     return $self->{'_model_description'};
 265 }
 266
 267 =head2 time_used
 268
 269  Title   : time_used
 270  Usage   : $obj->time_used($newval)
 271  Function: Get/Set the time it took to run this analysis
 272  Returns : value of time_used (a scalar)
 273  Args    : on set, new value (a scalar or undef, optional)
 274
 275
 276 =cut
 277
 278 sub time_used{
 279     my $self = shift;
 280     return $self->{'_time_used'} = shift if @_;
 281     return $self->{'_time_used'};
 282 }
 283
 284 =head2 kappa
 285
 286  Title   : kappa
 287  Usage   : $obj->kappa($newval)
 288  Function: Get/Set kappa (ts/tv)
 289  Returns : value of kappa (a scalar)
 290  Args    : on set, new value (a scalar or undef, optional)
 291
 292
 293 =cut
 294
 295 sub kappa{
 296     my $self = shift;
 297     return $self->{'_kappa'} = shift if @_;
 298     return $self->{'_kappa'};
 299 }
 300
 301 =head2 num_site_classes
 302
 303  Title   : num_site_classes
 304  Usage   : $obj->num_site_classes($newval)
 305  Function: Get/Set the number of site classes for this model
 306  Returns : value of num_site_classes (a scalar)
 307  Args    : on set, new value (a scalar or undef, optional)
 308
 309
 310 =cut
 311
 312 sub num_site_classes{
 313     my $self = shift;
 314     return $self->{'_num_site_classes'} = shift if @_;
 315     return $self->{'_num_site_classes'};
 316 }
 317
 318 =head2 dnds_site_classes
 319
 320  Title   : dnds_site_classes
 321  Usage   : $obj->dnds_site_classes($newval)
 322  Function: Get/Set dN/dS site classes, a hashref
 323            with 2 keys, 'p' and 'w' which point to arrays
 324            one slot for each site class.
 325  Returns : value of dnds_site_classes (a hashref)
 326  Args    : on set, new value (a scalar or undef, optional)
 327
 328
 329 =cut
 330
 331 sub dnds_site_classes{
 332     my $self = shift;
 333     return $self->{'_dnds_site_classes'} = shift if @_;
 334     return $self->{'_dnds_site_classes'};
 335 }
 336
 337 =head2 get_pos_selected_sites
 338
 339  Title   : get_pos_selected_sites
 340  Usage   : my @sites = $modelresult->get_pos_selected_sites();
 341  Function: Get the sites which PAML has identified as under positive
 342            selection (w > 1).  This returns an array with each slot
 343            being a site, 4 values,
 344            site location (in the original alignment)
 345            Amino acid    (I *think* in the first sequence)
 346            P             (P value)
 347            Significance  (** indicated > 99%, * indicates >=95%)
 348  Returns : Array
 349  Args    : none
 350
 351
 352 =cut
 353
 354 sub get_pos_selected_sites{
 355    return @{$_[0]->{'_posselsites'} || []};
 356 }
 357
 358 =head2 add_pos_selected_site
 359
 360  Title   : add_pos_selected_site
 361  Usage   : $result->add_pos_selected_site($site,$aa,$pvalue,$signif);
 362  Function: Add a site to the list of positively selected sites
 363  Returns : count of the number of sites stored
 364  Args    : $site   - site number (in the alignment)
 365            $aa     - amino acid under selection
 366            $pvalue - float from 0->1 represent probability site is under selection according to this model
 367            $signif - significance (coded as either empty, '*', or '**'
 368
 369 =cut
 370
 371 sub add_pos_selected_site{
 372    my ($self,$site,$aa,$pvalue,$signif) = @_;
 373    push @{$self->{'_posselsites'}}, [ $site,$aa,$pvalue,$signif ];
 374    return scalar @{$self->{'_posselsites'}};
 375 }
 376
 377 =head2 get_NEB_pos_selected_sites
 378
 379  Title   : get_NEB_pos_selected_sites
 380  Usage   : my @sites = $modelresult->get_NEB_pos_selected_sites();
 381  Function: Get the sites which PAML has identified as under positive
 382            selection (w > 1) using Naive Empirical Bayes.
 383            This returns an array with each slot being a site, 4 values,
 384            site location (in the original alignment)
 385            Amino acid    (I *think* in the first sequence)
 386            P             (P value)
 387            Significance  (** indicated > 99%, * indicates > 95%)
 388            post mean for w
 389  Returns : Array
 390  Args    : none
 391
 392
 393 =cut
 394
 395 sub get_NEB_pos_selected_sites{
 396    return @{$_[0]->{'_NEBposselsites'} || []};
 397 }
 398
 399 =head2 add_NEB_pos_selected_site
 400
 401  Title   : add_NEB_pos_selected_site
 402  Usage   : $result->add_NEB_pos_selected_site($site,$aa,$pvalue,$signif);
 403  Function: Add a site to the list of positively selected sites
 404  Returns : count of the number of sites stored
 405  Args    : $site   - site number (in the alignment)
 406            $aa     - amino acid under selection
 407            $pvalue - float from 0->1 represent probability site is under selection according to this model
 408            $signif - significance (coded as either empty, '*', or '**'
 409            $postmean - post mean for w
 410
 411 =cut
 412
 413 sub add_NEB_pos_selected_site{
 414    my ($self,@args) = @_;
 415    push @{$self->{'_NEBposselsites'}}, [ @args ];
 416    return scalar @{$self->{'_NEBposselsites'}};
 417 }
 418
 419
 420
 421 =head2 get_BEB_pos_selected_sites
 422
 423  Title   : get_BEB_pos_selected_sites
 424  Usage   : my @sites = $modelresult->get_BEB_pos_selected_sites();
 425  Function: Get the sites which PAML has identified as under positive
 426            selection (w > 1) using Bayes Empirical Bayes.
 427            This returns an array with each slot being a site, 6 values,
 428            site location (in the original alignment)
 429            Amino acid    (I *think* in the first sequence)
 430            P             (P value)
 431            Significance  (** indicated > 99%, * indicates > 95%)
 432            post mean for w (mean)
 433            Standard Error for w (SE)
 434  Returns : Array
 435  Args    : none
 436
 437 =cut
 438
 439 sub get_BEB_pos_selected_sites{
 440    return @{$_[0]->{'_BEBposselsites'} || []};
 441 }
 442
 443 =head2 add_BEB_pos_selected_site
 444
 445  Title   : add_BEB_pos_selected_site
 446  Usage   : $result->add_BEB_pos_selected_site($site,$aa,$pvalue,$signif);
 447  Function: Add a site to the list of positively selected sites
 448  Returns : count of the number of sites stored
 449  Args    : $site   - site number (in the alignment)
 450            $aa     - amino acid under selection
 451            $pvalue - float from 0->1 represent probability site is under selection according to this model
 452            $signif - significance (coded as either empty, '*', or '**'
 453            $postmean - post mean for w
 454            $SE       - Standard Error for w
 455
 456 =cut
 457
 458 sub add_BEB_pos_selected_site{
 459    my ($self,@args) = @_;
 460    push @{$self->{'_BEBposselsites'}}, [ @args ];
 461    return scalar @{$self->{'_BEBposselsites'}};
 462 }
 463
 464 =head2 next_tree
 465
 466  Title   : next_tree
 467  Usage   : my $tree = $factory->next_tree;
 468  Function: Get the next tree from the factory
 469  Returns : L<Bio::Tree::TreeI>
 470  Args    : none
 471
 472 =cut
 473
 474 sub next_tree{
 475    my ($self,@args) = @_;
 476    return $self->{'_trees'}->[$self->{'_treeiterator'}++] || undef;
 477 }
 478
 479 =head2 get_trees
 480
 481  Title   : get_trees
 482  Usage   : my @trees = $result->get_trees;
 483  Function: Get all the parsed trees as an array
 484  Returns : Array of trees
 485  Args    : none
 486
 487
 488 =cut
 489
 490 sub get_trees{
 491    my ($self) = @_;
 492    return @{$self->{'_trees'} || []};
 493 }
 494
 495 =head2 rewind_tree_iterator
 496
 497  Title   : rewind_tree_iterator
 498  Usage   : $result->rewind_tree_iterator()
 499  Function: Rewinds the tree iterator so that next_tree can be
 500            called again from the beginning
 501  Returns : none
 502  Args    : none
 503
 504 =cut
 505
 506 sub rewind_tree_iterator {
 507     shift->{'_treeiterator'} = 0;
 508 }
 509
 510 =head2 add_tree
 511
 512  Title   : add_tree
 513  Usage   : $result->add_tree($tree);
 514  Function: Adds a tree
 515  Returns : integer which is the number of trees stored
 516  Args    : L<Bio::Tree::TreeI>
 517
 518 =cut
 519
 520 sub add_tree{
 521    my ($self,$tree) = @_;
 522    if( $tree && ref($tree) && $tree->isa('Bio::Tree::TreeI') ) {
 523        push @{$self->{'_trees'}},$tree;
 524    }
 525    return scalar @{$self->{'_trees'}};
 526 }
 527
 528 =head2 shape_params
 529
 530  Title   : shape_params
 531  Usage   : $obj->shape_params($newval)
 532  Function: Get/Set shape params for the distribution, 'alpha', 'beta'
 533            which is a hashref
 534            with 1 keys, 'p' and 'q'
 535  Returns : value of shape_params (a scalar)
 536  Args    : on set, new value (a scalar or undef, optional)
 537
 538
 539 =cut
 540
 541 sub shape_params{
 542     my $self = shift;
 543     return $self->{'_shape_params'} = shift if @_;
 544     return $self->{'_shape_params'};
 545 }
 546
 547 =head2 likelihood
 548
 549  Title   : likelihood
 550  Usage   : $obj->likelihood($newval)
 551  Function: log likelihood
 552  Returns : value of likelihood (a scalar)
 553  Args    : on set, new value (a scalar or undef, optional)
 554
 555
 556 =cut
 557
 558 sub likelihood{
 559     my $self = shift;
 560     return $self->{'likelihood'} = shift if @_;
 561     return $self->{'likelihood'};
 562 }
 563
 564 1;