Bio::Tree Statistics
SummaryIncluded librariesPackage variablesSynopsisDescriptionGeneral documentationMethods
Toolbar
WebCvs
Summary
Bio::Tree::Statistics - Calculate certain statistics for a Tree
Package variables
No package variables defined.
Inherit
Bio::Root::Root
Synopsis
  use Bio::Tree::Statistics;
Description
This should be where Tree statistics are calculated. It was
previously where statistics from a Coalescent simulation.
It now contains several methods for calculating Tree-Trait statistics.
Methods
assess_bootstrapDescriptionCode
cherriesDescriptionCode
fitchDescriptionCode
psDescriptionCode
fitch_upDescriptionCode
fitch_downDescriptionCode
_persistence
No description
Code
persistenceDescriptionCode
_count_subclusters
No description
Code
count_subclustersDescriptionCode
_count_leaves
No description
Code
count_leavesDescriptionCode
_phylotype_length
No description
Code
phylotype_lengthDescriptionCode
_sum_of_leaf_distances
No description
Code
sum_of_leaf_distancesDescriptionCode
genetic_diversityDescriptionCode
statratioDescriptionCode
_node_ai
No description
Code
aiDescriptionCode
_node_mc
No description
Code
mcDescriptionCode
Methods description
assess_bootstrapcode    nextTop
 Title   : assess_bootstrap
Usage : my $tree_with_bs = $stats->assess_bootstrap(\@bs_trees);
Function: Calculates the bootstrap for internal nodes based on
Returns : Bio::Tree::TreeI
Args : Arrayref of Bio::Tree::TreeIs
cherriescodeprevnextTop
  Example    : cherries($tree, $node);
Description: Count number of paired leaf nodes
in a binary tree
Returns : integer
Exceptions :
Args : 1. Bio::Tree::TreeI object
2. Bio::Tree::NodeI object within the tree, optional
Commonly used statistics assume a binary tree, but this methods
returns a value even for trees with polytomies.
fitchcodeprevnextTop
  Example    : fitch($tree, $key, $node);
Description: Calculates Parsimony Score (PS) and internal trait
values using the Fitch 1971 parsimony algorithm for
the subtree a defined by the (internal) node.
Node defaults to the root.
Returns : true on success
Exceptions : leaf nodes have to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. trait name string
3. Bio::Tree::NodeI object within the tree, optional
Runs first fitch_up that calculates parsimony scores and then
fitch_down that should resolve most of the trait/character state
ambiguities.
Fitch, W.M., 1971. Toward defining the course of evolution: minimal
change for a specific tree topology. Syst. Zool. 20, 406-416.
You can access calculated parsimony values using:
  $score = $node->->get_tag_values('ps_score');
and the trait value with:
  $traitvalue = $node->->get_tag_values('ps_trait'); # only the first
@traitvalues = $node->->get_tag_values('ps_trait');
Note that there can be more that one trait value, especially for the
root node.
pscodeprevnextTop
  Example    : ps($tree, $key, $node);
Description: Calculates Parsimony Score (PS) from Fitch 1971
parsimony algorithm for the subtree as defined
by the (internal) node.
Node defaults to the root.
Returns : integer, 1 < PS < n, where n is number of branches
Exceptions : leaf nodes have to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. trait name string
3. Bio::Tree::NodeI object within the tree, optional
This is the first half of the Fitch algorithm that is enough for
calculating the resolved parsimony values. The trait/chararacter
states are commonly left in ambiguous state. To resolve them, run
fitch_down.
fitch_upcodeprevnextTop
  Example    : fitch_up($tree, $key, $node);
Description: Calculates Parsimony Score (PS) from the Fitch 1971
parsimony algorithm for the subtree as defined
by the (internal) node.
Node defaults to the root.
Returns : integer, 1< PS < n, where n is number of branches
Exceptions : leaf nodes have to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. trait name string
3. Bio::Tree::NodeI object within the tree, optional
This is a more generic name for ps and indicates that it performs
the first bottom-up tree traversal that calculates the parsimony score
but usually leaves trait/character states ambiguous. If you are
interested in internal trait states, running fitch_down should
resolve most of the ambiguities.
fitch_downcodeprevnextTop
  Example    : fitch_down($tree, $node);
