File Coverage

Bio/Tools/Phylo/Molphy/Result.pm
Criterion Covered Total %
statement 60 71 84.5
branch 24 36 66.6
condition 4 8 50.0
subroutine 10 11 90.9
pod 8 9 88.8
total 106 135 78.5


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::Tools::Phylo::Molphy::Result
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::Molphy::Result - container for data parsed from a ProtML run
17              
18             =head1 SYNOPSIS
19              
20             # do not use this object directly, you will get it back as part of a
21             # Molphy parser
22             use Bio::Tools::Phylo::Molphy;
23             my $parser = Bio::Tools::Phylo::Molphy->new(-file => 'output.protml');
24             while( my $r = $parser->next_result ) {
25             # r is a Bio::Tools::Phylo::Molphy::Result object
26              
27             # print the model name
28             print $r->model, "\n";
29              
30             # get the substitution matrix
31             # this is a hash of 3letter aa codes -> 3letter aa codes representing
32             # substitution rate
33             my $smat = $r->substitution_matrix;
34             print "Arg -> Gln substitution rate is %d\n",
35             $smat->{'Arg'}->{'Gln'}, "\n";
36              
37             # get the transition probablity matrix
38             # this is a hash of 3letter aa codes -> 3letter aa codes representing
39             # transition probabilty
40             my $tmat = $r->transition_probability_matrix;
41             print "Arg -> Gln transition probablity is %.2f\n",
42             $tmat->{'Arg'}->{'Gln'}, "\n";
43              
44             # get the frequency for each of the residues
45             my $rfreqs = $r->residue_frequencies;
46              
47             foreach my $residue ( keys %{$rfreqs} ) {
48             printf "residue %s expected freq: %.2f observed freq: %.2f\n",
49             $residue,$rfreqs->{$residue}->[0], $rfreqs->{$residue}->[1];
50             }
51              
52             my @trees;
53             while( my $t = $r->next_tree ) {
54             push @trees, $t;
55             }
56              
57             print "search space is ", $r->search_space, "\n",
58             "1st tree score is ", $trees[0]->score, "\n";
59              
60             # writing to STDOUT, use -file => '>filename' to specify a file
61             my $out = Bio::TreeIO->new(-format => "newick");
62             $out->write_tree($trees[0]); # writing only the 1st tree
63             }
64              
65              
66             =head1 DESCRIPTION
67              
68             A container for data parsed from a ProtML run.
69              
70              
71             =head1 FEEDBACK
72              
73             =head2 Mailing Lists
74              
75             User feedback is an integral part of the evolution of this and other
76             Bioperl modules. Send your comments and suggestions preferably to
77             the Bioperl mailing list. Your participation is much appreciated.
78              
79             bioperl-l@bioperl.org - General discussion
80             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
81              
82             =head2 Support
83              
84             Please direct usage questions or support issues to the mailing list:
85              
86             I
87              
88             rather than to the module maintainer directly. Many experienced and
89             reponsive experts will be able look at the problem and quickly
90             address it. Please include a thorough description of the problem
91             with code and data examples if at all possible.
92              
93             =head2 Reporting Bugs
94              
95             Report bugs to the Bioperl bug tracking system to help us keep track
96             of the bugs and their resolution. Bug reports can be submitted via the
97             web:
98              
99             https://github.com/bioperl/bioperl-live/issues
100              
101             =head1 AUTHOR - Jason Stajich
102              
103             Email jason-at-bioperl.org
104              
105             =head1 APPENDIX
106              
107             The rest of the documentation details each of the object methods.
108             Internal methods are usually preceded with a _
109              
110             =cut
111              
112              
113             # Let the code begin...
114              
115              
116             package Bio::Tools::Phylo::Molphy::Result;
117 1     1   3 use strict;
  1         2  
  1         25  
