File Coverage

blib/lib/Bio/Phylo/Models/Substitution/Dna.pm
Criterion Covered Total %
statement 86 280 30.7
branch 14 108 12.9
condition 1 42 2.3
subroutine 21 37 56.7
pod 24 24 100.0
total 146 491 29.7


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