[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