Description: Runs the second pass from Fitch 1971
parsimony algorithm to resolve ambiguous
trait states left by first pass.
by the (internal) node.
Node defaults to the root.
Returns : true
Exceptions : dies unless the trait is defined in all nodes
Args : 1. Bio::Tree::TreeI object
2. Bio::Tree::NodeI object within the tree, optional
Before running this method you should have ran fitch_up (alias to
ps ). Note that it is not guaranteed that all states are completely
resolved.
persistencecodeprevnextTop
  Example    : persistence($tree, $node);
Description: Calculates the persistence
for node in the subtree defined by the (internal)
node. Node defaults to the root.
Returns : int, number of generations trait value has to remain same
Exceptions : all the nodes need to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. Bio::Tree::NodeI object within the tree, optional
Persistence measures the stability that the trait value has in a
tree. It expresses the number of generations the trait value remains
the same. All the decendants of the root in the same generation have
to share the same value.
Depends on Fitch's parsimony score (PS).
count_subclusterscodeprevnextTop
  Example    : count_clusters($tree, $node);
Description: Calculates the number of sub-clusters
in the subtree defined by the (internal)
node. Node defaults to the root.
Returns : int, count
Exceptions : all the nodes need to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. Bio::Tree::NodeI object within the tree, optional
Depends on Fitch's parsimony score (PS).
count_leavescodeprevnextTop
  Example    : count_leaves($tree, $node);
Description: Calculates the number of leaves with same trait
value as root in the subtree defined by the (internal)
node. Requires an unbroken line of identical trait values.
Node defaults to the root.
Returns : int, number of leaves with this trait value
Exceptions : all the nodes need to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. Bio::Tree::NodeI object within the tree, optional
Depends on Fitch's parsimony score (PS).
phylotype_lengthcodeprevnextTop
  Example    : phylotype_length($tree, $node);
Description: Sums up the branch lengths within phylotype
exluding the subclusters where the trait values
are different
Returns : float, length
Exceptions : all the nodes need to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. Bio::Tree::NodeI object within the tree, optional
Depends on Fitch's parsimony score (PS).
sum_of_leaf_distancescodeprevnextTop
  Example    : sum_of_leaf_distances($tree, $node);
Description: Sums up the branch lengths from root to leaf
exluding the subclusters where the trait values
are different
Returns : float, length
Exceptions : all the nodes need to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. Bio::Tree::NodeI object within the tree, optional
Depends on Fitch's parsimony score (PS).
genetic_diversitycodeprevnextTop
  Example    : genetic_diversity($tree, $node);
Description: Diversity is the sum of root to leaf distances
within the phylotype normalised by number of leaf
nodes
Returns : float, value of genetic diversity
Exceptions : all the nodes need to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. Bio::Tree::NodeI object within the tree, optional
Depends on Fitch's parsimony score (PS).
statratiocodeprevnextTop
  Example    : statratio($tree, $node);
Description: Ratio of the stem length and the genetic diversity of the
phylotype genetic_diversity
Returns : float, separation score
Exceptions : all the nodes need to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. Bio::Tree::NodeI object within the tree, optional
Statratio gives a measure of separation and variability within the phylotype.
Larger values identify more rapidly evolving and recent phylotypes.
Depends on Fitch's parsimony score (PS).
aicodeprevnextTop
  Example    : ai($tree, $key, $node);
Description: Calculates the Association Index (AI) of Whang et
al. 2001 for the subtree defined by the (internal)
node. Node defaults to the root.
Returns : real
Exceptions : leaf nodes have to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. trait name string
3. Bio::Tree::NodeI object within the tree, optional
Association index (AI) gives a more fine grained results than PS since the result is a real number. ~0 <= AI.
Wang, T.H., Donaldson, Y.K., Brettle, R.P., Bell, J.E., Simmonds, P., 2001. Identification of shared populations of human immunodeficiency Virus Type 1 infecting microglia and tissue macrophages outside the central nervous system. J. Virol. 75 (23), 11686-11699.
mccodeprevnextTop
  Example    : mc($tree, $key, $node);
