Bio/Tools/Phylo/Molphy/Result.pm

   1 #
   2 # BioPerl module for Bio::Tools::Phylo::Molphy::Result
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Jason Stajich <jason@bioperl.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::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<bioperl-l@bioperl.org>
  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 use strict;
 118
 119 # Object preamble - inherits from Bio::Root::Root
 120
 121
 122
 123 use base qw(Bio::Root::Root);
 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   my($class,@args) = @_;
 138
 139   my $self = $class->SUPER::new(@args);
 140   my ($trees,$smat,$freq,
 141       $model, $sspace,
 142       ) = $self->_rearrange([qw(TREES SUBSTITUTION_MATRIX
 143                                 FREQUENCIES
 144                                 MODEL SEARCH_SPACE)], @args);
 145
 146   if( $trees ) {
 147       if(ref($trees) !~ /ARRAY/i ) {
 148           $self->warn("Must provide a valid array reference to initialize trees");
 149       } else {
 150           foreach my $t ( @$trees ) {
 151               $self->add_tree($t);
 152           }
 153       }
 154   }
 155   # initialize things through object methods to be a good
 156   # little OO programmer
 157   if( ref($smat) =~ /HASH/i ) {
 158       $self->substitution_matrix($smat);
 159   }
 160   if( ref($freq) =~ /HASH/i ) {
 161       $self->residue_frequencies($freq);
 162   }
 163
 164   $model && $self->model($model);
 165   $sspace && $self->search_space($sspace);
 166   $self->{'_treeiterator'} = 0;
 167
 168   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    my ($self,$value) = @_;
 184    if( defined $value) {
 185       $self->{'model'} = $value;
 186     }
 187     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    my ($self,$val) = @_;
 206    if(defined $val ) {
 207        if( ref($val) =~ /HASH/ ) {
 208            foreach my $v (values %{$val} ) {
 209                if( ref($v) !~ /HASH/i ) {
 210                    $self->warn("Must be a valid hashref of hashrefs for substition_matrix");
 211                    return;
 212                }
 213            }
 214            $self->{'_substitution_matrix'} = $val;
 215        } else {
 216            $self->warn("Must be a valid hashref of hashrefs for substition_matrix");
 217            return;
 218        }
 219    }
 220    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    my ($self,$type,$val) = @_;
 237    $type = '1PAM-1.0e7' unless defined $type;
 238    if(defined $val ) {
 239        if( ref($val) =~ /HASH/ ) {
 240            foreach my $v (values %{$val} ) {
 241                if( ref($v) !~ /HASH/i ) {
 242                    $self->warn("Must be a valid hashref of hashrefs for transition_probability_matrix");
 243                    return;
 244                }
 245            }
 246            $self->{'_TPM'}->{$type} = $val;
 247        } else {
 248            $self->warn("Must be a valid hashref of hashrefs for transition_probablity_matrix");
 249            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    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    my ($self,$val) = @_;
 278    if(defined $val ) {
 279        if( ref($val) =~ /HASH/ ) {
 280            $self->{'_residue_frequencies'} = $val;
 281        } else {
 282            $self->warn("Must be a valid hashref of hashrefs for residue_frequencies");
 283        }
 284    }
 285    return %{$self->{'_residue_frequencies'}};
 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<Bio::Tree::TreeI>
 294  Args    : none
 295
 296 =cut
 297
 298 sub next_tree{
 299    my ($self,@args) = @_;
 300    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     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<Bio::Tree::TreeI>
 325
 326 =cut
 327
 328 sub add_tree{
 329    my ($self,$tree) = @_;
 330    if( $tree && ref($tree) && $tree->isa('Bio::Tree::TreeI') ) {
 331        push @{$self->{'_trees'}},$tree;
 332    }
 333    return scalar @{$self->{'_trees'}};
 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    my ($self,$value) = @_;
 349    if( defined $value) {
 350       $self->{'search_space'} = $value;
 351     }
 352     return $self->{'search_space'};
 353 }
 354
 355 1;