source: branches/profile/PERL_SCRIPTS/ARBTOOLS/raxml2arb.pl

Last change on this file was 12553, checked in by westram, 5 years ago
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 9.8 KB
Line 
1#!/usr/bin/perl
2# =============================================================== #
3#                                                                 #
4#   File      : raxml2arb.pl                                      #
5#   Purpose   : Import XX best of YY calculated raxml trees into  #
6#               ARB and generate comment containing likelyhood    #
7#               and content of info file                          #
8#                                                                 #
9#   Coded by Ralf Westram (coder@reallysoft.de) in March 2008     #
10#   Institute of Microbiology (Technical University Munich)       #
11#   http://www.arb-home.de/                                       #
12#                                                                 #
13# =============================================================== #
14
15use strict;
16use warnings;
17
18my $MODE_NORMAL       = 0;
19my $MODE_BOOTSTRAPPED = 1;
20my $MODE_OPTIMIZED    = 2;
21
22sub arb_message($) {
23  my ($msg) = @_;
24  print "Message: $msg\n";
25  system("arb_message \"$msg\"");
26}
27
28sub error($) {
29  my ($msg) = @_;
30  if ($msg =~ /at.*\.pl/o) {
31    my ($pe,$loc) = ($`,$&.$');
32    chomp $loc;
33    $msg = $pe."\n(".$loc.")\nPlease check console for additional error reasons";
34  }
35  arb_message($msg);
36  die $msg;
37}
38
39my $loaded_log = undef;
40my @log = ();
41
42sub loadLog($) {
43  my ($log) = @_;
44  if (defined $loaded_log) {
45    if ($log ne $loaded_log) {
46      $loaded_log = undef;
47      @log = ();
48    }
49  }
50  if (not defined $loaded_log) {
51    open(LOG,'<'.$log) || die "can't load '$log' (Reason: $!)";
52    my $line;
53    while (defined($line=<LOG>)) {
54      chomp($line);
55      push @log, $line;
56    }
57    close(LOG);
58  }
59}
60
61sub firstLogLineMatching($$) {
62  my ($log,$regexp) = @_;
63  loadLog($log);
64  foreach my $line (@log) {
65    if ($line =~ $regexp) { return ($line,$1); }
66  }
67  return (undef,undef);
68}
69
70sub raxml_filename($$$) {
71  my ($type,$name,$run) = @_;
72  my $fname = 'RAxML_'.$type.'.'.$name;
73  if (defined $run) { $fname .= '.RUN.'.$run; }
74  return $fname;
75}
76
77sub someWhat($$) {
78  my ($count,$what) = @_;
79  if ($count==1) { return '1 '.$what; }
80  return $count.' '.$what.'s';
81}
82
83my @splitted_trees = ();
84
85sub split_treefile($$$) {
86  my ($in,$out,$treesExpected) = @_;
87  @splitted_trees = ();
88  my $outcount = 0;
89  open(IN,'<'.$in) || die "can't read '$in' (Reason: $!)";
90  my $line;
91  while (defined($line=<IN>)) {
92    if ($line =~ /\(/o) {
93      my $outname = "$out.$outcount"; $outcount++;
94      open(OUT,'>'.$outname) || die "can't write '$outname' (Reason: $!)";
95      print OUT $line;
96      close(OUT);
97      push @splitted_trees, $outname;
98    }
99    else {
100      print "Unexpected line in '$in':\n";
101      print $line;
102    }
103  }
104  close(IN);
105  my $treesFound = scalar(@splitted_trees);
106  if ($treesFound!=$treesExpected) {
107    die "Failed to split '$in' into single treefiles\n".
108    "(expected to find $treesExpected trees, found $treesFound)";
109  }
110}
111
112sub treeInfo_normal($$) {
113  my ($name,$run) = @_;
114
115  my $result       = raxml_filename('result',$name,$run);
116  my $bipartitions = raxml_filename('bipartitions',$name,$run);
117  my $parsimony    = raxml_filename('parsimonyTree',$name,$run);
118  my $log          = raxml_filename('log',$name,$run);
119
120  if (not -f $result) { # no result tree, try bipartitions or parsimonyTree
121    if (-f $bipartitions) { return ($bipartitions,'unknown'); }
122    if (-f $parsimony) { return ($parsimony,'unknown'); }
123  }
124
125  my ($line) = firstLogLineMatching($log,qr/./); # first non-empty line
126  if (not $line =~ / (.*)$/o) {
127    die "can't parse likelyhood from '$log'";
128  }
129
130  my $likelyhood = $1;
131  return ($result,$likelyhood);
132}
133
134sub treeInfo_bootstrapped($$) {
135  my ($name,$run) = @_;
136
137  my $info = raxml_filename('info',$name,undef);
138  if (defined $run) {
139    my $treeCount = scalar(@splitted_trees);
140    if ($run>=$treeCount) {
141      die "Invalid run number $run - has to be in [0 .. ".($treeCount-1)."]";
142    }
143    my ($line,$likelyhood) = firstLogLineMatching($info,qr/^Bootstrap\[$run\]:\s.*\slikelihood\s+([^\s,]+),/);
144    if (not defined $likelyhood) {
145      die "Failed to parse likelyhood for 'Bootstrap[$run]' from '$info'";
146    }
147    return ($splitted_trees[$run],$likelyhood);
148  }
149  else {
150    my $bestTree = raxml_filename('bipartitions',$name,undef);
151    my ($line,$likelyhood) = firstLogLineMatching($info,qr/^Final ML Optimization Likelihood:\s+(.*)$/);
152    if (not defined $likelyhood) {
153      arb_message("Failed to extract final likelyhood from '$info'");
154      $likelyhood = 'unknown';
155    }
156    return ($bestTree,$likelyhood);
157  }
158}
159
160
161
162sub treeInfo($$$) {
163  my ($mode,$name,$run) = @_;
164
165  if ($mode==$MODE_BOOTSTRAPPED)                     { return treeInfo_bootstrapped($name,$run); }
166  if ($mode==$MODE_NORMAL or $mode==$MODE_OPTIMIZED) { return treeInfo_normal($name,$run); }
167
168  die "Unhandled mode '$mode'";
169}
170
171sub findTrees($$$$\%) {
172  my ($name,$mode,$runs,$take,$likelyhood_r) = @_;
173
174  %$likelyhood_r = ();
175
176  my $singleTree = 0;
177  if    ($mode==$MODE_NORMAL)       { if ($runs==1) { $singleTree = 1; } }
178  elsif ($mode==$MODE_BOOTSTRAPPED) {
179    if ($take==1) { $singleTree = 1; }
180    else { arb_message("Note: to import raxml bootstrap tree, set 'Select ## best trees' to '1'"); }
181  }
182  elsif ($mode==$MODE_OPTIMIZED)    { $singleTree = 1; }
183  else { die; }
184
185  if ($singleTree==1) {
186    my ($tree,$likely) = treeInfo($mode,$name,undef);
187    $$likelyhood_r{$tree} = $likely;
188  }
189  else {
190    if ($mode==$MODE_BOOTSTRAPPED) {
191      my $raxml_out   = raxml_filename('bootstrap',$name,undef);
192      split_treefile($raxml_out, 'raxml2arb.tree',$runs);
193      print "Splitted '$raxml_out' into ".scalar(@splitted_trees)." treefiles\n";
194    }
195    for (my $r = 0; $r<$runs; $r++) {
196      my ($tree,$likely) = treeInfo($mode,$name,$r);
197      $$likelyhood_r{$tree} = $likely;
198    }
199  }
200}
201
202sub main() {
203  eval {
204    my $args = scalar(@ARGV);
205
206    if ($args != 6) {
207      die "Usage: raxml2arb.pl RUNNAME NUMBEROFRUNS [normal|bootstrapped|optimized] TAKETREES [import|consense] \"COMMENT\"\n".
208        "       import: Import the best TAKETREES trees of NUMBEROFRUNS generated trees into ARB\n".
209        "       consense: Create and import consensus tree of best TAKETREES trees of NUMBEROFRUNS generated trees\n";
210    }
211
212    my ($RUNNAME,$NUMBEROFRUNS,$MODE,$TAKETREES,$CONSENSE,$COMMENT) = @ARGV;
213
214    if ($NUMBEROFRUNS<1) { die "NUMBEROFRUNS has to be 1 or higher (NUMBEROFRUNS=$NUMBEROFRUNS)"; }
215
216    my %likelyhood  = (); # key=treefile, value=likelyhood
217    my @treesToTake = (); # treefiles
218
219    my $mode = $MODE_NORMAL;
220    if ($MODE eq 'bootstrapped') { $mode = $MODE_BOOTSTRAPPED; }
221    elsif ($MODE eq 'optimized') { $mode = $MODE_OPTIMIZED; }
222    elsif ($MODE ne 'normal') { die "Unexpected mode '$MODE' (accepted: 'normal', 'bootstrapped', 'optimized')"; }
223
224    my $calc_consense = 0;
225    if ($CONSENSE eq 'consense') { $calc_consense = 1; }
226    elsif ($CONSENSE ne 'import') { die "Unknown value '$CONSENSE' (expected 'import' or 'consense')"; }
227
228    findTrees($RUNNAME,$mode,$NUMBEROFRUNS,$TAKETREES,%likelyhood);
229
230    my $createdTrees = scalar(keys %likelyhood);
231    print "Found ".someWhat($createdTrees,'tree').":\n";
232
233    if ($TAKETREES > $createdTrees) {
234      arb_message("Cant take more trees (=$TAKETREES) than calculated (=$createdTrees) .. using all");
235      $TAKETREES = $createdTrees;
236    }
237
238    my @sortedTrees = sort { $likelyhood{$b} <=> $likelyhood{$a}; } keys %likelyhood;
239    foreach (@sortedTrees) { print "  $_ = ".$likelyhood{$_}."\n"; }
240
241    @treesToTake = splice(@sortedTrees,0,$TAKETREES);
242
243    my $treesToTake = scalar(@treesToTake);
244    if ($treesToTake==0) { die "failed to detect RAxML output trees"; }
245
246    if ($calc_consense and $treesToTake<2) {
247      arb_message("Need to take at least 2 trees to create a consensus tree - importing..");
248      $calc_consense = 0;
249    }
250
251    my $treename = 'tree_RAxML_'.$$;
252    my $infofile = 'RAxML_info.'.$RUNNAME;
253
254    if ($calc_consense==0) {
255      print 'Importing '.someWhat($treesToTake,'tree').":\n";
256      my $count = undef;
257      if ($treesToTake>1) { $count = 1; }
258
259      foreach (@treesToTake) {
260        print "  $_ = ".$likelyhood{$_}."\n";
261        my $currTreename = $treename;
262        if (defined $count) {
263          $currTreename .= '_'.$count;
264          $count++;
265        }
266        my $command = 'arb_read_tree '.$currTreename.' '.$_.' "'.$COMMENT."\n".'likelyhood='.$likelyhood{$_}.'"';
267        if (-f $infofile) { $command .= ' -commentFromFile '.$infofile; }
268
269        system($command)==0 || die "can't execute '$command' (Reason: $?)";
270      }
271    }
272    else { # calc consense tree
273      my $minLH  = undef;
274      my $maxLH  = undef;
275      my $meanLH = 0;
276
277      print 'Consensing '.someWhat($treesToTake,'tree').":\n";
278      open(INTREE,'>intree') || die "can't write 'intree' (Reason: $!)";
279      foreach (@treesToTake) {
280        my $LH = $likelyhood{$_};
281        print "  $_ = $LH\n";
282        if (not defined $minLH) {
283          $minLH = $maxLH = $LH;
284        }
285        else {
286          if ($LH < $minLH) { $minLH = $LH; }
287          if ($LH > $maxLH) { $maxLH = $LH; }
288        }
289        $meanLH += $LH;
290
291        open(RAXMLTREE,'<'.$_) || die "can't read '$_' (Reason: $!)";
292        foreach (<RAXMLTREE>) {
293          print INTREE $_;
294        }
295        close(RAXMLTREE);
296        print INTREE "\n";
297      }
298      close(INTREE);
299
300      $meanLH /= $treesToTake;
301
302      my $command = 'arb_echo y | consense';
303      system($command)==0 || die "can't execute '$command' (Reason: $?)";
304
305      if (not -f 'outtree') {
306        die "Consense failed (no 'outtree' generated)";
307      }
308
309      my $comment = "$COMMENT\nConsensus tree of $treesToTake trees\nLikelyhood: min=$minLH mean=$meanLH max=$maxLH";
310      $command = 'arb_read_tree -consense '.$treesToTake.' '.$treename.' outtree "'.$comment.'"';
311      if (-f $infofile) { $command .= ' -commentFromFile '.$infofile; }
312      system($command)==0 || die "can't execute '$command' (Reason: $?)";
313    }
314  };
315  if ($@) { error($@); }
316}
317main();
Note: See TracBrowser for help on using the repository browser.