1package Bio::Phylo::Models::Substitution::Dna; 2use Bio::Phylo::Util::CONSTANT qw'/looks_like/ :objecttypes'; 3use Bio::Phylo::Util::Exceptions qw'throw'; 4use Bio::Phylo::IO qw(parse unparse); 5use Bio::Phylo::Util::Logger':levels'; 6use File::Temp qw(tempfile cleanup); 7 8use strict; 9use warnings; 10 11sub _INDEX_OF_ { { A => 0, C => 1, G => 2, T => 3 } } 12sub _BASE_AT_ { [qw(A C G T)] } 13 14my $logger = Bio::Phylo::Util::Logger->new; 15 16=head1 NAME 17 18Bio::Phylo::Models::Substitution::Dna - DNA substitution model 19 20=head1 SYNOPSIS 21 22 use Bio::Phylo::Models::Substitution::Dna; 23 24 # create a DNA substitution model from scratch 25 my $model = Bio::Phylo::Models::Substitution::Dna->new( 26 '-type' => 'GTR', 27 '-pi' => [ 0.23, 0.27, 0.24, 0.26 ], 28 '-kappa' => 2, 29 '-alpha' => 0.9, 30 '-pinvar' => 0.5, 31 '-ncat' => 6, 32 '-median' => 1, 33 '-rate' => [ 34 [ 0.23, 0.23, 0.23, 0.23 ], 35 [ 0.23, 0.26, 0.26, 0.26 ], 36 [ 0.27, 0.26, 0.26, 0.26 ], 37 [ 0.24, 0.26, 0.26, 0.26 ] 38 ] 39 ); 40 41 # get substitution rate from A to C 42 my $rate = $model->get_rate('A', 'C'); 43 44 # get model representation that can be used by Garli 45 my $modelstr = $model->to_string( '-format' => 'garli' ) 46 47=head1 DESCRIPTION 48 49This is a superclass for models of DNA evolution. Classes that inherit from this 50class provide methods for retreiving general parameters such as substitution rates 51or the number of states as well as model-specific parameters. Currently most of the 52popular models are implemented. The static function C<modeltest> determines the 53substitution model from a L<Bio::Phylo::Matrices::Matrix> object and returns the 54appropriate instance of the subclass. This class also provides serialization 55of a model to standard phylogenetics file formats. 56 57=head1 METHODS 58 59=head2 CONSTRUCTOR 60 61=over 62 63=item new 64 65Dna model constructor. 66 67 Type : Constructor 68 Title : new 69 Usage : my $model = Bio::Phylo::Models::Substitution::Dna->new(%args); 70 Function: Instantiates a Bio::Phylo::Models::Substitution::Dna object. 71 Returns : A Bio::Phylo::Models::Substitution::Dna object. 72 Args : Optional: 73 -type => type of model, one of GTR, F81, HKY85, JC69, K80 74 -pi => base frequencies of bases A, C, G, T 75 -kappa => ratio transitions/transversions 76 -alpha => shape parameter (for models of GTR family) 77 -mu => overall mutation rate 78 -pinvar => proportion of invariant sites 79 -ncat => number of distinct rate categories 80 -median => median for gamma-modeled rate categories 81 -rate => Array of Arrays (4x4) giving substitution rates betwen A, C, T, G 82 -catweights => weights for rate categories 83=cut 84 85sub new { 86 my $class = shift; 87 my %args = looks_like_hash @_; 88 $class .= '::' . uc $args{'-type'} if $args{'-type'}; 89 delete $args{'-type'}; 90 my $self = {}; 91 bless $self, looks_like_class $class; 92 while ( my ( $key, $value ) = each %args ) { 93 $key =~ s/^-/set_/; 94 $self->$key($value); 95 } 96 return $self; 97} 98 99=item get_catrates 100 101Getter for rate categories, implemented by child classes. 102 103 Type : method 104 Title : get_catrates 105 Usage : $model->get_catrates; 106 Function: Getter for rate categories. 107 Returns : scalar or array 108 Args : None. 109 110=cut 111 112sub get_catrates { 113 throw 'NotImplemented' => 'FIXME'; 114} 115 116=item get_nst 117 118Getter for number of transition rate parameters. 119 120 Type : method 121 Title : get_nst 122 Usage : $model->get_nst; 123 Function: Getter for number of transition rate parameters. 124 Returns : scalar 125 Args : None. 126 127=cut 128 129sub get_nst { 6 } 130 131=item get_rate 132 133Getter for substitution rate. If bases are given as arguments, 134returns corresponding rate. If no arguments given, returns rate matrix or 135overall rate, dependent on model. 136 137 Type : method 138 Title : get_rate 139 Usage : $model->get_rate('A', 'C'); 140 Function: Getter for transition rate between nucleotides. 141 Returns : scalar or array 142 Args : Optional: 143 base1: scalar 144 base2: scalar 145=cut 146 147sub get_rate { 148 my $self = shift; 149 if (@_) { 150 my $src = _INDEX_OF_()->{ uc shift }; 151 my $target = _INDEX_OF_()->{ uc shift }; 152 $self->{'_rate'} = [] if not $self->{'_rate'}; 153 if ( not $self->{'_rate'}->[$src] ) { 154 $self->{'_rate'}->[$src] = []; 155 } 156 return $self->{'_rate'}->[$src]->[$target]; 157 } 158 else { 159 return $self->{'_rate'}; 160 } 161} 162 163=item get_nstates 164 165Getter for number of states (bases). 166 167 Type : method 168 Title : get_nstates 169 Usage : $model->get_nstates; 170 Function: Getter for transition rate between nucleotides. 171 Returns : scalar 172 Args : None 173 174=cut 175 176sub get_nstates { 177 my $states = _BASE_AT_; 178 return scalar @{ $states }; 179} 180 181=item get_ncat 182 183Getter for number of rate categories. 184 185 Type : method 186 Title : get_ncat 187 Usage : $model->get_ncat; 188 Function: Getter for number of rate categories. 189 Returns : scalar 190 Args : None 191 192=cut 193 194sub get_ncat { shift->{'_ncat'} } 195 196=item get_catweights 197 198Getter for weights on rate categories. 199 200 Type : method 201 Title : get_catweights 202 Usage : $model->get_catweights; 203 Function: Getter for number of rate categories. 204 Returns : array 205 Args : None 206 207=cut 208 209sub get_catweights { shift->{'_catweights'} } 210 211=item get_kappa 212 213Getter for transition/transversion ratio. 214 215 Type : method 216 Title : get_kappa 217 Usage : $model->get_kappa; 218 Function: Getter for transition/transversion ratio. 219 Returns : scalar 220 Args : None 221 222=cut 223 224sub get_kappa { shift->{'_kappa'} } 225 226=item get_alpha 227 228Getter for shape parameter. 229 230 Type : method 231 Title : get_alpha 232 Usage : $model->get_alpha; 233 Function: Getter for shape parameter. 234 Returns : scalar 235 Args : None 236 237=cut 238 239sub get_alpha { shift->{'_alpha'} } 240 241=item get_mu 242 243Getter for overall mutation rate. 244 245 Type : method 246 Title : get_mu 247 Usage : $model->get_mu; 248 Function: Getter for overall mutation rate. 249 Returns : scalar 250 Args : None 251 252=cut 253 254sub get_mu { shift->{'_mu'} } 255 256=item get_pinvar 257 258Getter for proportion of invariant sites. 259 260 Type : method 261 Title : get_pinvar 262 Usage : $model->get_pinvar; 263 Function: Getter for proportion of invariant sites. 264 Returns : scalar 265 Args : None 266 267=cut 268 269sub get_pinvar { shift->{'_pinvar'} } 270 271=item get_pi 272 273Getter for base frequencies. 274 275 Type : method 276 Title : get_pi 277 Usage : $model->get_pi; 278 Function: Getter for base frequencies. 279 Returns : array 280 Args : Optional: 281 Base (A, C, T or G) 282 283=cut 284 285sub get_pi { 286 my $self = shift; 287 $self->{'_pi'} = [] if not $self->{'_pi'}; 288 if (@_) { 289 my $base = uc shift; 290 return $self->{'_pi'}->[ _INDEX_OF_()->{$base} ]; 291 } 292 else { 293 return $self->{'_pi'}; 294 } 295} 296 297=item get_median 298 299Getter for median for gamma-modeled rate categories. 300 301 Type : method 302 Title : get_median 303 Usage : $model->get_median; 304 Function: Getter for median. 305 Returns : scalar 306 Args : None 307 308=cut 309 310sub get_median { shift->{'_median'} } 311 312=item set_rate 313 314Setter for substitution rate. 315 316 Type : method 317 Title : set_rate 318 Usage : $model->set_rate(1); 319 Function: Set nucleotide transition rates. 320 Returns : A Bio::Phylo::Models::Substitution::Dna object. 321 Args : scalar or array of arrays (4x4) 322 323=cut 324 325sub set_rate { 326 my ( $self, $q ) = @_; 327 ref $q eq 'ARRAY' or throw 'BadArgs' => 'Not an array ref!'; 328 scalar @{$q} == 4 or throw 'BadArgs' => 'Q matrix must be 4 x 4'; 329 for my $row ( @{$q} ) { 330 scalar @{$row} == 4 or throw 'BadArgs' => 'Q matrix must be 4 x 4'; 331 } 332 $self->{'_rate'} = $q; 333 return $self; 334} 335 336=item set_ncat 337 338Setter for number of rate categories. 339 340 Type : method 341 Title : set_ncat 342 Usage : $model->set_ncat(6); 343 Function: Set the number of rate categoeries. 344 Returns : A Bio::Phylo::Models::Substitution::Dna object. 345 Args : scalar 346 347=cut 348 349sub set_ncat { 350 my $self = shift; 351 $self->{'_ncat'} = shift; 352 return $self; 353} 354 355=item set_catweights 356 357Setter for weights on rate categories. 358 359 Type : method 360 Title : set_catweights 361 Usage : $model->get_catweights; 362 Function: Set number of rate categories. 363 Returns : A Bio::Phylo::Models::Substitution::Dna object. 364 Args : array 365 366=cut 367 368sub set_catweights { 369 my $self = shift; 370 $self->{'_catweights'} = shift; 371 return $self; 372} 373 374=item set_kappa 375 376Setter for weights on rate categories. 377 378 Type : method 379 Title : set_kappa 380 Usage : $model->set_kappa(2); 381 Function: Set transition/transversion ratio. 382 Returns : A Bio::Phylo::Models::Substitution::Dna object. 383 Args : scalar 384 385=cut 386 387sub set_kappa { 388 my $self = shift; 389 $self->{'_kappa'} = shift; 390 return $self; 391} 392 393=item set_alpha 394 395Setter for shape parameter. 396 397 Type : method 398 Title : set_alpha 399 Usage : $model->set_alpha(1); 400 Function: Set shape parameter. 401 Returns : A Bio::Phylo::Models::Substitution::Dna object. 402 Args : scalar 403 404=cut 405 406sub set_alpha { 407 my $self = shift; 408 $self->{'_alpha'} = shift; 409 return $self; 410} 411 412=item set_mu 413 414Setter for overall mutation rate. 415 416 Type : method 417 Title : set_mu 418 Usage : $model->set_mu(0.5); 419 Function: Set overall mutation rate. 420 Returns : A Bio::Phylo::Models::Substitution::Dna object. 421 Args : scalar 422 423=cut 424 425sub set_mu { 426 my $self = shift; 427 $self->{'_mu'} = shift; 428 return $self; 429} 430 431=item set_pinvar 432 433Set for proportion of invariant sites. 434 435 Type : method 436 Title : set_pinvar 437 Usage : $model->set_pinvar(0.1); 438 Function: Set proportion of invariant sites. 439 Returns : A Bio::Phylo::Models::Substitution::Dna object. 440 Args : scalar 441 442=cut 443 444sub set_pinvar { 445 my $self = shift; 446 my $pinvar = shift; 447 if ( $pinvar <= 0 || $pinvar >= 1 ) { 448 throw 'BadArgs' => "Pinvar not between 0 and 1"; 449 } 450 $self->{'_pinvar'} = $pinvar; 451 return $self; 452} 453 454=item set_pi 455 456Setter for base frequencies. 457 458 Type : method 459 Title : get_pi 460 Usage : $model->set_pi((0.2, 0.2, 0.3, 0.3)); 461 Function: Set base frequencies. 462 Returns : A Bio::Phylo::Models::Substitution::Dna object. 463 Args : array of four base frequencies (A, C, G, T) 464 Comments: Base frequencies must sum to one 465 466=cut 467 468sub set_pi { 469 my ( $self, $pi ) = @_; 470 ref $pi eq 'ARRAY' or throw 'BadArgs' => "Not an array ref!"; 471 my $total = 0; 472 $total += $_ for @{$pi}; 473 my $epsilon = 0.000001; 474 abs(1 - $total) < $epsilon or throw 'BadArgs' => 'Frequencies must sum to one'; 475 $self->{'_pi'} = $pi; 476 return $self; 477} 478 479=item set_median 480 481Setter for median for gamma-modeled rate categories. 482 483 Type : method 484 Title : set_median 485 Usage : $model->set_median(1); 486 Function: Setter for median. 487 Returns : A Bio::Phylo::Models::Substitution::Dna object. 488 Args : scalar 489 490=cut 491 492sub set_median { 493 my $self = shift; 494 $self->{'_median'} = !!shift; 495 return $self; 496} 497 498=item modeltest 499 500Performing a modeltest using the package 'phangorn' in 501R (Schliep, Bioinformatics (2011) 27 (4): 592-593) from an 502DNA alignment. If no tree is given as argument, a neighbor-joining 503tree is generated from the alignment to perform model testing. 504Selects the model with the minimum AIC. 505 506 Type : method 507 Title : modeltest 508 Usage : $model->modeltest(-matrix=>$matrix); 509 Function: Determine DNA substitution model from alignment. 510 Returns : An object which is subclass of Bio::Phylo::Models::Substitution::Dna. 511 Args : -matrix: A Bio::Phylo::Matrices::Matrix object 512 Optional: 513 -tree: A Bio::Phylo::Forest::Tree object 514 -timeout: Timeout in seconds to prevent getting stuck in an R process. 515 Comments: Prerequisites: Statistics::R, R, and the R package phangorn. 516 517=cut 518 519sub modeltest { 520 my ($self, %args) = @_; 521 522 my $matrix = $args{'-matrix'}; 523 my $tree = $args{'-tree'}; 524 my $timeout = $args{'-timeout'}; 525 526 my $model; 527 528 if ( looks_like_class 'Statistics::R' ) { 529 530 eval { 531 # phangorn needs files as input 532 my ($fasta_fh, $fasta) = tempfile(); 533 print $fasta_fh unparse('-phylo'=>$matrix, '-format'=>'fasta'); 534 close $fasta_fh; 535 536 # instanciate R and lcheck if phangorn is installed 537 my $R = Statistics::R->new; 538 $R->timeout($timeout) if $timeout; 539 $R->run(q[options(device=NULL)]); 540 $R->run(q[package <- require("phangorn")]); 541 542 if ( ! $R->get(q[package]) eq "TRUE") { 543 $logger->warn("R library phangorn must be installed to run modeltest"); 544 return $model; 545 } 546 547 # read data 548 $R->run(qq[data <- read.FASTA("$fasta")]); 549 550 # remove temp file 551 cleanup(); 552 553 if ( $tree ) { 554 # make copy of tree since it will be pruned 555 my $current_tree = parse('-format'=>'newick', '-string'=>$tree->to_newick)->first; 556 # prune out taxa from tree that are not present in the data 557 my @taxon_names = map {$_->get_name} @{ $matrix->get_entities }; 558 $logger->debug('pruning input tree'); 559 $current_tree->keep_tips(\@taxon_names); 560 $logger->debug('pruned input tree: ' . $current_tree->to_newick); 561 562 if ( ! $current_tree or scalar( @{ $current_tree->get_terminals } ) < 3 ) { 563 $logger->warn('pruned tree has too few tip labels, determining substitution model using NJ tree'); 564 $R->run(q[test <- modelTest(phyDat(data))]); 565 } 566 else { 567 my $newick = $current_tree->to_newick; 568 569 $R->run(qq[tree <- read.tree(text="$newick")]); 570 # call modelTest 571 $logger->debug("calling modelTest from R package phangorn"); 572 $R->run(q[test <- modelTest(phyDat(data), tree=tree)]); 573 } 574 } 575 else { 576 # modelTest will estimate tree 577 $R->run(q[test <- modelTest(phyDat(data))]); 578 } 579 580 # get model with lowest Aikaike information criterion 581 $R->run(q[model <- test[which(test$AIC==min(test$AIC)),]$Model]); 582 my $modeltype = $R->get(q[model]); 583 $logger->info("estimated DNA evolution model $modeltype"); 584 585 # determine model parameters 586 $R->run(q[env <- attr(test, "env")]); 587 $R->run(q[fit <- eval(get(model, env), env)]); 588 589 # get base freqs 590 my $pi = $R->get(q[fit$bf]); 591 592 # get overall mutation rate 593 my $mu = $R->get(q[fit$rate]); 594 595 # get lower triangle of rate matrix (column order ACGT) 596 # and fill whole matrix; set diagonal values to 1 597 my $q = $R->get(q[fit$Q]); 598 my $rate_matrix = [ [ 1, $q->[0], $q->[1], $q->[3] ], 599 [ $q->[0], 1, $q->[2], $q->[4] ], 600 [ $q->[1], $q->[2], 1, $q->[5] ], 601 [ $q->[3], $q->[4], $q->[5], 1 ] 602 ]; 603 604 # create model with specific parameters dependent on primary model type 605 if ( $modeltype =~ /JC/ ) { 606 require Bio::Phylo::Models::Substitution::Dna::JC69; 607 $model = Bio::Phylo::Models::Substitution::Dna::JC69->new(); 608 } 609 elsif ( $modeltype =~ /F81/ ) { 610 require Bio::Phylo::Models::Substitution::Dna::F81; 611 $model = Bio::Phylo::Models::Substitution::Dna::F81->new('-pi' => $pi); 612 } 613 elsif ( $modeltype =~ /GTR/ ) { 614 require Bio::Phylo::Models::Substitution::Dna::GTR; 615 $model = Bio::Phylo::Models::Substitution::Dna::GTR->new('-pi' => $pi); 616 } 617 elsif ( $modeltype =~ /HKY/ ) { 618 require Bio::Phylo::Models::Substitution::Dna::HKY85; 619 # transition/transversion ratio kappa determined by transiton A->G/A->C in Q matrix 620 my $kappa = $R->get(q[fit$Q[2]/fit$Q[1]]); 621 $model = Bio::Phylo::Models::Substitution::Dna::HKY85->new('-kappa' => $kappa, '-pi' => $pi ); 622 } 623 elsif ( $modeltype =~ /K80/ ) { 624 require Bio::Phylo::Models::Substitution::Dna::K80; 625 my $kappa = $R->get(q[fit$Q[2]]); 626 $model = Bio::Phylo::Models::Substitution::Dna::K80->new( 627 '-pi' => $pi, 628 '-kappa' => $kappa ); 629 } 630 # Model is unknown (e.g. phangorn's SYM ?) 631 else { 632 $logger->debug("unknown model type, setting to generic DNA substitution model"); 633 $model = Bio::Phylo::Models::Substitution::Dna->new( 634 '-pi' => $pi ); 635 } 636 637 # set gamma parameters 638 if ( $modeltype =~ /\+G/ ) { 639 $logger->debug("setting gamma parameters for $modeltype model"); 640 # shape of gamma distribution 641 my $alpha = $R->get(q[fit$shape]); 642 $model->set_alpha($alpha); 643 # number of categories for Gamma distribution 644 my $ncat = $R->get(q[fit$k]); 645 $model->set_ncat($ncat); 646 # weights for rate categories 647 my $catweights = $R->get(q[fit$w]); 648 $model->set_catweights($catweights); 649 } 650 651 # set invariant parameters 652 if ( $modeltype =~ /\+I/ ) { 653 $logger->debug("setting invariant site parameters for $modeltype model"); 654 # get proportion of invariant sites 655 my $pinvar = $R->get(q[fit$inv]); 656 $model->set_pinvar($pinvar); 657 } 658 # set universal parameters 659 $model->set_rate($rate_matrix); 660 $model->set_mu($mu); 661 }; 662 # catch possible R errors (e.g. timeout) 663 if ($@) { 664 $logger->warn("modeltest not successful : " . $@); 665 } 666 } 667 else { 668 $logger->warn("Statistics::R must be installed to run modeltest"); 669 } 670 671 return $model; 672} 673 674=item to_string 675 676Get string representation of model in specified format 677(paup, phyml, mrbayes or garli) 678 679 Type : method 680 Title : to_string 681 Usage : $model->to_string(-format=>'mrbayes'); 682 Function: Write model to string. 683 Returns : scalar 684 Args : scalar 685 Comments: format must be either paup, phyml, mrbayes or garli 686 687=cut 688 689sub to_string { 690 my $self = shift; 691 my %args = looks_like_hash @_; 692 if ( $args{'-format'} =~ m/paup/i ) { 693 return $self->_to_paup_string(@_); 694 } 695 if ( $args{'-format'} =~ m/phyml/i ) { 696 return $self->_to_phyml_string(@_); 697 } 698 if ( $args{'-format'} =~ m/mrbayes/i ) { 699 return $self->_to_mrbayes_string(@_); 700 } 701 if ( $args{'-format'} =~ m/garli/i ) { 702 return $self->_to_garli_string(@_); 703 } 704} 705 706sub _to_garli_string { 707 my $self = shift; 708 my $nst = $self->get_nst; 709 my $string = "ratematrix ${nst}\n"; 710 if ( my $pinvar = $self->get_pinvar ) { 711 $string .= "invariantsites fixed\n"; 712 } 713 if ( my $ncat = $self->get_ncat ) { 714 $string .= "numratecats ${ncat}\n"; 715 } 716 if ( my $alpha = $self->get_alpha ) { 717 $string .= "ratehetmodel gamma\n"; 718 } 719 return $string; 720} 721 722sub _to_mrbayes_string { 723 my $self = shift; 724 my $string = 'lset '; 725 $string .= ' nst=' . $self->get_nst; 726 if ( $self->get_pinvar && $self->get_alpha ) { 727 $string .= ' rates=invgamma'; 728 if ( $self->get_ncat ) { 729 $string .= ' ngammacat=' . $self->get_ncat; 730 } 731 } 732 elsif ( $self->get_pinvar ) { 733 $string .= ' rates=propinv'; 734 } 735 elsif ( $self->get_alpha ) { 736 $string .= ' rates=gamma'; 737 if ( $self->get_ncat ) { 738 $string .= ' ngammacat=' . $self->get_ncat; 739 } 740 } 741 $string .= ";\n"; 742 if ( $self->get_kappa && $self->get_nst == 2 ) { 743 $string .= 'prset tratiopr=fixed(' . $self->get_kappa . ");\n"; 744 } 745 my @rates; 746 push @rates, $self->get_rate( 'A' => 'C' ); 747 push @rates, $self->get_rate( 'A' => 'G' ); 748 push @rates, $self->get_rate( 'A' => 'T' ); 749 push @rates, $self->get_rate( 'C' => 'G' ); 750 push @rates, $self->get_rate( 'C' => 'T' ); 751 push @rates, $self->get_rate( 'G' => 'T' ); 752 $string .= 'prset revmatpr=fixed(' . join( ',', @rates ) . ");\n"; 753 754 if ( $self->get_pi('A') 755 && $self->get_pi('C') 756 && $self->get_pi('G') 757 && $self->get_pi('T') ) 758 { 759 my @freqs; 760 push @freqs, $self->get_pi('A'); 761 push @freqs, $self->get_pi('C'); 762 push @freqs, $self->get_pi('G'); 763 push @freqs, $self->get_pi('T'); 764 $string .= 'prset statefreqpr=fixed(' . join( ',', @freqs ) . ");\n"; 765 } 766 if ( $self->get_alpha ) { 767 $string .= 'prset shapepr=fixed(' . $self->get_alpha . ");\n"; 768 } 769 if ( $self->get_pinvar ) { 770 $string .= 'prset pinvarpr=fixed(' . $self->get_pinvar . ");\n"; 771 } 772} 773 774sub _to_phyml_string { 775 my $self = shift; 776 my $m = ref $self; 777 $m =~ s/.+://; 778 my $string = "--model $m"; 779 if ( $self->get_pi('A') 780 && $self->get_pi('C') 781 && $self->get_pi('G') 782 && $self->get_pi('T') ) 783 { 784 my @freqs; 785 push @freqs, $self->get_pi('A'); 786 push @freqs, $self->get_pi('C'); 787 push @freqs, $self->get_pi('G'); 788 push @freqs, $self->get_pi('T'); 789 $string .= ' -f ' . join ' ', @freqs; 790 } 791 if ( $self->get_nst == 2 and defined( my $kappa = $self->get_kappa ) ) { 792 $string .= ' --ts/tv ' . $kappa; 793 } 794 if ( $self->get_pinvar ) { 795 $string .= ' --pinv ' . $self->get_pinvar; 796 } 797 if ( $self->get_ncat ) { 798 $string .= ' --nclasses ' . $self->get_ncat; 799 $string .= ' --use_median' if $self->get_median; 800 } 801 if ( $self->get_alpha ) { 802 $string .= ' --alpha ' . $self->get_alpha; 803 } 804 return $string; 805} 806 807sub _to_paup_string { 808 my $self = shift; 809 my $nst = $self->get_nst; 810 my $string = 'lset nst=' . $nst; 811 if ( $nst == 2 and defined( my $kappa = $self->get_kappa ) ) { 812 $string .= ' tratio=' . $kappa; 813 } 814 if ( $nst == 6 ) { 815 my @rates; 816 push @rates, $self->get_rate( 'A' => 'C' ); 817 push @rates, $self->get_rate( 'A' => 'G' ); 818 push @rates, $self->get_rate( 'A' => 'T' ); 819 push @rates, $self->get_rate( 'C' => 'G' ); 820 push @rates, $self->get_rate( 'C' => 'T' ); 821 $string .= ' rmatrix=(' . join( ' ', @rates ) . ')'; 822 } 823 if ( $self->get_pi('A') && $self->get_pi('C') && $self->get_pi('G') ) { 824 my @freqs; 825 push @freqs, $self->get_pi('A'); 826 push @freqs, $self->get_pi('C'); 827 push @freqs, $self->get_pi('G'); 828 $string .= ' basefreq=(' . join( ' ', @freqs ) . ')'; 829 } 830 if ( $self->get_alpha ) { 831 $string .= ' rates=gamma shape=' . $self->get_alpha; 832 } 833 if ( $self->get_ncat ) { 834 $string .= ' ncat=' . $self->get_ncat; 835 $string .= ' reprate=' . ( $self->get_median ? 'median' : 'mean' ); 836 } 837 if ( $self->get_pinvar ) { 838 $string .= ' pinvar=' . $self->get_pinvar; 839 } 840 return $string . ';'; 841} 842 843sub _type { _MODEL_ } 844 845=back 846 847=head1 SEE ALSO 848 849There is a mailing list at L<https://groups.google.com/forum/#!forum/bio-phylo> 850for any user or developer questions and discussions. 851 852=over 853 854=item L<Bio::Phylo::Manual> 855 856Also see the manual: L<Bio::Phylo::Manual> and L<http://rutgervos.blogspot.com>. 857 858=back 859 860=head1 CITATION 861 862If you use Bio::Phylo in published research, please cite it: 863 864B<Rutger A Vos>, B<Jason Caravas>, B<Klaas Hartmann>, B<Mark A Jensen> 865and B<Chase Miller>, 2011. Bio::Phylo - phyloinformatic analysis using Perl. 866I<BMC Bioinformatics> B<12>:63. 867L<http://dx.doi.org/10.1186/1471-2105-12-63> 868 869=cut 870 8711; 872