Description: Calculates the Monophyletic Clade (MC) size statistics
for the subtree a defined by the (internal) node.
Node defaults to the root;
Returns : hashref with trait values as keys
Exceptions : leaf nodes have to have the trait defined
Args : 1. Bio::Tree::TreeI object
2. trait name string
3. Bio::Tree::NodeI object within the tree, optional
Monophyletic Clade (MC) size statistics by Salemi at al 2005. It is calculated for each trait value. 1 <= MC <= nx, where nx is the
number of tips with value x:
pick the internal node with maximim value for number of of tips with only trait x MC was defined by Parker et al 2008. Salemi, M., Lamers, S.L., Yu, S., de Oliveira, T., Fitch, W.M., McGrath, M.S., 2005. Phylodynamic analysis of Human Immunodeficiency Virus Type 1 in distinct brain compartments provides a model for the neuropathogenesis of AIDS. J. Virol. 79 (17), 11343-11352. Parker, J., Rambaut A., Pybus O., 2008. Correlating viral phenotypes with phylogeny: Accounting for phylogenetic uncertainty Infection, Genetics and Evolution 8 (2008), 239-246.
Methods code
assess_bootstrapdescriptionprevnextTop
sub assess_bootstrap {
   my ($self,$bs_trees,$guide_tree) = @_;
   my @consensus;

   # internal nodes are defined by their children
my (%lookup,%internal); my $i = 0; for my $tree ( $guide_tree, @$bs_trees ) { # Do this as a top down approach, can probably be
# improved by caching internal node states, but not going
# to worry about it right now.
my @allnodes = $tree->get_nodes; my @internalnodes = grep { ! $_->is_Leaf } @allnodes; for my $node ( @internalnodes ) { my @tips = sort map { $_->id } grep { $_->is_Leaf() } $node->get_all_Descendents; my $id = "(".join(",", @tips).")"; if( $i == 0 ) { $internal{$id} = $node->internal_id; } else { $lookup{$id}++; } } $i++; } my @save; for my $l ( keys %lookup ) { if( defined $internal{$l} ) {#&& $lookup{$l} > $min_seen ) {
my $intnode = $guide_tree->find_node(-internal_id => $internal{$l}); $intnode->bootstrap(sprintf("%d",100 * $lookup{$l} / $i));
} } return $guide_tree;
}
cherriesdescriptionprevnextTop
sub cherries {
($;$) {    my $self = shift;
    my $tree = shift;
    my $node = shift || $tree->get_root_node;

    my $cherries = 0;
    my @descs = $node->each_Descendent;

    if ($descs[0]->is_Leaf and $descs[1]->is_Leaf) {
        if ($descs[3]) { #polytomy at leaf level
$cherries = 0; } else { $cherries = 1; } } else { # recurse
foreach my $desc (@descs) { $cherries += $self->cherries($tree, $desc); } } return $cherries;
}
fitchdescriptionprevnextTop
sub fitch {
    my $self = shift;
    my $tree = shift;
    my $key = shift || $self->throw("Trait name is needed");
    my $node = shift || $tree->get_root_node;

    $self->fitch_up($tree, $key, $node);
    $self->fitch_down($tree, $node);
}
psdescriptionprevnextTop
sub ps {
 shift->fitch_up(@_)
}
fitch_updescriptionprevnextTop
sub fitch_up {
    my $self = shift;
    my $tree = shift;
    my $key = shift || $self->throw("Trait name is needed");
    my $node = shift || $tree->get_root_node;

    if ($node->is_Leaf) {
        $self->throw ("ERROR: ". $node->internal_id. " needs a value for trait $key")
            unless $node->has_tag($key);
        $node->set_tag_value('ps_trait', $node->get_tag_values($key) );
        $node->set_tag_value('ps_score', 0 );
        return; # end of recursion
} foreach my $child ($node->each_Descendent) { $self->fitch_up($tree, $key, $child); } my %intersection; my %union; my $score; foreach my $child ($node->each_Descendent) { foreach my $trait ($child->get_tag_values('ps_trait') ) { $intersection{$trait}++ if $union{$trait}; $union{$trait}++; } $score += $child->get_tag_values('ps_score'); } if (keys %intersection) { $node->set_tag_value('ps_trait', keys %intersection); $node->set_tag_value('ps_score', $score); } else { $node->set_tag_value('ps_trait', keys %union); $node->set_tag_value('ps_score', $score+1); } if ($self->verbose) { print "-- node --------------------------\n"; print "iID: ", $node->internal_id, " (", $node->id, ")\n"; print "Trait: ", join (', ', $node->get_tag_values('ps_trait') ), "\n"; print "length :", scalar($node->get_tag_values('ps_score')) , "\n"; } return scalar $node->get_tag_values('ps_score');
}
fitch_downdescriptionprevnextTop
sub fitch_down {
    my $self = shift;
    my $tree = shift;
    my $node = shift || $tree->get_root_node;

    my $key = 'ps_trait';
    $self->throw ("ERROR: ". $node->internal_id. " needs a value for $key")
        unless $node->has_tag($key);

    my $nodev;
    foreach my $trait ($node->get_tag_values($key) ) {
        $nodev->{$trait}++;
    }

    foreach my $child ($node->each_Descendent) {
        next if $child->is_Leaf;  # end of recursion
my $intersection; foreach my $trait ($child->get_tag_values($key) ) { $intersection->{$trait}++ if $nodev->{$trait}; } $self->fitch_down($tree, $child); $child->set_tag_value($key, keys %$intersection); } return 1; # success
}
_persistencedescriptionprevnextTop
sub _persistence {
    my $self = shift;
    my $tree = shift;
    my $node = shift;
    my $value = shift || $self->throw("Value is needed");


    my $key  = 'ps_trait';

    $self->throw("Node is needed") unless $node->isa('Bio::Tree::NodeI');

    return 0 unless $node->get_tag_values($key) eq $value; # wrong value
return 1 if $node->is_Leaf; # end of recursion
my $persistence = 10000000; # an arbitrarily large number
foreach my $child ($node->each_Descendent) { my $pers = $self->_persistence($tree, $child, $value); $persistence = $pers if $pers < $persistence; } return $persistence + 1;
}
persistencedescriptionprevnextTop
sub persistence {
    my $self = shift;
    my $tree = shift;
    my $node = shift  || $tree->get_root_node;
    $self->throw("Node is needed") unless $node->isa('Bio::Tree::NodeI');

    my $key  = 'ps_trait';
    my $value = $node->get_tag_values($key);

    #calculate
my $persistence = $self->_persistence($tree, $node, $value); $node->set_tag_value('persistance', $persistence); return $persistence;
}
_count_subclustersdescriptionprevnextTop
sub _count_subclusters {
    my $self = shift;
    my $tree = shift;
    my $node = shift;
    my $value = shift || $self->throw("Value is needed");

    my $key  = 'ps_trait';

    $self->throw ("ERROR: ". $node->internal_id. " needs a value for trait $key")
        unless $node->has_tag($key);

    if ($node->get_tag_values($key) eq $value) {
        if ($node->get_tag_values('ps_score') == 0) {
            return 0;
        } else {
            my $count = 0;
            foreach my $child ($node->each_Descendent) {
                $count += $self->_count_subclusters($tree, $child, $value);
            }
            return $count;
        }
    }
    return 1;
}
count_subclustersdescriptionprevnextTop
sub count_subclusters {
    my $self = shift;
    my $tree = shift;
    my $node = shift  || $tree->get_root_node;
    $self->throw("Node is needed") unless $node->isa('Bio::Tree::NodeI');

    my $key  = 'ps_trait';
    my $value = $node->get_tag_values($key);

    return $self->_count_subclusters($tree, $node, $value);
}
_count_leavesdescriptionprevnextTop
sub _count_leaves {
    my $self = shift;
    my $tree = shift;
    my $node = shift  || $tree->get_root_node;
    my $value = shift;

    my $key  = 'ps_trait';

    $self->throw ("ERROR: ". $node->internal_id. " needs a value for trait $key")
        unless $node->has_tag($key);

    if ($node->get_tag_values($key) eq $value) {
        #print $node->id, ": ", $node->get_tag_values($key), "\n";
return 1 if $node->is_Leaf; # end of recursion
my $count = 0; foreach my $child ($node->each_Descendent) { $count += $self->_count_leaves($tree, $child, $value); } return $count; } return 0;
}
count_leavesdescriptionprevnextTop
sub count_leaves {
    my $self = shift;
    my $tree = shift;
    my $node = shift  || $tree->get_root_node;
    $self->throw("Node is needed") unless $node->isa('Bio::Tree::NodeI');

    my $key  = 'ps_trait';
    my $value = $node->get_tag_values($key);

    return $self->_count_leaves($tree, $node, $value);
}
_phylotype_lengthdescriptionprevnextTop
sub _phylotype_length {
    my $self = shift;
    my $tree = shift;
    my $node = shift;
    my $value = shift;

    my $key  = 'ps_trait';

    $self->throw ("ERROR: ". $node->internal_id. " needs a value for trait $key")
        unless $node->has_tag($key);

    return 0 if $node->get_tag_values($key) ne $value;
    return $node->branch_length if $node->is_Leaf; # end of recursion
my $length = 0; foreach my $child ($node->each_Descendent) { my $sub_len = $self->_phylotype_length($tree, $child, $value); $length += $sub_len; $length += $child->branch_length if not $child->is_Leaf and $sub_len; } return $length;
}
phylotype_lengthdescriptionprevnextTop
sub phylotype_length {
    my $self = shift;
    my $tree = shift;
    my $node = shift  || $tree->get_root_node;

    my $key  = 'ps_trait';
    my $value = $node->get_tag_values($key);

    return $self->_phylotype_length($tree, $node, $value);
}
_sum_of_leaf_distancesdescriptionprevnextTop
sub _sum_of_leaf_distances {
    my $self = shift;
    my $tree = shift;
    my $node = shift;
    my $value = shift;

    my $key  = 'ps_trait';

    $self->throw ("ERROR: ". $node->internal_id. " needs a value for trait $key")
        unless $node->has_tag($key);
    return 0 if $node->get_tag_values($key) ne $value;
    #return $node->branch_length if $node->is_Leaf; # end of recursion
return 0 if $node->is_Leaf; # end of recursion
my $length = 0; foreach my $child ($node->each_Descendent) { $length += $self->_count_leaves($tree, $child, $value) * $child->branch_length + $self->_sum_of_leaf_distances($tree, $child, $value); } return $length;
}
sum_of_leaf_distancesdescriptionprevnextTop
sub sum_of_leaf_distances {
    my $self = shift;
    my $tree = shift;
    my $node = shift  || $tree->get_root_node;

    my $key  = 'ps_trait';
    my $value = $node->get_tag_values($key);

    return $self->_sum_of_leaf_distances($tree, $node, $value);
}
genetic_diversitydescriptionprevnextTop
sub genetic_diversity {
    my $self = shift;
    my $tree = shift;
    my $node = shift  || $tree->get_root_node;

    return $self->sum_of_leaf_distances($tree, $node) /
$self->count_leaves($tree, $node);
}
statratiodescriptionprevnextTop
sub statratio {
    my $self = shift;
    my $tree = shift;
    my $node = shift  || $tree->get_root_node;

    my $div = $self->genetic_diversity($tree, $node);
    return 0 if $div == 0;
    return $node->branch_length / $div;
}
_node_aidescriptionprevnextTop
sub _node_ai {
    my $self = shift;
    my $node = shift;
    my $key = shift;

    my $traits;
    my $leaf_count = 0;
    for my $desc ( $node->get_all_Descendents ) {
        next unless $desc->is_Leaf;
        $leaf_count++;
        $self->throw ("Node ". $desc->id. " needs a value for trait [$key]")
            unless $desc->has_tag($key);
        my $trait = $desc->get_tag_values($key);
        $traits->{$trait}++;
    }
    my $most_common = 0;
    foreach ( keys %$traits) {
        $most_common = $traits->{$_} if $traits->{$_} > $most_common;
    }
    return sprintf "%1.6f", (1 - ($most_common/$leaf_count) ) / (2**($leaf_count-1) );
}
aidescriptionprevnextTop
sub ai {
    my $self = shift;
    my $tree = shift;
    my $key = shift || $self->throw("Trait name is needed");
    my $start_node = shift || $tree->get_root_node;
    return unless $start_node;

    my $sum = 0;
    for my $node ( $start_node->get_all_Descendents ) {
        next if $node->is_Leaf;
        $sum += $self->_node_ai($node, $key);
    }
    return $sum;
}
_node_mcdescriptionprevnextTop
sub _node_mc {
    my $self = shift;
    my $node = shift;
    my $key = shift;

    my $traits;
    my $leaf_count = 0;
    for my $node2 ( $node->get_all_Descendents ) {
        next unless $node2->is_Leaf;
        $leaf_count++;
        my $trait = $node2->get_tag_values($key);
        $traits->{$trait}++;
    }
    return $traits;
}
mcdescriptionprevnextTop
sub mc {
    my $self = shift;
    my $tree = shift;
    my $key = shift || die "Trait name is needed";
    my $start_node = shift || $tree->get_root_node;
    return unless $start_node;

    my $sum; # hashref, keys are trait values
my $keys; # hashref, keys are trait values
foreach my $node ( $start_node->get_all_Descendents ) { next if $node->is_Leaf; my $traits = $self->_node_mc($node, $key); if (scalar keys %$traits == 1) { my ($value) = keys %$traits; no warnings; $sum->{$value} = $traits->{$value} if $sum->{$value} < $traits->{$value}; } else { map { $keys->{$_} = 1 } keys %$traits; } } # check for cases where there are no clusters
foreach my $value (keys %$keys) { $sum->{$value} = 1 unless defined $sum->{$value}; } return $sum; } 1;
}
General documentation
FEEDBACKTop
Mailing ListsTop
User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to
the Bioperl mailing list. Your participation is much appreciated.
  bioperl-l@bioperl.org                  - General discussion
