Bio::Tree DistanceFactory
SummaryIncluded librariesPackage variablesSynopsisDescriptionGeneral documentationMethods
Toolbar
WebCvs
Summary
Bio::Tree::DistanceFactory - Construct a tree using distance based methods
Package variables
No package variables defined.
Included modules
Bio::Tree::Node
Bio::Tree::Tree
Inherit
Bio::Root::Root
Synopsis
  use Bio::Tree::DistanceFactory;
use Bio::AlignIO;
use Bio::Align::DNAStatistics;
my $tfactory = Bio::Tree::DistanceFactory->new(-method => "NJ");
my $stats = Bio::Align::DNAStatistics->new();
my $alnin = Bio::AlignIO->new(-format => 'clustalw', -file => 'file.aln'); my $aln = $alnin->next_aln; # Of course matrix can come from a different place # like PHYLIP if you prefer, Bio::Matrix::IO should be able # to parse many things my $jcmatrix = $stats->distance(-align => $aln, -method => 'Jukes-Cantor'); my $tree = $tfactory->make_tree($jcmatrix);
Description
This is a factory which will construct a phylogenetic tree based on
the pairwise sequence distances for a set of sequences. Currently
UPGMA (Sokal and Michener 1958) and NJ (Saitou and Nei 1987) tree
construction methods are implemented.
Methods
newDescriptionCode
make_treeDescriptionCode
_njDescriptionCode
_upgmaDescriptionCode
_upgma_distance
No description
Code
methodDescriptionCode
check_additivityDescriptionCode
Methods description
newcode    nextTop
 Title   : new
Usage : my $obj = Bio::Tree::DistanceFactory->new();
Function: Builds a new Bio::Tree::DistanceFactory object
Returns : an instance of Bio::Tree::DistanceFactory
Args : -method => 'NJ' or 'UPGMA'
make_treecodeprevnextTop
 Title   : make_tree
Usage : my $tree = $disttreefact->make_tree($matrix);
Function: Build a Tree based on a distance matrix
Returns : Bio::Tree::TreeI
Args : Bio::Matrix::MatrixI object
_njcodeprevnextTop
 Title   : _nj
