[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