compactor work
[jalview.git] / forester / perl / phylo_pl.pl
1 #!/usr/bin/perl -W
2 #
3 # $Id: phylo_pl.pl,v 1.32 2010/12/13 19:00:22 cmzmasek Exp $
4 #
5 # FORESTER -- software libraries and applications
6 # for evolutionary biology research and applications.
7 #
8 # Copyright (C) 2008-2014 Christian M. Zmasek
9 # Copyright (C) 2008-2009 Burnham Institute for Medical Research
10 # All rights reserved
11
12 # This library is free software; you can redistribute it and/or
13 # modify it under the terms of the GNU Lesser General Public
14 # License as published by the Free Software Foundation; either
15 # version 2.1 of the License, or (at your option) any later version.
16 #
17 # This library is distributed in the hope that it will be useful,
18 # but WITHOUT ANY WARRANTY; without even the implied warranty of
19 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
20 # Lesser General Public License for more details.
21
22 # You should have received a copy of the GNU Lesser General Public
23 # License along with this library; if not, write to the Free Software
24 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
25 #
26 # Contact: phylosoft @ gmail . com
27 #     WWW: www.phylosoft.org/forester
28 #
29 #
30 #
31 #  Requirements  phylo_pl is part of the FORESTER libraries.
32 #  ------------  Many of its global variables are set via forester.pm.
33 #
34 #
35 #  Note. Use xt.pl (for Pfam alignments) or mt.pl (for other alignments)
36 #  to run phylo_pl.pl on whole directories of alignments files.       
37 #
38 #
39 #
40 #
41 # =========================
42 #
43 # METHOD ORDER (IMPORTANT!)
44 # 1. FastME
45 # 2. phylip NJ
46 # 3. phylip fitch FM
47 # 4. phylip fitch ME
48 # 5. BIONJ
49 # 6. Weighbor
50 # 7. Raxml
51 # 8. phyml
52 # 9. phylip proml
53 # 10. phylip protpars
54 # 11. all
55 #==========================
56
57 use strict;
58
59 use FindBin;
60 use lib $FindBin::Bin;
61 use forester;
62
63 my $VERSION                = "1.0.1";
64 my $LAST_MODIFIED          = "2009.10.02";
65
66 my $RAXML_MODEL_BASE       = "PROTGAMMA";
67 my $RAXML_ALGORITHM        = "a";
68
69 my $TEMP_DIR_DEFAULT       = "/tmp/phylo_pl_"; # Where all the infiles, outfiles, etc will be created.
70    
71 my $bootstraps             = 100; # 0,1: do not bootstrap. Default: 100
72 my $matrix                 = 5;   # 0 = JTT
73                                   # 1 = PAM 
74                                   # 2 = BLOSUM 62
75                                   # 3 = mtREV24
76                                   # 5 = VT - default 
77                                   # 6 = WAG 
78                                   # 7 = auto in puzzle
79                                   # 8 = DCMut in PHYML, VT in TREE-PUZZLE
80 my $rate_heterogeneity     = 0;   # 0 = Uniform rate (default)
81                                   # 1 = 8 Gamma distributed rates
82                                   # 2 = Two rates (1 invariable + 1 variable)
83                                   # 3 = Mixed (1 invariable + 8 Gamma rates)
84 my $seed                   = 9;   # Seed for random number generators. Default: 9
85 my $keep_multiple_trees    = 0;   # 0: delete multiple tree file
86                                   # 1: do not delete multiple tree file
87 my $exact_parameter_est    = 0;   # 0: no; 1: yes
88
89 my $phyml_rel_substitution_rate_cat = 4;
90
91 my $jumbles                = 2;
92 my $use_fastme             = 0;   # 0: no; 1: yes
93 my $use_phylip_nj          = 0;   # 0: no; 1: yes
94 my $use_phylip_fitch_fm    = 0;   # 0: no; 1: yes
95 my $use_phylip_fitch_me    = 0;   # 0: no; 1: yes
96 my $use_bionj              = 0;   # 0: no; 1: yes
97 my $use_weighbor           = 0;   # 0: no; 1: yes
98 my $use_raxml              = 0;   # 0: no; 1: yes
99 my $use_phyml              = 0;   # 0: no; 1: yes
100 my $use_proml              = 0;   # 0: no; 1: yes
101 my $use_protpars           = 0;   # 0: no; 1: yes
102 my $use_global_rearr       = 0;   # 0: no; 1: yes
103 my $estimate_invar_sites   = 0;   # 0: no; 1: yes
104
105 my $fastme_init_tree_opt   = "NJ";
106
107 my %seqnames       = ();           # number   =>  seqname 
108 my %numbers        = ();           # seqname  =>  number
109 my $options        = "";
110 my $infile         = "";
111 my $pwdfile        = "";
112 my $outfile        = "";
113 my $logfile        = "";
114 my $multipwdfile   = "";
115 my $distancefile   = "";
116 my $log            = "";
117 my $ii             = 0;
118 my $temp_dir       = "";
119 my $current_dir    = "";
120 my @out            = ();
121 my $number_of_seqs = 0;
122 my $number_of_aa   = 0;
123
124 my $use_pwd_based_methods = 0;
125
126 print( "\n");
127 print( "phylo_pl $VERSION ($LAST_MODIFIED)\n" );
128 print( "__________________________________\n");
129 print( "\n\n");
130
131
132
133 unless ( @ARGV == 2 || @ARGV == 3 || @ARGV == 4 || @ARGV == 5 ) {
134     &printUsage();
135     exit ( -1 ); 
136 }
137
138
139
140 # Analyzes the options:
141 # ---------------------
142
143 if ( $ARGV[ 0 ] =~ /^-.+/ ) {
144     
145     unless ( @ARGV > 2 ) {
146         &printUsage();
147         exit ( -1 ); 
148     }
149     $options = $ARGV[ 0 ];
150    
151     if ( @ARGV != 3 && @ARGV != 4 ) {
152         &printUsage();
153         exit ( -1 ); 
154     }
155     $infile  = $ARGV[ 1 ];
156     $outfile = $ARGV[ 2 ];
157     if ( @ARGV == 4 ) {
158         $temp_dir = $ARGV[ 3 ];
159     }
160     if ( $options =~ /B(\d+)/ ) {
161         $bootstraps = $1;
162         if ( $bootstraps <= 1 ) {
163             $bootstraps = 0;
164         }
165         elsif ( $bootstraps <= 9 ) {
166             $bootstraps = 0;
167             print "\n\nphylo_pl: WARNING: Bootstrap number must be devisable by 10,\nno bootstrapping.\n\n";
168         }   
169         elsif ( $bootstraps % 10 != 0 ) {
170             $bootstraps = $bootstraps - $bootstraps % 10; # to ensure $bootstraps % 10 == 0
171             print "\n\nphylo_pl: WARNING: Bootstrap number must be devisable by 10,\nset to $bootstraps.\n\n";
172         }    
173     }
174     if ( $options =~ /n/ ) {
175         $use_phylip_nj = 1; 
176     }
177     if ( $options =~ /q@(\d)/ ) {
178         $use_fastme = 1;
179         my $opt = $1;
180         if ( $opt == 1 ) {
181             $fastme_init_tree_opt = "GME";
182         }
183         elsif ( $opt == 2 ) {
184             $fastme_init_tree_opt = "BME";
185         }
186         elsif ( $opt == 3 ) {
187             $fastme_init_tree_opt = "NJ";
188         }
189         else {
190             &printUsage();
191             exit ( -1 );
192         }
193     }
194     if ( $options =~ /f/ ) {
195         $use_phylip_fitch_fm = 1; 
196     }
197     if ( $options =~ /e/ ) {
198         $use_phylip_fitch_me = 1; 
199     }
200     if ( $options =~ /b/ ) {
201         $use_bionj = 1; 
202     }
203     if ( $options =~ /w/ ) {
204         $use_weighbor = 1;
205     }   
206     if ( $options =~ /x/ ) {
207         $use_raxml = 1;
208     }    
209     if ( $options =~ /y/ ) {
210         $use_phyml = 1;
211     }    
212     if ( $options =~ /o/ ) {
213         $use_proml = 1;
214     }
215     if ( $options =~ /p/ ) {
216         $use_protpars = 1;
217     }
218     if ( $options =~ /G/ ) {
219         $use_global_rearr = 1;
220     } 
221     if ( $options =~ /I/ ) {
222         $estimate_invar_sites = 1;
223     }
224     if ( $options =~ /j(\d+)/ ) {
225         $jumbles = $1;
226         if ( $jumbles < 1 ) {
227             $jumbles = 0;
228         }
229     }
230     if ( $options =~ /r(\d+)/ ) {
231         $phyml_rel_substitution_rate_cat = $1;
232         if ( $phyml_rel_substitution_rate_cat < 1 ) {
233             $phyml_rel_substitution_rate_cat = 1;
234         }
235     }
236     if ( $options =~ /J/ ) {
237         $matrix = 0;      # JTT
238     }
239     if ( $options =~ /P/ ) {
240         $matrix = 1;      # PAM
241     }
242     if ( $options =~ /L/ ) {
243         $matrix = 2;      # Blosum 62
244     }
245     if ( $options =~ /M/ ) {
246         $matrix = 3;      # mtREV24
247     }
248     if ( $options =~ /W/ ) {
249         $matrix = 6;      # WAG
250     }
251     if ( $options =~ /A/ ) {
252         $matrix = 7;      # auto
253     }
254     if ( $options =~ /D/ ) {
255         $matrix = 8;      # DCMut in PHYML and RAXML, VT in PUZZLE
256     }
257     if ( $options =~ /S(\d+)/ ) {
258         $seed = $1; 
259     }
260     if ( $options =~ /X/ ) {
261         $keep_multiple_trees = 1; 
262     }
263     if ( $options =~ /E/ ) {
264         $exact_parameter_est = 1;
265     }
266     if ( $options =~ /g/ ) {
267         $rate_heterogeneity = 1; 
268     }
269     if ( $options =~ /t/ ) {
270         $rate_heterogeneity = 2; 
271     }
272     if ( $options =~ /m/ ) {
273         $rate_heterogeneity = 3; 
274     }
275 }
276 else {
277     unless ( @ARGV == 2 || @ARGV == 3 ) {
278         &printUsage();
279         exit ( -1 );  
280     }
281     $infile  = $ARGV[ 0 ];
282     $outfile = $ARGV[ 1 ];
283     if ( @ARGV == 3 ) {
284         $temp_dir = $ARGV[ 2 ];
285     } 
286 }
287
288 if ( $use_fastme    != 1 &&
289      $use_phylip_nj != 1 &&
290      $use_phylip_fitch_fm != 1 &&
291      $use_phylip_fitch_me != 1 &&
292      $use_bionj != 1 &&         
293      $use_weighbor != 1 &&
294      $use_raxml != 1 &&
295      $use_phyml != 1 &&
296      $use_proml != 1 &&         
297      $use_protpars != 1 ) {
298     
299      $use_fastme = 1;
300      $use_phylip_nj = 1;
301      $use_phylip_fitch_fm = 1;
302      $use_phylip_fitch_me = 1;
303      $use_bionj = 1;    
304      $use_raxml = 1;     
305      $use_weighbor = 1;
306      $use_phyml = 1;
307      $use_proml = 1;         
308      $use_protpars = 1; 
309 }
310
311
312 if ( $use_fastme    == 1 ||
313      $use_phylip_nj == 1 ||
314      $use_phylip_fitch_fm == 1 ||
315      $use_phylip_fitch_me == 1 ||
316      $use_bionj == 1          ||
317      $use_weighbor == 1 ) { 
318     $use_pwd_based_methods = 1;
319 }
320 else {
321     $use_pwd_based_methods = 0;
322 }    
323
324 $current_dir = `pwd`;
325 $current_dir =~ s/\s//;
326
327 if ( $outfile !~ /^\// ) {
328     # outfile is not absolute path.
329     $outfile = $current_dir."/".$outfile;
330 }
331
332
333
334 # TREE-PUZZLE sets the option in this way:
335 # If two rates or mixed, exact parameter estimates are used.
336 if ( $rate_heterogeneity == 2
337 ||   $rate_heterogeneity == 3 ) {
338     $exact_parameter_est = 1
339 }
340
341
342 if ( $outfile =~ /\.xml$/i ) {
343     $outfile =~ s/\.xml//i;
344
345 elsif ( $outfile =~ /\.aln$/i ) {
346     $outfile =~ s/\.aln//i;
347 }
348 elsif ( $outfile =~ /\.fasta$/i ) {
349     $outfile =~ s/\.fasta//i;
350 }
351 elsif ( $outfile =~ /\.fas$/i ) {
352     $outfile =~ s/\.fas//i;
353 }
354 elsif ( $outfile =~ /\.seqs$/i ) {
355     $outfile =~ s/\.seqs//i;
356 }  
357
358
359 $logfile          = $outfile.$LOG_FILE_SUFFIX;
360 $multipwdfile     = $outfile.$MULTIPLE_PWD_FILE_SUFFIX;
361 $distancefile     = $outfile.$SUFFIX_PWD_NOT_BOOTS;
362
363 &dieIfFileExists( $logfile );
364 &dieIfFileExists( $multipwdfile );
365 &dieIfFileExists( $distancefile );
366
367
368 my $fastme_outtree    = $outfile."_fme.xml";
369 my $phylip_nj_outtree = $outfile."_pnj.xml";
370 my $phylip_fm_outtree = $outfile."_pfm.xml";
371 my $phylip_me_outtree = $outfile."_pme.xml";
372 my $bionj_outtree     = $outfile."_bionj.xml";
373 my $weighbor_outtree  = $outfile."_weigh.xml";
374 my $raxml_outtree     = $outfile."_raxml.xml";
375 my $phyml_outtree     = $outfile."_phyml.xml";
376 my $proml_outtree     = $outfile."_proml.xml";
377 my $protpars_outtree  = $outfile."_ppp.xml";
378 my $all_outtree       = $outfile."_comb.xml";
379
380 my $multitreefile_fastme    = $outfile."_fme".$MULTIPLE_TREES_FILE_SUFFIX;
381 my $multitreefile_phylip_nj = $outfile."_pnj".$MULTIPLE_TREES_FILE_SUFFIX;
382 my $multitreefile_phylip_fm = $outfile."_pfm".$MULTIPLE_TREES_FILE_SUFFIX;
383 my $multitreefile_phylip_me = $outfile."_pme".$MULTIPLE_TREES_FILE_SUFFIX;
384 my $multitreefile_bionj     = $outfile."_bionj".$MULTIPLE_TREES_FILE_SUFFIX;
385 my $multitreefile_weighbor  = $outfile."_weigh".$MULTIPLE_TREES_FILE_SUFFIX;
386 my $multitreefile_raxml     = $outfile."_raxml".$MULTIPLE_TREES_FILE_SUFFIX;
387 my $multitreefile_phyml     = $outfile."_phyml".$MULTIPLE_TREES_FILE_SUFFIX;
388 my $multitreefile_proml     = $outfile."_proml".$MULTIPLE_TREES_FILE_SUFFIX;
389 my $multitreefile_protpars  = $outfile."_ppp".$MULTIPLE_TREES_FILE_SUFFIX;
390
391 if ( $use_fastme == 1 ) {
392     &dieIfFileExists( $fastme_outtree );
393     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
394         &dieIfFileExists( $multitreefile_fastme );
395     }
396 }
397 if( $use_phylip_nj == 1 ) {
398     &dieIfFileExists( $phylip_nj_outtree );
399     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
400         &dieIfFileExists( $multitreefile_phylip_nj );
401     }
402 }
403 if( $use_phylip_fitch_fm == 1 ) {
404     &dieIfFileExists( $phylip_fm_outtree );
405     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
406         &dieIfFileExists( $multitreefile_phylip_fm );
407     }
408 }
409 if( $use_phylip_fitch_me == 1 ) {
410     &dieIfFileExists( $phylip_me_outtree );
411     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
412         &dieIfFileExists( $multitreefile_phylip_me );
413     }
414 }
415 if( $use_bionj == 1 ) {
416     &dieIfFileExists( $bionj_outtree );
417     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
418         &dieIfFileExists( $multitreefile_bionj );
419     }
420 }
421 if( $use_weighbor == 1 ) {
422     &dieIfFileExists( $weighbor_outtree );
423     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
424         &dieIfFileExists( $multitreefile_weighbor );
425     }
426 }
427 if( $use_raxml == 1 ) {
428     &dieIfFileExists( $raxml_outtree );
429     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
430         &dieIfFileExists( $multitreefile_raxml );
431     }
432 }
433 if( $use_phyml == 1 ) {
434     &dieIfFileExists( $phyml_outtree );
435     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
436         &dieIfFileExists( $multitreefile_phyml );
437     }
438 }
439 if( $use_proml == 1 ) {
440     &dieIfFileExists( $proml_outtree );
441     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
442         &dieIfFileExists( $multitreefile_proml );
443     }
444
445 if( $use_protpars == 1 ) {
446     &dieIfFileExists( $protpars_outtree );
447     if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
448         &dieIfFileExists( $multitreefile_protpars );
449     }
450
451 if ( $bootstraps > 1 ) {
452      &dieIfFileExists( $all_outtree );
453 }
454 if ( $infile ne "" ) {
455     unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
456         die "\n\nphylo_pl: Input alignment file \"$infile\" does not exist, is empty, or is not a plain textfile.\n\n";
457     }
458 }
459
460
461
462
463 # Prints out the options:
464 # -----------------------
465
466
467 $log = "\n$0 logfile:\n";
468 $log = $log."Version: $VERSION\n\n";
469
470
471
472 if ( $infile ne ""  ) {
473     $log = $log."Input                               : $infile\n";
474 }
475
476 if ( $keep_multiple_trees == 1 && $bootstraps >= 2 ) {
477     $log = $log."Multiple distance matrices          : $multipwdfile\n";
478 }
479
480
481 $log = $log."Bootstraps                          : $bootstraps\n";
482
483 if ( $use_pwd_based_methods == 1 ) {    
484     $log = $log."Prgrm to calculate pairwise dist.   : TREE-PUZZLE (version: $PUZZLE_VERSION)\n";
485 }
486
487
488 if ( $use_fastme == 1 ) {
489     $log = $log."Program to calculate tree           : FastME (version: $FASTME_VERSION)\n";   
490     $log = $log."Method for intial tree in FastME    : $fastme_init_tree_opt\n";
491     $log = $log."Tree swapping (NNI) in FastME       : balanced (default)\n";
492 }
493 if ( $use_phylip_nj == 1 ) {
494     $log = $log."Program to calculate tree           : PHYLIP NEIGHBOR NJ (version: $PHYLIP_VERSION)\n";
495 }
496 if ( $use_phylip_fitch_fm == 1 ) {
497     $log = $log."Program to calculate tree           : PHYLIP FITCH Fitch-Margoliash (version: $PHYLIP_VERSION)\n";
498 }
499 if ( $use_phylip_fitch_me == 1 ) {
500     $log = $log."Program to calculate tree           : PHYLIP FITCH Minimal Evolution (version: $PHYLIP_VERSION)\n";
501 }
502 if ( $use_bionj == 1 ) {
503     $log = $log."Program to calculate tree           : BIONJ (version: $BIONJ_VERSION)\n";
504 }
505 if ( $use_weighbor == 1 ) {
506     $log = $log."Program to calculate tree           : Weighbor [no invariable sites, b=14] (version: $WEIGHBOR_VERSION)\n";
507 }
508 if ( $use_raxml == 1 ) {
509     $log = $log."Program to calculate tree           : RAxML [$RAXML_MODEL_BASE] (uses its own bootstraps, if bootstrapped: -f $RAXML_ALGORITHM) (version: $RAXML_VERSION)\n";
510 }
511 if ( $use_phyml == 1 ) {
512     $log = $log."Program to calculate tree           : PHYML (MLE for gamma distr param and proportion of inv sites) (version: $PHYML_VERSION)\n";
513     $log = $log."# of rel subst rate categories      : $phyml_rel_substitution_rate_cat\n";
514 }
515 if ( $use_proml == 1 ) {
516     $log = $log."Program to calculate tree           : PHYLIP PROML (uses PAM unless JTT selected) (version: $PHYLIP_VERSION)\n";
517 }
518 if ( $use_protpars == 1 ) {
519     $log = $log."Program to calculate tree           : PHYLIP PROTPARS (with global rearrangements) (version: $PHYLIP_VERSION)\n";
520 }
521 if ( $use_phylip_fitch_fm == 1 || $use_phylip_fitch_me == 1 || $use_protpars == 1 || $use_proml ) {
522     $log = $log."Number of jumbles (input order rand): $jumbles\n"; 
523     
524 }
525 if ( $use_phylip_fitch_fm == 1 || $use_phylip_fitch_me == 1 || $use_proml ) {
526     if ( $use_global_rearr == 1 ) {
527         $log = $log."Global rearrangements               : true\n";
528     }
529     else {
530         $log = $log."Global rearrangements               : false\n";
531         
532     }
533 }
534
535 if ( $bootstraps > 0 ) {
536     $log = $log."Prgrm to calculate ML branch lenghts: TREE-PUZZLE (version: $PUZZLE_VERSION)\n";
537 }
538
539 $log = $log."Model                               : ";
540 if ( $matrix == 0 ) { 
541     $log = $log."JTT (Jones et al. 1992)\n";
542 }
543 elsif ( $matrix == 1 ) {
544     $log = $log."PAM (Dayhoff et al. 1978)\n";
545 }
546 elsif ( $matrix == 2 ) {
547     $log = $log."BLOSUM 62 (Henikoff-Henikoff 92)\n";
548 }
549 elsif ( $matrix == 3 ) {
550     $log = $log."mtREV24 (Adachi-Hasegawa 1996)\n";
551 }
552 elsif ( $matrix == 5 ) {
553     $log = $log."VT (Mueller-Vingron 2000)\n";
554 }
555 elsif ( $matrix == 6 ) {
556     $log = $log."WAG (Whelan-Goldman 2000)\n";
557 }
558 elsif ( $matrix == 7 ) {
559     $log = $log."auto in TREE-PUZZLE\n";
560 }
561 elsif ( $matrix == 8 ) {
562     $log = $log."DCMut (Kosial and Goldman, 2005) in PHYML and RAxML, VT in TREE-PUZZLE\n";
563 }
564 else {
565     &dieWithUnexpectedError( "Unknown model: matrix=$matrix" );
566 }
567 if ( $use_raxml == 1 || $use_phyml == 1 ) {
568     if ( $estimate_invar_sites == 1 ) {
569         $log = $log."Estimate proportion of invariable sites in RAXML and/or PHYML: true\n";
570     }
571     else {
572         $log = $log."Estimate proportion of invariable sites in RAXML and/or PHYML: false (proportion \"0.0\" is used in PHYML)\n";
573     }
574 }
575
576 $log = $log."Model of rate heterogeneity (PUZZLE): ";
577 if ( $rate_heterogeneity == 1 ) { 
578     $log = $log."8 Gamma distributed rates\n";
579 }
580 elsif ( $rate_heterogeneity == 2 ) {
581     $log = $log."Two rates (1 invariable + 1 variable)\n";
582 }
583 elsif ( $rate_heterogeneity == 3 ) {
584     $log = $log."Mixed (1 invariable + 8 Gamma rates)\n";
585 }
586 else {
587     $log = $log."Uniform rate\n";
588 }
589 $log = $log."Seed for random number generators   : $seed\n";
590 if ( $exact_parameter_est == 1 ) {
591     $log = $log."Exact parameter estimates in TREE-PUZZLE\n";
592 }
593
594 $log = $log."Start time/date                     : ".`date`;
595
596
597
598
599 # That's where the mischief starts....
600 # ------------------------------------
601
602 $ii = 0;
603
604 srand();
605 my $range = 1000000;
606 my $random_number = int( rand( $range ) );
607
608 if ( $temp_dir eq "" ) {
609     $temp_dir = $TEMP_DIR_DEFAULT;
610 }
611
612 $temp_dir = $temp_dir.$random_number.$ii;
613
614 while ( -e $temp_dir ) {
615     $ii++;
616     $temp_dir = $temp_dir.$random_number.$ii;
617 }
618
619 mkdir( $temp_dir, 0700 )
620 || die "\n\n$0: Could not create <<$temp_dir>>: $!\n\n";
621
622 unless ( ( -e $temp_dir ) && ( -d $temp_dir ) ) {
623     die "\n\n$0: <<$temp_dir>> does not exist, or is not a directory.\n\n";
624 }
625
626
627 &cp( $infile, $temp_dir."/INFILE" );
628 unless ( chmod ( 0600, $temp_dir."/INFILE" ) ) {
629     warn "\n\n$0: Could not chmod. $!\n\n";
630 }
631 $infile = "INFILE";
632
633 chdir ( $temp_dir ) 
634 || die "\n\n$0: Could not chdir to <<$temp_dir>>: $!\n\n";
635
636 &cp( $infile, "infile" );
637 @out = &getNumberOfSeqsAndAas( $infile );
638 $number_of_seqs = $out[ 0 ];
639 $number_of_aa   = $out[ 1 ];
640
641 my $SEQBOOT_OUTFILE = "seqboot_outfile"; 
642
643 if (  $bootstraps > 1 && ( $use_pwd_based_methods == 1
644                                     || $use_phyml == 1
645                                     || $use_proml == 1
646                                     || $use_protpars == 1 ) ) {
647     &executeSeqboot( $seed, $bootstraps );
648     &mv( "outfile", $SEQBOOT_OUTFILE );
649     &rm( "infile" );
650 }
651
652 &cp( $infile, "align" ); 
653
654 if ( $use_pwd_based_methods == 1 ) {
655     # Calculating the pairwise distances (saved in file "infile"): "puzzle"
656     if ( $bootstraps > 1 ) {
657         &executePuzzleBootstrapped( $SEQBOOT_OUTFILE,
658                                     $matrix,
659                                     $number_of_seqs,
660                                     $exact_parameter_est,
661                                     $rate_heterogeneity );
662
663         $pwdfile = $SEQBOOT_OUTFILE.".dist";
664     }
665     else {
666         &executePuzzle( "infile",
667                         $matrix,
668                         $number_of_seqs,
669                         $exact_parameter_est,
670                         $rate_heterogeneity );
671         $pwdfile = "infile.dist";
672     }
673
674
675 &rm( "infile" );
676
677 # Methods based on alignment
678 # --------------------------
679 my $OUTTREE_RAXML    = "outtree_rax";
680 my $OUTTREE_PHYML    = "outtree_phyml";
681 my $OUTTREE_PROML    = "outtree_proml";
682 my $OUTTREE_PROTPARS = "outtree_protpars";
683
684 my $CONSENSUS_RAXML    = "consensus_raxml";
685 my $CONSENSUS_PHYML    = "consensus_phyml";
686 my $CONSENSUS_PROML    = "consensus_proml";
687 my $CONSENSUS_PROTPARS = "consensus_protpars";
688
689 my $OUTTREES_ALL = "outtrees_all";
690 my $all_count = 0;
691
692 if ( $use_raxml == 1 ) {
693    
694     my $model = "---";
695     if ( $matrix == 0 ) { 
696         $model = "JTT";
697     }
698     elsif ( $matrix == 1 ) {
699         $model = "DAYHOFF";
700     }
701     elsif ( $matrix == 2 ) {
702         $model = "BLOSUM62";
703     }
704     elsif ( $matrix == 3 ) {
705         $model = "MTREV";
706     }
707     elsif ( $matrix == 5 ) {
708         $model = "VT";
709     }
710     elsif ( $matrix == 6 ) {
711         $model = "WAG";
712     }
713     elsif ( $matrix == 7 ) {
714         $model = "VT";
715     }
716     elsif ( $matrix == 8 ) {
717         $model = "DCMUT";
718     }
719     else {
720         &dieWithUnexpectedError( "Unknown model: matrix=$matrix" );
721     }
722
723     print( "\n========== RAxML begin =========\n\n" );    
724     # Six arguments:
725     # 1. DNA or Amino-Acids sequence filename (PHYLIP format)
726     # 2. Model, eg. PROTGAMMAIVT
727     # 3. Replicates (bootstrap)
728     # 4. Seed for bootstrap
729     # 5. Output suffix
730     # 6. Algorithm (only for bootstrap, default otherwise)
731     my $invar = "";
732     if ( $estimate_invar_sites == 1 ) {
733         $invar = "I";
734     }
735     
736     # NOTE. RaxML does its own bootstrapping.
737     &executeRaxml( "align", $RAXML_MODEL_BASE.$invar.$model."F", $bootstraps, $seed, "xxx", $RAXML_ALGORITHM );
738     print( "\n========== RAxML end =========\n\n" );
739     
740     &rm( "RAxML_log.xxx" );
741     &rm( "RAxML_parsimonyTree.xxx" );
742     &mv( "RAxML_info.xxx", $outfile."_raxml_info" );
743     if ( $bootstraps > 1 ) {
744         &rm( "RAxML_bestTree.xxx" );
745         &mv( "RAxML_bipartitions.xxx", $CONSENSUS_RAXML );
746         &append( "RAxML_bootstrap.xxx", $OUTTREES_ALL );
747         if ( $keep_multiple_trees == 1 ) {
748             &mv( "RAxML_bootstrap.xxx", $multitreefile_raxml );
749         }
750         else {
751             &rm( "RAxML_bootstrap.xxx" );
752         }
753         $all_count++;
754     }
755     else {
756         &mv( "RAxML_result.xxx", $OUTTREE_RAXML );
757     }
758 }
759
760
761 if ( $use_phyml == 1 ) {
762    
763     my $model = "---";
764     if ( $matrix == 0 ) { 
765         $model = "JTT";
766     }
767     elsif ( $matrix == 1 ) {
768         $model = "Dayhoff";
769     }
770     elsif ( $matrix == 2 ) {
771         $model = "Blosum62";
772     }
773     elsif ( $matrix == 3 ) {
774         $model = "MtREV";
775     }
776     elsif ( $matrix == 5 ) {
777         $model = "VT";
778     }
779     elsif ( $matrix == 6 ) {
780         $model = "WAG";
781     }
782     elsif ( $matrix == 7 ) {
783         $model = "VT";
784     }
785     elsif ( $matrix == 8 ) {
786         $model = "DCMut";
787     }
788     else {
789         &dieWithUnexpectedError( "Unknown model: matrix=$matrix" );
790     }
791
792     my $input = "";
793     if ( $bootstraps > 1 ) {
794         $input = $SEQBOOT_OUTFILE;
795     }
796     else {
797         $input = "align";
798     } 
799     print( "\n========== PHYML begin =========\n\n" );    
800     # Six arguments:
801     # 1. DNA or Amino-Acids sequence filename (PHYLIP format)
802     # 2. number of data sets to analyse (ex:3)
803     # 3. Model: JTT | MtREV | Dayhoff | WAG | VT | DCMut | Blosum62 (Amino-Acids)
804     # 4. number of relative substitution rate categories (ex:4), positive integer
805     # 5. starting tree filename (Newick format), your tree filename | BIONJ for a distance-based tree
806     # 6. 1 to estimate proportion of invariable sites, otherwise, fixed proportion "0.0" is used
807     # PHYML produces several results files :
808     # <sequence file name>_phyml_lk.txt : likelihood value(s)
809     # <sequence file name>_phyml_tree.txt : inferred tree(s)
810     # <sequence file name>_phyml_stat.txt : detailed execution stats 
811     &executePhyml( $input, $bootstraps, $model, $phyml_rel_substitution_rate_cat, "BIONJ", $estimate_invar_sites );
812     print( "\n========== PHYML end =========\n\n" );
813     
814     &rm( $input."_phyml_lk.txt" );
815     &mv( $input."_phyml_tree.txt", $OUTTREE_PHYML );
816     if ( -e $outfile."_phyml_stat" ) {
817         &rm( $outfile."_phyml_stat" ); 
818     }    
819     &mv( $input."_phyml_stat.txt", $outfile."_phyml_stat" );
820     if ( $bootstraps > 1 ) {
821         &append( $OUTTREE_PHYML, $OUTTREES_ALL );
822         $all_count++;
823     }
824
825 }
826
827 if ( $use_proml == 1 ) {
828     my $input = "";
829     if ( $bootstraps > 1 ) {
830         $input = $SEQBOOT_OUTFILE;
831     }
832     else {
833         $input = "align";
834     }    
835     print( "\n========== PHYLIP PROML begin =========\n\n" );    
836     # Five arguments:
837     # 1. name of alignment file (in correct format!)
838     # 2. number of bootstraps
839     # 3. jumbles: 0: do not jumble; >=1 number of jumbles
840     # 4. seed for random number generator
841     # 5. 1 for PAM instead of JTT
842     my $use_pam = 1;
843     if ( $matrix == 0 ) {
844         $use_pam = 0;
845     }
846     &executeProml( $input, $bootstraps, $jumbles, $seed, $use_pam, $use_global_rearr );
847     print( "\n========== PHYLIP PROML end =========\n\n" );
848     &mv( "outtree", $OUTTREE_PROML );
849     &rm( "outfile" ); 
850     if ( $bootstraps > 1 ) {
851         &append( $OUTTREE_PROML, $OUTTREES_ALL );
852         $all_count++;
853     }
854 }
855
856
857 if ( $use_protpars == 1 ) {
858     my $input = "";
859     if ( $bootstraps > 1 ) {
860         $input = $SEQBOOT_OUTFILE;
861     }
862     else {
863        $input = "align";
864     }    
865     print( "\n========== PHYLIP PROTPARS begin =========\n\n" );    
866     &executeProtpars( $input, $bootstraps, $jumbles, $seed );
867     print( "\n========== PHYLIP PROTPARS end =========\n\n" );
868     &mv( "outtree", $OUTTREE_PROTPARS );
869     &rm( $outfile."_protpars_outfile" ); 
870     &mv( "outfile", $outfile."_protpars_outfile" );
871     if ( $bootstraps > 1 ) {
872         &append( $OUTTREE_PROTPARS, $OUTTREES_ALL );
873         $all_count++;
874     }
875 }
876
877
878
879 # Methods based on PWD
880 # --------------------
881 my $OUTTREE_FASTME    = "outtree_fastme";
882 my $OUTTREE_PHYLIP_NJ = "outtree_phylip_nj";
883 my $OUTTREE_PHYLIP_FM = "outtree_phylip_fm";
884 my $OUTTREE_PHYLIP_ME = "outtree_phylip_me";
885 my $OUTTREE_BIONJ     = "outtree_bionj";
886 my $OUTTREE_WEIGHBOR  = "outtree_weighbor";
887
888 my $CONSENSUS_FASTME    = "consensus_fastme";
889 my $CONSENSUS_PHYLIP_NJ = "consensus_phylip_nj";
890 my $CONSENSUS_PHYLIP_FM = "consensus_phylip_fm";
891 my $CONSENSUS_PHYLIP_ME = "consensus_phylip_me";
892 my $CONSENSUS_BIONJ     = "consensus_bionj";
893 my $CONSENSUS_WEIGHBOR  = "consensus_weighbor";
894 my $CONSENSUS_ALL       = "consensus_all";
895
896
897 if ( $use_fastme == 1 ) {
898     print( "\n========== FASTME begin =========\n\n" );
899     &executeFastme( $pwdfile, $bootstraps, $fastme_init_tree_opt );
900     print( "\n========== FASTME end ===========\n\n" );
901     &rm( "output.d" );
902     &mv( "output.tre", $OUTTREE_FASTME );
903     if ( $bootstraps > 1 ) {
904         &append( $OUTTREE_FASTME, $OUTTREES_ALL ); 
905         $all_count++;
906     }
907 }
908 if ( $use_phylip_nj ) {
909     print( "\n========== PHYLIP NEIGHBOR begin =========\n\n" );
910     &executeNeighbor( $pwdfile, $bootstraps, $seed, 0 );
911     print( "\n========== PHYLIP NEIGHBOR end =========\n\n" );
912     &mv( "outtree", $OUTTREE_PHYLIP_NJ );
913     &rm( "outfile" );
914     if ( $bootstraps > 1 ) {
915         &append( $OUTTREE_PHYLIP_NJ, $OUTTREES_ALL );
916         $all_count++;
917     }
918 }
919 if ( $use_phylip_fitch_fm ) {
920     print( "\n========== PHYLIP FITCH FM begin =========\n\n" );
921     &executeFitch( $pwdfile, $bootstraps, $seed, $jumbles, 0, "FM", $use_global_rearr );
922     print( "\n========== PHYLIP FITCH FM end =========\n\n" );
923     &mv( "outtree", $OUTTREE_PHYLIP_FM );
924     &rm( "outfile" );
925     if ( $bootstraps > 1 ) {
926         &append( $OUTTREE_PHYLIP_FM, $OUTTREES_ALL );
927         $all_count++;
928     }    
929 }
930 if (  $use_phylip_fitch_me ) {
931     print( "\n========== PHYLIP FITCH ME begin =========\n\n" );
932     &executeFitch( $pwdfile, $bootstraps, $seed, $jumbles, 0, "ME", $use_global_rearr );
933     print( "\n========== PHYLIP FITCH ME end =========\n\n" );
934     &mv( "outtree", $OUTTREE_PHYLIP_ME );
935     &rm( "outfile" );
936     if ( $bootstraps > 1 ) {
937         &append( $OUTTREE_PHYLIP_ME, $OUTTREES_ALL );
938         $all_count++;
939     } 
940 }
941 if ( $use_bionj ) {
942     print( "\n========== BIONJ begin =========\n\n" );
943     &executeBionj( $pwdfile, $OUTTREE_BIONJ );
944     print( "\n========== BIONJ end =========\n\n" );
945     if ( $bootstraps > 1 ) {
946         &append( $OUTTREE_BIONJ, $OUTTREES_ALL );
947         $all_count++;
948     }    
949 }
950 if ( $use_weighbor ) {
951     print( "\n========== WEIGHBOR begin =========\n\n" );
952     &executeWeighbor( $number_of_aa, 14, $pwdfile, $OUTTREE_WEIGHBOR );
953     print( "\n========== WEIGHBOR  end =========\n\n" );
954     if ( $bootstraps > 1 ) {
955         &append( $OUTTREE_WEIGHBOR, $OUTTREES_ALL );
956         $all_count++;
957     }
958 }
959
960
961
962 if ( $bootstraps > 1 ) {
963     # Consense:
964     if ( $use_fastme == 1 ) {
965         &consense( $OUTTREE_FASTME, $CONSENSUS_FASTME );
966     }
967     if ( $use_phylip_nj == 1 ) {
968         &consense( $OUTTREE_PHYLIP_NJ, $CONSENSUS_PHYLIP_NJ );
969     }
970     if ( $use_phylip_fitch_fm == 1 ) {
971         &consense( $OUTTREE_PHYLIP_FM, $CONSENSUS_PHYLIP_FM );
972     }
973     if ( $use_phylip_fitch_me == 1 ) {
974         &consense( $OUTTREE_PHYLIP_ME, $CONSENSUS_PHYLIP_ME );
975     }
976     if ( $use_bionj == 1 ) {
977         &consense( $OUTTREE_BIONJ, $CONSENSUS_BIONJ );
978     }   
979     if ( $use_weighbor == 1 ) {
980         &consense( $OUTTREE_WEIGHBOR, $CONSENSUS_WEIGHBOR );
981     }
982     if ( $use_phyml == 1 ) {
983         &consense( $OUTTREE_PHYML, $CONSENSUS_PHYML );
984     } 
985     if ( $use_proml == 1 ) {
986         &consense( $OUTTREE_PROML, $CONSENSUS_PROML );
987     } 
988     if ( $use_protpars == 1 ) {
989         &consense( $OUTTREE_PROTPARS, $CONSENSUS_PROTPARS );
990     } 
991     if ( $all_count > 1 ) {
992         &consense( $OUTTREES_ALL, $CONSENSUS_ALL );
993     }
994     else {
995         &rm( $OUTTREES_ALL );
996     }
997    
998     my $INTREE_FOR_PUZZLE = "intree"; #why so serious?
999     &rm( $INTREE_FOR_PUZZLE );
1000     system( "touch", $INTREE_FOR_PUZZLE )
1001     && die("\n\n$0: could not \"touch $INTREE_FOR_PUZZLE\": $!\n\n");
1002      
1003     if ( $use_fastme == 1 ) {
1004         &append( $CONSENSUS_FASTME, $INTREE_FOR_PUZZLE );
1005     }
1006     if ( $use_phylip_nj == 1 ) {
1007         &append( $CONSENSUS_PHYLIP_NJ, $INTREE_FOR_PUZZLE );
1008     }
1009     if ( $use_phylip_fitch_fm == 1 ) {
1010         &append( $CONSENSUS_PHYLIP_FM, $INTREE_FOR_PUZZLE );
1011     }
1012     if ( $use_phylip_fitch_me == 1 ) {
1013         &append( $CONSENSUS_PHYLIP_ME, $INTREE_FOR_PUZZLE );
1014     }
1015     if ( $use_bionj == 1 ) {
1016         &append( $CONSENSUS_BIONJ, $INTREE_FOR_PUZZLE );
1017     }
1018     if ( $use_weighbor == 1 ) {
1019         &append( $CONSENSUS_WEIGHBOR, $INTREE_FOR_PUZZLE );
1020     }
1021     if ( $use_raxml == 1 ) {
1022         # Needed, because TREE-PUZZLE adds internal labels for all subsequent trees
1023         # when evaluating given trees (this seems a strange behaviour).
1024         removeSupportValues( $CONSENSUS_RAXML, $CONSENSUS_RAXML."_support_removed" );
1025         &append( $CONSENSUS_RAXML."_support_removed", $INTREE_FOR_PUZZLE );
1026         &rm( $CONSENSUS_RAXML."_support_removed" );
1027     }
1028     if ( $use_phyml == 1 ) {
1029         &append( $CONSENSUS_PHYML, $INTREE_FOR_PUZZLE );
1030     }
1031     if ( $use_proml == 1 ) {
1032         &append( $CONSENSUS_PROML, $INTREE_FOR_PUZZLE );
1033     }
1034     if ( $use_protpars == 1 ) {
1035         &append( $CONSENSUS_PROTPARS, $INTREE_FOR_PUZZLE );
1036     }
1037     if ( $all_count > 1 ) {
1038         &append( $CONSENSUS_ALL, $INTREE_FOR_PUZZLE );
1039     }
1040     
1041
1042     # Puzzle for ML branch lenghts:
1043     # The alignment is read from infile by default.
1044     # The tree is read from intree by default.
1045     &rm( "infile" );
1046     &mv( "align", "infile" ); # align = original alignment in phylip interleaved.
1047     
1048     &executePuzzleToCalculateBranchLenghts( $matrix,
1049                                             $exact_parameter_est,
1050                                             $rate_heterogeneity );
1051
1052     my $OUTTREE_PUZZLE = "outtree_puzzle";
1053  
1054     &rm( $outfile."_puzzle_outfile" ); 
1055    
1056     &mv( "outfile", $outfile."_puzzle_outfile" );
1057     &mv( "outtree", $OUTTREE_PUZZLE );
1058     &rm( "outdist" );
1059     &rm( "intree" );
1060
1061
1062     # Transfer
1063     # --------
1064     my $counter = 0;
1065     if ( $use_fastme == 1 ) {
1066         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_FASTME, $fastme_outtree, $counter++ );
1067         &rm( $CONSENSUS_FASTME );
1068     }
1069     if ( $use_phylip_nj == 1 ) {
1070         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PHYLIP_NJ, $phylip_nj_outtree, $counter++ );
1071         &rm( $CONSENSUS_PHYLIP_NJ );
1072     }
1073     if ( $use_phylip_fitch_fm == 1 ) {
1074         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PHYLIP_FM, $phylip_fm_outtree, $counter++ );
1075         &rm( $CONSENSUS_PHYLIP_FM );
1076     }
1077     if ( $use_phylip_fitch_me == 1 ) {
1078         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PHYLIP_ME, $phylip_me_outtree, $counter++ );
1079         &rm( $CONSENSUS_PHYLIP_ME );
1080     }
1081     if ( $use_bionj == 1 ) {
1082         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_BIONJ, $bionj_outtree, $counter++ );
1083         &rm( $CONSENSUS_BIONJ );
1084     }
1085     if ( $use_weighbor == 1 ) {
1086         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_WEIGHBOR, $weighbor_outtree, $counter++ );
1087         &rm( $CONSENSUS_WEIGHBOR );
1088     }
1089     if ( $use_raxml == 1 ) {
1090         &to_phyloxml( $CONSENSUS_RAXML, $raxml_outtree, 1, 1 );
1091         $counter++;
1092     }
1093     if ( $use_phyml == 1 ) {
1094         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PHYML, $phyml_outtree, $counter++ );
1095         &rm( $CONSENSUS_PHYML );
1096     }
1097     if ( $use_proml == 1 ) {
1098         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PROML, $proml_outtree, $counter++ );
1099         &rm( $CONSENSUS_PROML );
1100     } 
1101     if ( $use_protpars == 1 ) {
1102         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PROTPARS, $protpars_outtree, $counter++ );
1103         &rm( $CONSENSUS_PROTPARS );
1104     } 
1105     if ( $all_count > 1 ) {
1106         &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_ALL, $all_outtree, $counter++ );
1107         &rm( $CONSENSUS_ALL );
1108     }
1109     
1110     # Clean up
1111     # --------
1112     &rm( $OUTTREE_PUZZLE );
1113     &rm( $SEQBOOT_OUTFILE );
1114     if ( $keep_multiple_trees == 1 ) {
1115         if ( $use_fastme == 1 ) {
1116             &mv( $OUTTREE_FASTME, $multitreefile_fastme );
1117         }
1118         if ( $use_phylip_nj == 1 ) {
1119             &mv( $OUTTREE_PHYLIP_NJ, $multitreefile_phylip_nj );
1120         }
1121         if ( $use_phylip_fitch_fm == 1 ) {
1122             &mv( $OUTTREE_PHYLIP_FM, $multitreefile_phylip_fm );
1123         }
1124         if ( $use_phylip_fitch_me == 1 ) {
1125             &mv( $OUTTREE_PHYLIP_ME, $multitreefile_phylip_me );
1126         }
1127         if ( $use_bionj == 1 ) {
1128             &mv( $OUTTREE_BIONJ, $multitreefile_bionj );
1129         }
1130         if ( $use_weighbor == 1 ) {
1131             &mv( $OUTTREE_WEIGHBOR, $multitreefile_weighbor );
1132         }
1133         if ( $use_phyml == 1 ) {
1134             &mv( $OUTTREE_PHYML, $multitreefile_phyml );
1135         }
1136         if ( $use_proml == 1 ) {
1137             &mv( $OUTTREE_PROML, $multitreefile_proml );
1138         }
1139         if ( $use_protpars == 1 ) {
1140             &mv( $OUTTREE_PROTPARS, $multitreefile_protpars );
1141         }
1142         &mv( $pwdfile, $multipwdfile );
1143     }
1144     else {
1145         if ( $use_fastme == 1 ) {
1146             &rm( $OUTTREE_FASTME );
1147         }
1148         if ( $use_phylip_nj == 1 ) {
1149             &rm( $OUTTREE_PHYLIP_NJ );
1150         }
1151         if ( $use_phylip_fitch_fm == 1 ) {
1152             &rm( $OUTTREE_PHYLIP_FM );
1153         }
1154         if ( $use_phylip_fitch_me == 1 ) {
1155             &rm( $OUTTREE_PHYLIP_ME );
1156         }
1157         if ( $use_bionj == 1 ) {
1158             &rm( $OUTTREE_BIONJ );
1159         }
1160         if ( $use_weighbor == 1 ) {
1161             &rm( $OUTTREE_WEIGHBOR );
1162         }
1163         if ( $use_phyml == 1 ) {
1164             &rm( $OUTTREE_PHYML );
1165         }
1166         if ( $use_proml == 1 ) {
1167             &rm( $OUTTREE_PROML );
1168         }
1169         if ( $use_protpars == 1 ) {
1170             &rm( $OUTTREE_PROTPARS );
1171         }
1172         &rm( $pwdfile );
1173     }
1174     if ( $all_count > 1 ) {
1175         &rm( $OUTTREES_ALL );
1176     }    
1177 } # if ( $bootstraps > 1 )
1178 else {
1179     &rm( "infile.dist" );
1180    
1181     &rm( "infile.puzzle" );
1182     if ( $use_fastme == 1 ) {
1183         &to_phyloxml( $OUTTREE_FASTME, $fastme_outtree, 0, 1 );
1184     }
1185     if ( $use_phylip_nj == 1 ) {
1186         &to_phyloxml( $OUTTREE_PHYLIP_NJ, $phylip_nj_outtree, 0, 1);
1187     }
1188     if ( $use_phylip_fitch_fm == 1 ) {
1189         &to_phyloxml( $OUTTREE_PHYLIP_FM, $phylip_fm_outtree, 0, 1 );
1190     }
1191     if ( $use_phylip_fitch_me == 1 ) {
1192         &to_phyloxml( $OUTTREE_PHYLIP_ME, $phylip_me_outtree, 0, 1 );
1193     }
1194     if ( $use_bionj == 1 ) {
1195         &to_phyloxml( $OUTTREE_BIONJ, $bionj_outtree, 0, 1 );
1196     }
1197     if ( $use_weighbor == 1 ) {
1198         &to_phyloxml( $OUTTREE_WEIGHBOR, $weighbor_outtree, 0, 1 );
1199     }
1200     if ( $use_raxml == 1 ) {
1201         &to_phyloxml( $OUTTREE_RAXML, $raxml_outtree, 0, 1 );
1202     }
1203     if ( $use_phyml == 1 ) {
1204         &to_phyloxml( $OUTTREE_PHYML, $phyml_outtree, 0, 1 );
1205     }
1206     if ( $use_proml == 1 ) {
1207         &to_phyloxml( $OUTTREE_PROML, $proml_outtree, 0, 1 );
1208     }
1209     if ( $use_protpars == 1 ) {
1210         &to_phyloxml( $OUTTREE_PROTPARS, $protpars_outtree, 0, 1 );
1211     }
1212 } # if ( $bootstraps > 1 )
1213
1214 &rm( $infile );
1215 &rm( "infile" );
1216 &rm( "align" );
1217 &rm( "align.reduced" );
1218
1219
1220 $log = $log."Finish time/date                    : ".`date`;
1221
1222 if ( $bootstraps > 1 ) {
1223     $log = $log."Puzzle output file                  : ".$outfile."_puzzle_outfile\n";
1224 }
1225 $log = $log."Columns in alignment                : $number_of_aa\n";
1226 $log = $log."Number of sequences in alignment    : $number_of_seqs\n";
1227 if ( $all_count > 1 ) {
1228     $log = $log."Combined consensus                  : $all_outtree\n";
1229
1230
1231
1232 if ( $bootstraps > 1 ) {
1233     $log = $log."\n\n";
1234     $log = $log."Simple support value statistics (trees are numbered the same as for TREE PUZZLE output)\n";
1235     $log = $log."------------------------------- \n";
1236     $log = $log."\n";
1237 }    
1238
1239 open( OUT, ">$logfile" ) || die "\n$0: Cannot create file <<$logfile>>: $!\n";
1240 print OUT $log;
1241 close( OUT );
1242
1243 if ( $bootstraps > 1 ) {
1244     # Simple support statistics
1245     # -------------------------
1246     my $SS_OUT = $temp_dir."/ss_out";
1247     my @phylos = ();
1248     my $ounter = 0;
1249     if ( $use_fastme == 1 ) {
1250         $phylos[ $ounter++ ] = $fastme_outtree;
1251     }
1252     if ( $use_phylip_nj == 1 ) {
1253         $phylos[ $ounter++ ] = $phylip_nj_outtree;
1254     }
1255     if ( $use_phylip_fitch_fm == 1 ) {
1256         $phylos[ $ounter++ ] = $phylip_fm_outtree;
1257     }
1258     if ( $use_phylip_fitch_me == 1 ) {
1259         $phylos[ $ounter++ ] = $phylip_me_outtree;
1260     }
1261     if ( $use_bionj == 1 ) {
1262         $phylos[ $ounter++ ] = $bionj_outtree;
1263     }
1264     if ( $use_weighbor == 1 ) {
1265         $phylos[ $ounter++ ] = $weighbor_outtree;
1266     }
1267     if ( $use_raxml == 1 ) {
1268         $phylos[ $ounter++ ] = $raxml_outtree;
1269     }
1270     if ( $use_phyml == 1 ) {
1271         $phylos[ $ounter++ ] = $phyml_outtree;
1272     }
1273     if ( $use_proml == 1 ) {
1274         $phylos[ $ounter++ ] = $proml_outtree;
1275     }
1276     if ( $use_protpars == 1 ) {
1277         $phylos[ $ounter++ ] = $protpars_outtree;
1278     }
1279     if ( $all_count > 1) {
1280         $phylos[ $ounter++ ] = $all_outtree;
1281     }
1282     &executeSupportStatistics( $SS_OUT, @phylos );
1283     &append( $SS_OUT, $logfile );
1284     &rm( $SS_OUT );
1285     
1286     # Append parts of puzzle output file
1287     # ----------------------------------
1288     if ( $all_count > 1 ) {
1289         &parsePuzzleOutfile( $outfile."_puzzle_outfile", $logfile );
1290     }
1291 }
1292
1293 chdir( $current_dir ) 
1294 || die "\n\n$0: Could not chdir to <<$current_dir>>: $!\n\n";
1295
1296 rmdir( $temp_dir )
1297 || print "\n\n$0: Warning: Could not remove <<$temp_dir>>: $!\n\n";
1298
1299 print "\n\n\n$0 successfully comleted.\n\n";
1300
1301 exit( 0 ); 
1302     
1303
1304
1305 # Methods:
1306 # --------
1307
1308
1309 # Six arguments:
1310 # 1. DNA or Amino-Acids sequence filename (PHYLIP format)
1311 # 2. Model, eg. PROTGAMMAIVT
1312 # 3. Replicates (bootstrap)
1313 # 4. Seed for bootstrap
1314 # 5. Output suffix
1315 # 6. Algorithm (only for bootstrap, default otherwise)
1316 # NOTE. RaxML does its own bootstrapping.
1317 sub executeRaxml {
1318     my $msa            = $_[ 0 ]; 
1319     my $model          = $_[ 1 ]; 
1320     my $replicates     = $_[ 2 ]; 
1321     my $seed           = $_[ 3 ];  
1322     my $outfile_suffix = $_[ 4 ];
1323     my $algo           = $_[ 5 ];
1324     
1325     &testForTextFilePresence( $msa );
1326     my $command = "$RAXML -m $model -s $msa -n $outfile_suffix";
1327       
1328     if ( $replicates > 1 ) {
1329         $command = $command . " -x $seed -N $replicates";
1330         if ( $algo ) {
1331             $command = $command . " -f $algo";
1332         }
1333     }
1334       
1335     print( "\n$command\n");  
1336       
1337     system( $command )
1338     && &dieWithUnexpectedError( $command );
1339     
1340
1341
1342
1343 sub to_phyloxml {
1344     my $from = $_[ 0 ];
1345     my $to   = $_[ 1 ];
1346     my $internal_names_are_boots = $_[ 2 ];
1347     my $extract_taxonomy = $_[ 3 ];
1348     &dieIfFileExists( $to );
1349     &dieIfFileNotExists( $from );
1350     my $command = "$NEWICK_TO_PHYLOXML -f=nn $from $to";
1351     if ( $internal_names_are_boots == 1 ) {
1352         $command = $command . " -i";
1353     }
1354     if ( $extract_taxonomy  == 1 ) {
1355         $command = $command . " -xt";
1356     }
1357     system( $command  )
1358     && die "$0: Could not execute \"$command \"";
1359     &rm( $from );
1360 }
1361
1362
1363 sub mv {
1364     my $from = $_[ 0 ];
1365     my $to   = $_[ 1 ];
1366     &dieIfFileExists( $to );
1367     &dieIfFileNotExists( $from );
1368     system( "mv", $from, $to )
1369     && die "\n\n$0: could not move \"$from\" to \"$to\": $!\n\n";
1370 }
1371
1372 sub cp {
1373     my $from = $_[ 0 ];
1374     my $to   = $_[ 1 ];
1375     &dieIfFileExists( $to );
1376     &dieIfFileNotExists( $from );
1377    
1378     system( "cp", $from, $to )
1379     && die "\n\n$0: could not copy \"$from\" to \"$to\": $!\n\n";
1380 }
1381
1382 sub rm {
1383     my $f = $_[ 0 ];
1384     unlink( $f );
1385 }
1386
1387 sub consense {
1388     my $multi_in     = $_[ 0 ]; 
1389     my $consense_out = $_[ 1 ];
1390     &executeConsense( $multi_in );
1391     &mv( "outtree", $consense_out );
1392     &rm( "outfile" );
1393     
1394 }    
1395
1396
1397
1398 # 1. file to be appended
1399 # 2. file to append to
1400 sub append {
1401     my $to_be_appended = $_[ 0 ];
1402     my $append_to      = $_[ 1 ];
1403     &dieIfFileNotExists( $to_be_appended );
1404     system( "cat $to_be_appended >> $append_to" )
1405     && die "\n\n$0: could not execute \"cat $to_be_appended >> $append_to\": $!\n\n";
1406     
1407 }
1408
1409 sub dieIfFileExists {
1410     my $file = $_[ 0 ]; 
1411     if ( -e $file ) {
1412         die "\n\n$0: \"$file\" already exists\n\n";
1413     }
1414
1415
1416 sub dieIfFileNotExists {
1417     my $file = $_[ 0 ]; 
1418     unless ( ( -s $file ) && ( -f $file ) ) {
1419         die( "\n\n$0: \"$file\" does not exist or is empty" );
1420     }
1421
1422
1423
1424
1425
1426 # Two arguments:
1427 # 1. seed for random number generator
1428 # 2. number of bootstraps
1429 # Reads in "infile" by default.
1430 sub executeSeqboot {
1431
1432     my $s    = $_[ 0 ];
1433     my $bs   = $_[ 1 ];
1434     my $verb = "";
1435     
1436     &testForTextFilePresence( $infile );
1437
1438    
1439     $verb = "
1440 2";
1441     
1442     system( "$SEQBOOT << !
1443 r
1444 $bs$verb
1445 Y
1446 $s
1447 !" )
1448     && die "$0: Could not execute \"$SEQBOOT\"";
1449   
1450 }
1451
1452
1453
1454 # One/two/three argument(s):
1455 # Reads in tree from "intree" by default. (Presence of "intree" automatically 
1456 # switches into "User defined trees" mode.)
1457 # 1. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
1458 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
1459 # 2. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
1460 # 3. Model of rate heterogeneity:
1461 #    1 for "8 Gamma distributed rates"
1462 #    2 for "Two rates (1 invariable + 1 variable)"
1463 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
1464 #    otherwise: Uniform rate
1465 # Last modified: 09/08/03 (added 2nd and 3rd parameter)
1466 sub executePuzzleToCalculateBranchLenghts {
1467     my $matrix_option              = $_[ 0 ];
1468     my $parameter_estimates_option = $_[ 1 ];
1469     my $rate_heterogeneity_option  = $_[ 2 ];
1470     my $i             = 0;
1471     my $mat           = "";
1472     my $est           = "";
1473     my $rate          = "";
1474     
1475     unless ( ( -s "infile" ) && ( -f "infile" ) && ( -T "infile" ) ) {
1476         die "\n$0: executePuzzleToCalculateBranchLenghts: <<infile>> does not exist, is empty, or is not a plain textfile.\n";
1477     }
1478     unless ( ( -s "intree" ) && ( -f "intree" ) && ( -T "intree" ) ) {
1479         die "\n$0: executePuzzleToCalculateBranchLenghts: <<intree>> does not exist, is empty, or is not a plain textfile.\n";
1480     }
1481
1482     $mat = setModelForPuzzle( $matrix_option );
1483     if ( $parameter_estimates_option ) {
1484         $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
1485     }
1486     if ( $rate_heterogeneity_option ) {
1487         $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
1488     }
1489   
1490     
1491     system( "$PUZZLE << !
1492 $mat$est$rate
1493 x
1494 y
1495 !" )
1496     && die "$0: Could not execute \"$PUZZLE\" (mat=$mat est=$est rate=$rate)";
1497     
1498 }
1499
1500 # Two arguments:
1501 # 1. puzzle outfile
1502 # 2. file to append to
1503 sub parsePuzzleOutfile {
1504     my $puzzle_outfile    = $_[ 0 ];
1505     my $file_to_append_to = $_[ 1 ];
1506     &testForTextFilePresence( $puzzle_outfile );
1507     open( OUT, ">>$file_to_append_to" ) || &dieWithUnexpectedError( "Cannot open \"$file_to_append_to\"" );
1508     open( IN, "$puzzle_outfile" ) || &dieWithUnexpectedError( "Cannot open file \"$puzzle_outfile\"" );
1509     my $return_line;
1510     my $read = 0;
1511     print OUT "\nTREE PUZZLE output\n";
1512     print OUT "------------------\n";
1513     while ( $return_line = <IN> ) {
1514         if ( $return_line =~/COMPARISON OF USER TREES/ ) {
1515             $read = 1;
1516         }                         
1517         elsif( $return_line =~/TIME STAMP/  ) {
1518             $read = 0;
1519         }
1520         elsif( $read ) {
1521             print OUT $return_line;
1522         }
1523     }
1524     close( IN );
1525     close( OUT );
1526 }   
1527
1528 # Three/four arguments:
1529 # 1. Name of file containing tree with correct branch lengths
1530 # 2. Name of file containing tree with correct bootstraps
1531 # 3. Outputfilename
1532 # 4. Index of tree with correct branch lengths, in case more than one in file
1533 # Last modified: 2007.11.27
1534 sub executeSupportTransfer {
1535     my $tree_with_bl = $_[ 0 ];
1536     my $tree_with_bs = $_[ 1 ];
1537     my $out          = $_[ 2 ];
1538     my $index        = $_[ 3 ];
1539   
1540     &testForTextFilePresence( $tree_with_bl );
1541     &testForTextFilePresence( $tree_with_bs );
1542     my $command = "$SUPPORT_TRANSFER $tree_with_bl $tree_with_bs $out $index";
1543     system( $command )
1544     && die "$0: Could not execute \"$command\"";
1545 }
1546
1547
1548 # Two or more arguments:
1549 # 1. outfile
1550 # 2. phylogeny 1 with support values 
1551 # 3. phylogeny 2 with support values 
1552 # 4. ...
1553 sub executeSupportStatistics {
1554     my $outfile      = $_[ 0 ];
1555     &dieIfFileExists( $outfile );
1556     my $phylos = "";
1557     for( my $i = 1; $i < scalar(@_); ++$i ) {
1558         &testForTextFilePresence( $_[ $i ] );
1559         $phylos .= $_[ $i ]." ";
1560     }    
1561     my $command = "$SUPPORT_STATISTICS -o=$outfile $phylos";
1562     system( "$command" )
1563     && die "$0: Could not execute \"$command\"";
1564 }
1565
1566
1567 sub getNumberOfSeqsAndAas { 
1568     my $infile = $_[ 0 ];
1569     my $seqs = 0;
1570     my $aa   = 0;
1571     open( IN, "$infile" ) || die "\n$0: Cannot open file <<$infile>>: $!\n";
1572     while( <IN> ) { 
1573         if ( $_ =~ /^\s*(\d+)\s+(\d+)\s*$/ ) { 
1574             $seqs = $1;
1575             $aa   = $2;
1576         } 
1577     }
1578     close( IN );
1579     
1580     if (  $seqs == 0 ||  $aa  == 0 ) {
1581         die( "\n$0: Could not get number of seqs and aa from: $infile" );
1582     }
1583     return $seqs, $aa;
1584 }
1585
1586
1587
1588 sub removeSupportValues {
1589     my $infile  = $_[ 0 ];
1590     my $outfile = $_[ 1 ];
1591     &testForTextFilePresence( $infile );
1592     open( OUT, ">$outfile" ) || &dieWithUnexpectedError( "Cannot create file \"$outfile\"" );
1593     open( IN, "$infile" ) || &dieWithUnexpectedError( "Cannot open file \"$infile\"" );
1594     while ( my $line = <IN> ) {
1595         $line =~ s/\)\d+\.?\d*:/\):/g;
1596         print OUT "$line";
1597     }
1598     close( OUT );
1599     close( IN );   
1600 }
1601
1602
1603
1604
1605 # Six arguments:
1606 # 1. name of alignment file (in correct format!)
1607 # 2. number of bootstraps
1608 # 3. jumbles: 0: do not jumble; >=1 number of jumbles
1609 # 4. seed for random number generator
1610 # 5. 1 for PAM instead of JTT
1611 # 6. 1 to use globale rearragements
1612 sub executeProml {
1613     my $align            = $_[ 0 ];
1614     my $bs               = $_[ 1 ];
1615     my $rand             = $_[ 2 ];
1616     my $s                = $_[ 3 ];
1617     my $use_pam          = $_[ 4 ];
1618     my $use_global_rearr = $_[ 5 ];
1619     my $jumble = "";
1620     my $multi  = "";
1621     my $pam    = ""; 
1622    
1623     &testForTextFilePresence( $align );
1624
1625     if ( $bs > 1 && $rand < 1 ) {
1626         $rand = 1;
1627     }
1628
1629     if ( $rand >= 1 ) {
1630         $jumble = "
1631 J
1632 $s
1633 $rand"; 
1634     }
1635    
1636     if (  $bs > 1 ) {
1637         $multi = "
1638 M
1639 D
1640 $bs";
1641     }
1642    
1643     
1644    if ( $use_pam == 1 ) {
1645         $pam = "
1646 P
1647 P";
1648     }
1649     
1650     my $global = "";
1651     if ( $use_global_rearr == 1 ) { 
1652         $global = "
1653 G"; 
1654     }
1655
1656     system( "$PROML  2>&1 << !
1657 $align$jumble$multi$pam$global
1658 3
1659 Y
1660 !" )
1661     && &dieWithUnexpectedError( "Could not execute \"$PROML $align$jumble$multi$pam$global\"" );
1662     # 3: Do NOT print out tree
1663       
1664     return;
1665
1666 } ## executeProml
1667
1668
1669 sub printUsage {
1670
1671     print <<END;
1672
1673 Copyright (C) 2002-2007 Christian M. Zmasek
1674 All rights reserved
1675
1676 Author: Christian M. Zmasek
1677 phylosoft\@gmail.com
1678 http://www.phylosoft.org
1679
1680   Requirements  phylo_pl is part of the FORESTER collection of programs.
1681   ------------  Many of its global variables are set via forester.pm.
1682
1683   Note. Use xt.pl (for Pfam alignments) or mt.pl (for other alignments) 
1684   to run phylo_pl.pl on whole directories of alignments files. 
1685
1686   Usage
1687   -----
1688
1689       phylo_pl.pl [-options] <input alignment in SELEX (Pfam), PHYLIP
1690       sequential format, or Clustal W output> <outputfile>
1691       [path/name for temporary directory to be created]
1692
1693      Example:
1694      "% phylo_pl.pl -B100q\@1nbS9X IL5.aln IL5_tree"
1695
1696   Options
1697   -------
1698  
1699   Bx : Number of bootstraps. B0: do not bootstrap. Default is 100 bootstrapps.
1700        The number of bootstrapps should be divisible by 10.
1701   J  : Use JTT matrix (Jones et al. 1992) in TREE-PUZZLE and/or PHYML, RAXML, default: VT (Mueller-Vingron 2000).
1702   L  : Use BLOSUM 62 matrix (Henikoff-Henikoff 92) in TREE-PUZZLE and/or PHYML, RAXML, default: VT.
1703   M  : Use mtREV24 matrix (Adachi-Hasegawa 1996) in TREE-PUZZLE and/or PHYML, default: VT.
1704   W  : Use WAG matrix (Whelan-Goldman 2000) in TREE-PUZZLE and/or PHYML, RAXML, default: VT.
1705   P  : Use PAM matrix (Dayhoff et al. 1978) in TREE-PUZZLE and/or PHYML, RAXML, default: VT.
1706   D  : Use DCMut matrix (Kosial and Goldman, 2005) in PHYML, RAXML, VT in TREE-PUZZLE.
1707   A  : Let TREE-PUZZLE choose which matrix to use, default: VT
1708   E  : Exact parameter estimates in TREE-PUZZLE, default: Approximate.
1709        Model of rate heterogeneity in TREE-PUZZLE (default: Uniform rate):
1710   g  : 8 Gamma distributed rates
1711   t  : Two rates (1 invariable + 1 variable)
1712   m  : Mixed (1 invariable + 8 Gamma rates)
1713   q\@x: Use FastME, x: 1: GME
1714                       2: BME
1715                       3: NJ
1716   n  : Use PHYLIP Neighbor (NJ).                    
1717   f  : Use PHYLIP Fitch.
1718   e  : Use PHYLIP Minimal Evolution.
1719   b  : Use BIONJ.
1720   w  : Use Weighbor.
1721   x  : Use RAxML.
1722   y  : Use PHYML. 
1723   o  : Use PHYLIP proml. 
1724   p  : Use PHYLIP protpars.
1725   rx : Number of relative substitution rate categories in PHYML (default is 4).
1726   jx : Number of jumbles (input order randomization) for PHYLIP FM, ME, PROTPARS, and PROML (default is 2) (random seed set with Sx).
1727   I  : Estimate proportion of invariable sites in RAXML and/or PHYML (otherwise, proportion "0.0" is used in PHYML)
1728   G  : to turn on global rearrangements in PHYLIP FM, ME, and PROML
1729   Sx : Seed for random number generator(s). Must be 4n+1. Default is 9.
1730   X  : To keep multiple tree file (=trees from bootstrap resampled alignments) and 
1731        pairwise distance matrix file (in case of bootstrap analysis).
1732   
1733 END
1734   
1735 } ## printUsage