Usage : my $tree = $disttreefact->_nj($matrix);
Function: Construct a tree based on distance matrix using the
Neighbor Joining algorithm (Saitou and Nei, 1987)
Implementation based on Kevin Howe's Quicktree implementation
and uses his tricks (some based on Bill Bruno's work) to eliminate
negative branch lengths
Returns : Bio::Tree::TreeI
Args : Bio::Matrix::MatrixI object
_upgmacodeprevnextTop
 Title   : _upgma
Usage : my $tree = $disttreefact->_upgma($matrix);
Function: Construct a tree based on alignment using UPGMA
Returns : Bio::Tree::TreeI
Args : Bio::Matrix::MatrixI object
methodcodeprevnextTop
 Title   : method
Usage : $obj->method($newval)
Function:
Example :
Returns : value of method (a scalar)
Args : on set, new value (a scalar or undef, optional)
check_additivitycodeprevnextTop
 Title     : check_additivity
Usage : if( $distance->check_additivity($matrix) ) {
}
Function : See if matrix obeys additivity principal
Returns : boolean
Args : Bio::Matrix::MatrixI
References: Based on a Java implementation by
Peter Sestoft, sestoft@dina.kvl.dk 1999-12-07 version 0.3
http://www.dina.kvl.dk/~sestoft/bsa.html
which in turn is based on algorithms described in
R. Durbin, S. Eddy, A. Krogh, G. Mitchison.
Biological Sequence Analysis CUP 1998, Chapter 7.
Methods code
newdescriptionprevnextTop
sub new {
  my($class,@args) = @_;  
  my $self = $class->SUPER::new(@args);

  my ($method) = $self->_rearrange([qw(METHOD)],
				   @args);
  $self->method($method || $DefaultMethod);
  return $self;
}
make_treedescriptionprevnextTop
sub make_tree {
   my ($self,$matrix) = @_;
   if( ! defined $matrix || !ref($matrix) || 
       ! $matrix->isa('Bio::Matrix::MatrixI') ) {
       $self->warn("Need to provide a valid Bio::Matrix::MatrixI object to make_tree");
       return;
   }

   my $method = uc ($self->method);
   if( $method =~ /NJ/i ) {
       return $self->_nj($matrix);
   } elsif( $method =~ /UPGMA/i ) {
       return $self->_upgma($matrix);
   } else { 
       $self->warn("Unknown tree construction method '$method'.  Cannot run.");
       return;
   }
}
_njdescriptionprevnextTop
sub _nj {
   my ($self,$distmat) = @_;

   # we assume type checking of $aln has already been done
# client shouldn't be calling this directly anyways, using the
# make_tree method is preferred
# so that we can trim the number of digits shown as the branch length
my $precisionstr = "%.$Precision"."f"; my @names = $distmat->column_names; my $N = scalar @names; my ($i,$j,$m,@nodes,$mat,@r); my $L = $N; if( $N < 2 ) { $self->warn("Can only perform NJ treebuilding on sets of 2 or more species\n"); return; } elsif( $N == 2 ) { $i = 0; my $d = sprintf($precisionstr, $distmat->get_entry($names[0],$names[1]) / 2);
my $root = Bio::Tree::Node->new(); for my $nm ( @names ) { $root->add_Descendents( Bio::Tree::Node->new(-id => $nm, -branch_length => $d)); } return Bio::Tree::Tree(-root => $root); } my $c = 0; for ( $i = 0; $i < $N; $i++ ) { push @nodes, Bio::Tree::Node->new(-id => $names[$i]); my $ri = 0; for( $j = 0; $j < $N; $j++ ) { $mat->[$i][$j] = $distmat->get_entry($names[$i],$names[$j]); $ri += $mat->[$i][$j]; } $r[$i] = $ri / ($L -2);
} for( my $nodecount = 0; $nodecount < $N-3; $nodecount++) { my ($mini,$minj,$min); for($i = 0; $i < $N; $i++ ) { next unless defined $nodes[$i]; for( $j = 0; $j < $i; $j++ ) { next unless defined $nodes[$j]; my $dist = $mat->[$i][$j] - ($r[$i] + $r[$j]); if( ! defined $min || $dist <= $min) { ($mini,$minj,$min) = ($i,$j,$dist); } } } my $dij = $mat->[$mini][$minj]; my $dist_i = ($dij + $r[$mini] - $r[$minj]) / 2;
my $dist_j = $dij - $dist_i; # deal with negative branch lengths
# per code in K.Howe's quicktree
if( $dist_i < 0 ) { $dist_i = 0; $dist_j = $dij; $dist_j = 0 if( $dist_j < 0 ); } elsif( $dist_j < 0 ) { $dist_j = 0; $dist_i = $dij; $dist_i = 0 if( $dist_i < 0 ); } $nodes[$mini]->branch_length(sprintf($precisionstr,$dist_i)); $nodes[$minj]->branch_length(sprintf($precisionstr,$dist_j)); my $newnode = Bio::Tree::Node->new(-descendents => [ $nodes[$mini], $nodes[$minj] ]); $nodes[$mini] = $newnode; delete $nodes[$minj]; # update the distance matrix
$r[$mini] = 0; my ($dmi,$dmj); for( $m = 0; $m < $N; $m++ ) { next unless defined $nodes[$m]; if( $m != $mini ) { $dmj = $mat->[$m][$minj]; my ($row,$col); ($row,$col) = ($m,$mini); $dmi = $mat->[$row][$col]; # from K.Howe's notes in quicktree
# we can actually adjust r[m] here, by using the form:
# rm = ((rm * numseqs) - dmi - dmj + dmk) / (numseqs-1)
# Note: in Bill Bruno's method for negative branch
# elimination, then if either dist_i is positive and
# dist_j is 0, or dist_i is zero and dist_j is positive
# (after adjustment) then the matrix entry is formed
# from the distance to the node in question (m) to the
# node with the zero branch length (whichever it was).
# I think my code already has the same effect; this is
# certainly true if dij is equal to dist_i + dist_j,
# which it should have been fixed to
my $dmk = $mat->[$row][$col] = $mat->[$col][$row] = ($dmi + $dmj - $dij) / 2;
# If we don't want to try and correct negative brlens
# this is essentially what is in Edddy et al, BSA book.
# $r[$m] = (($r[$m] * $L) - $dmi - $dmj + $dmk) / ($L-1);
#
$r[$m] = (($r[$m] * ($L - 2)) - $dmi - $dmj + $mat->[$row][$col]) / ( $L - 3);
$r[$mini] += $dmk; } } $L--; $r[$mini] /= $L - 2;
} # should be 3 nodes left
my (@leftovernodes,@leftovers); for( my $k = 0; $k < $N; $k++ ) { if( defined $nodes[$k] ) { push @leftovers, $k; push @leftovernodes, $nodes[$k]; } } my ($l_0,$l_1,$l_2) = @leftovers; my $dist_i = ( $mat->[$l_1][$l_0] + $mat->[$l_2][$l_0] - $mat->[$l_2][$l_1] ) / 2;
my $dist_j = ( $mat->[$l_1][$l_0] - $dist_i); my $dist_k = ( $mat->[$l_2][$l_0] - $dist_i); # This is Kev's code to get rid of negative branch lengths
if( $dist_i < 0 ) { $dist_i = 0; $dist_j = $mat->[$l_1][$l_0]; $dist_k = $mat->[$l_2][$l_0]; if( $dist_j < 0 ) { $dist_j = 0; $dist_k = ( $mat->[$l_2][$l_0] + $mat->[$l_2][$l_1] ) / 2;
$dist_k = 0 if( $dist_k < 0 ); } elsif( $dist_k < 0 ) { $dist_k = 0; $dist_j = ($mat->[$l_1][$l_0] + $mat->[$l_2][$l_1]) / 2;
$dist_j = 0 if( $dist_j < 0 ); } } elsif( $dist_j < 0 ) { $dist_j = 0; $dist_i = $mat->[$l_1][$l_0]; $dist_k = $mat->[$l_2][$l_1]; if( $dist_i < 0 ) { $dist_i = 0; $dist_k = ( $mat->[$l_2][$l_0] + $mat->[$l_2][$l_1]) / 2;
$dist_k = 0 if( $dist_k < 0 ); } elsif( $dist_k < 0 ) { $dist_k = 0; $dist_i = ( $mat->[$l_1][$l_0] + $mat->[$l_2][$l_0]) / 2;
$dist_i = 0 if( $dist_i < 0 ); } } elsif( $dist_k < 0 ) { $dist_k = 0; $dist_i = $mat->[$l_2][$l_0]; $dist_j = $mat->[$l_2][$l_1]; if( $dist_i < 0 ) { $dist_i = 0; $dist_j = ( $mat->[$l_1][$l_0] + $mat->[$l_2][$l_1] ) / 2;
$dist_j = 0 if $dist_j < 0; } elsif( $dist_j < 0 ) { $dist_j = 0; $dist_i = ($mat->[$l_1][$l_0] + $mat->[$l_2][$l_0]) / 2;
$dist_i = 0 if $dist_i < 0; } } $leftovernodes[0]->branch_length(sprintf($precisionstr,$dist_i)); $leftovernodes[1]->branch_length(sprintf($precisionstr,$dist_j)); $leftovernodes[2]->branch_length(sprintf($precisionstr,$dist_k)); Bio::Tree::Tree->new(-root => Bio::Tree::Node->new (-descendents =>\@ leftovernodes));
}
_upgmadescriptionprevnextTop
sub _upgma {
   my ($self,$distmat) = @_;
   # we assume type checking of $matrix has already been done
# client shouldn't be calling this directly anyways, using the
# make_tree method is preferred
# algorithm, from Eddy, Durbin, Krogh, Mitchison, 1998
# originally by Sokal and Michener 1956
my $precisionstr = "%.$Precision"."f"; my ($i,$j,$x,$y,@dmat,@orig,@nodes); my @names = $distmat->column_names; my $c = 0; my @clusters = map { my $r = { 'id' => $c, 'height' => 0, 'contains' => [$c], }; $c++; $r; } @names; my $K = scalar @clusters; my (@mins,$min); for ( $i = 0; $i < $K; $i++ ) { for( $j = $i+1; $j < $K; $j++ ) { my $d = $distmat->get_entry($names[$i],$names[$j]); # get Min here on first time around, save 1 cycle
$dmat[$j][$i] = $dmat[$i][$j] = $d; $orig[$i][$j] = $orig[$j][$i] = $d; if ( ! defined $min || $d <= $min ) { if( defined $min && $min == $d ) { push @mins, [$i,$j]; } else { @mins = [$i,$j]; $min = $d; } } } } # distance between each cluster is avg distance
# between pairs of sequences from each cluster
while( $K > 1 ) { # fencepost - we already have found the $min
# so very first time loop is executed we can skip checking
unless( defined $min ) { for($i = 0; $i < $K; $i++ ) { for( $j = $i+1; $j < $K; $j++ ) { my $dij = $dmat[$i][$j]; if( ! defined $min || $dij <= $min) { if( defined $min && $min == $dij ) { push @mins, [$i,$j]; } else { @mins = [ $i,$j ]; $min = $dij; } } } } } # randomly break ties
($x,$y) = @{ $mins[int(rand(scalar @mins))] }; # now we are going to join clusters x and y, make a new cluster
my $node = Bio::Tree::Node->new(); my @subids; for my $cid ( $x,$y ) { my $nid = $clusters[$cid]->{'id'}; if( ! defined $nodes[$nid] ) { $nodes[$nid] = Bio::Tree::Node->new(-id => $names[$nid]); } $nodes[$nid]->branch_length (sprintf($precisionstr,$min/2 - $clusters[$cid]->{'height'}));
$node->add_Descendent($nodes[$nid]); push @subids, @{ $clusters[$cid]->{'contains'} }; } my $cluster = { 'id' => $c++, 'height' => $min / 2,
'contains' => [
@subids],
};
$K--; # we are going to drop the last node so go ahead and decrement K
$nodes[$cluster->{'id'}] = $node; if ( $y != $K ) { $clusters[$y] = $clusters[$K]; $dmat[$y] = $dmat[$K]; for ( $i = 0; $i < $K; $i++ ) { $dmat[$i][$y] = $dmat[$y][$i]; } } delete $clusters[$K]; $clusters[$x] = $cluster; # now recalculate @dmat
for( $i = 0; $i < $K; $i++ ) { if( $i != $x) { $dmat[$i][$x] = $dmat[$x][$i] = &_upgma_distance($clusters[$i],$clusters[$x],\@orig); } else { $dmat[$i][$i] = 0; } } # reset so next loop iteration
# we will find minimum distance
@mins = (); $min = undef; } Bio::Tree::Tree->new(-root => $nodes[-1]); } # calculate avg distance between clusters - be they
# single sequences or the combination of multiple seqences
# $cluster_i and $cluster_j are the clusters to operate on
# and $distances is a matrix (arrayref of arrayrefs) of pairwise
# differences indexed on the sequence ids -
# so $distances->[0][1] is the distance between sequences 0 and 1
}
_upgma_distancedescriptionprevnextTop
sub _upgma_distance {
     my ($cluster_i, $cluster_j, $distances) = @_;
    my $ilen = scalar @{ $cluster_i->{'contains'} };
    my $jlen = scalar @{ $cluster_j->{'contains'} };
    my ($d,$count);
    for( my $i = 0; $i < $ilen; $i++ ) {
	my $i_id = $cluster_i->{'contains'}->[$i];
	for( my $j = 0; $j < $jlen; $j++) {	    
	    my $j_id = $cluster_j->{'contains'}->[$j];
	    if( ! defined $distances->[$i_id][$j_id] ) {
		warn("no value for $i_id $j_id\n");
	    } else { 
		$d += $distances->[$i_id][$j_id];
	    }
	    $count++;
	}
    }
    return $d / $count;
}
methoddescriptionprevnextTop
sub method {
    my $self = shift;
    return $self->{'_method'} = shift if @_;
    return $self->{'_method'};
}
check_additivitydescriptionprevnextTop
sub check_additivity {
   my ($self,$matrix) = @_;
   my @names = $matrix->column_names;
   my $len = scalar @names;
   return unless $len >= 4;
   # look at all sets of 4
for( my $i = 0; $i < $len; $i++ ) { for( my $j = $i+1; $j< $len; $j++) { for( my $k = $j+1; $k < $len; $k ++ ) { for( my $m = $k +1; $m < $len; $m++ ) { my $DijDkm = $matrix->get_entry($names[$i],$names[$j]) + $matrix->get_entry($names[$k],$names[$m]); my $DikDjm = $matrix->get_entry($names[$i],$names[$k]) + $matrix->get_entry($names[$j],$names[$m]); my $DimDjk = $matrix->get_entry($names[$i],$names[$m]) + $matrix->get_entry($names[$j],$names[$k]); if( !( ( $DijDkm == $DikDjm && $DijDkm >= $DimDjk) || ( $DijDkm == $DimDjk && $DijDkm >= $DikDjm) || ( $DikDjm == $DimDjk && $DikDjm >= $DijDkm) )) { return 0; } } } } } return 1; } 1;
}
General documentation
REFERENCESTop
Eddy SR, Durbin R, Krogh A, Mitchison G, (1998) "Biological Sequence Analysis",
Cambridge Univ Press, Cambridge, UK.
Howe K, Bateman A, Durbin R, (2002) "QuickTree: building huge
Neighbour-Joining trees of protein sequences." Bioinformatics
18(11):1546-1547.
Saitou N and Nei M, (1987) "The neighbor-joining method: a new method
for reconstructing phylogenetic trees." Mol Biol Evol 4(4):406-25.
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 the web:
  https://redmine.open-bio.org/projects/bioperl/
AUTHOR - Jason StajichTop
Email jason-at-bioperl.org
APPENDIXTop
The rest of the documentation details each of the object methods.
Internal methods are usually preceded with a _