http://bioperl.org/wiki/Mailing_lists - About the mailing lists
Support Top
Please direct usage questions or support issues to the mailing list:
bioperl-l@bioperl.org
rather than to the module maintainer directly. Many experienced and
reponsive experts will be able look at the problem and quickly
address it. Please include a thorough description of the problem
with code and data examples if at all possible.
Reporting BugsTop
Report bugs to the Bioperl bug tracking system to help us keep track
of the bugs and their resolution. Bug reports can be submitted via
the web:
  https://redmine.open-bio.org/projects/bioperl/
AUTHOR - Jason StajichTop
Email jason AT bioperl.org
CONTRIBUTORSTop
Heikki Lehvaslaiho, heikki at bioperl dot org
APPENDIXTop
The rest of the documentation details each of the object methods.
Internal methods are usually preceded with a _
newTop
 Title   : new
Usage : my $obj = Bio::Tree::Statistics->new();
Function: Builds a new Bio::Tree::Statistics object
Returns : Bio::Tree::Statistics
Args :
Tree-Trait statisticsTop
The following methods produce descriptors of trait distribution among
leaf nodes within the trees. They require that a trait has been set
for each leaf node. The tag methods of Bio::Tree::Node are used to
store them as key/value pairs. In this way, one tree can store more
than one trait.
Trees have method add_traits() to set trait values from a file. See the
add_trait() method in Bio::Tree::TreeFunctionsI.