[Bioperl-l] PAML nssites model result object

[Jason Stajich]

Resend with code pasted....

#!/usr/bin/perl -w
use strict;
use Bio::Tools::Phylo::PAML;

my $outcodeml = shift(@ARGV);

my $paml_parser = new Bio::Tools::Phylo::PAML(-file => "./$outcodeml",
					      -dir => "./");
my $result = $paml_parser->next_result();
my $MLmatrix = $result->get_MLmatrix(); # get MaxLikelihood Matrix
my @otus = $result->get_seqs;
if( $#{$MLmatrix} < 0 ) {
     for my $tree ($result->next_tree ) {
	for my $node ( $tree->get_nodes ) {
	    my $id;
	    if( $node->is_Leaf() ) {
		$id = $node->id;
	    } else {
		$id = "(".join(",", map { $_->id } grep { $_->is_Leaf }
			   $node->get_all_Descendents) .")";
	    }
	    if( ! $node->ancestor || ! $node->has_tag('t') ) {
		# skip when no values have been associated with this node
		# (like the root node)
		next;
	    }
             # I know this looks complicated
	    # but we use the get_tag_values method to pull out the annotations
	    # for each branch
	    # The ()[0] around the call is because get_tag_values returns a  
list
	    # if we want to just get the 1st item in the list we have
	    # to tell Perl we are treating it like an array.
	    # in the future get_tag_values needs to be smart and just
	    # return the 1st item in the array if called in scalar
	    # context
	
	    printf  
"%s\tt=%.3f\tS=%.1f\tN=%.1f\tdN/ 
dS=%.4f\tdN=%.4f\tdS=%.4f\tS*dS=%.1f\tN*dN=%.1f\n",
	    $id,
	    map { ($node->get_tag_values($_))[0] }
	    qw(t S N dN/dS dN dS), 'S*dS', 'N*dN';
	}
     }
} else {
     my $i =0;
     my @seqs = $result->get_seqs;
     for my $row ( @$MLmatrix ) {
	print $seqs[$i++]->display_id, join("\t",@$row), "\n";
     }
}

On Mar 9, 2005, at 5:41 PM, Jason Stajich wrote:

> I just updated things last week so this is brand-spanking-new.  I  
> don't know if I connected everything up for NSsites stuff quite yet   
> as that is handled in - the branch-specific parsing should work now.    
> I don't know if the synopsis code is really up to snuff either.  When  
> I get around to it I will try and see what still needs to be connected  
> in NSsites parsing.
>
> I don't think $node->param() is going to work -  
> $node->get_tag_values() is the way I've implemented it.
>
> <00parse_codeml.pl>
>
> -jason
> --
> Jason Stajich
> jason.stajich at duke.edu
> http://www.duke.edu/~jes12/
>
> On Mar 9, 2005, at 5:23 PM, Edward Chuong wrote:
>
>> Hi all,
>>
>> I'm trying to parse PAML results, and running into some trouble. I'm
>> using branch specific omega model, and I want to get the branch
>> specific ka/ks values out.
>> http://cvs.bioperl.org/cgi-bin/viewcvs/viewcvs.cgi/bioperl-live/Bio/ 
>> Tools/Phylo/PAML.pm?rev=HEAD&cvsroot=bioperl&content-type=text/ 
>> vnd.viewcvs-markup
>> says that $node->param('omega') should work, but Data::Dumper shows
>> that this value isn't stored in the node (only branch lengths and seq
>> IDs appear to be stored).
>>
>> I'm assuming that I can get these values out of the
>> get_NSSite_result() Bio::Tools::Phylo::PAML::ModelResult object, but
>> I'm not sure how to call it. The current synopsis uses
>> "get_model_params" but it seems to be out of date because it's not in
>> the current souce. The docs at
>> http://cvs.bioperl.org/cgi-bin/viewcvs/viewcvs.cgi/bioperl-live/Bio/ 
>> Tools/Phylo/PAML/Result.pm?rev=HEAD&cvsroot=bioperl&content- 
>> type=text/vnd.viewcvs-markup
>> say to use my
>> @results = @{$self->get_NSSite_results};
>> --that looks like a mistake, and I've tried
>> @result = $result->get_NSSite_results but that doesn't work either
>> (just get undefined objs).
>>
>> Am I doing something wrong, or is this functionality still being
>> worked on? I've tried using both 1.4 and the LIVE versions. Any help
>> is appreciated, thanks!
>>
>> -Ed
>> _______________________________________________
>> Bioperl-l mailing list
>> Bioperl-l at portal.open-bio.org
>> http://portal.open-bio.org/mailman/listinfo/bioperl-l
>>