File Coverage

Bio/Tools/Phylo/PAML/ModelResult.pm
Criterion Covered Total %
statement 95 111 85.5
branch 47 66 71.2
condition 6 13 46.1
subroutine 19 21 90.4
pod 19 19 100.0
total 186 230 80.8


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::Tools::Phylo::PAML::ModelResult
3             #
4             # Please direct questions and support issues to
5             #
6             # Cared for by Jason Stajich
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
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 1     1   4 use strict;
  1         2  
  1         25  
94              
95             # Object preamble - inherits from Bio::Root::Root
96              
97              
98 1     1   3 use base qw(Bio::Root::Root);
  1         1  
  1         1800  
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 9     9 1 51 my($class,@args) = @_;
136              
137 9         41 my $self = $class->SUPER::new(@args);
138 9         59 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 9 50       35 if( $trees ) {
156 9 50       46 if(ref($trees) !~ /ARRAY/i ) {
157 0         0 $self->warn("Must provide a valid array reference to initialize trees");
158             } else {
159 9         18 foreach my $t ( @$trees ) {
160 9         35 $self->add_tree($t);
161             }
162             }
163             }
164 9         17 $self->{'_treeiterator'} = 0;
165 9 100       24 if( $pos_sites ) {
166 6 50       21 if(ref($pos_sites) !~ /ARRAY/i ) {
167 0         0 $self->warn("Must provide a valid array reference to initialize pos_sites");
168             } else {
169 6         12 foreach my $s ( @$pos_sites ) {
170 8 50       16 if( ref($s) !~ /ARRAY/i ) {
171 0         0 $self->warn("Need an array reference for each entry in the pos_sites object");
172 0         0 next;
173             }
174 8         11 $self->add_pos_selected_site(@$s);
175             }
176             }
177             }
178 9 100       25 if( $beb_sites ) {
179 6 50       19 if(ref($beb_sites) !~ /ARRAY/i ) {
180 0         0 $self->warn("Must provide a valid array reference to initialize beb_sites");
181             } else {
182 6         11 foreach my $s ( @$beb_sites ) {
183 9 50       26 if( ref($s) !~ /ARRAY/i ) {
184 0         0 $self->warn("need an array ref for each entry in the beb_sites object");
185 0         0 next;
186             }
187 9         17 $self->add_BEB_pos_selected_site(@$s);
188             }
189             }
190             }
191 9 100       19 if( $neb_sites ) {
192 6 50       18 if(ref($neb_sites) !~ /ARRAY/i ) {
193 0         0 $self->warn("Must provide a valid array reference to initialize neb_sites");
194             } else {
195 6         11 foreach my $s ( @$neb_sites ) {
196 8 50       18 if( ref($s) !~ /ARRAY/i ) {
197 0         0 $self->warn("need an array ref for each entry in the neb_sites object");
198 0         0 next;
199             }
200 8         12 $self->add_NEB_pos_selected_site(@$s);
201             }
202             }
203             }
204              
205 9 100       72 defined $modelnum && $self->model_num($modelnum);
206 9 100       27 defined $modeldesc && $self->model_description($modeldesc);
207 9 50       37 defined $kappa && $self->kappa($kappa);
208 9 100       33 defined $timeused && $self->time_used($timeused);
209 9 50       36 defined $likelihood && $self->likelihood($likelihood);
210              
211 9   100     38 $self->num_site_classes($num_site_classes || 0);
212 9 100       19 if( defined $dnds_classes ) {
213 8 50 33     72 if( ref($dnds_classes) !~ /HASH/i ||
214             ! defined $dnds_classes->{'p'} ||
215             ! defined $dnds_classes->{'w'} ) {
216 0         0 $self->warn("-dnds_site_classes expects a hashref with keys p and w");
217             } else {
218 8         21 $self->dnds_site_classes($dnds_classes);
219             }
220             }
221 9 100       22 if( defined $shape_params ) {
222 3 50       11 if( ref($shape_params) !~ /HASH/i ) {
223 0         0 $self->warn("-shape_params expects a hashref not $shape_params\n");
224             } else {
225 3         9 $self->shape_params($shape_params);
226             }
227             }
228 9         83 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 14     14 1 2476 my $self = shift;
245 14 100       44 return $self->{'_num'} = shift if @_;
246 7         28 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 12     12 1 20 my $self = shift;
263 12 100       31 return $self->{'_model_description'} = shift if @_;
264 5         87 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 7     7 1 8 my $self = shift;
280 7 50       20 return $self->{'_time_used'} = shift if @_;
281 0         0 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 15     15 1 25 my $self = shift;
297 15 100       35 return $self->{'_kappa'} = shift if @_;
298 6         22 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 17     17 1 32 my $self = shift;
314 17 100       46 return $self->{'_num_site_classes'} = shift if @_;
315 8         36 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 12     12 1 19 my $self = shift;
333 12 100       34 return $self->{'_dnds_site_classes'} = shift if @_;
334 4         10 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 1 50   1 1 832 return @{$_[0]->{'_posselsites'} || []};
  1         8  
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 8     8 1 9 my ($self,$site,$aa,$pvalue,$signif) = @_;
373 8         9 push @{$self->{'_posselsites'}}, [ $site,$aa,$pvalue,$signif ];
  8         13  
374 8         7 return scalar @{$self->{'_posselsites'}};
  8         9  
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 1 50   1 1 841 return @{$_[0]->{'_NEBposselsites'} || []};
  1         7  
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 8     8 1 10 my ($self,@args) = @_;
415 8         4 push @{$self->{'_NEBposselsites'}}, [ @args ];
  8         22  
416 8         5 return scalar @{$self->{'_NEBposselsites'}};
  8         12  
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 0 0   0 1 0 return @{$_[0]->{'_BEBposselsites'} || []};
  0         0  
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 9     9 1 17 my ($self,@args) = @_;
460 9         7 push @{$self->{'_BEBposselsites'}}, [ @args ];
  9         30  
461 9         10 return scalar @{$self->{'_BEBposselsites'}};
  9         15  
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
470             Args : none
471              
472             =cut
473              
474             sub next_tree{
475 5     5 1 9 my ($self,@args) = @_;
476 5   50     22 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 2     2 1 3 my ($self) = @_;
492 2 50       3 return @{$self->{'_trees'} || []};
  2         12  
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 0     0 1 0 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
517              
518             =cut
519              
520             sub add_tree{
521 9     9 1 14 my ($self,$tree) = @_;
522 9 50 33     87 if( $tree && ref($tree) && $tree->isa('Bio::Tree::TreeI') ) {
      33        
523 9         15 push @{$self->{'_trees'}},$tree;
  9         29  
524             }
525 9         9 return scalar @{$self->{'_trees'}};
  9         34  
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 4     4 1 5 my $self = shift;
543 4 100       12 return $self->{'_shape_params'} = shift if @_;
544 1         3 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 15     15 1 20 my $self = shift;
560 15 100       45 return $self->{'likelihood'} = shift if @_;
561 6         20 return $self->{'likelihood'};
562             }
563              
564             1;