Molphy documentation.

  use Bio::Tools::Phylo::Molphy;
  my $parser = new Bio::Tools::Phylo::Molphy(-file => 'output.protml');
  while( my $r = $parser->next_result ) {
    # r is a Bio::Tools::Phylo::Molphy::Result object

    # print the model name
    print $r->model, "\n";

    # get the substitution matrix
    # this is a hash of 3letter aa codes -> 3letter aa codes representing
    # substitution rate
    my $smat = $r->substitution_matrix;
    print "Arg -> Gln substitution rate is %d\n", 
          $smat->{'Arg'}->{'Gln'}, "\n";

    # get the transition probablity matrix
    # this is a hash of 3letter aa codes -> 3letter aa codes representing
    # transition probabilty
    my $tmat = $r->transition_probability_matrix;
    print "Arg -> Gln transition probablity is %.2f\n", 
          $tmat->{'Arg'}->{'Gln'}, "\n";

    # get the frequency for each of the residues
    my $rfreqs = $r->residue_frequencies;

    foreach my $residue ( keys %{$rfreqs} ) {
       printf "residue %s  expected freq: %.2f observed freq: %.2f\n",
              $residue,$rfreqs->{$residue}->[0], $rfreqs->{$residue}->[1];     
    }

    my @trees;
    while( my $t = $r->next_tree ) {
        push @trees, $t;
    }

    print "search space is ", $r->search_space, "\n",
          "1st tree score is ", $trees[0]->score, "\n";

    # writing to STDOUT, use -file => '>filename' to specify a file
    my $out = new Bio::TreeIO(-format => "newick");
    $out->write_tree($trees[0]); # writing only the 1st tree
  }

sub next_result {

   my ($self) = @_;

   # A little statemachine for the parser here
   my ($state,$transition_ct,
       @transition_matrix, %transition_mat, @resloc,) = ( 0,0);
   my ( %subst_matrix, @treelines, @treedata, %frequencies);
   my ( $treenum,$possible_trees, $model);
   my ($trans_type,$trans_amount);
   my $parsed = 0;
   while( defined ( $_ = $self->_readline()) ) {
       $parsed = 1;
       if( /^Relative Substitution Rate Matrix/ ) {
	   if( %subst_matrix ) { 
	       $self->_pushback($_);
	       last;
	   }
	   $state = 0;
	   my ( @tempdata);
	   @resloc = ();
	   while( defined ($_ = $self->_readline) ) {
	       last if (/^\s+$/);
	       # remove leading/trailing spaces
	       s/^\s+//;
	       s/\s+$//;
	       my @data = split;
	       my $i = 0;
	       for my $l ( @data ) {
		   if( $l =~ /\D+/ ) { 
		       push @resloc, $l;
		   }
		   $i++;
	       }
	       push @tempdata,\@ data;
	   }
	   my $i = 0;
	   for my $row ( @tempdata ) {
	       my $j = 0;
	       for my $col ( @$row ) {
		   if( $i == $j ) {
		       # empty string for diagonals
		       $subst_matrix{$resloc[$i]}->{$resloc[$j]} = '';
		   } else {
		       $subst_matrix{$resloc[$i]}->{$resloc[$j]} = $col;
		   }
		   $j++;
	       }
	       $i++;
	   }
       } elsif( /^Transition Probability Matrix/ ) {	   
	   if( /(1\.0e(5|7))\)\s+(\S+)/ ) {
	       $state = 1;
	       my $newtrans_type = "$3-$1";
	       $trans_amount = $1;
	       if( defined $trans_type ) {
		   # finish processing the transition_matrix
		   my $i =0;
		   foreach my $row ( @transition_matrix ) {
		       my $j = 0;
		       foreach my $col ( @$row ) {
			   $transition_mat{$trans_type}->{$resloc[$i]}->{$resloc[$j]} = $col;
			   $j++;
		       }
		       $i++;
		   }
	       }
	       $trans_type = $newtrans_type;
	       $transition_ct = 0;
	       @transition_matrix = ();
	   }
       } elsif ( /Acid Frequencies/ ) {
	   $state = 0;
	   $self->_readline(); # skip the next line
	   while( defined( $_ = $self->_readline) ) {
	       unless( /^\s+/) {
		   $self->_pushback($_);
		   last;
	       }
	       s/^\s+//;
	       s/\s+$//;
	       my ($index,$res,$model,$data) = split;
	       $frequencies{$res} = [ $model,$data];
	   }
       } elsif( /^(\d+)\s*\/\s*(\d+)\s+(.+)\s+model/ ) {
	   my @save = ($1,$2,$3);	   
	   # finish processing the transition_matrix
	   my $i =0;
	   foreach my $row ( @transition_matrix ) {
	       my $j = 0;
	       foreach my $col ( @$row ) {
		   $transition_mat{$trans_type}->{$resloc[$i]}->{$resloc[$j]} = $col;
		   $j++;
	       }
	       $i++;
	   }	   
	   if( defined $treenum ) { 	       
	       $self->_pushback($_);
	       last;
	   }
	   
	   $state = 2;	   
	   ($treenum,$possible_trees, $model) = @save;
	   $model =~ s/\s+/ /g;
       } elsif( $state == 1 ) {
	   next if( /^\s+$/ || /^\s+Ala/);
	   s/^\s+//;
	   s/\s+$//;
	   if( $trans_type eq '1PAM-1.0e7' ) {
	       # because the matrix is split up into 2-10 column sets 
	       push @{$transition_matrix[$transition_ct++]}, split ;	   
	       $transition_ct = 0 if $transition_ct % 20 == 0;
	   } elsif( $trans_type eq '1PAM-1.0e5' ) {
	       # because the matrix is split up into 2-10 column sets 
	       my ($res,@row) = split;
	       next if $transition_ct >= 20; # skip last 
	       push @{$transition_matrix[$transition_ct++]}, @row;	   	       
	   }
       } elsif( $state == 2 ) {
	   if( s/^(\d+)\s+(\-?\d+(\.\d+)?)\s+// ) {
	       push @treedata, [ $1,$2];
	   }
	   # save this for the end so that we can 
	   # be efficient and only open one tree parser
	   push @treelines, $_;
       }
   }
   # waiting till the end to do this, is it better
   my @trees;
   if( @treelines ) {
       my $strdat = IO::String->new(join('',@treelines));
       my $treeio = new Bio::TreeIO(-fh => $strdat,
				    -format => 'newick');
       while( my $tree = $treeio->next_tree ) {
	   if( @treedata ) {
	       my $dat = shift @treedata;
	       # set the associated information
	       $tree->id($dat->[0]);
	       $tree->score($dat->[1]);
	   }
	   push @trees, $tree;
       }
   }
   return undef unless( $parsed );
   my $result = new Bio::Tools::Phylo::Molphy::Result
       (-trees =>\@ trees,
	-substitution_matrix =>\% subst_matrix,
	-frequencies         =>\% frequencies,
	-model               => $model,
	-search_space        => $possible_trees,
	);
   while( my ($type,$mat) = each %transition_mat ) {
       $result->transition_probability_matrix( $type,$mat);
   }
   $result;

}

new	Description	Code
next_result	Description	Code