Bio::Tree::DistanceFactory - Construct a tree using distance based methods

  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);

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  https://redmine.open-bio.org/projects/bioperl/

 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'

 Title   : make_tree
 Usage   : my $tree = $disttreefact->make_tree($matrix);
 Function: Build a Tree based on a distance matrix
 Returns : L<Bio::Tree::TreeI>
 Args    : L<Bio::Matrix::MatrixI> object

 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 : L<Bio::Tree::TreeI>
 Args    : L<Bio::Matrix::MatrixI> object

 Title   : _upgma
 Usage   : my $tree = $disttreefact->_upgma($matrix);
 Function: Construct a tree based on alignment using UPGMA
 Returns : L<Bio::Tree::TreeI>
 Args    : L<Bio::Matrix::MatrixI> object

 Title   : method
 Usage   : $obj->method($newval)
 Function: 
 Example : 
 Returns : value of method (a scalar)
 Args    : on set, new value (a scalar or undef, optional)

 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.

#
# BioPerl module for Bio::Tree::DistanceFactory
#
# Please direct questions and support issues to <bioperl-l@bioperl.org> 
#
# Cared for by Jason Stajich <jason-at-bioperl.org>
#
# Copyright Jason Stajich
#
# You may distribute this module under the same terms as perl itself

# POD documentation - main docs before the code

package Bio::Tree::DistanceFactory;
use vars qw($DefaultMethod $Precision);
use strict;

# some defaults
$DefaultMethod = 'UPGMA';
$Precision = 5;

use Bio::Tree::Node;
use Bio::Tree::Tree;

use base qw(Bio::Root::Root);

sub new {
  my($class,@args) = @_;  
  my $self = $class->SUPER::new(@args);

  my ($method) = $self->_rearrange([qw(METHOD)],
				   @args);
  $self->method($method || $DefaultMethod);
  return $self;
}

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;
   }
   
}

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));
}

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

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;
}

sub method{
    my $self = shift;
    return $self->{'_method'} = shift if @_;
    return $self->{'_method'};
}

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;