[Bioperl-l] PAML nssites model result object

Edward Chuong echuong at gmail.com
Thu Mar 10 15:47:17 EST 2005


Hi,

I think the problem is that the $result object, according to Dumper,
doesn't store any ModelResult (NSSite_results), so the for loop
condition in this code ($result->get_NSSite_results) is never true. Is
this working on a mlc file that you have, and if so, can you send it
so I can see if it's a problem on my side?

Thanks
-Ed


On Thu, 10 Mar 2005 09:48:02 -0500, Jason Stajich
<jason.stajich at duke.edu> wrote:
> 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
> >
> 
> 


-- 
Edward Chuong
(949) 939-2732
AIM: edawad85


More information about the Bioperl-l mailing list