branch-1-6/Bio/Tree/RandomFactory.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::Tree::RandomFactory
   4 #
   5 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   6 #
   7 # Cared for by Jason Stajich <jason@bioperl.org>
   8 #
   9 # Copyright Jason Stajich
  10 #
  11 # You may distribute this module under the same terms as perl itself
  12
  13 # POD documentation - main docs before the code
  14
  15 =head1 NAME
  16
  17 Bio::Tree::RandomFactory - TreeFactory for generating Random Trees
  18
  19 =head1 SYNOPSIS
  20
  21   use Bio::Tree::RandomFactory
  22   my @taxonnames;
  23   my $factory = Bio::Tree::RandomFactory->new( -taxa => \@taxonnames,
  24                                               -maxcount => 10);
  25
  26   # or for anonymous samples
  27
  28   my $factory = Bio::Tree::RandomFactory->new( -num_taxa => 6,
  29                                               -maxcount => 50);
  30
  31
  32   my $tree = $factory->next_tree;
  33
  34 =head1 DESCRIPTION
  35
  36 Builds a random tree every time next_tree is called or up to -maxcount times.
  37
  38 This module was originally written for Coalescent simulations see
  39 L<Bio::PopGen::Simulation::Coalescent>.  I've left the next_tree
  40 method intact although it is not generating random trees in the
  41 phylogenetic sense.  I would be happy for someone to provide
  42 alternative implementations which can be used here.  As written it
  43 will generate random topologies but the branch lengths are built from
  44 assumptions in the coalescent and are not appropriate for phylogenetic
  45 analyses.
  46
  47 This algorithm is based on the make_tree algorithm from Richard Hudson 1990.
  48
  49 Hudson, R. R. 1990. Gene genealogies and the coalescent
  50        process. Pp. 1-44 in D. Futuyma and J.  Antonovics, eds. Oxford
  51        surveys in evolutionary biology. Vol. 7. Oxford University
  52        Press, New York
  53
  54 Sanderson, M ...
  55
  56 =head1 FEEDBACK
  57
  58 =head2 Mailing Lists
  59
  60 User feedback is an integral part of the evolution of this and other
  61 Bioperl modules. Send your comments and suggestions preferably to
  62 the Bioperl mailing list.  Your participation is much appreciated.
  63
  64   bioperl-l@bioperl.org                  - General discussion
  65   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  66
  67 =head2 Support
  68
  69 Please direct usage questions or support issues to the mailing list:
  70
  71 I<bioperl-l@bioperl.org>
  72
  73 rather than to the module maintainer directly. Many experienced and
  74 reponsive experts will be able look at the problem and quickly
  75 address it. Please include a thorough description of the problem
  76 with code and data examples if at all possible.
  77
  78 =head2 Reporting Bugs
  79
  80 Report bugs to the Bioperl bug tracking system to help us keep track
  81 of the bugs and their resolution. Bug reports can be submitted via
  82 the web:
  83
  84   http://bugzilla.open-bio.org/
  85
  86 =head1 AUTHOR - Jason Stajich
  87
  88 Email jason-AT-bioperl.org
  89
  90 =head1 CONTRIBUTORS
  91
  92 Matthew Hahn, E<lt>matthew.hahn@duke.eduE<gt>
  93 Mike Sanderson
  94
  95 =head1 APPENDIX
  96
  97 The rest of the documentation details each of the object methods.
  98 Internal methods are usually preceded with a _
  99
 100 =cut
 101
 102
 103 # Let the code begin...
 104
 105
 106 package Bio::Tree::RandomFactory;
 107 use vars qw($PRECISION_DIGITS $DefaultNodeType %Defaults);
 108 use strict;
 109
 110 $PRECISION_DIGITS = 3; # Precision for the branchlength
 111 $DefaultNodeType = 'Bio::Tree::Node';
 112 %Defaults = ('YuleRate'          => 1.0, # as set by Sanderson in Rates
 113              'Speciation'        => 1.0, #
 114              'DefaultTreeMethod' => 'yule',
 115              );
 116
 117 use Bio::Tools::RandomDistFunctions;
 118 use Bio::Tree::Tree;
 119
 120 use base qw(Bio::Root::Root Bio::Factory::TreeFactoryI);
 121
 122 =head2 new
 123
 124  Title   : new
 125  Usage   : my $factory = Bio::Tree::RandomFactory->new(-samples => \@samples,
 126                                                       -maxcount=> $N);
 127  Function: Initializes a Bio::Tree::RandomFactory object
 128  Returns : Bio::Tree::RandomFactory
 129  Args    : -nodetype => Type of Nodes to create [default Bio::Tree::Node]
 130            -maxcount => [optional] Maximum num trees to create
 131            -randtype => Type of random trees so far support
 132                - yule/backward_yule/BY [default]
 133                - forward_yule/FY
 134                - birthdeath_forward/BDF
 135                - birthdeath_backwards/BDB
 136
 137
 138           ONE of the following must be specified
 139            -taxa     => $arrayref of taxa names
 140            -num_taxa => integer indicating number of taxa in the tree
 141
 142 =cut
 143
 144 sub new{
 145    my ($class,@args) = @_;
 146    my $self = $class->SUPER::new(@args);
 147
 148    $self->{'_treecounter'} = 0;
 149    $self->{'_maxcount'} = 0;
 150    my ($nodetype,$randtype,
 151        $maxcount, $samps,$samplesize,
 152        $taxa, $num_taxa) = $self->_rearrange([qw(NODETYPE
 153                                                  RANDTYPE
 154                                                  MAXCOUNT
 155                                                  SAMPLES
 156                                                  SAMPLE_SIZE
 157                                                  TAXA
 158                                                  NUM_TAXA)],
 159                                              @args);
 160    my @taxa;
 161    $nodetype ||= $DefaultNodeType;
 162    $self->nodetype($nodetype);
 163    $taxa = $samps if defined $samps && ! defined $taxa;
 164    $num_taxa = $samplesize if $samplesize && ! $num_taxa;
 165    if( ! defined $taxa ) {
 166        if( ! defined $num_taxa || $num_taxa <= 0 ) {
 167            $self->throw("Must specify a valid num_taxa if parameter -TAXA is not specified");
 168        }
 169        foreach ( 1..$num_taxa ) { push @taxa, "Taxon$_"; }
 170    } else {
 171        if( ref($taxa) !~ /ARRAY/i ) {
 172            $self->throw("Must specify a valid ARRAY reference to the parameter -TAXA, did you forget a leading '\\'? for $taxa");
 173        }
 174        @taxa = @$taxa;
 175    }
 176
 177    $self->taxa(\@taxa);
 178    defined $maxcount && $self->maxcount($maxcount);
 179    $self->{'_count'} = 0;
 180    return $self;
 181 }
 182
 183 =head2 next_tree
 184
 185  Title   : next_tree
 186  Usage   : my $tree = $factory->next_tree
 187  Function: Returns a random tree based on the initialized number of nodes
 188            NOTE: if maxcount is not specified on initialization or
 189                  set to a valid integer, subsequent calls to next_tree will
 190                  continue to return random trees and never return undef
 191
 192  Returns : Bio::Tree::TreeI object
 193  Args    : none
 194
 195 =cut
 196
 197
 198 sub next_tree{
 199    my ($self,%options) = @_;
 200    return if $self->maxcount &&
 201        $self->{'_count'}++ >= $self->maxcount;
 202    my $rand_type = $options{'randtype'} || $self->random_tree_method;
 203    my $nodetype = $self->nodetype;
 204    my $treearray;
 205
 206    if( $rand_type =~ /(birthdeath_forward|birth|BDF)/i ) {
 207
 208    } elsif ( $rand_type =~ /(birthdeath_backward|BDB)/i ) {
 209        $treearray = $self->rand_birthdeath_backwards_tree;
 210    } elsif( $rand_type =~ /(BY|backwards_yule)/i ||
 211             $rand_type =~ /^yule/i ) {
 212        my $speciation = $options{'speciation'}; # can be undef
 213        $treearray = $self->rand_yule_c_tree($speciation);
 214    } else {
 215        $self->warn("unrecognized random type $rand_type");
 216    }
 217
 218    my @nodes = ();
 219    foreach my $n ( @$treearray ) {
 220        for my $k ( qw(desc1 desc2) ) {
 221            next unless defined $n->{$k};
 222            push @{$n->{'descendents'}}, $nodes[$n->{$k}];
 223        }
 224        push @nodes,
 225        $nodetype->new(-id            => $n->{'nodenum'},
 226                       -branch_length => $n->{'time'},
 227                       -descendents   => $n->{'descendents'},
 228                       );
 229    }
 230    my $T = Bio::Tree::Tree->new(-root => pop @nodes );
 231    return $T;
 232 }
 233
 234
 235 =head2 maxcount
 236
 237  Title   : maxcount
 238  Usage   : $obj->maxcount($newval)
 239  Function:
 240  Returns : Maxcount value
 241  Args    : newvalue (optional)
 242
 243
 244 =cut
 245
 246 sub maxcount{
 247    my ($self,$value) = @_;
 248    if( defined $value) {
 249        if( $value =~ /^(\d+)/ ) {
 250            $self->{'_maxcount'} = $1;
 251        } else {
 252            $self->warn("Must specify a valid Positive integer to maxcount");
 253            $self->{'_maxcount'} = 0;
 254        }
 255   }
 256    return $self->{'_maxcount'};
 257 }
 258
 259
 260 =head2 reset_tree_count
 261
 262  Title   : reset_tree_count
 263  Usage   : $factory->reset_tree_count;
 264  Function: Reset the tree counter
 265  Returns : none
 266  Args    : none
 267
 268
 269 =cut
 270
 271 sub reset_count{
 272     shift->{'_count'} = 0;
 273 }
 274
 275 =head2 taxa
 276
 277  Title   : taxa
 278  Usage   : $obj->taxa($newval)
 279  Function: Set the leaf node names
 280  Returns : value of taxa
 281  Args    : Arrayref of Taxon names
 282
 283
 284 =cut
 285
 286 sub taxa {
 287     my ($self,$value) = @_;
 288     if( defined $value) {
 289         if( ref($value) !~ /ARRAY/i ) {
 290             $self->warn("Must specify a valid array ref to the method 'taxa'");
 291             $value = [];
 292         }
 293         $self->{'_taxa'} = $value;
 294         $self->{'_num_taxa'} = scalar @$value;
 295     }
 296     return $self->{'_taxa'};
 297
 298 }
 299
 300 =head2 num_taxa
 301
 302  Title   : num_taxa
 303  Usage   : $obj->num_taxa($newval)
 304  Function: Get the number of Taxa
 305  Returns : value of num_taxa
 306  Args    : none
 307
 308
 309 =cut
 310
 311 sub num_taxa {
 312     my ($self) = @_;
 313     return  $self->{'_num_taxa'};
 314 }
 315
 316 # alias old methods
 317 *num_samples = \&num_taxa;
 318 *samples = \&taxa;
 319
 320 =head2 random
 321
 322  Title   : random
 323  Usage   : my $rfloat = $node->random($size)
 324  Function: Generates a random number between 0 and $size
 325            This is abstracted so that someone can override and provide their
 326            own special RNG.  This is expected to be a uniform RNG.
 327  Returns : Floating point random
 328  Args    : $maximum size for random number (defaults to 1)
 329
 330
 331 =cut
 332
 333 sub random{
 334    my ($self,$max) = @_;
 335    return rand($max);
 336 }
 337
 338
 339 =head2 random_tree_method
 340
 341  Title   : random_tree_method
 342  Usage   : $obj->random_tree_method($newval)
 343  Function:
 344  Example :
 345  Returns : value of random_tree_method (a scalar)
 346  Args    : on set, new value (a scalar or undef, optional)
 347
 348
 349 =cut
 350
 351 sub random_tree_method{
 352     my $self = shift;
 353
 354     return $self->{'random_tree_method'} = shift if @_;
 355     return $self->{'random_tree_method'} || $Defaults{'DefaultTreeMethod'};
 356 }
 357
 358 =head2 nodetype
 359
 360  Title   : nodetype
 361  Usage   : $obj->nodetype($newval)
 362  Function:
 363  Example :
 364  Returns : value of nodetype (a scalar)
 365  Args    : on set, new value (a scalar or undef, optional)
 366
 367
 368 =cut
 369
 370 sub nodetype{
 371    my ($self,$value) = @_;
 372    if( defined $value) {
 373        eval "require $value";
 374        if( $@ ) { $self->throw("$@: Unrecognized Node type for ".ref($self).
 375                                "'$value'");}
 376
 377        my $a = bless {},$value;
 378        unless( $a->isa('Bio::Tree::NodeI')  ) {
 379            $self->throw("Must provide a valid Bio::Tree::NodeI or child class to SeqFactory Not $value");
 380        }
 381       $self->{'nodetype'} = $value;
 382     }
 383     return $self->{'nodetype'};
 384 }
 385
 386
 387 # The assignment of times are based on Mike Sanderson's r8s code
 388 # The topology assignment code is based on Richard Hudson's
 389 # make_trees
 390
 391
 392 sub rand_yule_c_tree {
 393     my ($self,$speciation) = @_;
 394     $speciation ||= $Defaults{'Speciation'};
 395     my $n_taxa = $self->num_taxa;
 396     my $taxa = $self->taxa || [];
 397     my $nodetype = $self->nodetype;
 398
 399     my $randfuncs = Bio::Tools::RandomDistFunctions->new();
 400     my $rate = $Defaults{'YuleRate'};
 401     my (@tree,@list,@times,$i,$in);
 402     my $max = 2 * $n_taxa - 1;
 403     for($in=0;$in < $max; $in++ ) {
 404         push @tree, { 'nodenum' => "Node$in" };
 405     }
 406     # setup leaf nodes
 407     for($in=0;$in < $n_taxa;$in++)  {
 408         $tree[$in]->{'time'} = 0;
 409         $tree[$in]->{'desc1'} = undef;
 410         $tree[$in]->{'desc2'} = undef;
 411         if( my $r = $taxa->[$in] ) {
 412             $tree[$in]->{'nodenum'} = $r;
 413         }
 414         push @list, $in;
 415     }
 416
 417     for( $i = 0; $i < $n_taxa - 1; $i++ ) {
 418         # draw random interval times
 419         push @times, $randfuncs->rand_birth_distribution($speciation);
 420     }
 421     # sort smallest to largest
 422     @times = sort {$a <=> $b} @times;
 423      # topology generation
 424     for ($in = $n_taxa; $in > 1; $in-- ) {
 425         my $time = shift @times;
 426
 427         my $pick = int $self->random($in);
 428         my $nodeindex = $list[$pick];
 429         $tree[$list[$pick]]->{'time'} = $time;
 430         my $swap = 2 * $n_taxa - $in;
 431         $tree[$swap]->{'desc1'} = $nodeindex;
 432         $list[$pick] = $list[$in-1];
 433
 434         $pick = int rand($in - 1);
 435         $nodeindex = $list[$pick];
 436         $tree[$list[$pick]]->{'time'} = $time;
 437         $tree[$swap]->{'desc2'} = $nodeindex;
 438         $list[$pick] = $swap;
 439     }
 440     $tree[-1]->{'time'} = shift @times;
 441     return \@tree;
 442 }
 443
 444
 445
 446 sub rand_birthdeath_backwards_tree {
 447     my ($self) = @_;
 448     my $n_taxa = $self->num_taxa;
 449     my $taxa = $self->taxa || [];
 450
 451     my $randfuncs = Bio::Tools::RandomDistFunctions->new();
 452     my $rate = $Defaults{'YuleRate'};
 453     my (@tree,@list,@times,$i,$in);
 454     my $max = 2 * $n_taxa - 1;
 455     for($in=0;$in < $max; $in++ ) {
 456         push @tree, { 'nodenum' => "Node$in" };
 457     }
 458     # setup leaf nodes
 459     for($in=0;$in < $n_taxa;$in++)  {
 460         $tree[$in]->{'time'} = 0;
 461         $tree[$in]->{'desc1'} = undef;
 462         $tree[$in]->{'desc2'} = undef;
 463         if( my $r = $taxa->[$in] ) {
 464             # deal with pre-labeled nodes
 465             $tree[$in]->{'nodenum'} = $r;
 466         }
 467         push @list, $in;
 468     }
 469     my ($time) = (0);
 470
 471      # topology generation
 472     for ($in = $n_taxa; $in > 1; $in-- ) {
 473         my $pick = int $self->random($in);
 474         my $nodeindex = $list[$pick];
 475         my $swap = 2 * $n_taxa - $in;
 476         $time += $randfuncs->rand_geometric_distribution($n_taxa * $rate);;
 477         $tree[$list[$pick]]->{'time'} = $time;
 478         $tree[$swap]->{'desc1'} = $nodeindex;
 479         $list[$pick] = $list[$in-1];
 480
 481         $pick = int rand($in - 1);
 482         $nodeindex = $list[$pick];
 483         $tree[$list[$pick]]->{'time'} = $time;
 484         $tree[$swap]->{'desc2'} = $nodeindex;
 485         $list[$pick] = $swap;
 486     }
 487     my $root = $tree[-1];
 488     $time += $randfuncs->rand_geometric_distribution($n_taxa * $rate);;
 489     $root->{'time'} = $time;
 490
 491     # Normalize times by the root node...
 492     for my $node ( @tree ) {
 493         $node->{'time'} /= $root->{'time'};
 494     }
 495     return \@tree;
 496 }
 497
 498
 499 # The assignment of times are based on Mike Sanderson's r8s code
 500 # The topology assignment code is based on Richard Hudson's
 501 # make_trees
 502
 503 sub rand_birth_death_tree {
 504 # Still need to finish
 505 #     my ($self,$spec_rate,$extinct_rate,$char_rate) = @_;
 506 #     my $n_taxa =  $self->num_taxa;
 507 #     my $dt = 0.1 / $n_taxa;
 508 #     my @tree;
 509 #     my $max = 3 * $n_taxa - 1;
 510 #     # setup leaf nodes
 511
 512 #     for($in=0;$in < $size;$in++)  {
 513 #       push @tree, { 'nodenum' => $taxa->[$in] || "Node$in",
 514 #                     'time'    => 0,
 515 #                     'desc1'   => undef,
 516 #                     'desc2'   => undef,
 517 #                 };
 518 #     }
 519 #     my $time = $dt;
 520 #     my $idx = 0;
 521 #     while( $n_taxa > 1 ) {
 522 #       if ( event($dt * $spec_rate, $n_taxa) ) {
 523 #           my $pick = int $self->random($n_taxa);
 524 #           my $pick2 = int $self->random($n_taxa);
 525 #           while( $pick2 == $pick ) {
 526 #               $pick2 = int $self->random($n_taxa);
 527 #           }
 528             # to finish....
 529
 530 #           $tree[$swap]->{'desc1'} = $nodeindex;
 531 #       }
 532 #     }
 533
 534
 535
 536 #       $list[$pick] = $list[$in-1];
 537
 538 #       $pick = int rand($in - 1);
 539 #       $nodeindex = $list[$pick];
 540 #       $tree[$swap]->{'desc2'} = $nodeindex;
 541 #       $list[$pick] = $swap;
 542 #       $tree[$swap]->{'time'} = $times[$ix++];
 543 #     }
 544 }
 545
 546
 547 1;