118              
119             # Object preamble - inherits from Bio::Root::Root
120              
121              
122              
123 1     1   17 use base qw(Bio::Root::Root);
  1         1  
  1         366  
124              
125             =head2 new
126              
127             Title : new
128             Usage : my $obj = Bio::Tools::Phylo::Molphy::Result->new();
129             Function: Builds a new Bio::Tools::Phylo::Molphy::Result object
130             Returns : Bio::Tools::Phylo::Molphy::Result
131             Args :
132              
133              
134             =cut
135              
136             sub new {
137 2     2 1 8 my($class,@args) = @_;
138              
139 2         16 my $self = $class->SUPER::new(@args);
140 2         11 my ($trees,$smat,$freq,
141             $model, $sspace,
142             ) = $self->_rearrange([qw(TREES SUBSTITUTION_MATRIX
143             FREQUENCIES
144             MODEL SEARCH_SPACE)], @args);
145              
146 2 50       6 if( $trees ) {
147 2 50       9 if(ref($trees) !~ /ARRAY/i ) {
148 0         0 $self->warn("Must provide a valid array reference to initialize trees");
149             } else {
150 2         3 foreach my $t ( @$trees ) {
151 10         10 $self->add_tree($t);
152             }
153             }
154             }
155             # initialize things through object methods to be a good
156             # little OO programmer
157 2 50       6 if( ref($smat) =~ /HASH/i ) {
158 2         4 $self->substitution_matrix($smat);
159             }
160 2 50       5 if( ref($freq) =~ /HASH/i ) {
161 2         5 $self->residue_frequencies($freq);
162             }
163            
164 2 50       6 $model && $self->model($model);
165 2 50       5 $sspace && $self->search_space($sspace);
166 2         3 $self->{'_treeiterator'} = 0;
167              
168 2         5 return $self;
169             }
170              
171             =head2 model
172              
173             Title : model
174             Usage : $obj->model($newval)
175             Function:
176             Returns : value of model
177             Args : newvalue (optional)
178              
179              
180             =cut
181              
182             sub model{
183 4     4 1 188 my ($self,$value) = @_;
184 4 100       7 if( defined $value) {
185 2         3 $self->{'model'} = $value;
186             }
187 4         10 return $self->{'model'};
188              
189             }
190              
191             =head2 substitution_matrix
192              
193             Title : substitution_matrix
194             Usage : my $smat = $result->subsitution_matrix;
195             Function: Get the relative substitution matrix calculated in the ML procedure
196             Returns : reference to hash of hashes where key is the aa/nt name and value
197             is another hash ref which contains keys for all the aa/nt
198             possibilities
199             Args : none
200              
201              
202             =cut
203              
204             sub substitution_matrix{
205 3     3 1 627 my ($self,$val) = @_;
206 3 100       12 if(defined $val ) {
207 2 50       7 if( ref($val) =~ /HASH/ ) {
208 2         2 foreach my $v (values %{$val} ) {
  2         6  
209 20 50       35 if( ref($v) !~ /HASH/i ) {
210 0         0 $self->warn("Must be a valid hashref of hashrefs for substition_matrix");
211 0         0 return;
212             }
213             }
214 2         3 $self->{'_substitution_matrix'} = $val;
215             } else {
216 0         0 $self->warn("Must be a valid hashref of hashrefs for substition_matrix");
217 0         0 return;
218             }
219             }
220 3         8 return $self->{'_substitution_matrix'};
221             }
222              
223             =head2 transition_probability_matrix
224              
225             Title : transition_probability_matrix
226             Usage : my $matrixref = $molphy->transition_probablity_matrix();
227             Function: Gets the observed transition probability matrix
228             Returns : hash of hashes of aa/nt transition to each other aa/nt
229             Args : Transition matrix type, typically
230             '1PAM-1.0e05' or '1PAM-1.0e07'
231              
232              
233             =cut
234              
235             sub transition_probability_matrix {
236 3     3 1 4 my ($self,$type,$val) = @_;
237 3 100       7 $type = '1PAM-1.0e7' unless defined $type;
238 3 100       5 if(defined $val ) {
239 2 50       5 if( ref($val) =~ /HASH/ ) {
240 2         2 foreach my $v (values %{$val} ) {
  2         6  
241 40 50       65 if( ref($v) !~ /HASH/i ) {
242 0         0 $self->warn("Must be a valid hashref of hashrefs for transition_probability_matrix");
243 0         0 return;
244             }
245             }
246 2         4 $self->{'_TPM'}->{$type} = $val;
247             } else {
248 0         0 $self->warn("Must be a valid hashref of hashrefs for transition_probablity_matrix");
249 0         0 return;
250             }
251             }
252              
253             # fix this for nucml where there are 2 values (one is just a transformation
254             # of the either, but how to represent?)
255 3         9 return $self->{'_TPM'}->{$type};
256             }
257              
258             =head2 residue_frequencies
259              
260             Title : residue_frequencies
261             Usage : my %data = $molphy->residue_frequencies()
262             Function: Get the modeled and expected frequencies for
263             each of the residues in the sequence
264             Returns : hash of either aa (protml) or nt (nucml) frequencies
265             each key will point to an array reference where
266             1st slot is model's expected frequency
267             2nd slot is observed frequency in the data
268             $hash{'A'}->[0] =
269             Args : none
270              
271              
272             =cut
273              
274             #'
275              
276             sub residue_frequencies {
277 3     3 1 907 my ($self,$val) = @_;
278 3 100       5 if(defined $val ) {
279 2 50       5 if( ref($val) =~ /HASH/ ) {
280 2         3 $self->{'_residue_frequencies'} = $val;
281             } else {
282 0         0 $self->warn("Must be a valid hashref of hashrefs for residue_frequencies");
283             }
284             }
285 3         3 return %{$self->{'_residue_frequencies'}};
  3         14  
286             }
287              
288             =head2 next_tree
289              
290             Title : next_tree
291             Usage : my $tree = $factory->next_tree;
292             Function: Get the next tree from the factory
293             Returns : L
294             Args : none
295              
296             =cut
297              
298             sub next_tree{
299 12     12 1 25 my ($self,@args) = @_;
300 12   100     28 return $self->{'_trees'}->[$self->{'_treeiterator'}++] || undef;
301             }
302              
303             =head2 rewind_tree
304              
305             Title : rewind_tree_iterator
306             Usage : $result->rewind_tree()
307             Function: Rewinds the tree iterator so that next_tree can be
308             called again from the beginning
309             Returns : none
310             Args : none
311              
312             =cut
313              
314             sub rewind_tree_iterator {
315 0     0 0 0 shift->{'_treeiterator'} = 0;
316             }
317              
318             =head2 add_tree
319              
320             Title : add_tree
321             Usage : $result->add_tree($tree);
322             Function: Adds a tree
323             Returns : integer which is the number of trees stored
324             Args : L
325              
326             =cut
327              
328             sub add_tree{
329 10     10 1 8 my ($self,$tree) = @_;
330 10 50 33     51 if( $tree && ref($tree) && $tree->isa('Bio::Tree::TreeI') ) {
      33        
331 10         7 push @{$self->{'_trees'}},$tree;
  10         12  
332             }
333 10         5 return scalar @{$self->{'_trees'}};
  10         14  
334             }
335              
336             =head2 search_space
337              
338             Title : search_space
339             Usage : $obj->search_space($newval)
340             Function:
341             Returns : value of search_space
342             Args : newvalue (optional)
343              
344              
345             =cut
346              
347             sub search_space{
348 4     4 1 5 my ($self,$value) = @_;
349 4 100       7 if( defined $value) {
350 2         3 $self->{'search_space'} = $value;
351             }
352 4         8 return $self->{'search_space'};
353             }
354              
355             1;