[Bioperl-l] PAML nssites model result object

Jason Stajich jason.stajich at duke.edu
Thu Mar 10 09:48:02 EST 2005


The script needs to be adjusted for NSsites because their are trees are  
associated with each model result so you need one more loop on the  
get_NSSite_results. I added some code to the script to print out the  
positively selected sites as well.

#!/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;
# process the NSsites results
for my $ns_result ( $result->get_NSSite_results ) {
     print "model ", $ns_result->model_num, " ",  
$ns_result->model_description, "\n";
     while ( my $tree = $ns_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;
	    }
		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';
	}
     }
     print "positively selected sites:\n";
     #  get the positively select sites
     for my $site ( $ns_result->get_pos_selected_sites ) {
	print join(" ", @$site, "\n");
     }
     print "\n";
}


--
Jason Stajich
jason.stajich at duke.edu
http://www.duke.edu/~jes12/

On Mar 9, 2005, at 6:37 PM, Edward Chuong wrote:

> Hi Jason,
>
> Thanks for the help.
>
> The code seems to get stuck at
>
>  if( ! $node->ancestor || ! $node->has_tag('t') ) {
> (this condition turns out true for every node, not just root, so it
> always hits "next")
>
> I used Data::Dumper to check on the node and I've pasted the
> results--it seems like those tags aren't being sent in?
>
>
> Thanks!
> -Ed
>
> '_root_cleanup_methods' => [
> 		  sub { "DUMMY" }
> 		],
> '_creation_id' => 0,
> '_branch_length' => '0.613722',
> '_desc' => {},
> '_id' => 'NP_033437.2_mus',
> '_ancestor' => bless( {
> 	 '_root_cleanup_methods' => [
> 					  $VAR1->{'_root_cleanup_methods'}[0]
> 					],
> 	 '_creation_id' => 3,
> 	 '_desc' => {
> 				  '2' => bless( {
> 					  '_root_cleanup_methods' => [
> 						   $VAR1->{'_root_cleanup_methods'}[0]
> 						 ],
> 					  '_creation_id' => 2,
> 					  '_branch_length' => '0.768322',
> 					  '_desc' => {},
> 					  '_id' => 'PM_BWp0001H02f',
> 					  '_ancestor' => $VAR1->{'_ancestor'},
> 					  '_root_verbose' => 0
> 					}, 'Bio::Tree::Node' ),
> 				  '0' => $VAR1,
> 				  '1' => bless( {
> 					  '_root_cleanup_methods' => [
> 												   $VAR1->{'_root_cleanup_methods'}[0]
> 												 ],
> 					  '_creation_id' => 1,
> 					  '_branch_length' => '0.366319',
> 					  '_desc' => {},
> 					  '_id' => 'NP_742070.1_rat',
> 					  '_ancestor' => $VAR1->{'_ancestor'},
> 					  '_root_verbose' => 0
> 					}, 'Bio::Tree::Node' )
> 				},
> 	 '_id' => '',
> 	 '_height' => undef,
> 	 '_root_verbose' => 0
>    }, 'Bio::Tree::Node' ),
> '_root_verbose' => 0
> }, 'Bio::Tree::Node' );
>
>
> On Wed, 9 Mar 2005 18:01:34 -0500, Jason Stajich  
> <jason.stajich at duke.edu> wrote:
>> 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
>>>>
>>
>>
>
>
> -- 
> Edward Chuong
> (949) 939-2732
> AIM: edawad85
>



More information about the Bioperl-l mailing list