[Bio] / FigKernelPackages / gjonewicklib.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/gjonewicklib.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.23 - (view) (download) (as text)

1 : olson 1.17 # This is a SAS component.
2 :    
3 : olson 1.2 #
4 : golsen 1.20 # Copyright (c) 2003-2010 University of Chicago and Fellowship
5 : olson 1.2 # for Interpretations of Genomes. All Rights Reserved.
6 :     #
7 :     # This file is part of the SEED Toolkit.
8 : golsen 1.9 #
9 : olson 1.2 # The SEED Toolkit is free software. You can redistribute
10 :     # it and/or modify it under the terms of the SEED Toolkit
11 : golsen 1.9 # Public License.
12 : olson 1.2 #
13 :     # You should have received a copy of the SEED Toolkit Public License
14 :     # along with this program; if not write to the University of Chicago
15 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
16 :     # Genomes at veronika@thefig.info or download a copy from
17 :     # http://www.theseed.org/LICENSE.TXT.
18 :     #
19 :    
20 : golsen 1.1 package gjonewicklib;
21 :    
22 :     #===============================================================================
23 :     # perl functions for dealing with trees
24 :     #
25 :     # Usage:
26 :     # use gjonewicklib
27 :     #
28 :     #
29 :     #===============================================================================
30 :     # Tree data structures:
31 :     #===============================================================================
32 :     #
33 :     # Elements in newick text file are:
34 :     #
35 :     # [c1] ( desc_1, desc_2, ... ) [c2] label [c3] : [c4] x [c5]
36 :     #
37 :     # Note that:
38 : golsen 1.9 #
39 : golsen 1.1 # Comment list 1 can exist on any subtree, but its association with
40 :     # tree components can be upset by rerooting
41 :     # Comment list 2 cannot exist without a descendant list
42 :     # Comment list 3 cannot exist without a label
43 :     # Comment lists 4 and 5 cannot exist without a branch length
44 :     #
45 :     # Elements in perl representation are:
46 :     #
47 :     # $tree = \@rootnode;
48 :     #
49 : overbeek 1.7 # $node = [ \@desc, # reference to list of descendants
50 : golsen 1.1 # $label, # node label
51 :     # $x, # branch length
52 :     # \@c1, # reference to comment list 1
53 :     # \@c2, # reference to comment list 2
54 :     # \@c3, # reference to comment list 3
55 :     # \@c4, # reference to comment list 4
56 :     # \@c5 # reference to comment list 5
57 : overbeek 1.7 # ]
58 : golsen 1.1 #
59 :     # At present, no routine tests or enforces the length of the list (a single
60 :     # element list could be a valid internal node).
61 :     #
62 :     # All empty lists can be [] or undef
63 :     #
64 :     # Putting the comments at the end allows a shorter list nearly all the
65 :     # time, but is different from the prolog representation.
66 :     #
67 :     #
68 : overbeek 1.7 # Ross Overbeek has a different tree node structure:
69 :     #
70 :     # $node = [ Label,
71 :     # DistanceToParent,
72 :     # [ ParentPointer, ChildPointer1, ... ],
73 :     # [ Name1\tVal1, Name2\tVal2, ... ]
74 :     # ]
75 :     #
76 :     # So:
77 :     #
78 :     #===============================================================================
79 :     # Tree format interconversion:
80 :     #===============================================================================
81 :     #
82 : golsen 1.9 # $bool = is_overbeek_tree( $tree )
83 :     # $bool = is_gjonewick_tree( $tree )
84 :     #
85 : overbeek 1.7 # $gjonewick = overbeek_to_gjonewick( $overbeek )
86 :     # $overbeek = gjonewick_to_overbeek( $gjonewick )
87 :     #
88 : golsen 1.1 #===============================================================================
89 :     # Tree data extraction:
90 :     #===============================================================================
91 :     #
92 :     # $listref = newick_desc_ref( $noderef )
93 :     # $label = newick_lbl( $noderef )
94 :     # $x = newick_x( $noderef )
95 :     # $listref = newick_c1( $noderef )
96 :     # $listref = newick_c2( $noderef )
97 :     # $listref = newick_c3( $noderef )
98 :     # $listref = newick_c4( $noderef )
99 :     # $listref = newick_c5( $noderef )
100 :     #
101 :     # @desclist = newick_desc_list( $noderef )
102 :     # $n = newick_n_desc( $noderef )
103 :     # $descref = newick_desc_i( $noderef, $i ) # 1-based numbering
104 :     # $bool = newick_is_tip( $noderef )
105 :     #
106 :     # set_newick_desc_ref( $noderef, $listref )
107 :     # set_newick_lbl( $noderef, $label )
108 :     # set_newick_x( $noderef, $x )
109 :     # set_newick_c1( $noderef, $listref )
110 :     # set_newick_c2( $noderef, $listref )
111 :     # set_newick_c3( $noderef, $listref )
112 :     # set_newick_c4( $noderef, $listref )
113 :     # set_newick_c5( $noderef, $listref )
114 :     # set_newick_desc_list( $noderef, @desclist )
115 : golsen 1.9 # set_newick_desc_i( $noderef1, $i, $noderef2 ) # 1-based numbering
116 : golsen 1.8 #
117 :     # $bool = newick_is_valid( $noderef ) # verify that tree is valid
118 : golsen 1.1 #
119 :     # $bool = newick_is_rooted( $noderef ) # 2 branches from root
120 :     # $bool = newick_is_unrooted( $noderef ) # 3 or more branches from root
121 :     # $bool = newick_is_tip_rooted( $noderef ) # 1 branch from root
122 :     # $bool = newick_is_bifurcating( $noderef )
123 :     #
124 :     # $n = newick_tip_count( $noderef )
125 :     # @tiprefs = newick_tip_ref_list( $noderef )
126 : golsen 1.9 # \@tiprefs = newick_tip_ref_list( $noderef )
127 : golsen 1.1 # @tips = newick_tip_list( $noderef )
128 : golsen 1.9 # \@tips = newick_tip_list( $noderef )
129 :     #
130 : golsen 1.1 # $tipref = newick_first_tip_ref( $noderef )
131 :     # $tip = newick_first_tip( $noderef )
132 : golsen 1.9 #
133 : golsen 1.1 # @tips = newick_duplicated_tips( $noderef )
134 : golsen 1.9 # \@tips = newick_duplicated_tips( $noderef )
135 :     #
136 : golsen 1.1 # $bool = newick_tip_in_tree( $noderef, $tipname )
137 : golsen 1.9 #
138 : golsen 1.1 # @tips = newick_shared_tips( $tree1, $tree2 )
139 : golsen 1.9 # \@tips = newick_shared_tips( $tree1, $tree2 )
140 : golsen 1.1 #
141 :     # $length = newick_tree_length( $noderef )
142 : golsen 1.9 #
143 :     # %tip_distances = newick_tip_distances( $noderef )
144 :     # \%tip_distances = newick_tip_distances( $noderef )
145 :     #
146 : golsen 1.1 # $xmax = newick_max_X( $noderef )
147 :     # ( $tipref, $xmax ) = newick_most_distant_tip_ref( $noderef )
148 : golsen 1.9 # ( $tipname, $xmax ) = newick_most_distant_tip_name( $noderef )
149 : golsen 1.1 #
150 : golsen 1.20 # Provide a standard name by which two trees can be compared for same topology
151 :     #
152 :     # $stdname = std_tree_name( $tree )
153 :     #
154 : golsen 1.8 # Tree tip insertion point (tip is on branch of length x that
155 :     # is inserted into branch connecting node1 and node2, a distance
156 :     # x1 from node1 and x2 from node2):
157 :     #
158 : golsen 1.9 # [ $node1, $x1, $node2, $x2, $x ] = newick_tip_insertion_point( $tree, $tip )
159 : golsen 1.8 #
160 :     # Standardized label for a node in terms of intersection of 3 lowest sorting
161 :     # tips (sort is lower case):
162 :     #
163 : golsen 1.9 # @TipOrTips = std_node_name( $tree, $node )
164 : golsen 1.8 #
165 : golsen 1.1 #-------------------------------------------------------------------------------
166 :     # Paths from root of tree:
167 :     #-------------------------------------------------------------------------------
168 :     #
169 :     # Path descriptions are of form:
170 :     # ( $node0, $i0, $node1, $i1, $node2, $i2, ..., $nodeN )
171 :     # () is returned upon failure
172 :     #
173 :     # @path = path_to_tip( $treenode, $tipname )
174 :     # @path = path_to_named_node( $treenode, $nodename )
175 :     # @path = path_to_node_ref( $treenode, $noderef )
176 :     #
177 :     # @path = path_to_node( $node, $tip1, $tip2, $tip3 ) # 3 tip names
178 :     # @path = path_to_node( $node, [ $tip1, $tip2, $tip3 ] ) # Array of tips
179 :     # @path = path_to_node( $node, $tip1 ) # Use path_to_tip
180 :     # @path = path_to_node( $node, [ $tip1 ] ) # Use path_to_tip
181 :     #
182 :     # $distance = newick_path_length( @path )
183 :     # $distance = tip_to_tip_distance( $tree, $tip1, $tip2 )
184 :     # $distance = node_to_node_distance( $tree, $node1, $node2 )
185 :     #
186 :     #
187 :     #===============================================================================
188 :     # Tree manipulations:
189 :     #===============================================================================
190 :     #
191 : overbeek 1.4 # $treecopy = copy_newick_tree( $tree )
192 : golsen 1.1 #
193 :     #-------------------------------------------------------------------------------
194 : golsen 1.8 # The following modify the existing tree, and possibly any components of that
195 : golsen 1.1 # tree that are reached by reference. If the old version is still needed, copy
196 :     # before modifying.
197 :     #-------------------------------------------------------------------------------
198 :     #
199 :     # Modify labels:
200 :     #
201 :     # $newtree = newick_relabel_nodes( $node, \%new_name )
202 :     # $newtree = newick_relabel_nodes_i( $node, \%new_name )
203 :     # $newtree = newick_relabel_tips( $node, \%new_name )
204 :     # $newtree = newick_relabel_tips_i( $node, \%new_name )
205 :     #
206 :     # Modify branches:
207 :     #
208 :     # $n_changed = newick_set_undefined_branches( $node, $x )
209 :     # $n_changed = newick_set_all_branches( $node, $x )
210 : golsen 1.5 # $n_changed = newick_fix_negative_branches( $tree )
211 : overbeek 1.7 # $node = newick_rescale_branches( $node, $factor )
212 : golsen 1.15 # $node = newick_modify_branches( $node, \&function )
213 :     # $node = newick_modify_branches( $node, \&function, \@func_parms )
214 : golsen 1.1 #
215 : golsen 1.8 # Modify comments:
216 :     #
217 :     # $node = newick_strip_comments( $node )
218 :     #
219 : golsen 1.1 # Modify rooting and/or order:
220 :     #
221 :     # $nrmtree = normalize_newick_tree( $tree )
222 :     # $revtree = reverse_newick_tree( $tree )
223 :     # $stdtree = std_unrooted_newick( $tree )
224 :     # $newtree = aesthetic_newick_tree( $tree, $direction )
225 :     # $rndtree = random_order_newick_tree( $tree )
226 :     # $newtree = reroot_newick_by_path( @path )
227 :     # $newtree = reroot_newick_to_tip( $tree, $tip )
228 :     # $newtree = reroot_newick_next_to_tip( $tree, $tip )
229 :     # $newtree = reroot_newick_to_node( $tree, @node )
230 :     # $newtree = reroot_newick_to_node_ref( $tree, $noderef )
231 : golsen 1.9 # $newtree = reroot_newick_between_nodes( $tree, $node1, $node2, $fraction )
232 : golsen 1.12 # $newtree = reroot_newick_to_midpoint( $tree ) # unweighted
233 :     # $newtree = reroot_newick_to_midpoint_w( $tree ) # weight by tips
234 : golsen 1.5 # $newtree = reroot_newick_to_approx_midpoint( $tree ) # unweighted
235 :     # $newtree = reroot_newick_to_approx_midpoint_w( $tree ) # weight by tips
236 : golsen 1.1 # $newtree = uproot_tip_rooted_newick( $tree )
237 :     # $newtree = uproot_newick( $tree )
238 :     #
239 :     # $newtree = prune_from_newick( $tree, $tip )
240 : golsen 1.12 # $newtree = rooted_newick_subtree( $tree, @tips )
241 :     # $newtree = rooted_newick_subtree( $tree, \@tips )
242 : golsen 1.1 # $newtree = newick_subtree( $tree, @tips )
243 :     # $newtree = newick_subtree( $tree, \@tips )
244 : golsen 1.12 # $newtree = newick_covering_subtree( $tree, @tips )
245 :     # $newtree = newick_covering_subtree( $tree, \@tips )
246 : golsen 1.1 #
247 : golsen 1.5 # $newtree = collapse_zero_length_branches( $tree )
248 :     #
249 : golsen 1.8 # $node = newick_insert_at_node( $node, $subtree )
250 :     # $tree = newick_insert_between_nodes( $tree, $subtree, $node1, $node2, $fraction )
251 :     #
252 : golsen 1.1 #===============================================================================
253 : golsen 1.9 # Tree neighborhood: subtree of n tips to represent a larger tree.
254 :     #===============================================================================
255 :     #
256 :     # Focus around root:
257 :     #
258 :     # $subtree = root_neighborhood_representative_tree( $tree, $n, \%tip_priority )
259 :     # $subtree = root_neighborhood_representative_tree( $tree, $n )
260 :     # @tips = root_neighborhood_representative_tips( $tree, $n, \%tip_priority )
261 :     # @tips = root_neighborhood_representative_tips( $tree, $n )
262 :     # \@tips = root_neighborhood_representative_tips( $tree, $n, \%tip_priority )
263 :     # \@tips = root_neighborhood_representative_tips( $tree, $n )
264 :     #
265 :     # Focus around a tip insertion point (the tip is not in the subtree):
266 :     #
267 :     # $subtree = tip_neighborhood_representative_tree( $tree, $tip, $n, \%tip_priority )
268 :     # $subtree = tip_neighborhood_representative_tree( $tree, $tip, $n )
269 :     # @tips = tip_neighborhood_representative_tips( $tree, $tip, $n, \%tip_priority )
270 :     # @tips = tip_neighborhood_representative_tips( $tree, $tip, $n )
271 :     # \@tips = tip_neighborhood_representative_tips( $tree, $tip, $n, \%tip_priority )
272 :     # \@tips = tip_neighborhood_representative_tips( $tree, $tip, $n )
273 :     #
274 :     #===============================================================================
275 : golsen 1.22 # Random trees
276 :     #===============================================================================
277 :     #
278 :     # $tree = random_equibranch_tree( @tips, \%options )
279 :     # $tree = random_equibranch_tree( \@tips, \%options )
280 :     # $tree = random_equibranch_tree( @tips )
281 :     # $tree = random_equibranch_tree( \@tips )
282 :     #
283 :     # $tree = random_ultrametric_tree( @tips, \%options )
284 :     # $tree = random_ultrametric_tree( \@tips, \%options )
285 :     # $tree = random_ultrametric_tree( @tips )
286 :     # $tree = random_ultrametric_tree( \@tips )
287 :     #
288 :     #===============================================================================
289 : golsen 1.1 # Tree reading and writing:
290 :     #===============================================================================
291 : golsen 1.9 # Write machine-readable trees:
292 : golsen 1.1 #
293 :     # writeNewickTree( $tree )
294 :     # writeNewickTree( $tree, $file )
295 : overbeek 1.7 # writeNewickTree( $tree, \*FH )
296 :     # fwriteNewickTree( $file, $tree ) # Matches the C arg list for f... I/O
297 : golsen 1.1 # $treestring = swriteNewickTree( $tree )
298 :     # $treestring = formatNewickTree( $tree )
299 : golsen 1.9 #
300 :     # Write human-readable trees:
301 :     #
302 : golsen 1.1 # @textlines = text_plot_newick( $node, $width, $min_dx, $dy )
303 :     # printer_plot_newick( $node, $file, $width, $min_dx, $dy )
304 :     #
305 : golsen 1.9 # Read trees:
306 :     #
307 : golsen 1.23 # $tree = read_newick_tree( )
308 :     # $tree = read_newick_tree( \*FH )
309 :     # $tree = read_newick_tree( $file )
310 :     #
311 :     # @trees = read_newick_trees( )
312 :     # @trees = read_newick_trees( \*FH )
313 :     # @trees = read_newick_trees( $file )
314 :     #
315 : overbeek 1.7 # $tree = parse_newick_tree_str( $string )
316 : golsen 1.1 #
317 :     #===============================================================================
318 :    
319 :    
320 : golsen 1.8 use Carp;
321 :     use Data::Dumper;
322 : golsen 1.9 use strict;
323 : golsen 1.8
324 : golsen 1.1 require Exporter;
325 :    
326 :     our @ISA = qw(Exporter);
327 :     our @EXPORT = qw(
328 : golsen 1.9 is_overbeek_tree
329 :     is_gjonewick_tree
330 : overbeek 1.7 overbeek_to_gjonewick
331 :     gjonewick_to_overbeek
332 : golsen 1.8 newick_is_valid
333 : golsen 1.1 newick_is_rooted
334 :     newick_is_unrooted
335 :     tree_rooted_on_tip
336 :     newick_is_bifurcating
337 :     newick_tip_count
338 : golsen 1.9 newick_tip_ref_list
339 : golsen 1.1 newick_tip_list
340 : golsen 1.9
341 : golsen 1.1 newick_first_tip
342 :     newick_duplicated_tips
343 :     newick_tip_in_tree
344 :     newick_shared_tips
345 :    
346 :     newick_tree_length
347 : golsen 1.9 newick_tip_distances
348 : golsen 1.1 newick_max_X
349 :     newick_most_distant_tip_ref
350 :     newick_most_distant_tip_name
351 : golsen 1.9
352 :     newick_tip_insertion_point
353 :    
354 : golsen 1.20 std_tree_name
355 : golsen 1.1
356 :     path_to_tip
357 :     path_to_named_node
358 :     path_to_node_ref
359 :     path_to_node
360 :    
361 :     newick_path_length
362 :     tip_to_tip_distance
363 :     node_to_node_distance
364 :    
365 : overbeek 1.4 copy_newick_tree
366 : golsen 1.1
367 :     newick_relabel_nodes
368 :     newick_relabel_nodes_i
369 :     newick_relabel_tips
370 :     newick_relabel_tips_i
371 :    
372 :     newick_set_undefined_branches
373 :     newick_set_all_branches
374 : golsen 1.5 newick_fix_negative_branches
375 : overbeek 1.7 newick_rescale_branches
376 : golsen 1.15 newick_modify_branches
377 : golsen 1.1
378 : golsen 1.8 newick_strip_comments
379 :    
380 : golsen 1.1 normalize_newick_tree
381 :     reverse_newick_tree
382 :     std_unrooted_newick
383 :     aesthetic_newick_tree
384 :     unaesthetic_newick_tree
385 :     random_order_newick_tree
386 :    
387 :     reroot_newick_by_path
388 :     reroot_newick_to_tip
389 :     reroot_newick_next_to_tip
390 :     reroot_newick_to_node
391 :     reroot_newick_to_node_ref
392 : golsen 1.9 reroot_newick_between_nodes
393 : golsen 1.12 reroot_newick_to_midpoint
394 :     reroot_newick_to_midpoint_w
395 : golsen 1.1 reroot_newick_to_approx_midpoint
396 : golsen 1.5 reroot_newick_to_approx_midpoint_w
397 : golsen 1.1 uproot_tip_rooted_newick
398 :     uproot_newick
399 :    
400 :     prune_from_newick
401 : golsen 1.12 rooted_newick_subtree
402 : golsen 1.1 newick_subtree
403 : golsen 1.12 newick_covering_subtree
404 : golsen 1.5 collapse_zero_length_branches
405 : golsen 1.1
406 : golsen 1.8 newick_insert_at_node
407 :     newick_insert_between_nodes
408 :    
409 : golsen 1.9 root_neighborhood_representative_tree
410 :     root_neighborhood_representative_tips
411 :     tip_neighborhood_representative_tree
412 :     tip_neighborhood_representative_tips
413 :    
414 : golsen 1.22 random_equibranch_tree
415 :     random_ultrametric_tree
416 :    
417 : golsen 1.1 writeNewickTree
418 :     fwriteNewickTree
419 :     strNewickTree
420 :     formatNewickTree
421 : overbeek 1.7
422 :     read_newick_tree
423 :     read_newick_trees
424 : golsen 1.1 parse_newick_tree_str
425 :    
426 :     printer_plot_newick
427 :     text_plot_newick
428 :     );
429 :    
430 :     our @EXPORT_OK = qw(
431 :     newick_desc_ref
432 :     newick_lbl
433 :     newick_x
434 :     newick_c1
435 :     newick_c2
436 :     newick_c3
437 :     newick_c4
438 :     newick_c5
439 :     newick_desc_list
440 :     newick_n_desc
441 :     newick_desc_i
442 :     newick_is_tip
443 :     newick_is_valid
444 :    
445 :     set_newick_desc_ref
446 :     set_newick_lbl
447 :     set_newick_x
448 :     set_newick_c1
449 :     set_newick_c2
450 :     set_newick_c3
451 :     set_newick_c4
452 :     set_newick_c5
453 :    
454 :     set_newick_desc_list
455 :     set_newick_desc_i
456 :    
457 :     add_to_newick_branch
458 :     dump_tree
459 :     );
460 :    
461 :    
462 :     #-------------------------------------------------------------------------------
463 :     # Internally used definitions
464 :     #-------------------------------------------------------------------------------
465 :    
466 : golsen 1.21 sub array_ref { $_[0] && ref( $_[0] ) eq 'ARRAY' }
467 :     sub hash_ref { $_[0] && ref( $_[0] ) eq 'HASH' }
468 : golsen 1.1
469 :    
470 :     #===============================================================================
471 : golsen 1.9 # Interconvert overbeek and gjonewick trees:
472 : overbeek 1.7 #===============================================================================
473 :    
474 : golsen 1.9 sub is_overbeek_tree { array_ref( $_[0] ) && array_ref( $_[0]->[2] ) }
475 :    
476 :     sub is_gjonewick_tree { array_ref( $_[0] ) && array_ref( $_[0]->[0] ) }
477 :    
478 : overbeek 1.7 sub overbeek_to_gjonewick
479 :     {
480 :     return () unless ref( $_[0] ) eq 'ARRAY';
481 :     my ( $lbl, $x, $desc ) = @{ $_[0] };
482 :     my ( undef, @desc ) = ( $desc && ref( $desc ) eq 'ARRAY' ) ? @$desc : ();
483 :     [ [ map { overbeek_to_gjonewick( $_ ) } @desc ], $lbl, $x ]
484 :     }
485 :    
486 :     sub gjonewick_to_overbeek
487 :     {
488 :     return () unless ref( $_[0] ) eq 'ARRAY';
489 :     my ( $desc, $lbl, $x ) = @{ $_[0] };
490 :     my @desc = ( $desc && ref( $desc ) eq 'ARRAY' ) ? @$desc : ();
491 :     my $parent = $_[1];
492 :     my $node = [ $lbl, $x, undef, [] ];
493 :     $node->[2] = [ $parent, map { gjonewick_to_overbeek( $_, $node ) } @desc ];
494 :     return $node;
495 :     }
496 :    
497 : golsen 1.9
498 : overbeek 1.7 #===============================================================================
499 : golsen 1.1 # Extract tree structure values:
500 :     #===============================================================================
501 :     #
502 :     # $listref = newick_desc_ref( $noderef )
503 :     # $string = newick_lbl( $noderef )
504 :     # $real = newick_x( $noderef )
505 :     # $listref = newick_c1( $noderef )
506 :     # $listref = newick_c2( $noderef )
507 :     # $listref = newick_c3( $noderef )
508 :     # $listref = newick_c4( $noderef )
509 :     # $listref = newick_c5( $noderef )
510 :     # @list = newick_desc_list( $noderef )
511 :     # $int = newick_n_desc( $noderef )
512 :     # $listref = newick_desc_i( $noderef )
513 :     # $bool = node_is_tip( $noderef )
514 :     # $bool = node_is_valid( $noderef )
515 :     #
516 :     #-------------------------------------------------------------------------------
517 :    
518 : golsen 1.12 sub newick_desc_ref { ref($_[0]) ? $_[0]->[0] : Carp::confess() } # = ${$_[0]}[0]
519 : golsen 1.8 sub newick_lbl { ref($_[0]) ? $_[0]->[1] : Carp::confess() }
520 : golsen 1.12 sub newick_x { ref($_[0]) ? $_[0]->[2] : Carp::confess() }
521 : golsen 1.10 sub newick_c1 { ref($_[0]) ? $_[0]->[3] : Carp::confess() }
522 : golsen 1.12 sub newick_c2 { ref($_[0]) ? $_[0]->[4] : Carp::confess() }
523 :     sub newick_c3 { ref($_[0]) ? $_[0]->[5] : Carp::confess() }
524 :     sub newick_c4 { ref($_[0]) ? $_[0]->[6] : Carp::confess() }
525 :     sub newick_c5 { ref($_[0]) ? $_[0]->[7] : Carp::confess() }
526 : golsen 1.1
527 :     sub newick_desc_list {
528 :     my $node = $_[0];
529 : golsen 1.21 array_ref( $node ) && array_ref( $node->[0] ) ? @{ $node->[0] } : ();
530 : golsen 1.1 }
531 :    
532 :     sub newick_n_desc {
533 :     my $node = $_[0];
534 : golsen 1.21 array_ref( $node ) && array_ref( $node->[0] ) ? scalar @{ $node->[0] } : 0;
535 : golsen 1.1 }
536 :    
537 :     sub newick_desc_i {
538 :     my ( $node, $i ) = @_;
539 : golsen 1.21 array_ref( $node ) && $i && array_ref( $node->[0] ) ? $node->[0]->[$i-1] : undef;
540 : golsen 1.1 }
541 :    
542 :     sub node_is_tip {
543 :     my $node = $_[0];
544 :     ! array_ref( $node ) ? undef : # Not a node ref
545 :     array_ref( $node->[0] ) ? @{ $node->[0] } == 0 : # Empty descend list?
546 :     1 ; # No descend list
547 :     }
548 :    
549 :     sub node_is_valid { # An array ref with nonempty descend list or a label
550 :     my $node = $_[0];
551 :     array_ref( $node ) && ( array_ref( $node->[0] ) && @{ $node->[0] }
552 :     || defined( $node->[1] )
553 :     )
554 :     }
555 :    
556 :    
557 :     #-------------------------------------------------------------------------------
558 :     # Set tree structure values
559 :     #-------------------------------------------------------------------------------
560 :    
561 :     sub set_newick_desc_ref { $_[0]->[0] = $_[1] }
562 :     sub set_newick_lbl { $_[0]->[1] = $_[1] }
563 :     sub set_newick_x { $_[0]->[2] = $_[1] }
564 :     sub set_newick_c1 { $_[0]->[3] = $_[1] }
565 :     sub set_newick_c2 { $_[0]->[4] = $_[1] }
566 :     sub set_newick_c3 { $_[0]->[5] = $_[1] }
567 :     sub set_newick_c4 { $_[0]->[6] = $_[1] }
568 :     sub set_newick_c5 { $_[0]->[7] = $_[1] }
569 :    
570 :     sub set_newick_desc_list {
571 :     my $node = shift;
572 :     array_ref( $node ) || return;
573 :     if ( array_ref( $node->[0] ) ) { @{ $node->[0] } = @_ }
574 :     else { $node->[0] = [ @_ ] }
575 :     }
576 :    
577 :     sub set_newick_desc_i {
578 :     my ( $node1, $i, $node2 ) = @_;
579 :     array_ref( $node1 ) && array_ref( $node2 ) || return;
580 :     if ( array_ref( $node1->[0] ) ) { $node1->[0]->[$i-1] = $node2 }
581 :     else { $node1->[0] = [ $node2 ] }
582 :     }
583 :    
584 :    
585 :     #===============================================================================
586 :     # Some tree property tests:
587 :     #===============================================================================
588 : golsen 1.8 # Tree is valid?
589 :     #
590 :     # $bool = newick_is_valid( $node, $verbose )
591 :     #-------------------------------------------------------------------------------
592 :     sub newick_is_valid
593 :     {
594 :     my $node = shift;
595 :    
596 :     if ( ! array_ref( $node ) )
597 :     {
598 :     print STDERR "Node is not array reference\n" if $_[0];
599 :     return 0;
600 :     }
601 :    
602 :     my @node = @$node;
603 :     if ( ! @node )
604 :     {
605 :     print STDERR "Node is empty array reference\n" if $_[0];
606 :     return 0;
607 :     }
608 :    
609 :     # Must have descendant or label:
610 :    
611 :     if ( ! ( array_ref( $node[0] ) && @{ $node[0] } ) && ! $node[2] )
612 :     {
613 :     print STDERR "Node has neither descendant nor label\n" if $_[0];
614 :     return 0;
615 :     }
616 :    
617 :     # If comments are present, they must be array references
618 :    
619 :     foreach ( ( @node > 3 ) ? @node[ 3 .. $#node ] : () )
620 :     {
621 :     if ( defined( $_ ) && ! array_ref( $_ ) )
622 :     {
623 :     print STDERR "Node has neither descendant or label\n" if $_[0];
624 :     return 0;
625 :     }
626 :     }
627 :    
628 :     # Inspect the descendants:
629 :    
630 :     foreach ( array_ref( $node[0] ) ? @{ $node[0] } : () )
631 :     {
632 :     newick_is_valid( $_, @_ ) || return 0
633 :     }
634 :    
635 :     return 1;
636 :     }
637 :    
638 :    
639 :     #-------------------------------------------------------------------------------
640 : golsen 1.1 # Tree is rooted (2 branches at root node)?
641 :     #
642 :     # $bool = newick_is_rooted( $node )
643 :     #-------------------------------------------------------------------------------
644 :     sub newick_is_rooted {
645 :     my $node = $_[0];
646 :     ! array_ref( $node ) ? undef : # Not a node ref
647 :     array_ref( $node->[0] ) ? @{ $node->[0] } == 2 : # 2 branches
648 :     0 ; # No descend list
649 :     }
650 :    
651 :    
652 :     #-------------------------------------------------------------------------------
653 :     # Tree is unrooted (> 2 branches at root node)?
654 :     #
655 :     # $bool = newick_is_unrooted( $node )
656 :     #-------------------------------------------------------------------------------
657 :     sub newick_is_unrooted {
658 :     my $node = $_[0];
659 :     ! array_ref( $node ) ? undef : # Not a node ref
660 :     array_ref( $node->[0] ) ? @{ $node->[0] } >= 3 : # Over 2 branches
661 :     0 ; # No descend list
662 :     }
663 :    
664 :    
665 :     #-------------------------------------------------------------------------------
666 :     # Tree is rooted on a tip (1 branch at root node)?
667 :     #
668 :     # $bool = newick_is_tip_rooted( $node )
669 :     #-------------------------------------------------------------------------------
670 :     sub newick_is_tip_rooted {
671 :     my $node = $_[0];
672 :     ! array_ref( $node ) ? undef : # Not a node ref
673 :     array_ref( $node->[0] ) ? @{ $node->[0] } == 1 : # 1 branch
674 :     0 ; # No descend list
675 :     }
676 :    
677 :     #===============================================================================
678 :     # Everything below this point refers to parts of the tree structure using
679 :     # only the routines above.
680 :     #===============================================================================
681 :     # Tree is bifurcating? If so, return number of descendents of root node.
682 :     #
683 :     # $n_desc = newick_is_bifurcating( $node )
684 :     #-------------------------------------------------------------------------------
685 :     sub newick_is_bifurcating {
686 :     my ( $node, $not_root ) = @_;
687 :     if ( ! array_ref( $node ) ) { return undef } # Bad arg
688 :    
689 :     my $n = newick_n_desc( $node );
690 :     $n == 0 && ! $not_root ? 0 :
691 :     $n == 1 && $not_root ? 0 :
692 :     $n == 3 && $not_root ? 0 :
693 :     $n > 3 ? 0 :
694 :     $n > 2 && ! newick_is_bifurcating( newick_desc_i( $node, 3, 1 ) ) ? 0 :
695 :     $n > 1 && ! newick_is_bifurcating( newick_desc_i( $node, 2, 1 ) ) ? 0 :
696 :     $n > 0 && ! newick_is_bifurcating( newick_desc_i( $node, 1, 1 ) ) ? 0 :
697 :     $n
698 :     }
699 :    
700 :    
701 :     #-------------------------------------------------------------------------------
702 :     # Number of tips:
703 :     #
704 :     # $n = newick_tip_count( $node )
705 :     #-------------------------------------------------------------------------------
706 :     sub newick_tip_count {
707 :     my ( $node, $not_root ) = @_;
708 :    
709 :     my $imax = newick_n_desc( $node );
710 :     if ( $imax < 1 ) { return 1 }
711 :    
712 :     # Special case for tree rooted on tip
713 :    
714 :     my $n = ( $imax == 1 && ( ! $not_root ) ) ? 1 : 0;
715 :    
716 :     foreach ( newick_desc_list( $node ) ) { $n += newick_tip_count( $_, 1 ) }
717 :    
718 :     $n;
719 :     }
720 :    
721 :    
722 :     #-------------------------------------------------------------------------------
723 :     # List of tip nodes:
724 :     #
725 : golsen 1.9 # @tips = newick_tip_ref_list( $noderef )
726 :     # \@tips = newick_tip_ref_list( $noderef )
727 : golsen 1.1 #-------------------------------------------------------------------------------
728 :     sub newick_tip_ref_list {
729 :     my ( $node, $not_root ) = @_;
730 :    
731 :     my $imax = newick_n_desc( $node );
732 :     if ( $imax < 1 ) { return $node }
733 :    
734 :     my @list = ();
735 :    
736 :     # Tree rooted on tip?
737 :     if ( $imax == 1 && ! $not_root && newick_lbl( $node ) ) { push @list, $node }
738 :    
739 :     foreach ( newick_desc_list( $node ) ) {
740 :     push @list, newick_tip_ref_list( $_, 1 );
741 :     }
742 :    
743 : golsen 1.9 wantarray ? @list : \@list;
744 : golsen 1.1 }
745 :    
746 :    
747 :     #-------------------------------------------------------------------------------
748 :     # List of tips:
749 :     #
750 :     # @tips = newick_tip_list( $node )
751 : golsen 1.9 # \@tips = newick_tip_list( $node )
752 : golsen 1.1 #-------------------------------------------------------------------------------
753 :     sub newick_tip_list {
754 : golsen 1.9 my @tips = map { newick_lbl( $_ ) } newick_tip_ref_list( $_[0] );
755 :     wantarray ? @tips : \@tips;
756 : golsen 1.1 }
757 :    
758 :    
759 :     #-------------------------------------------------------------------------------
760 :     # First tip node in tree:
761 :     #
762 :     # $tipref = newick_first_tip_ref( $node )
763 :     #-------------------------------------------------------------------------------
764 :     sub newick_first_tip_ref {
765 :     my ( $node, $not_root ) = @_;
766 :     valid_node( $node ) || return undef;
767 :    
768 :     # Arrived at tip, or start of a tip-rooted tree?
769 :     my $n = newick_n_desc( $node );
770 :     if ( ( $n < 1 ) || ( $n == 1 && ! $not_root ) ) { return $node }
771 :    
772 :     newick_first_tip_ref( newick_desc_i( $node, 1 ), 1 );
773 :     }
774 :    
775 :    
776 :     #-------------------------------------------------------------------------------
777 :     # First tip name in tree:
778 :     #
779 :     # $tip = newick_first_tip( $node )
780 :     #-------------------------------------------------------------------------------
781 :     sub newick_first_tip {
782 :     my ( $noderef ) = @_;
783 :    
784 :     my $tipref;
785 :     array_ref( $tipref = newick_first_tip_ref( $noderef ) ) ? newick_lbl( $tipref )
786 : golsen 1.9 : undef;
787 : golsen 1.1 }
788 :    
789 :    
790 :     #-------------------------------------------------------------------------------
791 :     # List of duplicated tip labels.
792 :     #
793 :     # @tips = newick_duplicated_tips( $node )
794 : golsen 1.9 # \@tips = newick_duplicated_tips( $node )
795 : golsen 1.1 #-------------------------------------------------------------------------------
796 :     sub newick_duplicated_tips {
797 : golsen 1.9 my @tips = &duplicates( newick_tip_list( $_[0] ) );
798 :     wantarray ? @tips : \@tips;
799 : golsen 1.1 }
800 :    
801 :    
802 :     #-------------------------------------------------------------------------------
803 :     # Tip in tree?
804 :     #
805 :     # $bool = newick_tip_in_tree( $node, $tipname )
806 :     #-------------------------------------------------------------------------------
807 :     sub newick_tip_in_tree {
808 :     my ( $node, $tip, $not_root ) = @_;
809 :    
810 :     my $n = newick_n_desc( $node );
811 :     if ( $n < 1 ) { return ( newick_lbl( $node ) eq $tip) ? 1 : 0 }
812 :    
813 :     # Special case for tree rooted on tip
814 :    
815 :     if ( $n == 1 && ( ! $not_root ) && newick_lbl( $node ) eq $tip ) { return 1 }
816 :    
817 :     foreach ( newick_desc_list( $node ) ) {
818 :     if ( newick_tip_in_tree( $_, $tip, 1 ) ) { return 1 }
819 :     }
820 :    
821 :     0; # Fall through means not found
822 :     }
823 :    
824 :    
825 :     #-------------------------------------------------------------------------------
826 :     # Tips shared between 2 trees.
827 :     #
828 :     # @tips = newick_shared_tips( $tree1, $tree2 )
829 : golsen 1.9 # \@tips = newick_shared_tips( $tree1, $tree2 )
830 : golsen 1.1 #-------------------------------------------------------------------------------
831 :     sub newick_shared_tips {
832 : golsen 1.9 my ( $tree1, $tree2 ) = @_;
833 :     my $tips1 = newick_tip_list( $tree1 );
834 :     my $tips2 = newick_tip_list( $tree2 );
835 :     my @tips = &intersection( $tips1, $tips2 );
836 :     wantarray ? @tips : \@tips;
837 : golsen 1.1 }
838 :    
839 :    
840 :     #-------------------------------------------------------------------------------
841 :     # Tree length.
842 :     #
843 :     # $length = newick_tree_length( $node )
844 :     #-------------------------------------------------------------------------------
845 :     sub newick_tree_length {
846 :     my ( $node, $not_root ) = @_;
847 :    
848 :     my $x = $not_root ? newick_x( $node ) : 0;
849 :     defined( $x ) || ( $x = 1 ); # Convert undefined to 1
850 :    
851 :     foreach ( newick_desc_list( $node ) ) { $x += newick_tree_length( $_, 1 ) }
852 :    
853 :     $x;
854 :     }
855 :    
856 :    
857 :     #-------------------------------------------------------------------------------
858 : golsen 1.9 # Hash of tip nodes and corresponding distances from root:
859 :     #
860 :     # %tip_distances = newick_tip_distances( $node )
861 :     # \%tip_distances = newick_tip_distances( $node )
862 :     #-------------------------------------------------------------------------------
863 :     sub newick_tip_distances
864 :     {
865 :     my ( $node, $x, $hash ) = @_;
866 :     my $root = ! $hash;
867 :     ref( $hash ) eq 'HASH' or $hash = {};
868 :    
869 :     $x ||= 0;
870 :     $x += newick_x( $node ) || 0;
871 :    
872 :     # Is it a tip?
873 :    
874 :     my $n_desc = newick_n_desc( $node );
875 :     if ( ! $n_desc )
876 :     {
877 :     $hash->{ newick_lbl( $node ) } = $x;
878 :     return $hash;
879 :     }
880 :    
881 :     # Tree rooted on tip?
882 :    
883 :     if ( ( $n_desc == 1 ) && $root && ( newick_lbl( $node ) ) )
884 :     {
885 :     $hash->{ newick_lbl( $node ) } = 0; # Distance to root is zero
886 :     }
887 :    
888 :     foreach ( newick_desc_list( $node ) )
889 :     {
890 :     newick_tip_distances( $_, $x, $hash );
891 :     }
892 :    
893 :     wantarray ? %$hash : $hash;
894 :     }
895 :    
896 :    
897 :     #-------------------------------------------------------------------------------
898 : golsen 1.1 # Tree max X.
899 :     #
900 :     # $xmax = newick_max_X( $node )
901 :     #-------------------------------------------------------------------------------
902 :     sub newick_max_X {
903 :     my ( $node, $not_root ) = @_;
904 :    
905 :     my $xmax = 0;
906 :     foreach ( newick_desc_list( $node ) ) {
907 :     my $x = newick_max_X( $_, 1 );
908 :     if ( $x > $xmax ) { $xmax = $x }
909 :     }
910 :    
911 :     my $x = $not_root ? newick_x( $node ) : 0;
912 :     $xmax + ( defined( $x ) ? $x : 1 ); # Convert undefined to 1
913 :     }
914 :    
915 :    
916 :     #-------------------------------------------------------------------------------
917 :     # Most distant tip from root: distance and path.
918 :     #
919 :     # ( $xmax, @path ) = newick_most_distant_tip_path( $tree )
920 :     #-------------------------------------------------------------------------------
921 :     sub newick_most_distant_tip_path {
922 :     my ( $node, $not_root ) = @_;
923 :    
924 :     my $imax = newick_n_desc( $node );
925 :     my $xmax = ( $imax > 0 ) ? -1 : 0;
926 :     my @pmax = ();
927 :     for ( my $i = 1; $i <= $imax; $i++ ) {
928 :     my ( $x, @path ) = newick_most_distant_tip_path( newick_desc_i( $node, $i ), 1 );
929 :     if ( $x > $xmax ) { $xmax = $x; @pmax = ( $i, @path ) }
930 :     }
931 :    
932 :     my $x = $not_root ? newick_x( $node ) : 0;
933 :     $xmax += defined( $x ) ? $x : 0; # Convert undefined to 1
934 :     ( $xmax, $node, @pmax );
935 :     }
936 :    
937 :    
938 :     #-------------------------------------------------------------------------------
939 :     # Most distant tip from root, and its distance.
940 :     #
941 :     # ( $tipref, $xmax ) = newick_most_distant_tip_ref( $tree )
942 :     #-------------------------------------------------------------------------------
943 :     sub newick_most_distant_tip_ref {
944 :     my ( $node, $not_root ) = @_;
945 :    
946 :     my $imax = newick_n_desc( $node );
947 :     my $xmax = ( $imax > 0 ) ? -1 : 0;
948 :     my $tmax = $node;
949 :     foreach ( newick_desc_list( $node ) ) {
950 :     my ( $t, $x ) = newick_most_distant_tip_ref( $_, 1 );
951 :     if ( $x > $xmax ) { $xmax = $x; $tmax = $t }
952 :     }
953 :    
954 :     my $x = $not_root ? newick_x( $node ) : 0;
955 :     $xmax += defined( $x ) ? $x : 1; # Convert undefined to 1
956 :     ( $tmax, $xmax );
957 :     }
958 :    
959 :    
960 :     #-------------------------------------------------------------------------------
961 :     # Name of most distant tip from root, and its distance.
962 :     #
963 :     # ( $tipname, $xmax ) = newick_most_distant_tip_name( $tree )
964 :     #-------------------------------------------------------------------------------
965 :     sub newick_most_distant_tip_name {
966 :     my ( $tipref, $xmax ) = newick_most_distant_tip_ref( $_[0] );
967 :     ( newick_lbl( $tipref ), $xmax )
968 :     }
969 :    
970 :    
971 :     #-------------------------------------------------------------------------------
972 : golsen 1.8 # Tree tip insertion point (with standard node labels):
973 :     #
974 :     # [ $node1, $x1, $node2, $x2, $x ]
975 :     # = newick_tip_insertion_point( $tree, $tip )
976 :     #
977 :     # Which means: tip is on a branch of length x that is inserted into the branch
978 :     # connecting node1 and node2, at distance x1 from node1 and x2 from node2.
979 :     #
980 :     # x1 +------ n1a (lowest sorting tip of this subtree)
981 :     # +--------n1
982 :     # | +------n1b (lowest sorting tip of this subtree)
983 :     # tip-------n
984 :     # x | +------------- n2a (lowest sorting tip of this subtree)
985 :     # +------n2
986 :     # x2 +-------- n2b (lowest sorting tip of this subtree)
987 :     #
988 :     # The designations of 1 vs 2, and a vs b are chosen such that:
989 :     # n1a < n1b, and n2a < n2b, and n1a < n2a
990 :     #
991 :     # Then the statandard description becomes:
992 :     #
993 :     # [ [ $n1a, min(n1b,n2a), max(n1b,n2a) ], x1,
994 :     # [ $n2a, min(n2b,n1a), max(n2b,n1a) ], x2,
995 :     # x
996 :     # ]
997 :     #
998 :     #-------------------------------------------------------------------------------
999 :     sub newick_tip_insertion_point
1000 :     {
1001 :     my ( $tree, $tip ) = @_;
1002 :     $tree && $tip && ref( $tree ) eq 'ARRAY' or return undef;
1003 :     $tree = copy_newick_tree( $tree ) or return undef;
1004 :     $tree = reroot_newick_to_tip( $tree, $tip ) or return undef;
1005 :     my $node = $tree;
1006 :    
1007 :     my $x = 0; # Distance to node
1008 :     my $dl = newick_desc_ref( $node ); # Descendent list of tip node;
1009 :     $node = $dl->[0]; # Node adjacent to tip
1010 :     $dl = newick_desc_ref( $node );
1011 :     while ( $dl && ( @$dl == 1 ) ) # Traverse unbranched nodes
1012 :     {
1013 :     $node = $dl->[0];
1014 :     $x += newick_x( $node );
1015 :     $dl = newick_desc_ref( $node );
1016 :     }
1017 :     $x += newick_x( $node );
1018 :    
1019 :     # We are now at the node that is the insertion point.
1020 :     # Is it a tip?
1021 :    
1022 :     my @description;
1023 :    
1024 :     if ( ( ! $dl ) || @$dl == 0 )
1025 :     {
1026 :     @description = ( [ newick_lbl( $node ) ], 0, undef, 0, $x );
1027 :     }
1028 :    
1029 :     # Is it a trifurcation or greater, in which case it does not go
1030 :     # away with tip deletion?
1031 :    
1032 :     elsif ( @$dl > 2 )
1033 :     {
1034 :     @description = ( [ std_node_name( $node, $node ) ], 0, undef, 0, $x );
1035 :     }
1036 :    
1037 :     # The node is bifurcating. We need to describe it.
1038 :    
1039 :     else
1040 :     {
1041 : golsen 1.9 my ( $n1, $x1 ) = describe_descendant( $dl->[0] );
1042 :     my ( $n2, $x2 ) = describe_descendant( $dl->[1] );
1043 : golsen 1.8
1044 :     if ( @$n1 == 2 ) { push @$n1, $n2->[0] }
1045 :     if ( @$n2 == 2 )
1046 :     {
1047 :     @$n2 = sort { lc $a cmp lc $b } ( @$n2, $n1->[0] );
1048 :     }
1049 :     if ( @$n1 == 3 ) { @$n2 = sort { lc $a cmp lc $b } @$n2 }
1050 :     @description = ( $n1, $x1, $n2, $x2, $x );
1051 :     }
1052 :    
1053 :     return wantarray ? @description : \@description;
1054 :     }
1055 :    
1056 :    
1057 : golsen 1.9 sub describe_descendant
1058 : golsen 1.8 {
1059 :     my $node = shift;
1060 :    
1061 :     my $x = 0; # Distance to node
1062 :     my $dl = newick_desc_ref( $node ); # Descendent list of tip node;
1063 :     while ( $dl && ( @$dl == 1 ) ) # Traverse unbranched nodes
1064 :     {
1065 :     $node = $dl->[0];
1066 : golsen 1.9 $x += newick_x( $node );
1067 : golsen 1.8 $dl = newick_desc_ref( $node );
1068 :     }
1069 :     $x += newick_x( $node );
1070 :    
1071 :     # Is it a tip? Return list of one tip;
1072 :    
1073 :     if ( ( ! $dl ) || @$dl == 0 )
1074 :     {
1075 :     return ( [ newick_lbl( $node ) ], $x );
1076 :     }
1077 :    
1078 :     # Get tips of each descendent, keeping lowest sorting from each.
1079 :     # Return the two lowest of those (the third will come from the
1080 :     # other side of the original node).
1081 :    
1082 : golsen 1.9 my @rep_tips = sort { lc $a cmp lc $b }
1083 :     map { ( sort { lc $a cmp lc $b } newick_tip_list( $_ ) )[0] }
1084 :     @$dl;
1085 :     return ( [ @rep_tips[0,1] ], $x );
1086 : golsen 1.8 }
1087 :    
1088 :    
1089 :     #-------------------------------------------------------------------------------
1090 : golsen 1.1 # Standard node name:
1091 :     # Tip label if at a tip
1092 :     # Three sorted tip labels intersecting at node, each being smallest
1093 :     # of all the tips of their subtrees
1094 :     #
1095 : golsen 1.9 # @TipOrTips = std_node_name( $tree, $node )
1096 : golsen 1.1 #-------------------------------------------------------------------------------
1097 :     sub std_node_name {
1098 :     my $tree = $_[0];
1099 :    
1100 :     # Node reference is last element of path to node
1101 :    
1102 :     my $noderef = ( path_to_node( @_ ) )[-1];
1103 :     defined( $noderef ) || return ();
1104 :    
1105 :     if ( node_is_tip( $noderef ) || $noderef eq $tree ) { # Is it a tip?
1106 :     return newick_lbl( $noderef );
1107 :     }
1108 :    
1109 :     # Work through lists of tips in descendant subtrees, removing them from
1110 :     # @rest, and keeping the best tip for each subtree.
1111 :    
1112 : golsen 1.8 my @rest = newick_tip_list( $tree );
1113 : golsen 1.9 my @best = map
1114 :     {
1115 : golsen 1.8 my @tips = sort { lc $a cmp lc $b } newick_tip_list( $_ );
1116 : golsen 1.9 @rest = &set_difference( \@rest, \@tips );
1117 : golsen 1.1 $tips[0];
1118 : golsen 1.9 } newick_desc_list( $noderef );
1119 : golsen 1.1
1120 :     # Best of the rest of the tree
1121 :     push @best, ( sort { lc $a cmp lc $b } @rest )[0];
1122 :    
1123 :     # Take the top 3, in order:
1124 :    
1125 :     ( @best >= 3 ) ? ( sort { lc $a cmp lc $b } @best )[0 .. 2] : ();
1126 :     }
1127 :    
1128 :    
1129 :     #===============================================================================
1130 :     # Functions to find paths in trees.
1131 :     #
1132 :     # Path descriptions are of form:
1133 :     # ( $node0, $i0, $node1, $i1, $node2, $i2, ..., $nodeN ) # Always odd
1134 :     # () is returned upon failure
1135 :     #
1136 :     # Numbering of descendants is 1-based.
1137 :     #===============================================================================
1138 :     # Path to tip:
1139 :     #
1140 :     # @path = path_to_tip( $treenode, $tipname )
1141 :     #-------------------------------------------------------------------------------
1142 :     sub path_to_tip {
1143 :     my ( $node, $tip, @path0 ) = @_;
1144 :    
1145 :     push @path0, $node;
1146 :     my $imax = newick_n_desc( $node );
1147 :     if ( $imax < 1 ) { return ( newick_lbl( $node ) eq $tip ) ? @path0 : () }
1148 :    
1149 :     # Special case for tree rooted on tip
1150 :    
1151 :     if ( ( $imax == 1 ) # One descendant
1152 :     && ( @path0 == 1 ) # First step in path
1153 :     && ( newick_lbl( $node ) eq $tip ) # Label matches
1154 :     ) { return @path0 }
1155 :    
1156 :     my @path;
1157 :     for (my $i = 1; $i <= $imax; $i++ ) {
1158 :     @path = path_to_tip( newick_desc_i( $node, $i ), $tip, ( @path0, $i ) );
1159 :     if ( @path ) { return @path }
1160 :     }
1161 :    
1162 :     (); # Not found
1163 :     }
1164 :    
1165 :    
1166 :     #-------------------------------------------------------------------------------
1167 :     # Path to named node.
1168 :     # Like path to tip, but will find named internal nodes as well.
1169 :     #
1170 :     # @path = path_to_named_node( $treenode, $name )
1171 :     #-------------------------------------------------------------------------------
1172 :     sub path_to_named_node {
1173 :     my ( $node, $name, @path0 ) = @_;
1174 :    
1175 :     push @path0, $node;
1176 :     if ( newick_lbl( $node ) eq $name ) { return @path0 }
1177 :    
1178 :     my @path;
1179 :     my $imax = newick_n_desc( $node );
1180 :     for ( my $i = 1; $i <= $imax; $i++ ) {
1181 :     @path = path_to_named_node( newick_desc_i( $node, $i ), $name, ( @path0, $i ) );
1182 :     if ( @path ) { return @path }
1183 :     }
1184 :    
1185 :     (); # Not found
1186 :     }
1187 :    
1188 :    
1189 :     #-------------------------------------------------------------------------------
1190 :     # Path to node reference.
1191 :     #
1192 :     # @path = path_to_node_ref( $treenode, $noderef )
1193 :     #-------------------------------------------------------------------------------
1194 :     sub path_to_node_ref {
1195 :     my ( $node, $noderef, @path0 ) = @_;
1196 :    
1197 :     push @path0, $node;
1198 :     if ( $node eq $noderef ) { return @path0 }
1199 :    
1200 :     my @path;
1201 :     my $imax = newick_n_desc( $node );
1202 :     for ( my $i = 1; $i <= $imax; $i++ ) {
1203 : golsen 1.12 @path = path_to_node_ref( newick_desc_i( $node, $i ), $noderef, ( @path0, $i ) );
1204 :     return @path if @path;
1205 : golsen 1.1 }
1206 :    
1207 :     (); # Not found
1208 :     }
1209 :    
1210 :    
1211 :     #-------------------------------------------------------------------------------
1212 :     # Path to node, as defined by 1 or 3 tips.
1213 :     #
1214 :     # @path = path_to_node( $node, $tip1, $tip2, $tip3 ) # 3 tip names
1215 :     # @path = path_to_node( $node, [ $tip1, $tip2, $tip3 ] ) # Allow array ref
1216 :     # @path = path_to_node( $node, $tip1 ) # Use path_to_tip
1217 :     # @path = path_to_node( $node, [ $tip1 ] ) # Use path_to_tip
1218 :     #-------------------------------------------------------------------------------
1219 :     sub path_to_node {
1220 :     my ( $node, $tip1, $tip2, $tip3 ) = @_;
1221 :     array_ref( $node ) && defined( $tip1 ) || return ();
1222 :    
1223 :     # Allow arg 2 to be an array reference
1224 :     if ( array_ref( $tip1 ) ) { ( $tip1, $tip2, $tip3 ) = @$tip1 }
1225 :    
1226 :     my @p1 = path_to_tip( $node, $tip1 ); # Path to first tip
1227 :     @p1 || return (); # Was the tip found?
1228 :     defined( $tip2 ) && defined( $tip3 ) || return @p1; # Two more defined?
1229 :    
1230 :     my @p2 = path_to_tip( $node, $tip2 );
1231 :     my @p3 = path_to_tip( $node, $tip3 );
1232 :     @p2 && @p3 || return (); # Were they found?
1233 :    
1234 :     # Find the common prefix for each pair of paths
1235 : golsen 1.9 my @p12 = &common_prefix( \@p1, \@p2 );
1236 :     my @p13 = &common_prefix( \@p1, \@p3 );
1237 :     my @p23 = &common_prefix( \@p2, \@p3 );
1238 : golsen 1.1
1239 :     # Return the longest common prefix of any two paths
1240 :     ( @p12 >= @p13 && @p12 >= @p23 ) ? @p12 :
1241 :     ( @p13 >= @p23 ) ? @p13 :
1242 :     @p23 ;
1243 :     }
1244 :    
1245 :    
1246 :     #-------------------------------------------------------------------------------
1247 :     # Distance along path.
1248 :     #
1249 :     # $distance = newick_path_length( @path )
1250 :     #-------------------------------------------------------------------------------
1251 :     sub newick_path_length {
1252 :     my $node = shift; # Discard the first node
1253 :     array_ref( $node ) || return undef;
1254 :     @_ ? distance_along_path_2( @_ ) : 0;
1255 :     }
1256 :    
1257 :    
1258 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1259 :     # This expects to get path minus root node:
1260 :     #
1261 :     # $distance = distance_along_path_2( @path )
1262 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1263 :     sub distance_along_path_2 {
1264 :     shift; # Discard descendant number
1265 :     my $node = shift;
1266 :     array_ref( $node ) || return undef;
1267 :     my $d1 = newick_x( $node );
1268 : golsen 1.8 my $d2 = @_ ? distance_along_path_2( @_ ) : 0;
1269 : golsen 1.1 defined( $d1 ) && defined( $d2 ) ? $d1 + $d2 : undef;
1270 :     }
1271 :    
1272 :    
1273 :     #-------------------------------------------------------------------------------
1274 :     # Tip-to-tip distance.
1275 :     #
1276 :     # $distance = tip_to_tip_distance( $tree, $tip1, $tip2 )
1277 :     #-------------------------------------------------------------------------------
1278 :     sub tip_to_tip_distance {
1279 :     my ( $node, $tip1, $tip2 ) = @_;
1280 :    
1281 :     array_ref( $node ) && defined( $tip1 )
1282 :     && defined( $tip2 ) || return undef;
1283 :     my @p1 = path_to_tip( $node, $tip1 );
1284 :     my @p2 = path_to_tip( $node, $tip2 );
1285 :     @p1 && @p2 || return undef; # Were they found?
1286 :    
1287 :     # Find the unique suffixes of the two paths
1288 : golsen 1.9 my ( $suf1, $suf2 ) = &unique_suffixes( \@p1, \@p2 ); # Common node is lost
1289 : golsen 1.1 my $d1 = @$suf1 ? distance_along_path_2( @$suf1 ) : 0;
1290 :     my $d2 = @$suf2 ? distance_along_path_2( @$suf2 ) : 0;
1291 :    
1292 :     defined( $d1 ) && defined( $d2 ) ? $d1 + $d2 : undef;
1293 :     }
1294 :    
1295 :    
1296 :     #-------------------------------------------------------------------------------
1297 :     # Node-to-node distance.
1298 :     # Nodes can be: $tipname
1299 :     # [ $tipname ]
1300 :     # [ $tipname1, $tipname2, $tipname3 ]
1301 :     #
1302 :     # $distance = node_to_node_distance( $tree, $node1, $node2 )
1303 :     #-------------------------------------------------------------------------------
1304 :     sub node_to_node_distance {
1305 :     my ( $node, $node1, $node2 ) = @_;
1306 :    
1307 :     array_ref( $node ) && defined( $node1 )
1308 :     && defined( $node2 ) || return undef;
1309 : golsen 1.8 my @p1 = path_to_node( $node, $node1 ) or return undef;
1310 :     my @p2 = path_to_node( $node, $node2 ) or return undef;
1311 : golsen 1.1
1312 :     # Find the unique suffixes of the two paths
1313 : golsen 1.9 my ( $suf1, $suf2 ) = &unique_suffixes( \@p1, \@p2 ); # Common node is lost
1314 : golsen 1.1 my $d1 = @$suf1 ? distance_along_path_2( @$suf1 ) : 0;
1315 :     my $d2 = @$suf2 ? distance_along_path_2( @$suf2 ) : 0;
1316 :    
1317 :     defined( $d1 ) && defined( $d2 ) ? $d1 + $d2 : undef;
1318 :     }
1319 :    
1320 :    
1321 :     #===============================================================================
1322 :     # Tree manipulations:
1323 :     #===============================================================================
1324 :     # Copy tree.
1325 :     # Lists are copied, except that references to empty lists go to undef.
1326 :     # Only defined fields are added, so tree list may be shorter than 8 fields.
1327 :     #
1328 : overbeek 1.4 # $treecopy = copy_newick_tree( $tree )
1329 : golsen 1.1 #-------------------------------------------------------------------------------
1330 : overbeek 1.4 sub copy_newick_tree {
1331 : golsen 1.1 my ( $node ) = @_;
1332 :     array_ref( $node ) || return undef;
1333 :    
1334 :     my $nn = []; # Reference to a new node structure
1335 :     # Build a new descendant list, if not empty
1336 :     my @dl = newick_desc_list( $node );
1337 : overbeek 1.4 set_newick_desc_ref( $nn, @dl ? [ map { copy_newick_tree( $_ ) } @dl ]
1338 : golsen 1.1 : undef
1339 :     );
1340 :    
1341 :     # Copy label and x, if defined
1342 :     my ( $l, $x );
1343 :     if ( defined( $l = newick_lbl( $node ) ) ) { set_newick_lbl( $nn, $l ) }
1344 :     if ( defined( $x = newick_x( $node ) ) ) { set_newick_x( $nn, $x ) }
1345 :    
1346 :     # Build new comment lists, when not empty ( does not extend array unless
1347 :     # necessary)
1348 :     my $c;
1349 :     if ( $c = newick_c1( $node ) and @$c ) { set_newick_c1( $nn, [ @$c ] ) }
1350 :     if ( $c = newick_c2( $node ) and @$c ) { set_newick_c2( $nn, [ @$c ] ) }
1351 :     if ( $c = newick_c3( $node ) and @$c ) { set_newick_c3( $nn, [ @$c ] ) }
1352 :     if ( $c = newick_c4( $node ) and @$c ) { set_newick_c4( $nn, [ @$c ] ) }
1353 :     if ( $c = newick_c5( $node ) and @$c ) { set_newick_c5( $nn, [ @$c ] ) }
1354 :    
1355 :     $nn;
1356 :     }
1357 :    
1358 :    
1359 :     #-------------------------------------------------------------------------------
1360 :     # Use a hash to relabel the nodes in a newick tree.
1361 :     #
1362 :     # $newtree = newick_relabel_nodes( $node, \%new_name )
1363 :     #-------------------------------------------------------------------------------
1364 :     sub newick_relabel_nodes {
1365 :     my ( $node, $new_name ) = @_;
1366 :    
1367 :     my ( $lbl, $new );
1368 :     if ( defined( $lbl = newick_lbl( $node ) )
1369 :     && ( $lbl ne "" )
1370 :     && defined( $new = $new_name->{ $lbl } )
1371 :     ) {
1372 :     set_newick_lbl( $node, $new );
1373 :     }
1374 :    
1375 :     foreach ( newick_desc_list( $node ) ) {
1376 :     newick_relabel_nodes( $_, $new_name );
1377 :     }
1378 :    
1379 :     $node;
1380 :     }
1381 :    
1382 :    
1383 :     #-------------------------------------------------------------------------------
1384 :     # Use a hash to relabel the nodes in a newick tree (case insensitive).
1385 :     #
1386 :     # $newtree = newick_relabel_nodes_i( $node, \%new_name )
1387 :     #-------------------------------------------------------------------------------
1388 :     sub newick_relabel_nodes_i {
1389 :     my ( $node, $new_name ) = @_;
1390 :    
1391 :     # Add any necessary lowercase keys to the hash:
1392 :    
1393 :     my $lc_lbl;
1394 :     foreach ( keys %$new_name ) {
1395 :     $lc_lbl = lc $_;
1396 :     ( $lc_lbl eq $_ ) or ( $new_name->{ $lc_lbl } = $new_name->{ $_ } );
1397 :     }
1398 :    
1399 :     newick_relabel_nodes_i2( $node, $new_name );
1400 :     }
1401 :    
1402 :    
1403 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1404 :     # Do the actual relabeling
1405 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1406 :     sub newick_relabel_nodes_i2 {
1407 :     my ( $node, $new_name ) = @_;
1408 :    
1409 :     my ( $lbl, $new );
1410 :     if ( defined( $lbl = newick_lbl( $node ) )
1411 :     && ( $lbl ne "" )
1412 :     && defined( $new = $new_name->{ lc $lbl } )
1413 :     ) {
1414 :     set_newick_lbl( $node, $new );
1415 :     }
1416 :    
1417 :     foreach ( newick_desc_list( $node ) ) {
1418 :     newick_relabel_nodes_i2( $_, $new_name );
1419 :     }
1420 :    
1421 :     $node;
1422 :     }
1423 :    
1424 :    
1425 :     #-------------------------------------------------------------------------------
1426 :     # Use a hash to relabel the tips in a newick tree.
1427 :     #
1428 :     # $newtree = newick_relabel_tips( $node, \%new_name )
1429 :     #-------------------------------------------------------------------------------
1430 :     sub newick_relabel_tips {
1431 :     my ( $node, $new_name ) = @_;
1432 :    
1433 :     my @desc = newick_desc_list( $node );
1434 :    
1435 :     if ( @desc ) {
1436 :     foreach ( @desc ) { newick_relabel_tips( $_, $new_name ) }
1437 :     }
1438 :     else {
1439 :     my ( $lbl, $new );
1440 :     if ( defined( $lbl = newick_lbl( $node ) )
1441 :     && ( $lbl ne "" )
1442 :     && defined( $new = $new_name->{ $lbl } )
1443 :     ) {
1444 :     set_newick_lbl( $node, $new );
1445 :     }
1446 :     }
1447 :    
1448 :     $node;
1449 :     }
1450 :    
1451 :    
1452 :     #-------------------------------------------------------------------------------
1453 :     # Use a hash to relabel the tips in a newick tree (case insensitive).
1454 :     #
1455 :     # $newtree = newick_relabel_tips_i( $node, \%new_name )
1456 :     #-------------------------------------------------------------------------------
1457 :     sub newick_relabel_tips_i {
1458 :     my ( $node, $new_name ) = @_;
1459 :    
1460 :     # Add any necessary lowercase keys to the hash:
1461 :    
1462 :     my $lc_lbl;
1463 :     foreach ( keys %$new_name ) {
1464 :     $lc_lbl = lc $_;
1465 :     ( $lc_lbl eq $_ ) or ( $new_name->{ $lc_lbl } = $new_name->{ $_ } );
1466 :     }
1467 :    
1468 :     newick_relabel_tips_i2( $node, $new_name );
1469 :     }
1470 :    
1471 :    
1472 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1473 :     # Do the actual relabeling
1474 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1475 :     sub newick_relabel_tips_i2 {
1476 :     my ( $node, $new_name ) = @_;
1477 :    
1478 :     my @desc = newick_desc_list( $node );
1479 :    
1480 :     if ( @desc ) {
1481 :     foreach ( @desc ) { newick_relabel_tips_i2( $_, $new_name ) }
1482 :     }
1483 :     else {
1484 :     my ( $lbl, $new );
1485 :     if ( defined( $lbl = newick_lbl( $node ) )
1486 :     && ( $lbl ne "" )
1487 :     && defined( $new = $new_name->{ lc $lbl } )
1488 :     ) {
1489 :     set_newick_lbl( $node, $new );
1490 :     }
1491 :     }
1492 :    
1493 :     $node;
1494 :     }
1495 :    
1496 :    
1497 :     #-------------------------------------------------------------------------------
1498 :     # Set undefined branch lenghts (except root) to length x.
1499 :     #
1500 :     # $n_changed = newick_set_undefined_branches( $node, $x )
1501 :     #-------------------------------------------------------------------------------
1502 :     sub newick_set_undefined_branches {
1503 :     my ( $node, $x, $not_root ) = @_;
1504 :    
1505 :     my $n = 0;
1506 :     if ( $not_root && ! defined( newick_x( $node ) ) ) {
1507 :     set_newick_x( $node, $x );
1508 :     $n++;
1509 :     }
1510 :    
1511 :     foreach ( newick_desc_list( $node ) ) {
1512 :     $n += newick_set_undefined_branches( $_, $x, 1 );
1513 :     }
1514 :    
1515 :     $n;
1516 :     }
1517 :    
1518 :    
1519 :     #-------------------------------------------------------------------------------
1520 :     # Set all branch lenghts (except root) to length x.
1521 :     #
1522 :     # $n_changed = newick_set_all_branches( $node, $x )
1523 :     #-------------------------------------------------------------------------------
1524 :     sub newick_set_all_branches {
1525 :     my ( $node, $x, $not_root ) = @_;
1526 :    
1527 :     my $n = 0;
1528 : overbeek 1.7 if ( $not_root )
1529 :     {
1530 : golsen 1.1 set_newick_x( $node, $x );
1531 :     $n++;
1532 :     }
1533 :    
1534 : overbeek 1.7 foreach ( newick_desc_list( $node ) )
1535 :     {
1536 : golsen 1.1 $n += newick_set_all_branches( $_, $x, 1 );
1537 :     }
1538 :    
1539 :     $n;
1540 :     }
1541 :    
1542 :    
1543 :     #-------------------------------------------------------------------------------
1544 : overbeek 1.7 # Rescale all branch lenghts by factor.
1545 :     #
1546 :     # $node = newick_rescale_branches( $node, $factor )
1547 :     #-------------------------------------------------------------------------------
1548 :     sub newick_rescale_branches {
1549 :     my ( $node, $factor ) = @_;
1550 :    
1551 :     my $x = newick_x( $node );
1552 :     set_newick_x( $node, $factor * $x ) if $x;
1553 :    
1554 :     foreach ( newick_desc_list( $node ) )
1555 :     {
1556 :     newick_rescale_branches( $_, $factor );
1557 :     }
1558 :    
1559 :     $node;
1560 :     }
1561 :    
1562 :    
1563 :     #-------------------------------------------------------------------------------
1564 : golsen 1.15 # Modify all branch lengths by a function.
1565 :     #
1566 :     # $node = newick_modify_branches( $node, \&function )
1567 :     # $node = newick_modify_branches( $node, \&function, \@func_parms )
1568 :     #
1569 :     # Function must have form
1570 :     #
1571 :     # $x2 = &$function( $x1 )
1572 :     # $x2 = &$function( $x1, @$func_parms )
1573 :     #
1574 :     #-------------------------------------------------------------------------------
1575 :     sub newick_modify_branches {
1576 :     my ( $node, $func, $parm ) = @_;
1577 :    
1578 :     set_newick_x( $node, &$func( newick_x( $node ), ( $parm ? @$parm : () ) ) );
1579 :     foreach ( newick_desc_list( $node ) )
1580 :     {
1581 :     newick_modify_branches( $_, $func, $parm )
1582 :     }
1583 :    
1584 :     $node;
1585 :     }
1586 :    
1587 :    
1588 :     #-------------------------------------------------------------------------------
1589 : golsen 1.5 # Set negative branches to zero. The original tree is modfied.
1590 :     #
1591 :     # $n_changed = newick_fix_negative_branches( $tree )
1592 :     #-------------------------------------------------------------------------------
1593 :     sub newick_fix_negative_branches {
1594 :     my ( $tree ) = @_;
1595 :     array_ref( $tree ) or return undef;
1596 :     my $n_changed = 0;
1597 :     my $x = newick_x( $tree );
1598 :     if ( defined( $x ) and $x < 0 )
1599 :     {
1600 :     set_newick_x( $tree, 0 );
1601 :     $n_changed++;
1602 :     }
1603 :    
1604 :     foreach ( newick_desc_list( $tree ) )
1605 :     {
1606 :     $n_changed += newick_fix_negative_branches( $_ );
1607 :     }
1608 :    
1609 :     $n_changed;
1610 :     }
1611 :    
1612 :    
1613 :     #-------------------------------------------------------------------------------
1614 : golsen 1.8 # Remove comments from a newick tree (e.g., before writing for phylip).
1615 :     #
1616 :     # $node = newick_strip_comments( $node )
1617 :     #-------------------------------------------------------------------------------
1618 :     sub newick_strip_comments {
1619 :     my ( $node ) = @_;
1620 :    
1621 :     @$node = @$node[ 0 .. 2 ];
1622 :     foreach ( newick_desc_list( $node ) ) { newick_strip_comments( $_ ) }
1623 :     $node;
1624 :     }
1625 :    
1626 :    
1627 :     #-------------------------------------------------------------------------------
1628 : golsen 1.1 # Normalize tree order (in place).
1629 :     #
1630 :     # ( $tree, $label1 ) = normalize_newick_tree( $tree )
1631 :     #-------------------------------------------------------------------------------
1632 :     sub normalize_newick_tree {
1633 :     my ( $node ) = @_;
1634 :    
1635 :     my @descends = newick_desc_list( $node );
1636 :     if ( @descends == 0 ) { return ( $node, lc newick_lbl( $node ) ) }
1637 :    
1638 : golsen 1.9 my %hash = map { ( normalize_newick_tree($_) )[1] => $_ } @descends;
1639 : golsen 1.1 my @keylist = sort { $a cmp $b } keys %hash;
1640 :     set_newick_desc_list( $node, map { $hash{$_} } @keylist );
1641 :    
1642 :     ( $node, $keylist[0] );
1643 :     }
1644 :    
1645 :    
1646 :     #-------------------------------------------------------------------------------
1647 :     # Reverse tree order (in place).
1648 :     #
1649 :     # $tree = reverse_newick_tree( $tree )
1650 :     #-------------------------------------------------------------------------------
1651 :     sub reverse_newick_tree {
1652 :     my ( $node ) = @_;
1653 :    
1654 :     my @descends = newick_desc_list( $node );
1655 :     if ( @descends ) {
1656 :     set_newick_desc_list( $node, reverse @descends );
1657 :     foreach ( @descends ) { reverse_newick_tree( $_ ) }
1658 :     }
1659 :     $node;
1660 :     }
1661 :    
1662 :    
1663 :     #-------------------------------------------------------------------------------
1664 :     # Standard unrooted tree (in place).
1665 :     #
1666 :     # $stdtree = std_unrooted_newick( $tree )
1667 :     #-------------------------------------------------------------------------------
1668 :     sub std_unrooted_newick {
1669 :     my ( $tree ) = @_;
1670 :    
1671 :     my ( $mintip ) = sort { lc $a cmp lc $b } newick_tip_list( $tree );
1672 :     ( normalize_newick_tree( reroot_newick_next_to_tip( $tree, $mintip ) ) )[0];
1673 :     }
1674 :    
1675 :    
1676 :     #-------------------------------------------------------------------------------
1677 : golsen 1.20 # Standard name for a Newick tree topology
1678 :     #
1679 :     # $stdname = std_tree_name( $tree )
1680 :     #
1681 :     #-------------------------------------------------------------------------------
1682 :     sub std_tree_name
1683 :     {
1684 :     my ( $tree ) = @_;
1685 :     my ( $mintip ) = sort { lc $a cmp lc $b } newick_tip_list( $tree );
1686 :     ( std_tree_name_2( reroot_newick_next_to_tip( copy_newick_tree( $tree ), $mintip ) ) )[0];
1687 :     }
1688 :    
1689 :    
1690 :     #
1691 :     # ( $name, $mintip ) = std_tree_name_2( $node )
1692 :     #
1693 :     sub std_tree_name_2
1694 :     {
1695 :     my ( $node ) = @_;
1696 :    
1697 :     my @descends = newick_desc_list( $node );
1698 :     if ( @descends == 0 )
1699 :     {
1700 :     my $lbl = newick_lbl( $node );
1701 :     return ( $lbl, $lbl );
1702 :     }
1703 :    
1704 :     my @list = sort { lc $a->[1] cmp lc $b->[1] || $a->[1] cmp $b->[1] }
1705 :     map { [ std_tree_name_2( $_ ) ] }
1706 :     @descends;
1707 :     my $mintip = $list[0]->[1];
1708 :     my $name = '(' . join( "\t", map { $_->[0] } @list ) . ')';
1709 :    
1710 :     return ( $name, $mintip );
1711 :     }
1712 :    
1713 :    
1714 :     #-------------------------------------------------------------------------------
1715 : golsen 1.1 # Move largest groups to periphery of tree (in place).
1716 :     #
1717 :     # dir <= -2 for up-sweeping tree (big groups always first),
1718 :     # = -1 for big group first, balanced tree,
1719 :     # = 0 for balanced tree,
1720 :     # = 1 for small group first, balanced tree, and
1721 :     # >= 2 for down-sweeping tree (small groups always top)
1722 :     #
1723 :     # $tree = aesthetic_newick_tree( $treeref, $dir )
1724 :     #-------------------------------------------------------------------------------
1725 :     sub aesthetic_newick_tree {
1726 :     my ( $tree, $dir ) = @_;
1727 :     my %cnt;
1728 :    
1729 :     $dir = ! $dir ? 0 : # Undefined or zero
1730 :     $dir <= -2 ? -1000000 :
1731 :     $dir < 0 ? -1 :
1732 :     $dir >= 2 ? 1000000 :
1733 :     1 ;
1734 :     build_tip_count_hash( $tree, \%cnt );
1735 :     reorder_by_tip_count( $tree, \%cnt, $dir );
1736 :     }
1737 :    
1738 :    
1739 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1740 :     # Build a hash to look up the number of descendants of each node.
1741 :     # Access count with $cntref->{$noderef}
1742 :     #
1743 :     # $count = build_tip_count_hash( $node, $cnt_hash_ref )
1744 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1745 :     sub build_tip_count_hash {
1746 :     my ( $node, $cntref ) = @_;
1747 :     my ( $i, $cnt );
1748 :    
1749 :     if ( newick_n_desc( $node ) < 1 ) { $cnt = 1 }
1750 :     else {
1751 :     $cnt = 0;
1752 :     foreach ( newick_desc_list( $node ) ) {
1753 :     $cnt += build_tip_count_hash( $_, $cntref );
1754 :     }
1755 :     }
1756 :    
1757 :     $cntref->{$node} = $cnt;
1758 :     $cnt;
1759 :     }
1760 :    
1761 :    
1762 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1763 :     # $node = reorder_by_tip_count( $node, $cntref, $dir )
1764 :     # dir < 0 for upward branch (big group first),
1765 :     # = 0 for no change, and
1766 :     # > 0 for downward branch (small group first).
1767 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1768 :     sub reorder_by_tip_count {
1769 :     my ( $node, $cntref, $dir ) = @_;
1770 :    
1771 :     my $nd = newick_n_desc( $node );
1772 :     if ( $nd < 1 ) { return $node } # Do nothing to a tip
1773 :    
1774 : golsen 1.16 my $dl_ref = newick_desc_ref( $node );
1775 :    
1776 :     # Reorder this subtree (biggest subtrees to outside)
1777 :    
1778 :     if ( $dir )
1779 :     {
1780 :     # Big group first
1781 :     my @dl = sort { $cntref->{$b} <=> $cntref->{$a} } @$dl_ref;
1782 :    
1783 :     my ( @dl1, @dl2 );
1784 :     for ( my $i = 0; $i < $nd; $i++ ) {
1785 :     if ( $i & 1 ) { push @dl2, $dl[$i] } else { push @dl1, $dl[$i] }
1786 :     }
1787 : golsen 1.1
1788 : golsen 1.16 @$dl_ref = ( $dir < 0 ) ? ( @dl1, reverse @dl2 )
1789 :     : ( @dl2, reverse @dl1 );
1790 : golsen 1.1 }
1791 :    
1792 :     # Reorder within descendant subtrees:
1793 :    
1794 :     my $step = 0;
1795 :     if ( abs( $dir ) < 1e5 ) {
1796 :     $dir = 1 - $nd; # Midgroup => as is
1797 :     # $dir = 1 - $nd + ( $dir < 0 ? -0.5 : 0.5 ); # Midgroup => outward
1798 :     $step = 2;
1799 :     }
1800 :    
1801 :     for ( my $i = 0; $i < $nd; $i++ ) {
1802 :     reorder_by_tip_count( $dl_ref->[$i], $cntref, $dir );
1803 :     $dir += $step;
1804 :     }
1805 :    
1806 :     $node;
1807 :     }
1808 :    
1809 :    
1810 :     #-------------------------------------------------------------------------------
1811 :     # Move smallest groups to periphery of tree (in place).
1812 :     #
1813 :     # dir <= -2 for up-sweeping tree (big groups always first),
1814 :     # = -1 for big group first, balanced tree,
1815 :     # = 0 for balanced tree,
1816 :     # = 1 for small group first, balanced tree, and
1817 :     # >= 2 for down-sweeping tree (small groups always top)
1818 :     #
1819 :     # $tree = unaesthetic_newick_tree( $treeref, $dir )
1820 :     #-------------------------------------------------------------------------------
1821 : golsen 1.9 sub unaesthetic_newick_tree
1822 :     {
1823 : golsen 1.1 my ( $tree, $dir ) = @_;
1824 :     my %cnt;
1825 :    
1826 :     $dir = ! $dir ? 0 : # Undefined or zero
1827 :     $dir <= -2 ? -1000000 :
1828 :     $dir < 0 ? -1 :
1829 :     $dir >= 2 ? 1000000 :
1830 :     1 ;
1831 :     build_tip_count_hash( $tree, \%cnt );
1832 :     reorder_against_tip_count( $tree, \%cnt, $dir );
1833 :     }
1834 :    
1835 :    
1836 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1837 :     # $node = reorder_by_tip_count( $node, $cntref, $dir )
1838 :     # dir < 0 for upward branch (big group first),
1839 :     # = 0 for no change, and
1840 :     # > 0 for downward branch (small group first).
1841 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1842 : golsen 1.9 sub reorder_against_tip_count
1843 :     {
1844 : golsen 1.1 my ( $node, $cntref, $dir ) = @_;
1845 :    
1846 :     my $nd = newick_n_desc( $node );
1847 :     if ( $nd < 1 ) { return $node } # Do nothing to a tip
1848 :    
1849 :     # Reorder this subtree:
1850 :    
1851 :     my $dl_ref = newick_desc_ref( $node );
1852 :     if ( $dir > 0 ) { # Big group first
1853 :     @$dl_ref = sort { $cntref->{$b} <=> $cntref->{$a} } @$dl_ref;
1854 :     }
1855 :     elsif ( $dir < 0 ) { # Small group first
1856 :     @$dl_ref = sort { $cntref->{$a} <=> $cntref->{$b} } @$dl_ref;
1857 :     }
1858 :    
1859 :     # Reorder within descendant subtrees:
1860 :    
1861 :     my $step = 0;
1862 :     if (abs( $dir ) < 1e5) {
1863 :     $dir = 1 - $nd; # Midgroup => as is
1864 :     # $dir = 1 - $nd + ( $dir < 0 ? -0.5 : 0.5 ); # Midgroup => outward
1865 :     $step = 2;
1866 :     }
1867 :    
1868 :     for ( my $i = 0; $i < $nd; $i++ ) {
1869 :     reorder_by_tip_count( $dl_ref->[$i], $cntref, $dir );
1870 :     $dir += $step;
1871 :     }
1872 :    
1873 :     $node;
1874 :     }
1875 :    
1876 :    
1877 :     #-------------------------------------------------------------------------------
1878 :     # Randomize descendant order at each node (in place).
1879 :     #
1880 :     # $tree = random_order_newick_tree( $tree )
1881 :     #-------------------------------------------------------------------------------
1882 : golsen 1.9 sub random_order_newick_tree
1883 :     {
1884 : golsen 1.1 my ( $node ) = @_;
1885 :    
1886 :     my $nd = newick_n_desc( $node );
1887 :     if ( $nd < 1 ) { return $node } # Do nothing to a tip
1888 :    
1889 :     # Reorder this subtree:
1890 :    
1891 :     my $dl_ref = newick_desc_ref( $node );
1892 : golsen 1.9 @$dl_ref = &random_order( @$dl_ref );
1893 : golsen 1.1
1894 :     # Reorder descendants:
1895 :    
1896 :     foreach ( @$dl_ref ) { random_order_newick_tree( $_ ) }
1897 :    
1898 :     $node;
1899 :     }
1900 :    
1901 :    
1902 :     #-------------------------------------------------------------------------------
1903 :     # Reroot a tree to the node that lies at the end of a path.
1904 :     #
1905 :     # $newtree = reroot_newick_by_path( @path )
1906 :     #-------------------------------------------------------------------------------
1907 : golsen 1.9 sub reroot_newick_by_path
1908 :     {
1909 : golsen 1.1 my ( $node1, $path1, @rest ) = @_;
1910 :     array_ref( $node1 ) || return undef; # Always expect a node
1911 :    
1912 :     defined( $path1 ) && @rest || return $node1; # If no path, we're done
1913 :    
1914 :     my $node2 = $rest[0]; # Next element in path is node 2
1915 :     newick_desc_i( $node1, $path1 ) eq $node2 || return undef; # Check link
1916 :    
1917 :     # Remove node 2 from node 1 descendant list. Could use a simple splice:
1918 :     #
1919 :     # splice( @$dl1, $path1-1, 1 );
1920 :     #
1921 : golsen 1.8 # But the following maintains the cyclic order of the nodes:
1922 : golsen 1.1
1923 :     my $dl1 = newick_desc_ref( $node1 );
1924 :     my $nd1 = @$dl1;
1925 :     if ( $path1 == 1 ) { shift @$dl1 }
1926 :     elsif ( $path1 == $nd1 ) { pop @$dl1 }
1927 :     else { @$dl1 = ( @$dl1[ $path1 .. $nd1-1 ]
1928 :     , @$dl1[ 0 .. $path1-2 ]
1929 :     )
1930 : golsen 1.8 }
1931 : golsen 1.1
1932 :     # Append node 1 to node 2 descendant list (does not alter numbering):
1933 :    
1934 :     my $dl2 = newick_desc_ref( $node2 );
1935 :     if ( array_ref( $dl2 ) ) { push @$dl2, $node1 }
1936 : golsen 1.8 else { set_newick_desc_list( $node2, $node1 ) }
1937 : golsen 1.1
1938 :     # Move c1 comments from node 1 to node 2:
1939 :    
1940 :     my $C11 = newick_c1( $node1 );
1941 :     my $C12 = newick_c1( $node2 );
1942 :     ! defined( $C11 ) || set_newick_c1( $node1, undef ); # Remove them from node 1
1943 : golsen 1.8 if ( $C12 && @$C12 ) { # If node 2 comments and
1944 :     if ( $C11 && @$C11 ) { unshift @$C12, @$C11 } # Node 1, prefix 1 to 2
1945 : golsen 1.1 }
1946 : golsen 1.8 elsif ( $C11 && @$C11 ) { set_newick_c1( $node2, $C11 ) } # Otherwise move node 1 link
1947 : golsen 1.1
1948 :     # Swap branch lengths and comments for reversal of link direction:
1949 :    
1950 :     my $x1 = newick_x( $node1 );
1951 :     my $x2 = newick_x( $node2 );
1952 :     ! defined( $x1 ) && ! defined ( $x2 ) || set_newick_x( $node1, $x2 );
1953 :     ! defined( $x1 ) && ! defined ( $x2 ) || set_newick_x( $node2, $x1 );
1954 :    
1955 :     my $c41 = newick_c4( $node1 );
1956 :     my $c42 = newick_c4( $node2 );
1957 :     ! defined( $c42 ) || ! @$c42 || set_newick_c4( $node1, $c42 );
1958 :     ! defined( $c41 ) || ! @$c41 || set_newick_c4( $node2, $c41 );
1959 :    
1960 :     my $c51 = newick_c5( $node1 );
1961 :     my $c52 = newick_c5( $node2 );
1962 :     ! defined( $c52 ) || ! @$c52 || set_newick_c5( $node1, $c52 );
1963 :     ! defined( $c51 ) || ! @$c51 || set_newick_c5( $node2, $c51 );
1964 :    
1965 :     reroot_newick_by_path( @rest ); # Node 2 is first element of rest
1966 :     }
1967 :    
1968 :    
1969 :     #-------------------------------------------------------------------------------
1970 :     # Move root of tree to named tip.
1971 :     #
1972 :     # $newtree = reroot_newick_to_tip( $tree, $tip )
1973 :     #-------------------------------------------------------------------------------
1974 :     sub reroot_newick_to_tip {
1975 :     my ( $tree, $tipname ) = @_;
1976 :     reroot_newick_by_path( path_to_tip( $tree, $tipname ) );
1977 :     }
1978 :    
1979 :    
1980 :     #-------------------------------------------------------------------------------
1981 :     # Move root of tree to be node adjacent to a named tip.
1982 :     #
1983 :     # $newtree = reroot_newick_next_to_tip( $tree, $tip )
1984 :     #-------------------------------------------------------------------------------
1985 :     sub reroot_newick_next_to_tip {
1986 :     my ( $tree, $tipname ) = @_;
1987 :     my @path = path_to_tip( $tree, $tipname );
1988 :     @path || return undef;
1989 :     @path == 1 ? reroot_newick_by_path( $tree, 1, newick_desc_i( $tree, 1 ) )
1990 :     : reroot_newick_by_path( @path[0 .. @path-3] );
1991 :     }
1992 :    
1993 :    
1994 :     #-------------------------------------------------------------------------------
1995 :     # Move root of tree to a node, defined by 1 or 3 tip names.
1996 :     #
1997 :     # $newtree = reroot_newick_to_node( $tree, @node )
1998 :     #-------------------------------------------------------------------------------
1999 :     sub reroot_newick_to_node {
2000 :     reroot_newick_by_path( path_to_node( @_ ) );
2001 :     }
2002 :    
2003 :    
2004 :     #-------------------------------------------------------------------------------
2005 :     # Move root of tree to a node, defined by reference.
2006 :     #
2007 :     # $newtree = reroot_newick_to_node_ref( $tree, $noderef )
2008 :     #-------------------------------------------------------------------------------
2009 :     sub reroot_newick_to_node_ref {
2010 :     my ( $tree, $node ) = @_;
2011 :     reroot_newick_by_path( path_to_node_ref( $tree, $node ) );
2012 :     }
2013 :    
2014 :    
2015 :     #-------------------------------------------------------------------------------
2016 : golsen 1.9 # Reroot a newick tree along the path between 2 nodes:
2017 :     #
2018 :     # $tree = reroot_newick_between_nodes( $tree, $node1, $node2, $fraction )
2019 :     #-------------------------------------------------------------------------------
2020 :     sub reroot_newick_between_nodes
2021 :     {
2022 :     my ( $tree, $node1, $node2, $fraction ) = @_;
2023 :     array_ref( $tree ) or return undef;
2024 :     $fraction >= 0 && $fraction <= 1 or return undef;
2025 :    
2026 :     # Find the paths to the nodes:
2027 :    
2028 :     my @path1 = path_to_node( $tree, $node1 ) or return undef;
2029 :     my @path2 = path_to_node( $tree, $node2 ) or return undef;
2030 :    
2031 : golsen 1.12 reroot_newick_between_nodes_by_path( $tree, \@path1, \@path2, $fraction )
2032 :     }
2033 :    
2034 :    
2035 :     #-------------------------------------------------------------------------------
2036 :     # Reroot a newick tree along the path between 2 nodes:
2037 :     #
2038 :     # $tree = reroot_newick_between_node_refs( $tree, $node1, $node2, $fraction )
2039 :     #-------------------------------------------------------------------------------
2040 :     sub reroot_newick_between_node_refs
2041 :     {
2042 :     my ( $tree, $node1, $node2, $fraction ) = @_;
2043 :     array_ref( $tree ) or return undef;
2044 :    
2045 :     # Find the paths to the nodes:
2046 :    
2047 :     my @path1 = path_to_node_ref( $tree, $node1 ) or return undef;
2048 :     my @path2 = path_to_node_ref( $tree, $node2 ) or return undef;;
2049 :    
2050 :     reroot_newick_between_nodes_by_path( $tree, \@path1, \@path2, $fraction )
2051 :     }
2052 :    
2053 :    
2054 :     #-------------------------------------------------------------------------------
2055 :     # Reroot a newick tree along the path between 2 nodes defined by paths:
2056 :     #
2057 :     # $tree = reroot_newick_between_nodes_by_path( $tree, $path1, $path2, $fraction )
2058 :     #-------------------------------------------------------------------------------
2059 :     sub reroot_newick_between_nodes_by_path
2060 :     {
2061 :     my ( $tree, $path1, $path2, $fraction ) = @_;
2062 :     array_ref( $tree ) and array_ref( $path1 ) and array_ref( $path2 )
2063 :     or return undef;
2064 :     $fraction >= 0 && $fraction <= 1 or return undef;
2065 :    
2066 :     my @path1 = @$path1;
2067 :     my @path2 = @$path2;
2068 :    
2069 : golsen 1.9 # Trim the common prefix, saving it:
2070 :    
2071 :     my @prefix = ();
2072 : golsen 1.15 while ( defined( $path1[1] ) && defined( $path2[1] ) && ( $path1[1] == $path2[1] ) )
2073 : golsen 1.9 {
2074 :     push @prefix, splice( @path1, 0, 2 );
2075 :     splice( @path2, 0, 2 );
2076 :     }
2077 :    
2078 :     my ( @path, $dist );
2079 :     if ( @path1 < 3 )
2080 :     {
2081 :     @path2 >= 3 or return undef; # node1 = node2
2082 :     $dist = $fraction * newick_path_length( @path2 );
2083 :     @path = @path2;
2084 :     }
2085 :     elsif ( @path2 < 3 )
2086 :     {
2087 :     $dist = ( 1 - $fraction ) * newick_path_length( @path1 );
2088 :     @path = @path1;
2089 :     }
2090 :     else
2091 :     {
2092 :     my $dist1 = newick_path_length( @path1 );
2093 :     my $dist2 = newick_path_length( @path2 );
2094 :     $dist = $fraction * ( $dist1 + $dist2 ) - $dist1;
2095 :     @path = ( $dist <= 0 ) ? @path1 : @path2;
2096 :     $dist = abs( $dist );
2097 :     }
2098 :    
2099 :     # Descend tree until we reach the insertion branch:
2100 :    
2101 :     my $x;
2102 :     while ( ( $dist > ( $x = newick_x( $path[2] ) ) ) && ( @path > 3 ) )
2103 :     {
2104 :     $dist -= $x;
2105 :     push @prefix, splice( @path, 0, 2 );
2106 :     }
2107 :    
2108 :     # Insert the new node:
2109 :    
2110 :     my $newnode = [ [ $path[2] ], undef, $dist ];
2111 :     set_newick_desc_i( $path[0], $path[1], $newnode );
2112 :     set_newick_x( $path[2], ( ( $x > $dist ) ? ( $x - $dist ) : 0 ) );
2113 :    
2114 :     # We can now build the path from root to the new node
2115 :    
2116 :     reroot_newick_by_path( @prefix, @path[0,1], $newnode );
2117 :     }
2118 :    
2119 :    
2120 :     #-------------------------------------------------------------------------------
2121 : golsen 1.1 # Move root of tree to an approximate midpoint.
2122 :     #
2123 :     # $newtree = reroot_newick_to_approx_midpoint( $tree )
2124 :     #-------------------------------------------------------------------------------
2125 :     sub reroot_newick_to_approx_midpoint {
2126 :     my ( $tree ) = @_;
2127 : golsen 1.5
2128 : golsen 1.1 # Compile average tip to node distances assending
2129 :    
2130 :     my $dists1 = average_to_tips_1( $tree );
2131 :    
2132 : golsen 1.12 # Compile average tip to node distances descending, returning midpoint
2133 :     # cadidates as a list of [ $node1, $node2, $fraction ]
2134 :    
2135 :     my @mids = average_to_tips_2( $dists1, undef, undef );
2136 :    
2137 :     # Reroot to first midpoint candidate
2138 :    
2139 :     return $tree if ! @mids;
2140 :     my ( $node1, $node2, $fraction ) = @{ $mids[0] };
2141 :     reroot_newick_to_node_ref( $tree, $fraction >= 0.5 ? $node2 : $node1 );
2142 :     }
2143 :    
2144 :    
2145 :     #-------------------------------------------------------------------------------
2146 :     # Move root of tree to a midpoint.
2147 :     #
2148 :     # $newtree = reroot_newick_to_midpoint( $tree )
2149 :     #-------------------------------------------------------------------------------
2150 :     sub reroot_newick_to_midpoint {
2151 :     my ( $tree ) = @_;
2152 :    
2153 :     # Compile average tip to node distances assending
2154 :    
2155 :     my $dists1 = average_to_tips_1( $tree );
2156 : golsen 1.1
2157 : golsen 1.12 # Compile average tip to node distances descending, returning midpoint
2158 :     # [ $node1, $node2, $fraction ]
2159 : golsen 1.1
2160 : golsen 1.12 my @mids = average_to_tips_2( $dists1, undef, undef );
2161 : golsen 1.1
2162 : golsen 1.12 @mids ? reroot_newick_between_node_refs( $tree, @{ $mids[0] } ) : $tree;
2163 : golsen 1.1 }
2164 :    
2165 :    
2166 : golsen 1.12 #-------------------------------------------------------------------------------
2167 :     # Compile average tip to node distances assending
2168 :     #-------------------------------------------------------------------------------
2169 : golsen 1.1 sub average_to_tips_1 {
2170 :     my ( $node ) = @_;
2171 :    
2172 :     my @desc_dists = map { average_to_tips_1( $_ ) } newick_desc_list( $node );
2173 :     my $x_below = 0;
2174 :     if ( @desc_dists )
2175 :     {
2176 :     foreach ( @desc_dists ) { $x_below += $_->[0] }
2177 :     $x_below /= @desc_dists;
2178 :     }
2179 : golsen 1.12
2180 : golsen 1.1 my $x = newick_x( $node ) || 0;
2181 :     my $x_net = $x_below + $x;
2182 :    
2183 :     [ $x_net, $x, $x_below, [ @desc_dists ], $node ]
2184 :     }
2185 :    
2186 : golsen 1.9
2187 : golsen 1.12 #-------------------------------------------------------------------------------
2188 :     # Compile average tip to node distances descending, returning midpoint as
2189 :     # [ $node1, $node2, $fraction_of_dist_between ]
2190 :     #-------------------------------------------------------------------------------
2191 : golsen 1.1 sub average_to_tips_2 {
2192 :     my ( $dists1, $x_above, $anc_node ) = @_;
2193 :     my ( undef, $x, $x_below, $desc_list, $node ) = @$dists1;
2194 :    
2195 :     # Are we done? Root is in this node's branch, or "above"?
2196 :    
2197 : golsen 1.12 my @mids = ();
2198 : golsen 1.1 if ( defined( $x_above ) && ( ( $x_above + $x ) >= $x_below ) )
2199 :     {
2200 :     # At this point the root can only be in this node's branch,
2201 :     # or "above" it in the current rooting of the tree (which
2202 :     # would mean that the midpoint is actually down a different
2203 :     # path from the root of the current tree).
2204 :     #
2205 :     # Is the root in the current branch?
2206 :    
2207 :     if ( ( $x_below + $x ) >= $x_above )
2208 :     {
2209 : golsen 1.12 # We will need to make a new node for the root, $fract of
2210 :     # the way from $node to $anc_node:
2211 :     my $fract = ( $x > 0 ) ? 0.5 * ( ( $x_above - $x_below ) / $x + 1 )
2212 :     : 0.5;
2213 :     push @mids, [ $node, $anc_node, $fract ];
2214 : golsen 1.1 }
2215 :     }
2216 :    
2217 : golsen 1.12 # The root might be somewhere below this node:
2218 : golsen 1.1
2219 : golsen 1.5 my $n_1 = @$desc_list - ( $anc_node ? 0 : 1 );
2220 : golsen 1.1 my $ttl_dist = ( @$desc_list * $x_below ) + ( defined( $x_above ) ? ( $x_above + $x ) : 0 );
2221 :    
2222 :     foreach ( @$desc_list )
2223 :     {
2224 :     # If input tree is tip_rooted, $n-1 can be 0, so:
2225 :    
2226 :     my $above2 = $n_1 ? ( ( $ttl_dist - $_->[0] ) / $n_1 ) : 0;
2227 : golsen 1.12 push @mids, average_to_tips_2( $_, $above2, $node );
2228 : golsen 1.1 }
2229 :    
2230 : golsen 1.12 return @mids;
2231 : golsen 1.1 }
2232 :    
2233 : golsen 1.9
2234 : golsen 1.1 #-------------------------------------------------------------------------------
2235 : golsen 1.5 # Move root of tree to an approximate midpoint. Weight by tips.
2236 :     #
2237 :     # $newtree = reroot_newick_to_approx_midpoint_w( $tree )
2238 :     #-------------------------------------------------------------------------------
2239 :     sub reroot_newick_to_approx_midpoint_w {
2240 :     my ( $tree ) = @_;
2241 : golsen 1.21 array_ref( $tree ) or return undef;
2242 : golsen 1.5
2243 : golsen 1.12 # Compile average tip to node distances assending from tips
2244 :    
2245 :     my $dists1 = average_to_tips_1_w( $tree );
2246 :    
2247 :     # Compile average tip to node distances descending, returning midpoints
2248 :    
2249 :     my @mids = average_to_tips_2_w( $dists1, undef, undef, undef );
2250 :    
2251 :     # Reroot to first midpoint candidate
2252 :    
2253 :     return $tree if ! @mids;
2254 :     my ( $node1, $node2, $fraction ) = @{ $mids[0] };
2255 :     reroot_newick_to_node_ref( $tree, $fraction >= 0.5 ? $node2 : $node1 );
2256 :     }
2257 :    
2258 :    
2259 :     #-------------------------------------------------------------------------------
2260 :     # Move root of tree to an approximate midpoint. Weight by tips.
2261 :     #
2262 :     # $newtree = reroot_newick_to_midpoint_w( $tree )
2263 :     #-------------------------------------------------------------------------------
2264 :     sub reroot_newick_to_midpoint_w {
2265 :     my ( $tree ) = @_;
2266 : golsen 1.21 array_ref( $tree ) or return ();
2267 : golsen 1.12
2268 : golsen 1.5 # Compile average tip to node distances assending
2269 :    
2270 :     my $dists1 = average_to_tips_1_w( $tree );
2271 :    
2272 :     # Compile average tip to node distances descending, returning midpoint node
2273 :    
2274 : golsen 1.12 my @mids = average_to_tips_2_w( $dists1, undef, undef, undef );
2275 : golsen 1.5
2276 : golsen 1.12 # Reroot at first candidate midpoint
2277 : golsen 1.5
2278 : golsen 1.12 @mids ? reroot_newick_between_node_refs( $tree, @{ $mids[0] } ) : $tree;
2279 : golsen 1.5 }
2280 :    
2281 :    
2282 : golsen 1.21 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2283 : golsen 1.5 sub average_to_tips_1_w {
2284 :     my ( $node ) = @_;
2285 :    
2286 :     my @desc_dists = map { average_to_tips_1_w( $_ ) } newick_desc_list( $node );
2287 :     my $x_below = 0;
2288 :     my $n_below = 1;
2289 :     if ( @desc_dists )
2290 :     {
2291 :     $n_below = 0;
2292 :     my $n;
2293 :     foreach ( @desc_dists )
2294 :     {
2295 :     $n_below += $n = $_->[1];
2296 :     $x_below += $n * $_->[0];
2297 :     }
2298 :     $x_below /= $n_below;
2299 :     }
2300 : golsen 1.12
2301 : golsen 1.5 my $x = newick_x( $node ) || 0;
2302 :     my $x_net = $x_below + $x;
2303 :    
2304 :     [ $x_net, $n_below, $x, $x_below, [ @desc_dists ], $node ]
2305 :     }
2306 :    
2307 :    
2308 : golsen 1.21 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2309 : golsen 1.5 sub average_to_tips_2_w {
2310 :     my ( $dists1, $x_above, $n_above, $anc_node ) = @_;
2311 :     my ( undef, $n_below, $x, $x_below, $desc_list, $node ) = @$dists1;
2312 :    
2313 :     # Are we done? Root is in this node's branch, or "above"?
2314 :    
2315 : golsen 1.12 my @mids = ();
2316 : golsen 1.5 if ( defined( $x_above ) && ( ( $x_above + $x ) >= $x_below ) )
2317 :     {
2318 :     # At this point the root can only be in this node's branch,
2319 :     # or "above" it in the current rooting of the tree (which
2320 :     # would mean that the midpoint is actually down a different
2321 :     # path from the root of the current tree).
2322 :     #
2323 : golsen 1.12 # Is their a root in the current branch?
2324 : golsen 1.5
2325 :     if ( ( $x_below + $x ) >= $x_above )
2326 :     {
2327 : golsen 1.12 # We will need to make a new node for the root, $fract of
2328 :     # the way from $node to $anc_node:
2329 :     my $fract = ( $x > 0 ) ? 0.5 * ( ( $x_above - $x_below ) / $x + 1 )
2330 :     : 0.5;
2331 :     push @mids, [ $node, $anc_node, $fract ];
2332 : golsen 1.5 }
2333 :     }
2334 :    
2335 :     # The root must be some where below this node:
2336 :    
2337 :     $n_above ||= 0;
2338 :     my $n = $n_above + $n_below;
2339 :     my $ttl_w_dist = ( $n_below * $x_below )
2340 :     + ( defined( $x_above ) ? $n_above * ( $x_above + $x ) : 0 );
2341 :    
2342 :     foreach ( @$desc_list )
2343 :     {
2344 :     my $n_2 = $_->[1]; # n in subtree
2345 :     my $n_above2 = $n - $n_2; # tip rooted has 1 above
2346 :    
2347 :     # If input tree is tip_rooted, $n_above2 can be 0, so:
2348 :    
2349 :     my $x_above2 = $n_above2 ? ( ( $ttl_w_dist - $n_2 * $_->[0] ) / $n_above2 )
2350 :     : 0;
2351 : golsen 1.12 push @mids, average_to_tips_2_w( $_, $x_above2, $n_above2 || 1, $node );
2352 : golsen 1.5 }
2353 :    
2354 : golsen 1.12 return @mids;
2355 : golsen 1.5 }
2356 :    
2357 : golsen 1.9
2358 : golsen 1.5 #-------------------------------------------------------------------------------
2359 : golsen 1.1 # Move root of tree from tip to adjacent node.
2360 :     #
2361 :     # $newtree = uproot_tip_rooted_newick( $tree )
2362 :     #-------------------------------------------------------------------------------
2363 :     sub uproot_tip_rooted_newick {
2364 :     my ( $node ) = @_;
2365 :     newick_is_tip_rooted( $node ) || return $node;
2366 :    
2367 :     # Path to the sole descendant:
2368 :    
2369 :     reroot_newick_by_path( $node, 1, newick_desc_i( $node, 1 ) );
2370 :     }
2371 :    
2372 :    
2373 :     #-------------------------------------------------------------------------------
2374 :     # Remove root bifurcation.
2375 :     #
2376 :     # Root node label, label comment and descendant list comment are discarded.
2377 :     #
2378 :     # $newtree = uproot_newick( $tree )
2379 :     #-------------------------------------------------------------------------------
2380 :     sub uproot_newick {
2381 :     my ( $node0 ) = @_;
2382 :     newick_is_rooted( $node0 ) || return $node0;
2383 :    
2384 :     my ( $node1, $node2 ) = newick_desc_list( $node0 );
2385 :    
2386 :     # Ensure that node1 has at least 1 descendant
2387 :    
2388 :     if ( newick_n_desc( $node1 ) ) {
2389 :     push @{ newick_desc_ref( $node1 ) }, $node2; # Add node2 to descend list
2390 :     }
2391 :    
2392 :     # Or node2 has at least 1 descendant
2393 :    
2394 :     elsif ( newick_n_desc( $node2 ) ) {
2395 :     unshift @{ newick_desc_ref( $node2 ) }, $node1; # Add node1 to descend list
2396 :     ( $node1, $node2 ) = ( $node2, $node1 ); # And reverse labels
2397 :     }
2398 :    
2399 :     # We could make this into a tip rooted tree, but for now:
2400 :    
2401 :     else { return $node0 }
2402 :    
2403 :     # Prefix node1 branch to that of node2:
2404 :    
2405 :     add_to_newick_branch( $node2, $node1 );
2406 :     set_newick_x( $node1, undef );
2407 :    
2408 :     # Tree prefix comment lists (as references):
2409 :    
2410 :     my $C10 = newick_c1( $node0 );
2411 :     my $C11 = newick_c1( $node1 );
2412 : golsen 1.9 if ( $C11 && @$C11 ) {
2413 : golsen 1.1 if ( $C10 && @$C10 ) { unshift @$C11, @$C10 } # Prefix to node1 comments
2414 :     }
2415 :     elsif ( $C10 && @$C10 ) {
2416 :     set_newick_c1( $node1, $C10 ) # Or move node0 comments to node1
2417 :     }
2418 :    
2419 :     $node1;
2420 :     }
2421 :    
2422 :    
2423 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2424 :     # Prefix branch of node2 to that of node1:
2425 :     #
2426 :     # $node1 = add_to_newick_branch( $node1, $node2 )
2427 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2428 :     sub add_to_newick_branch {
2429 :     my ( $node1, $node2 ) = @_;
2430 :     array_ref( $node1 ) || die "add_to_newick_branch: arg 1 not array ref\n";
2431 :     array_ref( $node2 ) || die "add_to_newick_branch: arg 2 not array ref\n";
2432 :    
2433 :     # Node structure template:
2434 :     # my ( $DL, $L, $X, $C1, $C2, $C3, $C4, $C5 ) = @$node;
2435 :    
2436 :     # Fix branch lengths for joining of two branches:
2437 :    
2438 :     set_newick_x( $node1, newick_x( $node1 ) + newick_x( $node2 ) );
2439 :    
2440 :     # Merge branch length comments:
2441 :    
2442 :     my $C41 = newick_c4( $node1 ); # Ref to node1 C4
2443 :     my $C42 = newick_c4( $node2 ); # Ref to node2 C4
2444 :     if ( $C41 && @$C41 ) {
2445 :     if ( $C42 && @$C42 ) { unshift @$C41, @$C42 } # Add node2 comment
2446 :     }
2447 :     elsif ( $C42 && @$C42 ) { set_newick_c4( $node1, $C42 ) } # Or move node1 comment
2448 :    
2449 :     my $C51 = newick_c5( $node1 ); # Ref to node1 C5
2450 :     my $C52 = newick_c5( $node2 ); # Ref to node2 C5
2451 :     if ( $C51 && @$C51 ) {
2452 :     if ( $C52 && @$C52 ) { unshift @$C51, @$C52 } # Add node2 comment
2453 :     }
2454 :     elsif ( $C52 && @$C52 ) { set_newick_c5( $node1, $C52 ) } # Or move node1 comment
2455 :    
2456 :     $node1;
2457 :     }
2458 :    
2459 :    
2460 :     #-------------------------------------------------------------------------------
2461 : golsen 1.5 # Collapse zero-length branches to make multifurcation. The original tree
2462 :     # is modified.
2463 :     #
2464 :     # $tree = collapse_zero_length_branches( $tree )
2465 :     # $tree = collapse_zero_length_branches( $tree, $not_root )
2466 :     #-------------------------------------------------------------------------------
2467 :     sub collapse_zero_length_branches {
2468 :     my ( $tree, $not_root ) = @_;
2469 :     array_ref( $tree ) || return undef;
2470 :    
2471 :     my @desc = newick_desc_list( $tree );
2472 :     @desc or return ( $tree ); # Cannot collapse terminal branch
2473 :    
2474 :     # Analyze descendants:
2475 :    
2476 :     $not_root ||= 0;
2477 :     my @new_desc = ();
2478 :     my $changed = 0;
2479 :     foreach ( @desc )
2480 :     {
2481 :     my ( undef, @to_add ) = collapse_zero_length_branches( $_, $not_root+1 );
2482 :     if ( @to_add )
2483 :     {
2484 :     push @new_desc, @to_add;
2485 :     $changed = 1;
2486 :     }
2487 :     else
2488 :     {
2489 :     push @new_desc, $_;
2490 :     }
2491 :     }
2492 :     set_newick_desc_ref( $tree, [ @new_desc ] ) if $changed;
2493 :    
2494 :     # Collapse if not root, not tip and zero (or negative) branch:
2495 :    
2496 :     my $collapse = $not_root && @new_desc && ( newick_x( $tree ) <= 0 ) ? 1 : 0;
2497 :     ( $tree, ( $collapse ? @new_desc : () ) );
2498 :     }
2499 :    
2500 : golsen 1.8 #-------------------------------------------------------------------------------
2501 :     # Add a subtree to a newick tree node:
2502 :     #
2503 :     # $node = newick_insert_at_node( $node, $subtree )
2504 :     #-------------------------------------------------------------------------------
2505 :     sub newick_insert_at_node
2506 :     {
2507 :     my ( $node, $subtree ) = @_;
2508 :     array_ref( $node ) && array_ref( $subtree ) or return undef;
2509 :    
2510 :     # We could check validity of trees, but ....
2511 :    
2512 :     my $dl = newick_desc_ref( $node );
2513 :     if ( array_ref( $dl ) )
2514 :     {
2515 :     push @$dl, $subtree;
2516 :     }
2517 :     else
2518 :     {
2519 :     set_newick_desc_ref( $node, [ $subtree ] );
2520 :     }
2521 :     return $node;
2522 :     }
2523 :    
2524 :    
2525 :     #-------------------------------------------------------------------------------
2526 :     # Insert a subtree into a newick tree along the path between 2 nodes:
2527 :     #
2528 :     # $tree = newick_insert_between_nodes( $tree, $subtree, $node1, $node2, $fraction )
2529 :     #-------------------------------------------------------------------------------
2530 :     sub newick_insert_between_nodes
2531 :     {
2532 :     my ( $tree, $subtree, $node1, $node2, $fraction ) = @_;
2533 :     array_ref( $tree ) && array_ref( $subtree ) or return undef;
2534 :     $fraction >= 0 && $fraction <= 1 or return undef;
2535 :    
2536 :     # Find the paths to the nodes:
2537 :    
2538 :     my @path1 = path_to_node( $tree, $node1 ) or return undef;
2539 :     my @path2 = path_to_node( $tree, $node2 ) or return undef;
2540 :    
2541 :     # Trim the common prefix:
2542 :    
2543 :     while ( $path1[1] == $path2[1] )
2544 :     {
2545 :     splice( @path1, 0, 2 );
2546 :     splice( @path2, 0, 2 );
2547 :     }
2548 :    
2549 :     my ( @path, $dist );
2550 :     if ( @path1 < 3 )
2551 :     {
2552 :     @path2 >= 3 or return undef; # node1 = node2
2553 : golsen 1.9 $dist = $fraction * newick_path_length( @path2 );
2554 : golsen 1.8 @path = @path2;
2555 :     }
2556 :     elsif ( @path2 < 3 )
2557 :     {
2558 :     $dist = ( 1 - $fraction ) * newick_path_length( @path1 );
2559 :     @path = @path1;
2560 :     }
2561 :     else
2562 :     {
2563 :     my $dist1 = newick_path_length( @path1 );
2564 :     my $dist2 = newick_path_length( @path2 );
2565 :     $dist = $fraction * ( $dist1 + $dist2 ) - $dist1;
2566 :     @path = ( $dist <= 0 ) ? @path1 : @path2;
2567 :     $dist = abs( $dist );
2568 :     }
2569 :    
2570 :     # Descend tree until we reach the insertion branch:
2571 :    
2572 :     my $x;
2573 :     while ( ( $dist > ( $x = newick_x( $path[2] ) ) ) && ( @path > 3 ) )
2574 :     {
2575 :     $dist -= $x;
2576 :     splice( @path, 0, 2 );
2577 :     }
2578 :    
2579 :     # Insert the new node:
2580 :    
2581 :     set_newick_desc_i( $path[0], $path[1], [ [ $path[2], $subtree ], undef, $dist ] );
2582 :     set_newick_x( $path[2], ( ( $x > $dist ) ? ( $x - $dist ) : 0 ) );
2583 :    
2584 :     return $tree;
2585 :     }
2586 :    
2587 : golsen 1.5
2588 :     #-------------------------------------------------------------------------------
2589 : golsen 1.1 # Prune one or more tips from a tree:
2590 :     # Caveat: if one tip is listed, the original tree is modified.
2591 : golsen 1.5 # if more than one tip is listed, a copy of the tree is returned
2592 : golsen 1.1 # (even if it is just listing the same tip twice!).
2593 :     #
2594 :     # $newtree = prune_from_newick( $tree, $tip )
2595 :     # $newtree = prune_from_newick( $tree, @tips )
2596 :     # $newtree = prune_from_newick( $tree, \@tips )
2597 :     #-------------------------------------------------------------------------------
2598 :     sub prune_from_newick {
2599 :     my ( $tr, @tips ) = @_;
2600 :     if ( @tips == 1 && ref( $tips[0] ) eq "ARRAY" ) { @tips = @{ $tips[0] } }
2601 : golsen 1.9
2602 : golsen 1.1 if ( @tips == 0 ) { return $tr }
2603 :     if ( @tips == 1 ) { return prune_1_from_newick( $tr, @tips ) }
2604 :    
2605 :     my %del = map { ( $_, 1 ) } @tips;
2606 :     my @keep = grep { ! $del{ $_ } } newick_tip_list( $tr );
2607 :     newick_subtree( $tr, @keep );
2608 :     }
2609 :    
2610 :    
2611 :     #-------------------------------------------------------------------------------
2612 :     # Prune a tip from a tree:
2613 :     #
2614 :     # $newtree = prune_1_from_newick( $tree, $tip )
2615 :     #-------------------------------------------------------------------------------
2616 :     sub prune_1_from_newick {
2617 :     my ( $tr, $tip ) = @_;
2618 :     my @path = path_to_tip( $tr, $tip );
2619 :     if ( @path < 3 ) { return $tr }
2620 :    
2621 :     my $node = $path[-1]; # Node with the tip
2622 :     my $i1 = $path[-2]; # Descendant number of node in ancestor desc list
2623 :     my $anc1 = $path[-3]; # Ancestor of node
2624 :     my $nd1 = newick_n_desc( $anc1 ); # Number of descendants of ancestor
2625 :     my $anc2 = ( @path >= 5 ) ? $path[-5] : undef; # Ancestor of anc1
2626 :    
2627 :     # dump_tree( $node );
2628 :     # print STDERR "i1 = $i1\n";
2629 :     # dump_tree( $anc1 );
2630 :     # print STDERR "nd1 = $nd1\n";
2631 :     # defined( $anc2 ) && dump_tree( $anc2 );
2632 :    
2633 :     if ( $nd1 > 3 || ( $anc2 && $nd1 > 2 ) ) { # Tip joins at multifurcation
2634 :     splice( @{ $anc1->[0] }, $i1-1, 1 ); # delete the descendant
2635 :     }
2636 :    
2637 :     elsif ( $anc2 ) { # Tip joins at internal bifurcation
2638 :     my $sis = newick_desc_i( $anc1, 3-$i1 ); # find sister node
2639 :     add_to_newick_branch( $sis, $anc1 ); # combine internal branches
2640 :     set_newick_desc_i( $anc2, $path[-4], $sis ); # remove $anc1
2641 :     }
2642 :    
2643 :     elsif ( $nd1 == 2) { # Tip joins bifurcating root node
2644 :     my $sis = newick_desc_i( $anc1, 3-$i1 ); # find sister node
2645 :     $sis->[1] = $anc1->[1] if ! $sis->[1] && $anc1->[1]; # root label
2646 :     $sis->[2] = undef; # root branch len
2647 :     $sis->[3] = $anc1->[3] if ! $sis->[3] && $anc1->[3]; # tree comment
2648 :     $sis->[4] = $anc1->[4] if ! $sis->[4] && $anc1->[4]; # desc list comment
2649 : golsen 1.9 $sis->[5] = $anc1->[5] if ! $sis->[5] && $anc1->[5]; # label comment
2650 : golsen 1.1 $sis->[6] = undef if $sis->[6]; # root branch comment
2651 :     $sis->[7] = undef if $sis->[7]; # root branch comment
2652 :     $tr = $sis; # sister is new root
2653 :     }
2654 :    
2655 :     elsif ( $nd1 == 3 ) { # Tip joins trifurcating root:
2656 :     splice( @{ $anc1->[0] }, $i1-1, 1 ); # delete the descendant, and
2657 :     $tr = uproot_newick( $tr ); # fix the rooting
2658 :     }
2659 :    
2660 :     else {
2661 :     return undef;
2662 :     }
2663 :    
2664 :     return $tr;
2665 :     }
2666 :    
2667 :    
2668 :     #-------------------------------------------------------------------------------
2669 : golsen 1.12 # Produce a potentially rooted subtree with the desired tips:
2670 :     #
2671 :     # Except for (some) tip nodes, the tree produced is a copy.
2672 :     # There is no check that requested tips exist.
2673 :     #
2674 :     # $newtree = rooted_newick_subtree( $tree, @tips )
2675 :     # $newtree = rooted_newick_subtree( $tree, \@tips )
2676 :     #-------------------------------------------------------------------------------
2677 :     sub rooted_newick_subtree {
2678 :     my ( $tr, @tips ) = @_;
2679 :     if ( @tips == 1 && ref( $tips[0] ) eq "ARRAY" ) { @tips = @{ $tips[0] } }
2680 :    
2681 :     if ( @tips < 2 ) { return undef }
2682 :     my $keephash = { map { ( $_, 1 ) } @tips };
2683 :     my $tr2 = subtree1( $tr, $keephash );
2684 :     $tr2->[2] = undef if $tr2; # undef root branch length
2685 :     $tr2;
2686 :     }
2687 :    
2688 :    
2689 :     #-------------------------------------------------------------------------------
2690 : golsen 1.1 # Produce a subtree with the desired tips:
2691 :     #
2692 :     # Except for (some) tip nodes, the tree produced is a copy.
2693 :     # There is no check that requested tips exist.
2694 :     #
2695 :     # $newtree = newick_subtree( $tree, @tips )
2696 :     # $newtree = newick_subtree( $tree, \@tips )
2697 :     #-------------------------------------------------------------------------------
2698 :     sub newick_subtree {
2699 :     my ( $tr, @tips ) = @_;
2700 :     if ( @tips == 1 && ref( $tips[0] ) eq "ARRAY" ) { @tips = @{ $tips[0] } }
2701 :    
2702 :     if ( @tips < 2 ) { return undef }
2703 :     my $was_rooted = newick_is_rooted( $tr );
2704 :     my $keephash = { map { ( $_, 1 ) } @tips };
2705 :     my $tr2 = subtree1( $tr, $keephash );
2706 :     $tr2 = uproot_newick( $tr2 ) if ! $was_rooted && newick_is_rooted( $tr2 );
2707 :     $tr2->[2] = undef if $tr2; # undef root branch length
2708 :     $tr2;
2709 :     }
2710 :    
2711 :    
2712 :     sub subtree1 {
2713 :     my ( $tr, $keep ) = @_;
2714 :     my @desc1 = newick_desc_list( $tr );
2715 :    
2716 :     # Is this a tip, and is it in the keep list?
2717 :    
2718 :     if ( @desc1 < 1 ) {
2719 :     return ( $keep->{ newick_lbl( $tr ) } ) ? $tr : undef;
2720 :     }
2721 :    
2722 :     # Internal node: analyze the descendants:
2723 :    
2724 :     my @desc2 = ();
2725 :     foreach ( @desc1 ) {
2726 :     my $desc = subtree1( $_, $keep );
2727 :     if ( $desc && @$desc ) { push @desc2, $desc }
2728 :     }
2729 :    
2730 :     if ( @desc2 == 0 ) { return undef }
2731 :     if ( @desc2 > 1 ) { return [ \@desc2, @$tr[ 1 .. @$tr - 1 ] ] }
2732 :    
2733 :     # Exactly 1 descendant
2734 :    
2735 :     my $desc = $desc2[ 0 ];
2736 :     my @nn = ( $desc->[0],
2737 :     $desc->[1] ? $desc->[1] : $tr->[1],
2738 :     defined( $tr->[2] ) ? $desc->[2] + $tr->[2] : undef
2739 :     );
2740 :    
2741 :     # Merge comments (only recreating the ones that existed):
2742 :    
2743 :     if ( $tr->[3] && @{$tr->[3]} || $desc->[3] && @{$desc->[3]} ) {
2744 :     $nn[3] = [ $tr->[3] ? @{$tr->[3]} : (), $desc->[3] ? @{$desc->[3]} : () ];
2745 :     }
2746 :     if ( $tr->[4] && @{$tr->[4]} || $desc->[4] && @{$desc->[4]} ) {
2747 :     $nn[4] = [ $tr->[4] ? @{$tr->[4]} : (), $desc->[4] ? @{$desc->[4]} : () ];
2748 :     }
2749 :     if ( $tr->[5] && @{$tr->[5]} || $desc->[5] && @{$desc->[5]} ) {
2750 :     $nn[5] = [ $tr->[5] ? @{$tr->[5]} : (), $desc->[5] ? @{$desc->[5]} : () ];
2751 :     }
2752 :     if ( $tr->[6] && @{$tr->[6]} || $desc->[6] && @{$desc->[6]} ) {
2753 :     $nn[6] = [ $tr->[6] ? @{$tr->[6]} : (), $desc->[6] ? @{$desc->[6]} : () ];
2754 :     }
2755 :     if ( $tr->[7] && @{$tr->[7]} || $desc->[7] && @{$desc->[7]} ) {
2756 :     $nn[7] = [ $tr->[7] ? @{$tr->[7]} : (), $desc->[7] ? @{$desc->[7]} : () ];
2757 :     }
2758 :    
2759 :     return \@nn;
2760 :     }
2761 :    
2762 :    
2763 : golsen 1.12 #-------------------------------------------------------------------------------
2764 :     # The smallest subtree of rooted tree that includes @tips:
2765 :     #
2766 :     # $node = newick_covering_subtree( $tree, @tips )
2767 :     # $node = newick_covering_subtree( $tree, \@tips )
2768 :     #-------------------------------------------------------------------------------
2769 :    
2770 :     sub newick_covering_subtree {
2771 :     my $tree = shift;
2772 :     my %tips = map { $_ => 1 } ( ( ref( $_[0] ) eq 'ARRAY' ) ? @{ $_[0] } : @_ );
2773 :    
2774 :     # Return smallest covering node, if any:
2775 :    
2776 :     ( newick_covering_subtree( $tree, \%tips ) )[ 0 ];
2777 :     }
2778 :    
2779 :    
2780 :     sub newick_covering_subtree_1 {
2781 :     my ( $node, $tips ) = @_;
2782 :     my $n_cover = 0;
2783 :     my @desc = newick_desc_list( $node );
2784 :     if ( @desc )
2785 :     {
2786 :     foreach ( @desc )
2787 :     {
2788 :     my ( $subtree, $n ) = newick_covering_subtree_1( $_, $tips );
2789 :     return ( $subtree, $n ) if $subtree;
2790 :     $n_cover += $n;
2791 :     }
2792 :     }
2793 :     elsif ( $tips->{ newick_lbl( $node ) } )
2794 :     {
2795 :     $n_cover++;
2796 :     }
2797 :    
2798 :     # If all tips are covered, return node
2799 :    
2800 :     ( $n_cover == keys %$tips ) ? ( $node, $n_cover ) : ( undef, $n_cover );
2801 :     }
2802 :    
2803 :    
2804 : golsen 1.1 #===============================================================================
2805 :     #
2806 : golsen 1.9 # Representative subtrees
2807 :     #
2808 :     #===============================================================================
2809 :     # Find subtree of size n representating vicinity of the root:
2810 :     #
2811 :     # $subtree = root_neighborhood_representative_tree( $tree, $n, \%tip_priority )
2812 :     # $subtree = root_neighborhood_representative_tree( $tree, $n )
2813 :     #
2814 :     # Note that if $tree is rooted, then the subtree will also be. This can have
2815 :     # consequences on downstream programs.
2816 :     #-------------------------------------------------------------------------------
2817 :     sub root_neighborhood_representative_tree
2818 :     {
2819 :     my ( $tree, $n, $tip_priority ) = @_;
2820 :     array_ref( $tree ) && ( $n >= 2 ) or return undef;
2821 :     if ( newick_tip_count( $tree ) <= $n ) { return $tree }
2822 :    
2823 :     $tip_priority ||= default_tip_priority( $tree );
2824 :     my @tips = map { representative_tip_of_newick_node( $_, $tip_priority ) }
2825 :     root_proximal_newick_subtrees( $tree, $n );
2826 :    
2827 :     newick_subtree( copy_newick_tree( $tree ), \@tips );
2828 :     }
2829 :    
2830 :    
2831 :     #-------------------------------------------------------------------------------
2832 :     # Find n tips to represent tree lineages in vicinity of another tip.
2833 :     # Default tip priority is short total branch length.
2834 :     #
2835 :     # \@tips = root_neighborhood_representative_tips( $tree, $n, \%tip_priority )
2836 :     # @tips = root_neighborhood_representative_tips( $tree, $n, \%tip_priority )
2837 :     # \@tips = root_neighborhood_representative_tips( $tree, $n )
2838 :     # @tips = root_neighborhood_representative_tips( $tree, $n )
2839 :     #-------------------------------------------------------------------------------
2840 :     sub root_neighborhood_representative_tips
2841 :     {
2842 :     my ( $tree, $n, $tip_priority ) = @_;
2843 :     array_ref( $tree ) && ( $n >= 2 ) or return undef;
2844 :    
2845 :     my @tips;
2846 :     if ( newick_tip_count( $tree ) <= $n )
2847 :     {
2848 :     @tips = newick_tip_list( $tree );
2849 :     }
2850 :     else
2851 :     {
2852 :     $tip_priority ||= default_tip_priority( $tree );
2853 :     @tips = map { representative_tip_of_newick_node( $_, $tip_priority ) }
2854 :     root_proximal_newick_subtrees( $tree, $n );
2855 :     }
2856 :    
2857 :     wantarray ? @tips : \@tips;
2858 :     }
2859 :    
2860 :    
2861 :     #-------------------------------------------------------------------------------
2862 :     # Find subtree of size n representating vicinity of a tip:
2863 :     #
2864 :     # $subtree = tip_neighborhood_representative_tree( $tree, $tip, $n, \%tip_priority )
2865 :     # $subtree = tip_neighborhood_representative_tree( $tree, $tip, $n )
2866 :     #-------------------------------------------------------------------------------
2867 :     sub tip_neighborhood_representative_tree
2868 :     {
2869 :     my ( $tree, $tip, $n, $tip_priority ) = @_;
2870 :     array_ref( $tree ) && $tip && ( $n >= 2 ) or return undef;
2871 :     newick_tip_in_tree( $tree, $tip ) or return undef;
2872 :    
2873 :     my $tree1 = copy_newick_tree( $tree );
2874 :     if ( newick_tip_count( $tree1 ) - 1 <= $n )
2875 :     {
2876 :     return prune_from_newick( $tree1, $tip )
2877 :     }
2878 :    
2879 :     $tree1 = reroot_newick_to_tip( $tree1, $tip );
2880 :     $tree1 = newick_desc_i( $tree1, 1 ); # Node immediately below tip
2881 :     my @tips = root_neighborhood_representative_tips( $tree1, $n, $tip_priority );
2882 :     newick_subtree( copy_newick_tree( $tree ), \@tips );
2883 :     }
2884 :    
2885 :    
2886 :     #-------------------------------------------------------------------------------
2887 :     # Find n tips to represent tree lineages in vicinity of another tip.
2888 :     # Default tip priority is short total branch length.
2889 :     #
2890 :     # \@tips = tip_neighborhood_representative_tips( $tree, $tip, $n, \%tip_priority )
2891 :     # @tips = tip_neighborhood_representative_tips( $tree, $tip, $n, \%tip_priority )
2892 :     # \@tips = tip_neighborhood_representative_tips( $tree, $tip, $n )
2893 :     # @tips = tip_neighborhood_representative_tips( $tree, $tip, $n )
2894 :     #-------------------------------------------------------------------------------
2895 :     sub tip_neighborhood_representative_tips
2896 :     {
2897 :     my ( $tree, $tip, $n, $tip_priority ) = @_;
2898 :     array_ref( $tree ) && $tip && ( $n >= 2 ) or return undef;
2899 :     newick_tip_in_tree( $tree, $tip ) or return undef;
2900 :    
2901 :     my @tips = newick_tip_list( $tree );
2902 :     if ( newick_tip_count( $tree ) - 1 <= $n )
2903 :     {
2904 :     @tips = grep { $_ ne $tip } @tips;
2905 :     }
2906 :     else
2907 :     {
2908 :     my $tree1 = copy_newick_tree( $tree );
2909 :     $tree1 = reroot_newick_to_tip( $tree1, $tip );
2910 :     $tree1 = newick_desc_i( $tree1, 1 ); # Node immediately below tip
2911 :     @tips = root_neighborhood_representative_tips( $tree1, $n, $tip_priority );
2912 :     }
2913 :    
2914 :     wantarray ? @tips : \@tips;
2915 :     }
2916 :    
2917 :    
2918 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2919 :     # Anonymous hash of the negative distance from root to each tip:
2920 :     #
2921 :     # \%tip_priority = default_tip_priority( $tree )
2922 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2923 :     sub default_tip_priority
2924 :     {
2925 :     my ( $tree ) = @_;
2926 :     my $tip_distances = newick_tip_distances( $tree ) || {};
2927 :     return { map { $_ => -$tip_distances->{$_} } keys %$tip_distances };
2928 :     }
2929 :    
2930 :    
2931 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2932 :     # Select a tip from a subtree base on a priority value:
2933 :     #
2934 :     # $tip = representative_tip_of_newick_node( $node, \%tip_priority )
2935 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2936 :     sub representative_tip_of_newick_node
2937 :     {
2938 :     my ( $node, $tip_priority ) = @_;
2939 :     my ( $tip ) = sort { $b->[1] <=> $a->[1] } # The best
2940 :     map { [ $_, $tip_priority->{ $_ } ] }
2941 :     newick_tip_list( $node );
2942 :     $tip->[0]; # Label from label-priority pair
2943 :     }
2944 :    
2945 :    
2946 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2947 :     # Find n subtrees focused around the root of a tree. Typically each will
2948 :     # then be reduced to a single tip to make a representative tree:
2949 :     #
2950 :     # @subtrees = root_proximal_newick_subtrees( $tree, $n )
2951 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2952 :     sub root_proximal_newick_subtrees
2953 :     {
2954 :     my ( $tree, $n ) = @_;
2955 :     my $node_start_end = newick_branch_intervals( $tree );
2956 :     n_representative_branches( $n, $node_start_end );
2957 :     }
2958 :    
2959 :    
2960 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2961 :     # @node_start_end = newick_branch_intervals( $tree )
2962 :     # \@node_start_end = newick_branch_intervals( $tree )
2963 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2964 :     sub newick_branch_intervals
2965 :     {
2966 :     my ( $node, $parent_x ) = @_;
2967 :     $parent_x ||= 0;
2968 :     my ( $desc, undef, $dx ) = @$node;
2969 :     my $x = $parent_x + $dx;
2970 :     my $interval = [ $node, $parent_x, $desc && @$desc ? $x : 1e100 ];
2971 :     my @intervals = ( $interval,
2972 :     map { &newick_branch_intervals( $_, $x ) } @$desc
2973 :     );
2974 :     return wantarray ? @intervals : \@intervals;
2975 :     }
2976 :    
2977 :    
2978 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2979 :     # @ids = n_representative_branches( $n, @id_start_end )
2980 :     # @ids = n_representative_branches( $n, \@id_start_end )
2981 :     # \@ids = n_representative_branches( $n, @id_start_end )
2982 :     # \@ids = n_representative_branches( $n, \@id_start_end )
2983 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2984 :     sub n_representative_branches
2985 :     {
2986 :     my $n = shift;
2987 :     # Sort intervals by start point:
2988 :     my @unprocessed = sort { $a->[1] <=> $b->[1] }
2989 :     ( @_ == 1 ) ? @{ $_[0] } : @_;
2990 :     my @active = ();
2991 :     my ( $interval, $current_point );
2992 :     foreach $interval ( @unprocessed )
2993 :     {
2994 :     $current_point = $interval->[1];
2995 :     # Filter out intervals that have ended. This is N**2 in the number
2996 :     # of representatives. Fixing this would require maintaining a sorted
2997 :     # active list.
2998 :     @active = grep { $_->[2] > $current_point } @active;
2999 :     push @active, $interval;
3000 :     last if ( @active >= $n );
3001 :     }
3002 :    
3003 :     my @ids = map { $_->[0] } @active;
3004 :     return wantarray() ? @ids : \@ids;
3005 :     }
3006 :    
3007 :    
3008 :     #===============================================================================
3009 : golsen 1.22 # Random trees
3010 :     #===============================================================================
3011 :     #
3012 :     # $tree = random_equibranch_tree( @tips, \%options )
3013 :     # $tree = random_equibranch_tree( \@tips, \%options )
3014 :     # $tree = random_equibranch_tree( @tips )
3015 :     # $tree = random_equibranch_tree( \@tips )
3016 :     #
3017 :     # Options:
3018 :     #
3019 :     # length => $branch_length # D = 1
3020 :     #
3021 :     #-------------------------------------------------------------------------------
3022 :     sub random_equibranch_tree
3023 :     {
3024 :     my $opts = $_[ 0] && ref $_[ 0] eq 'HASH' ? shift
3025 :     : $_[-1] && ref $_[-1] eq 'HASH' ? pop
3026 :     : {};
3027 :     return undef if ! defined $_[0];
3028 :    
3029 :     my @tips = ref $_[0] ? @{ $_[0] } : @_;
3030 :     return undef if @tips < 2;
3031 :    
3032 :     my $len = $opts->{ length } ||= 1;
3033 :    
3034 :     if ( @tips == 2 )
3035 :     {
3036 :     return [ [ map { [ [], $_, $len ] } @tips ], undef, 0 ];
3037 :     }
3038 :    
3039 :     my $tree = [ [ ], undef, 0 ];
3040 :    
3041 :     my @links; # \$anc_dl[i], i.e. a reference to an element in a descendent list
3042 :    
3043 :     my $anc_dl = $tree->[0];
3044 :     foreach my $tip ( splice( @tips, 0, 3 ) )
3045 :     {
3046 :     my $node = [ [], $tip, $len ];
3047 :     push @$anc_dl, $node;
3048 :     push @links, \$anc_dl->[-1]; # Ref to the just added descendent list entry
3049 :     }
3050 :    
3051 :     foreach my $tip ( @tips )
3052 :     {
3053 :     my $link = $links[ int( rand( scalar @links ) ) ];
3054 :     my $newtip = [ [], $tip, $len ];
3055 :     my $new_dl = [ $$link, $newtip ];
3056 :     my $newnode = [ $new_dl, undef, $len ];
3057 :     $$link = $newnode;
3058 :     push @links, \$new_dl->[0], \$new_dl->[1]
3059 :     }
3060 :    
3061 :     return $tree;
3062 :     }
3063 :    
3064 :    
3065 :     #-------------------------------------------------------------------------------
3066 :     #
3067 :     # $tree = random_ultrametric_tree( @tips, \%options )
3068 :     # $tree = random_ultrametric_tree( \@tips, \%options )
3069 :     # $tree = random_ultrametric_tree( @tips )
3070 :     # $tree = random_ultrametric_tree( \@tips )
3071 :     #
3072 :     # Options:
3073 :     #
3074 :     # depth => $root_to_tip_dist # D = 1
3075 :     #
3076 :     #-------------------------------------------------------------------------------
3077 :     sub random_ultrametric_tree
3078 :     {
3079 :     my $opts = $_[ 0] && ref $_[ 0] eq 'HASH' ? shift
3080 :     : $_[-1] && ref $_[-1] eq 'HASH' ? pop
3081 :     : {};
3082 :     return undef if ! defined $_[0];
3083 :    
3084 :     my @tips = ref $_[0] ? @{ $_[0] } : @_;
3085 :     return undef if @tips < 2;
3086 :    
3087 :     my $d2tip = $opts->{ depth } ||= 1;
3088 :    
3089 :     # Random tip addition order (for rooted tree it matters):
3090 :    
3091 :     @tips = sort { rand() <=> 0.5 } @tips;
3092 :     my $tree = [ [ ], undef, 0 ];
3093 :    
3094 :     my $subtree_size = { $tree => 0 }; # total branch length of each subtree
3095 :    
3096 :     # We start with root bifurcation:
3097 :    
3098 :     foreach my $tip ( splice( @tips, 0, 2 ) )
3099 :     {
3100 :     my $node = [ [], $tip, $d2tip ];
3101 :     push @{ $tree->[0] }, $node;
3102 :     $subtree_size->{ $node } = $d2tip;
3103 :     $subtree_size->{ $tree } += $d2tip;
3104 :     }
3105 :    
3106 :     # Add each remaining tip at $pos, measured along the contour length
3107 :     # of the tree (with no retracing along branches).
3108 :    
3109 :     foreach my $tip ( @tips )
3110 :     {
3111 :     my $pos = rand( $subtree_size->{ $tree } );
3112 :     random_add_to_ultrametric_tree( $tree, $tip, $subtree_size, $pos, $d2tip );
3113 :     }
3114 :    
3115 :     return $tree;
3116 :     }
3117 :    
3118 :    
3119 :     sub random_add_to_ultrametric_tree
3120 :     {
3121 :     my ( $node, $tip, $subtree_size, $pos, $d2tip ) = @_;
3122 :     $node && $node->[0] && ref $node->[0] eq 'ARRAY' or die "Bad tree node passed to random_add_to_ultrametric_tree().\n";
3123 :    
3124 :     # Find the descendent line that it goes in:
3125 :    
3126 :     my $i;
3127 :     my $dl = $node->[0];
3128 :     my $size0 = $subtree_size->{ $dl->[0] };
3129 :     if ( $size0 > $pos ) { $i = 0 } else { $i = 1; $pos -= $size0 }
3130 :     my $desc = $dl->[$i];
3131 :    
3132 :     # Does it go within the subtree, or the branch to the subtree?
3133 :    
3134 :     my $len;
3135 :     my $added;
3136 :     if ( ( $len = $desc->[2] ) <= $pos )
3137 :     {
3138 :     $added = random_add_to_ultrametric_tree( $desc, $tip, $subtree_size, $pos - $len, $d2tip - $len );
3139 :     }
3140 :     else
3141 :     {
3142 :     # If not in subtree, then it goes in the branch to the descendent node
3143 :     #
3144 :     # ----- node ------------ node
3145 :     # ^ / \ ^ / \
3146 :     # | \ | pos \l1
3147 :     # | \ v \
3148 :     # | len\ ---------- newnode
3149 :     # | \ / \ l2
3150 :     # d2tip | \ / \
3151 :     # | desc / desc
3152 :     # | / \ l3/ / \
3153 :     # | . . / . .
3154 :     # v . . / . .
3155 :     # ----- . . newtip . .
3156 :    
3157 :     my $l1 = $pos;
3158 :     my $l2 = $len - $pos;
3159 :     my $l3 = $d2tip - $pos;
3160 :     my $newtip = [ [], $tip, $l3 ];
3161 :     my $newnode = [ [ $desc, $newtip ], undef, $l1 ];
3162 :     $dl->[$i] = $newnode;
3163 :     $subtree_size->{ $newtip } = $l3;
3164 :     $subtree_size->{ $newnode } = $subtree_size->{ $desc } + $l3;
3165 :     $desc->[2] = $l2;
3166 :     $subtree_size->{ $desc } -= $l1;
3167 :     $added = $l3;
3168 :     }
3169 :    
3170 :     # New branch was inserted below this point:
3171 :    
3172 :     $subtree_size->{ $node } += $added;
3173 :     return $added;
3174 :     }
3175 :    
3176 :    
3177 :    
3178 :     #===============================================================================
3179 : golsen 1.9 #
3180 : golsen 1.1 # Tree writing and reading
3181 :     #
3182 :     #===============================================================================
3183 : overbeek 1.7 # writeNewickTree( $tree )
3184 :     # writeNewickTree( $tree, $file )
3185 :     # writeNewickTree( $tree, \*FH )
3186 : golsen 1.1 #-------------------------------------------------------------------------------
3187 :     sub writeNewickTree {
3188 :     my ( $tree, $file ) = @_;
3189 : overbeek 1.7 my ( $fh, $close ) = open_output( $file );
3190 :     $fh or return;
3191 :     print $fh ( strNewickTree( $tree ), "\n" );
3192 :     close $fh if $close;
3193 : golsen 1.1 }
3194 :    
3195 :    
3196 :     #-------------------------------------------------------------------------------
3197 :     # fwriteNewickTree( $file, $tree ) # Args reversed to writeNewickTree
3198 :     #-------------------------------------------------------------------------------
3199 :     sub fwriteNewickTree { writeNewickTree( $_[1], $_[0] ) }
3200 :    
3201 :    
3202 :     #-------------------------------------------------------------------------------
3203 :     # $treestring = strNewickTree( $tree )
3204 :     #-------------------------------------------------------------------------------
3205 :     sub strNewickTree {
3206 :     my $node = shift @_;
3207 :     strNewickSubtree( $node, "" ) . ";";
3208 :     }
3209 :    
3210 :    
3211 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3212 :     # $string = strNewickSubtree( $node, $prefix )
3213 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3214 :     sub strNewickSubtree {
3215 :     my ( $node, $prefix ) = @_;
3216 :     my $s;
3217 :    
3218 :     $s = strNewickComments( newick_c1( $node ), $prefix );
3219 :     if ( $s ) { $prefix = " " }
3220 :    
3221 :     my $ndesc;
3222 :     if ( $ndesc = newick_n_desc( $node ) ) {
3223 :     for (my $d = 1; $d <= $ndesc; $d++) {
3224 :     $s .= ( ( $d == 1 ) ? $prefix . "(" : "," )
3225 :     . strNewickSubtree( newick_desc_i( $node, $d ), " " );
3226 :     }
3227 :    
3228 :     $s .= ")" . strNewickComments( newick_c2( $node ), " " );
3229 :     $prefix = " ";
3230 :     }
3231 :    
3232 :     if ( defined( newick_lbl( $node ) ) && newick_lbl( $node ) ) {
3233 :     $s .= $prefix
3234 :     . q_newick_lbl( $node )
3235 :     . strNewickComments( newick_c3( $node ), " " );
3236 :     }
3237 :    
3238 :     if ( defined( newick_x( $node ) ) ) {
3239 :     $s .= ":"
3240 :     . strNewickComments( newick_c4( $node ), " " )
3241 :     . sprintf( " %.6f", newick_x( $node ) )
3242 :     . strNewickComments( newick_c5( $node ), " " );
3243 :     }
3244 :    
3245 :     $s;
3246 :     }
3247 :    
3248 :    
3249 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3250 :     # $string = strNewickComments( $clist, $prefix )
3251 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3252 :     sub strNewickComments {
3253 :     my ( $clist, $prefix ) = @_;
3254 :     array_ref( $clist ) && ( @$clist > 0 ) || return "";
3255 :     $prefix . "[" . join( "] [", @$clist ) . "]";
3256 :     }
3257 :    
3258 :    
3259 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3260 :     # $quoted_label = q_newick_lbl( $label )
3261 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3262 :     sub q_newick_lbl {
3263 :     my $lbl = newick_lbl( $_[0] );
3264 :     defined( $lbl ) && ( $lbl ne "" ) || return undef;
3265 :    
3266 :     if ( $lbl =~ m/^[^][()_:;,]+$/ # Anything but []()_:;,
3267 :     && $lbl !~ m/^'/ ) { # and does not start with '
3268 :     $lbl =~ s/ /_/g; # Recode blanks as _
3269 :     return $lbl;
3270 :     }
3271 :    
3272 :     else {
3273 :     $lbl =~ s/'/''/g; # Double existing single quote marks
3274 :     return q(') . $lbl . q('); # Wrap in single quote marks
3275 :     }
3276 :     }
3277 :    
3278 :    
3279 :     #===============================================================================
3280 :     # $treestring = formatNewickTree( $tree )
3281 :     #===============================================================================
3282 :     sub formatNewickTree {
3283 :     formatNewickSubtree( $_[0], "", "" ) . ";";
3284 :     }
3285 :    
3286 :    
3287 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3288 :     # $string = formatNewickSubtree( $node, $prefix, $indent )
3289 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3290 :     sub formatNewickSubtree {
3291 :     my ( $node, $prefix, $indent ) = @_;
3292 :     my $s;
3293 :    
3294 :     $s = formatNewickComments( newick_c1( $node ), $prefix, $indent );
3295 :     if ( $s ) { $prefix = "\n$indent" }
3296 :    
3297 :     if ( my $ndesc = newick_n_desc( $node ) ) {
3298 :     for (my $d = 1; $d <= $ndesc; $d++) {
3299 :     $s .= ( ( $d == 1 ) ? $prefix . "(" : ",\n$indent " )
3300 :     . formatNewickSubtree( newick_desc_i( $node, $d ), " ", $indent . " " );
3301 :     }
3302 :    
3303 :     $s .= "\n$indent)" . formatNewickComments( newick_c2( $node ), " ", $indent );
3304 :     $prefix = " ";
3305 :     }
3306 :    
3307 :     if ( defined( newick_lbl( $node ) ) && newick_lbl( $node ) ) {
3308 :     $s .= $prefix
3309 :     . q_newick_lbl( $node )
3310 :     . formatNewickComments( newick_c3( $node ), " ", $indent );
3311 :     }
3312 :    
3313 :     if ( defined( newick_x( $node ) ) ) {
3314 :     $s .= ":"
3315 :     . formatNewickComments( newick_c4( $node ), " ", $indent )
3316 :     . sprintf(" %.6f", newick_x( $node ) )
3317 :     . formatNewickComments( newick_c5( $node ), " ", $indent );
3318 :     }
3319 :    
3320 :     $s;
3321 :     }
3322 :    
3323 :    
3324 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3325 :     # $string = formatNewickComments( $clist, $prefix, $indent )
3326 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3327 :     sub formatNewickComments {
3328 :     my ( $clist, $prefix, $indent ) = @_;
3329 :     array_ref( $clist ) && @$clist || return "";
3330 :     $prefix . "[" . join( "] [", @$clist ) . "]";
3331 :     }
3332 :    
3333 :    
3334 :     #===============================================================================
3335 : golsen 1.23 # Read to a semicolon
3336 :     #
3337 :     # $tree = read_newick_tree( )
3338 :     # $tree = read_newick_tree( \*FH )
3339 :     # $tree = read_newick_tree( $file )
3340 :     #
3341 :     # Read to end of file:
3342 :     # @trees = read_newick_trees( )
3343 :     # @trees = read_newick_trees( \*FH )
3344 :     # @trees = read_newick_trees( $file )
3345 : overbeek 1.7 #===============================================================================
3346 :    
3347 :     sub read_newick_tree
3348 :     {
3349 :     my $file = shift;
3350 :     my ( $fh, $close ) = open_input( $file );
3351 :     my $tree;
3352 :     my @lines = ();
3353 : golsen 1.9 foreach ( <$fh> )
3354 : overbeek 1.7 {
3355 :     chomp;
3356 :     push @lines, $_;
3357 :     if ( /;/ )
3358 :     {
3359 :     $tree = parse_newick_tree_str( join( ' ', @lines ) );
3360 :     last;
3361 :     }
3362 :     }
3363 :     close $fh if $close;
3364 :    
3365 :     $tree;
3366 :     }
3367 :    
3368 :    
3369 :     sub read_newick_trees
3370 :     {
3371 :     my $file = shift;
3372 :     my ( $fh, $close ) = open_input( $file );
3373 :     my @trees = ();
3374 :     my @lines = ();
3375 : golsen 1.9 foreach ( <$fh> )
3376 : overbeek 1.7 {
3377 :     chomp;
3378 :     push @lines, $_;
3379 :     if ( /;/ )
3380 :     {
3381 :     push @trees, parse_newick_tree_str( join( ' ', @lines ) );
3382 :     @lines = ()
3383 :     }
3384 :     }
3385 :     close $fh if $close;
3386 :    
3387 :     @trees;
3388 :     }
3389 :    
3390 :    
3391 :     #===============================================================================
3392 : golsen 1.1 # Tree reader adapted from the C language reader in fastDNAml
3393 :     #
3394 :     # $tree = parse_newick_tree_str( $string )
3395 :     #===============================================================================
3396 :     sub parse_newick_tree_str {
3397 :     my $s = shift @_;
3398 :    
3399 :     my ( $ind, $rootnode ) = parse_newick_subtree( $s, 0 );
3400 :     if ( substr( $s, $ind, 1 ) ne ";") { warn "warning: tree missing ';'\n" }
3401 :     $rootnode;
3402 :     }
3403 :    
3404 :    
3405 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3406 :     # Read a subtrees recursively (everything of tree but a semicolon)
3407 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3408 :     sub parse_newick_subtree {
3409 :     my ( $s, $ind ) = @_;
3410 :    
3411 :     my $newnode = [];
3412 :     my @dlist = ();
3413 :     my ( $lbl, $x, $c1, $c2, $c3, $c4, $c5 );
3414 : golsen 1.9
3415 : golsen 1.1 ( $ind, $c1 ) = getNextTreeChar( $s, $ind ); # Comment 1
3416 :     if ( ! defined( $ind ) ) { treeParseError( "missing subtree" ) }
3417 :     if ( $c1 && @$c1 ) { set_newick_c1( $newnode, $c1 ) }
3418 :    
3419 :     if ( substr( $s, $ind, 1 ) eq "(" ) { # New internal node
3420 :     while ( ! @dlist || ( substr( $s, $ind, 1 ) eq "," ) ) {
3421 :     my $desc;
3422 :     ( $ind, $desc ) = parse_newick_subtree( $s, $ind+1 );
3423 :     if (! $ind) { return () }
3424 :     push @dlist, $desc;
3425 :     }
3426 :     if ( substr( $s, $ind, 1 ) ne ")" ) { treeParseError( "missing ')'" ) }
3427 :    
3428 :     ( $ind, $c2 ) = getNextTreeChar( $s, $ind+1 ); # Comment 2
3429 :     if ( $c2 && @$c2 ) { set_newick_c2( $newnode, $c2 ) }
3430 :     ( $ind, $lbl ) = parseTreeNodeLabel( $s, $ind ); # Node label
3431 :     }
3432 :    
3433 :     elsif ( substr( $s, $ind, 1 ) =~ /[^][(,):;]/ ) { # New tip
3434 :     ( $ind, $lbl ) = parseTreeNodeLabel( $s, $ind ); # Tip label
3435 :     if (! $ind) { return () }
3436 :     }
3437 :    
3438 :     @dlist || $lbl || treeParseError( "no descendant list or label" );
3439 :    
3440 :     if ( @dlist ) { set_newick_desc_ref( $newnode, \@dlist ) }
3441 :     if ( $lbl ) { set_newick_lbl( $newnode, $lbl ) }
3442 :    
3443 :     ( $ind, $c3 ) = getNextTreeChar( $s, $ind ); # Comment 3
3444 :     if ( $c3 && @$c3 ) { set_newick_c3( $newnode, $c3 ) }
3445 :    
3446 :     if (substr( $s, $ind, 1 ) eq ":") { # Branch length
3447 :     ( $ind, $c4 ) = getNextTreeChar( $s, $ind+1 ); # Comment 4
3448 :     if ( $c4 && @$c4 ) { set_newick_c4( $newnode, $c4 ) }
3449 :     ( $ind, $x ) = parseBranchLength( $s, $ind );
3450 :     if ( defined( $x ) ) { set_newick_x( $newnode, $x ) }
3451 :     ( $ind, $c5 ) = getNextTreeChar( $s, $ind ); # Comment 5
3452 :     if ( $c5 && @$c5 ) { set_newick_c5( $newnode, $c5 ) }
3453 :     }
3454 :    
3455 :     ( $ind, $newnode );
3456 :     }
3457 :    
3458 :    
3459 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3460 :     # Read a Newick tree label
3461 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3462 :     sub parseTreeNodeLabel { # Empty string is permitted
3463 :     my ( $s, $ind ) = @_;
3464 :     my ( $lbl, $c );
3465 :    
3466 :     if ( substr( $s, $ind, 1 ) eq "'") {
3467 :     my $ind1 = ++$ind;
3468 :    
3469 :     while ( ) {
3470 :     if ( ! defined( $c = substr( $s, $ind, 1 ) ) || $c eq "" ) {
3471 :     treeParseError( "missing close quote on label '" . substr( $s, $ind1 ) . "'" )
3472 :     }
3473 :     elsif ( $c ne "'" ) { $ind++ }
3474 :     elsif ( substr( $s, $ind, 2 ) eq "''" ) { $ind += 2 }
3475 :     else { last }
3476 :     }
3477 :    
3478 :     $lbl = substr( $s, $ind1, $ind-$ind1 );
3479 :     $lbl =~ s/''/'/g;
3480 :     $ind++;
3481 :     }
3482 :    
3483 :     else {
3484 :     my $ind1 = $ind;
3485 :     while ( defined( $c = substr($s, $ind, 1) ) && $c ne "" && $c !~ /[][\s(,):;]/ ) { $ind++ }
3486 :     $lbl = substr( $s, $ind1, $ind-$ind1 );
3487 :     $lbl =~ s/_/ /g;
3488 :     }
3489 :    
3490 :     ( $ind, $lbl );
3491 :     }
3492 :    
3493 :    
3494 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3495 :     # Read a Newick tree branch length
3496 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3497 :     sub parseBranchLength {
3498 :     my ( $s, $ind ) = @_;
3499 :    
3500 :     my $c = substr( $s, $ind, 1 );
3501 :    
3502 :     my $sign = ( $c ne "-" ) ? 1 : -1; # Sign
3503 :     if ( $c =~ /[-+]/ ) { $c = substr( $s, ++$ind, 1 ) }
3504 :    
3505 :     if ( $c !~ /^[.0-9]$/ ) { # Allows starting with decimal
3506 :     treeParseError( "invalid branch length character '$c'" )
3507 :     }
3508 :    
3509 :     my $v = 0;
3510 :     while ( $c =~ /[0-9]/ ) { # Whole number
3511 :     $v = 10 * $v + $c;
3512 :     $c = substr( $s, ++$ind, 1 );
3513 :     }
3514 :    
3515 :     if ( $c eq "." ) { # Fraction
3516 :     my $f = 0.1;
3517 :     $c = substr( $s, ++$ind, 1 );
3518 :     while ( $c =~ /[0-9]/ ) {
3519 :     $v += $f * $c;
3520 :     $f *= 0.1;
3521 :     $c = substr( $s, ++$ind, 1 );
3522 :     }
3523 :     }
3524 :    
3525 :     $v *= $sign;
3526 :    
3527 :     if ( $c =~ /[dDeEgG]/ ) { # Exponent
3528 :     $c = substr( $s, ++$ind, 1 );
3529 :     my $esign = ( $c ne "-" ) ? 1 : -1;
3530 :     if ( $c =~ /^[-+]$/ ) { $c = substr( $s, ++$ind, 1 ) }
3531 :     if ( $c !~ /^[0-9]$/ ) {
3532 :     treeParseError( "missing branch length exponent '$c'" )
3533 :     }
3534 :    
3535 :     my $e = 0;
3536 :     while ( $c =~ /[0-9]/ ) {
3537 :     $e = 10 * $e + $c;
3538 :     $c = substr( $s, ++$ind, 1 );
3539 :     }
3540 :     $e *= $esign;
3541 :     $v *= 10**$e;
3542 :     }
3543 :    
3544 :     ( $ind, $v );
3545 :     }
3546 :    
3547 :    
3548 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3549 :     # ( $index, /@commentlist ) = getNextTreeChar( $string, $index )
3550 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3551 :     sub getNextTreeChar { # Move to next nonblank, noncomment character
3552 :     my ( $s, $ind ) = @_;
3553 :    
3554 :     my @clist = ();
3555 :    
3556 :     # Skip white space
3557 :     if ( substr( $s, $ind ) =~ /^(\s+)/ ) { $ind += length( $1 ) }
3558 :    
3559 :     # Loop while it is a comment:
3560 :     while ( substr( $s, $ind, 1 ) eq "[" ) {
3561 :     $ind++;
3562 : golsen 1.11 my $depth = 1;
3563 :     my $ind2 = $ind;
3564 : golsen 1.1
3565 :     # Find end
3566 : golsen 1.11 while ( $depth > 0 )
3567 :     {
3568 :     if ( substr( $s, $ind2 ) =~ /^([^][]*\[)/ ) # nested [ ... ]
3569 :     {
3570 :     $ind2 += length( $1 ); # Points at char just past [
3571 :     $depth++; # If nested comments are allowed
3572 :     }
3573 :     elsif ( substr( $s, $ind2 ) =~ /^([^][]*\])/ ) # close bracket
3574 :     {
3575 :     $ind2 += length( $1 ); # Points at char just past ]
3576 :     $depth--;
3577 :     }
3578 :     else
3579 :     {
3580 :     treeParseError( "comment missing closing bracket '["
3581 :     . substr( $s, $ind ) . "'" )
3582 :     }
3583 : golsen 1.1 }
3584 :    
3585 : golsen 1.11 my $comment = substr( $s, $ind, $ind2-$ind-1 );
3586 : golsen 1.1 if ( $comment =~ m/\S/ ) { push @clist, $comment }
3587 :    
3588 : golsen 1.11 $ind = $ind2;
3589 : golsen 1.1
3590 :     # Skip white space
3591 :     if ( substr( $s, $ind ) =~ /^(\s+)/ ) { $ind += length( $1 ) }
3592 :     }
3593 :    
3594 :     ( $ind, @clist ? \@clist : undef )
3595 :     }
3596 :    
3597 :    
3598 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3599 :     # treeParseError( $message )
3600 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3601 :     sub treeParseError { die "Error: parse_newick_subtree: " . $_[0] . "\n" }
3602 :    
3603 :    
3604 :     #===============================================================================
3605 :     # Make a printer plot of a tree:
3606 :     #
3607 : golsen 1.13 # printer_plot_newick( $node, $file, $width, $min_dx, $dy )
3608 :     # printer_plot_newick( $node, $file, \%options )
3609 :     #
3610 :     # $node # newick tree root node
3611 :     # $file # undef = \*STDOUT, \*FH, or a file name.
3612 :     # $width # the approximate characters for the tree without labels (D = 68)
3613 :     # $min_dx # the minimum horizontal branch length (D = 2)
3614 :     # $dy # the vertical space per taxon (D = 1, most compressed)
3615 :     #
3616 :     # Options:
3617 :     #
3618 :     # dy => nat_number # the vertical space per taxon
3619 :     # chars => key # line drawing character set:
3620 :     # # html_unicode
3621 :     # # text (default)
3622 :     # min_dx => whole_number # the minimum horizontal branch length
3623 :     # width => whole_number # approximate tree width without labels
3624 : golsen 1.1 #
3625 : golsen 1.19 #===============================================================================
3626 : golsen 1.13 sub printer_plot_newick
3627 :     {
3628 :     my ( $node, $file, @opts ) = @_;
3629 : golsen 1.1
3630 : overbeek 1.7 my ( $fh, $close ) = open_output( $file );
3631 :     $fh or return;
3632 : golsen 1.1
3633 : golsen 1.13 my $html = $opts[0] && ref($opts[0]) eq 'HASH'
3634 :     && $opts[0]->{ chars }
3635 :     && $opts[0]->{ chars } =~ /html/;
3636 :     print $fh '<PRE>' if $html;
3637 :     print $fh join( "\n", text_plot_newick( $node, @opts ) ), "\n";
3638 :     print $fh "</PRE>\n" if $html;
3639 :    
3640 : golsen 1.1 if ( $close ) { close $fh }
3641 :     }
3642 :    
3643 :    
3644 : golsen 1.14 #===============================================================================
3645 :     # Character sets for printer plot trees:
3646 : golsen 1.13 #-------------------------------------------------------------------------------
3647 : golsen 1.14
3648 : golsen 1.13 my %char_set =
3649 :     ( text1 => { space => ' ',
3650 :     horiz => '-',
3651 :     vert => '|',
3652 :     el_d_r => '/',
3653 :     el_u_r => '\\',
3654 :     el_d_l => '\\',
3655 :     el_u_l => '/',
3656 :     tee_l => '+',
3657 :     tee_r => '+',
3658 :     tee_u => '+',
3659 :     tee_d => '+',
3660 :     half_l => '-',
3661 :     half_r => '-',
3662 :     half_u => '|',
3663 :     half_d => '|',
3664 :     cross => '+',
3665 :     },
3666 :     text2 => { space => ' ',
3667 :     horiz => '-',
3668 :     vert => '|',
3669 :     el_d_r => '+',
3670 :     el_u_r => '+',
3671 :     el_d_l => '+',
3672 :     el_u_l => '+',
3673 :     tee_l => '+',
3674 :     tee_r => '+',
3675 :     tee_u => '+',
3676 :     tee_d => '+',
3677 :     half_l => '-',
3678 :     half_r => '-',
3679 :     half_u => '|',
3680 :     half_d => '|',
3681 :     cross => '+',
3682 :     },
3683 : golsen 1.14 html_box => { space => '&nbsp;',
3684 :     horiz => '&#9472;',
3685 :     vert => '&#9474;',
3686 :     el_d_r => '&#9484;',
3687 :     el_u_r => '&#9492;',
3688 :     el_d_l => '&#9488;',
3689 :     el_u_l => '&#9496;',
3690 :     tee_l => '&#9508;',
3691 :     tee_r => '&#9500;',
3692 :     tee_u => '&#9524;',
3693 :     tee_d => '&#9516;',
3694 :     half_l => '&#9588;',
3695 :     half_r => '&#9590;',
3696 :     half_u => '&#9589;',
3697 :     half_d => '&#9591;',
3698 :     cross => '&#9532;',
3699 : golsen 1.13 },
3700 : golsen 1.14 utf8_box => { space => ' ',
3701 : golsen 1.13 horiz => chr(226) . chr(148) . chr(128),
3702 :     vert => chr(226) . chr(148) . chr(130),
3703 :     el_d_r => chr(226) . chr(148) . chr(140),
3704 :     el_u_r => chr(226) . chr(148) . chr(148),
3705 :     el_d_l => chr(226) . chr(148) . chr(144),
3706 :     el_u_l => chr(226) . chr(148) . chr(152),
3707 :     tee_l => chr(226) . chr(148) . chr(164),
3708 :     tee_r => chr(226) . chr(148) . chr(156),
3709 :     tee_u => chr(226) . chr(148) . chr(180),
3710 :     tee_d => chr(226) . chr(148) . chr(172),
3711 :     half_l => chr(226) . chr(149) . chr(180),
3712 :     half_r => chr(226) . chr(149) . chr(182),
3713 :     half_u => chr(226) . chr(149) . chr(181),
3714 :     half_d => chr(226) . chr(149) . chr(183),
3715 :     cross => chr(226) . chr(148) . chr(188),
3716 :     },
3717 :     );
3718 :    
3719 : golsen 1.14 %{ $char_set{ html1 } } = %{ $char_set{ text1 } };
3720 :     $char_set{ html1 }->{ space } = '&nbsp;';
3721 :    
3722 :     %{ $char_set{ html2 } } = %{ $char_set{ text2 } };
3723 :     $char_set{ html2 }->{ space } = '&nbsp;';
3724 :    
3725 : golsen 1.13 # Define some synonyms
3726 : golsen 1.14
3727 :     $char_set{ html } = $char_set{ html_box };
3728 :     $char_set{ line } = $char_set{ utf8_box };
3729 :     $char_set{ symb } = $char_set{ utf8_box };
3730 : golsen 1.13 $char_set{ text } = $char_set{ text1 };
3731 : golsen 1.14 $char_set{ utf8 } = $char_set{ utf8_box };
3732 : golsen 1.13
3733 : golsen 1.14 # Define tree formats and synonyms
3734 :    
3735 :     my %tree_format =
3736 :     ( text => 'text',
3737 :     tree_tab_lbl => 'tree_tab_lbl',
3738 :     tree_lbl => 'tree_lbl',
3739 :     chrlist_lbl => 'chrlist_lbl',
3740 :     raw => 'chrlist_lbl',
3741 :     );
3742 : golsen 1.13
3743 : golsen 1.1 #===============================================================================
3744 :     # Make a text plot of a tree:
3745 :     #
3746 : golsen 1.13 # @lines = text_plot_newick( $node, $width, $min_dx, $dy )
3747 :     # @lines = text_plot_newick( $node, \%options )
3748 :     #
3749 :     # $node # newick tree root node
3750 :     # $width # the approximate characters for the tree without labels (D = 68)
3751 :     # $min_dx # the minimum horizontal branch length (D = 2)
3752 :     # $dy # the vertical space per taxon (D = 1, most compressed)
3753 :     #
3754 :     # Options:
3755 :     #
3756 : golsen 1.14 # chars => keyword # the output character set for the tree
3757 : golsen 1.13 # dy => nat_number # the vertical space per taxon
3758 : golsen 1.14 # format => keyword # output format of each line
3759 : golsen 1.13 # min_dx => whole_number # the minimum horizontal branch length
3760 :     # width => whole_number # approximate tree width without labels
3761 : golsen 1.1 #
3762 : golsen 1.14 # Character sets:
3763 :     #
3764 :     # html # synonym of html1
3765 :     # html_box # html encoding of unicode box drawing characters
3766 :     # html1 # text1 with nonbreaking spaces
3767 :     # html2 # text2 with nonbreaking spaces
3768 :     # line # synonym of utf8_box
3769 :     # raw # pass out the internal representation
3770 :     # symb # synonym of utf8_box
3771 :     # text # synonym of text1 (Default)
3772 :     # text1 # ascii characters: - + | / \ and space
3773 :     # text2 # ascii characters: - + | + + and space
3774 :     # utf8 # synonym of utf8_box
3775 :     # utf8_box # utf8 encoding of unicode box drawing characters
3776 :     #
3777 :     # Formats for row lines:
3778 :     #
3779 :     # text # $textstring # Default
3780 :     # tree_tab_lbl # $treestr \t $labelstr
3781 :     # tree_lbl # [ $treestr, $labelstr ]
3782 :     # chrlist_lbl # [ \@treechar, $labelstr ] # Forced with raw chars
3783 :     # raw # synonym of chrlist_lbl
3784 :     #
3785 : golsen 1.1 #===============================================================================
3786 : golsen 1.13 sub text_plot_newick
3787 :     {
3788 :     my $node = shift @_;
3789 : golsen 1.1 array_ref( $node ) || die "Bad node passed to text_plot_newick\n";
3790 : golsen 1.13
3791 : golsen 1.14 my ( $opts, $width, $min_dx, $dy, $chars, $fmt );
3792 : golsen 1.13 if ( $_[0] && ref $_[0] eq 'HASH' )
3793 :     {
3794 :     $opts = shift;
3795 :     }
3796 :     else
3797 :     {
3798 : golsen 1.14 ( $width, $min_dx, $dy ) = @_;
3799 : golsen 1.13 $opts = {};
3800 :     }
3801 :    
3802 : golsen 1.14 $chars = $opts->{ chars } || '';
3803 :     my $charH;
3804 :     $charH = $char_set{ $chars } || $char_set{ 'text1' } if ( $chars ne 'raw' );
3805 :     my $is_box = $charH eq $char_set{ html_box }
3806 :     || $charH eq $char_set{ utf8_box }
3807 :     || $chars eq 'raw';
3808 :    
3809 :     $fmt = ( $chars eq 'raw' ) ? 'chrlist_lbl' : $opts->{ format };
3810 :     $fmt = $tree_format{ $fmt || '' } || 'text';
3811 :    
3812 :     $dy ||= $opts->{ dy } || 1;
3813 :     $width ||= $opts->{ width } || 68;
3814 :     $min_dx = $opts->{ min_dx } if ( ! defined $min_dx || $min_dx < 0 );
3815 :     $min_dx = $is_box ? 1 : 2 if ( ! defined $min_dx || $min_dx < 0 );
3816 :    
3817 :     # Layout the tree:
3818 : golsen 1.1
3819 :     $min_dx = int( $min_dx );
3820 :     $dy = int( $dy );
3821 : golsen 1.5 my $x_scale = $width / ( newick_max_X( $node ) || 1 ); # Div by zero caught by RAE
3822 : golsen 1.1
3823 :     my $hash = {};
3824 :     layout_printer_plot( $node, $hash, 0, -0.5 * $dy, $x_scale, $min_dx, $dy );
3825 :    
3826 : golsen 1.14 # Generate the lines of the tree-one by-one:
3827 : golsen 1.1
3828 :     my ( $y1, $y2 ) = @{ $hash->{ $node } };
3829 : golsen 1.14 my @lines;
3830 :     foreach ( ( $y1 .. $y2 ) )
3831 : golsen 1.13 {
3832 : golsen 1.16 my $line = text_tree_row( $node, $hash, $_, [], 'tee_l', $dy >= 2 );
3833 : golsen 1.14 my $lbl = '';
3834 :     if ( @$line )
3835 : golsen 1.13 {
3836 : golsen 1.14 if ( $line->[-1] eq '' ) { pop @$line; $lbl = pop @$line }
3837 :     # Translate tree characters
3838 :     @$line = map { $charH->{ $_ } } @$line if $chars ne 'raw';
3839 : golsen 1.13 }
3840 :    
3841 : golsen 1.14 # Convert to requested output format:
3842 :    
3843 :     push @lines, $fmt eq 'text' ? join( '', @$line, ( $lbl ? " $lbl" : () ) )
3844 :     : $fmt eq 'text_tab_lbl' ? join( '', @$line, "\t", $lbl )
3845 :     : $fmt eq 'tree_lbl' ? [ join( '', @$line ), $lbl ]
3846 :     : $fmt eq 'chrlist_lbl' ? [ $line, $lbl ]
3847 :     : ();
3848 :     }
3849 :    
3850 :     # if ( $cells )
3851 :     # {
3852 :     # my $nmax = 0;
3853 :     # foreach ( @lines ) { $nmax = @$_ if @$_ > $nmax }
3854 :     # foreach ( @lines )
3855 :     # {
3856 :     # @$_ = map { "<TD>$_</TD>" } @$_;
3857 :     # my $span = $nmax - @$_ + 1;
3858 :     # $_->[-1] =~ s/^<TD>/<TD NoWrap ColSpan=$span>/;
3859 :     # }
3860 :     # }
3861 :     # elsif ( $tables )
3862 :     # {
3863 :     # my $nmax = 0;
3864 :     # foreach ( @lines ) { $nmax = @$_ if @$_ > $nmax }
3865 :     # foreach ( @lines )
3866 :     # {
3867 :     # @$_ = map { "<TD>$_</TD>" } @$_;
3868 :     # my $span = $nmax - @$_ + 1;
3869 :     # $_->[-1] =~ s/^<TD>/<TD NoWrap ColSpan=$span>/;
3870 :     # }
3871 :     # }
3872 :    
3873 :     wantarray ? @lines : \@lines;
3874 : golsen 1.1 }
3875 :    
3876 : golsen 1.14
3877 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3878 : golsen 1.14 # ( $xmax, $ymax, $root_y ) = layout_printer_plot( $node, $hash, $x0, $y0, $x_scale, $min_dx, $dy, $yrnd )
3879 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3880 : golsen 1.13 sub layout_printer_plot
3881 :     {
3882 : golsen 1.14 my ( $node, $hash, $x0, $y0, $x_scale, $min_dx, $dy, $yrnd ) = @_;
3883 : golsen 1.1 array_ref( $node ) || die "Bad node ref passed to layout_printer_plot\n";
3884 :     hash_ref( $hash ) || die "Bad hash ref passed to layout_printer_plot\n";
3885 :    
3886 :     my $dx = newick_x( $node );
3887 :     if ( defined( $dx ) ) {
3888 :     $dx *= $x_scale;
3889 : golsen 1.14 $dx = $min_dx if $dx < $min_dx;
3890 : golsen 1.1 }
3891 :     else {
3892 :     $dx = ( $x0 > 0 ) ? $min_dx : 0;
3893 :     }
3894 :     $dx = int( $dx + 0.4999 );
3895 :    
3896 :     my ( $x, $xmax, $y, $ymax, $y1, $y2, $yn1, $yn2 );
3897 :    
3898 :     $x = $x0 + $dx;
3899 :     $y1 = int( $y0 + 0.5 * $dy + 0.4999 );
3900 :     my @dl = newick_desc_list( $node );
3901 :    
3902 :     if ( ! @dl ) { # A tip
3903 :     $xmax = $x;
3904 :     $y = $yn1 = $yn2 = $y2 = $y1;
3905 :     $ymax = $y + 0.5 * $dy;
3906 :     }
3907 :    
3908 :     else { # A subtree
3909 :     $xmax = -1;
3910 :     my $xmaxi;
3911 :     my $yi;
3912 :     my @ylist = ();
3913 :     $ymax = $y0;
3914 :    
3915 :     foreach ( @dl ) {
3916 : golsen 1.14 ( $xmaxi, $ymax, $yi ) = layout_printer_plot( $_, $hash, $x, $ymax, $x_scale, $min_dx, $dy,
3917 :     ( 2*@ylist < @dl ? 0.5001 : 0.4999 )
3918 :     );
3919 : golsen 1.1 push @ylist, $yi;
3920 :     if ( $xmaxi > $xmax ) { $xmax = $xmaxi }
3921 :     }
3922 :    
3923 :     # Use of y-list is overkill for saving first and last values,
3924 :     # but eases implimentation of alternative y-value calculations.
3925 :    
3926 :     $yn1 = $ylist[ 0];
3927 :     $yn2 = $ylist[-1];
3928 : golsen 1.16 $y = int( 0.5 * ( $yn1 + $yn2 ) + ( $yrnd || 0.4999 ) );
3929 :    
3930 :     # Handle special case of internal node label. Put it between subtrees.
3931 :    
3932 :     if ( ( $dy >= 2 ) && newick_lbl( $node ) && ( @dl > 1 ) ) {
3933 :     # Find the descendents $i1 and $i2 to put the branch between
3934 :     my $i2 = 1;
3935 :     while ( ( $i2+1 < @ylist ) && ( $ylist[$i2] < $y ) ) { $i2++ }
3936 :     my $i1 = $i2 - 1;
3937 :     # Get bottom of subtree1 and top of subtree2:
3938 :     my $ymax1 = $hash->{ $dl[ $i1 ] }->[ 1 ];
3939 :     my $ymin2 = $hash->{ $dl[ $i2 ] }->[ 0 ];
3940 :     # Midway between bottom of subtree1 and top of subtree2, with
3941 :     # preferred rounding direction
3942 :     $y = int( 0.5 * ( $ymax1 + $ymin2 ) + ( $yrnd || 0.4999 ) );
3943 :     }
3944 : golsen 1.1 }
3945 :    
3946 :     $y2 = int( $ymax - 0.5 * $dy + 0.4999 );
3947 :    
3948 :     $hash->{ $node } = [ $y1, $y2, $x0, $x, $y, $yn1, $yn2 ];
3949 :     ( $xmax, $ymax, $y );
3950 :     }
3951 :    
3952 :    
3953 : golsen 1.14 # What symbol do we get if we add a leftward line to some other symbol?
3954 :    
3955 :     my %with_left_line = ( space => 'half_l',
3956 :     horiz => 'horiz',
3957 :     vert => 'tee_l',
3958 :     el_d_r => 'tee_d',
3959 :     el_u_r => 'tee_u',
3960 :     el_d_l => 'el_d_l',
3961 :     el_u_l => 'el_u_l',
3962 :     tee_l => 'tee_l',
3963 :     tee_r => 'cross',
3964 :     tee_u => 'tee_u',
3965 :     tee_d => 'tee_d',
3966 :     half_l => 'half_l',
3967 :     half_r => 'horiz',
3968 :     half_u => 'el_u_l',
3969 :     half_d => 'el_d_l',
3970 :     cross => 'cross',
3971 :     );
3972 :    
3973 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3974 : golsen 1.14 # Produce a description of one line of a printer plot tree.
3975 :     #
3976 : golsen 1.16 # \@line = text_tree_row( $node, $hash, $row, \@line, $symb, $ilbl )
3977 : golsen 1.14 #
3978 :     # \@line is the character descriptions accumulated so far, one per array
3979 :     # element, except for a label, which can be any number of characters.
3980 :     # Labels are followed by an empty string, so if $line->[-1] eq '',
3981 :     # then $line->[-2] is a label. The calling program translates the
3982 :     # symbol names to output characters.
3983 :     #
3984 :     # \@node is a newick tree node
3985 :     # \%hash contains tree layout information
3986 :     # $row is the row number (y value) that we are building
3987 :     # $symb is the plot symbol proposed for the current x and y position
3988 : golsen 1.16 # $ilbl is true if internal node labels are allowed
3989 : golsen 1.14 #
3990 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3991 : golsen 1.13 sub text_tree_row
3992 :     {
3993 : golsen 1.16 my ( $node, $hash, $row, $line, $symb, $ilbl ) = @_;
3994 : golsen 1.1
3995 :     my ( $y1, $y2, $x0, $x, $y, $yn1, $yn2 ) = @{ $hash->{ $node } };
3996 :     if ( $row < $y1 || $row > $y2 ) { return $line }
3997 :    
3998 : golsen 1.14 if ( @$line < $x0 ) { push @$line, ('space') x ( $x0 - @$line ) }
3999 : golsen 1.1
4000 :     if ( $row == $y ) {
4001 : golsen 1.14 while ( @$line > $x0 ) { pop @$line } # Actually 0-1 times
4002 :     push @$line, $symb,
4003 :     ( ( $x > $x0 ) ? ('horiz') x ($x - $x0) : () );
4004 : golsen 1.1 }
4005 :    
4006 :     elsif ( $row > $yn1 && $row < $yn2 ) {
4007 : golsen 1.14 if ( @$line < $x ) { push @$line, ('space') x ( $x - @$line ), 'vert' }
4008 :     else { $line->[$x] = 'vert' }
4009 : golsen 1.1 }
4010 :    
4011 :     my @dl = newick_desc_list( $node );
4012 :    
4013 : golsen 1.14 if ( @dl < 1 ) { push @$line, ( newick_lbl( $node ) || '' ), '' }
4014 : golsen 1.1
4015 :     else {
4016 : golsen 1.14 my @list = map { [ $_, 'tee_r' ] } @dl; # Line to the right
4017 :     if ( @list > 1 ) { # Fix top and bottom sympbols
4018 :     $list[ 0]->[1] = 'el_d_r';
4019 :     $list[-1]->[1] = 'el_u_r';
4020 :     }
4021 :     elsif ( @list ) { # Only one descendent
4022 :     $list[ 0]->[1] = 'half_r';
4023 :     }
4024 : golsen 1.1 foreach ( @list ) {
4025 :     my ( $n, $s ) = @$_;
4026 :     if ( $row >= $hash->{ $n }->[0] && $row <= $hash->{ $n }->[1] ) {
4027 : golsen 1.16 $line = text_tree_row( $n, $hash, $row, $line, $s, $ilbl );
4028 : golsen 1.1 }
4029 : golsen 1.16 }
4030 : golsen 1.1
4031 : golsen 1.14 if ( $row == $y ) {
4032 :     $line->[$x] = ( $line->[$x] eq 'horiz' ) ? 'tee_l'
4033 :     : $with_left_line{ $line->[$x] };
4034 : golsen 1.16 push( @$line, newick_lbl( $node ), '' ) if $ilbl && newick_lbl( $node );
4035 : golsen 1.13 }
4036 : golsen 1.1 }
4037 :    
4038 :     return $line;
4039 :     }
4040 :    
4041 :    
4042 : golsen 1.13 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4043 :     # Debug routine
4044 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4045 :     sub dump_tree {
4046 :     my ( $node, $prefix ) = @_;
4047 :     defined( $prefix ) or $prefix = "";
4048 :     print STDERR $prefix, join(", ", @$node), "\n";
4049 :     my @dl = $node->[0] ? @{$node->[0]} : ();
4050 :     foreach ( @dl ) { dump_tree( $_, $prefix . " " ) }
4051 :     $prefix or print STDERR "\n";
4052 :     }
4053 :    
4054 :    
4055 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4056 :     # Debug routine
4057 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4058 :     sub dump_tree_hash {
4059 :     my ( $node, $hash, $prefix ) = @_;
4060 :     defined( $prefix ) or print STDERR "node; [ y1, y2, x0, x, y, yn1, yn2 ]\n" and $prefix = "";
4061 :     print STDERR $prefix, join(", ", @$node), "; ", join(", ", @{ $hash->{ $node } } ), "\n";
4062 :     my @dl = $node->[0] ? @{$node->[0]} : ();
4063 :     foreach (@dl) { dump_tree_hash( $_, $hash, $prefix . " " ) }
4064 :     }
4065 :    
4066 :    
4067 : overbeek 1.7 #===============================================================================
4068 :     # Open an input file stream:
4069 :     #
4070 :     # ( $handle, undef ) = open_input( ); # \*STDIN
4071 :     # ( $handle, undef ) = open_input( \*FH );
4072 :     # ( $handle, 1 ) = open_input( $file ); # need to close $handle
4073 :     #
4074 :     #===============================================================================
4075 :     sub open_input
4076 :     {
4077 :     my $file = shift;
4078 :     my $fh;
4079 : golsen 1.21 if ( ! defined $file || $file eq '' ) { return ( \*STDIN ) }
4080 :     elsif ( ref( $file ) eq 'GLOB' ) { return ( $file ) }
4081 :     elsif ( open( $fh, "<$file" ) ) { return ( $fh, 1 ) } # Need to close
4082 : overbeek 1.7
4083 :     print STDERR "gjonewick::open_input could not open '$file' for reading\n";
4084 :     return undef;
4085 :     }
4086 :    
4087 :    
4088 :     #===============================================================================
4089 :     # Open an output file stream:
4090 :     #
4091 :     # ( $handle, undef ) = open_output( ); # \*STDOUT
4092 :     # ( $handle, undef ) = open_output( \*FH );
4093 :     # ( $handle, 1 ) = open_output( $file ); # need to close $handle
4094 :     #
4095 :     #===============================================================================
4096 :     sub open_output
4097 :     {
4098 :     my $file = shift;
4099 :     my $fh;
4100 : golsen 1.21 if ( ! defined $file || $file eq '' ) { return ( \*STDOUT ) }
4101 :     elsif ( ref( $file ) eq 'GLOB' ) { return ( $file ) }
4102 :     elsif ( ( open $fh, ">$file" ) ) { return ( $fh, 1 ) } # Need to close
4103 : overbeek 1.7
4104 :     print STDERR "gjonewick::open_output could not open '$file' for writing\n";
4105 :     return undef;
4106 :     }
4107 :    
4108 : golsen 1.9
4109 :     #===============================================================================
4110 :     # Some subroutines copied from gjolists
4111 :     #===============================================================================
4112 :     # Return the common prefix of two lists:
4113 :     #
4114 :     # @common = common_prefix( \@list1, \@list2 )
4115 :     #-----------------------------------------------------------------------------
4116 :     sub common_prefix
4117 :     {
4118 :     my ($l1, $l2) = @_;
4119 :     ref($l1) eq "ARRAY" || die "common_prefix: arg 1 is not an array ref\n";
4120 :     ref($l2) eq "ARRAY" || die "common_prefix: arg 2 is not an array ref\n";
4121 :    
4122 :     my $i = 0;
4123 :     my $l1_i;
4124 :     while ( defined( $l1_i = $l1->[$i] ) && $l1_i eq $l2->[$i] ) { $i++ }
4125 :    
4126 :     return @$l1[ 0 .. ($i-1) ]; # perl handles negative range
4127 :     }
4128 :    
4129 :    
4130 :     #-----------------------------------------------------------------------------
4131 :     # Return the unique suffixes of each of two lists:
4132 :     #
4133 :     # ( \@suffix1, \@suffix2 ) = unique_suffixes( \@list1, \@list2 )
4134 :     #-----------------------------------------------------------------------------
4135 :     sub unique_suffixes
4136 :     {
4137 :     my ($l1, $l2) = @_;
4138 :     ref($l1) eq "ARRAY" || die "common_prefix: arg 1 is not an array ref\n";
4139 :     ref($l2) eq "ARRAY" || die "common_prefix: arg 2 is not an array ref\n";
4140 :    
4141 :     my $i = 0;
4142 :     my @l1 = @$l1;
4143 :     my @l2 = @$l2;
4144 :     my $l1_i;
4145 :     while ( defined( $l1_i = $l1[$i] ) && $l1_i eq $l2[$i] ) { $i++ }
4146 :    
4147 :     splice @l1, 0, $i;
4148 :     splice @l2, 0, $i;
4149 :     return ( \@l1, \@l2 );
4150 :     }
4151 :    
4152 :    
4153 :     #-------------------------------------------------------------------------------
4154 :     # List of values duplicated in a list (stable in order by second occurance):
4155 :     #
4156 :     # @dups = duplicates( @list )
4157 :     #-------------------------------------------------------------------------------
4158 :     sub duplicates
4159 :     {
4160 :     my %cnt = ();
4161 :     grep { ++$cnt{$_} == 2 } @_;
4162 :     }
4163 :    
4164 :    
4165 :     #-------------------------------------------------------------------------------
4166 :     # Randomize the order of a list:
4167 :     #
4168 :     # @random = random_order( @list )
4169 :     #-------------------------------------------------------------------------------
4170 :     sub random_order
4171 :     {
4172 :     my ( $i, $j );
4173 :     for ( $i = @_ - 1; $i > 0; $i-- )
4174 :     {
4175 :     $j = int( ($i+1) * rand() );
4176 :     ( $_[$i], $_[$j] ) = ( $_[$j], $_[$i] ); # Interchange i and j
4177 :     }
4178 :    
4179 :     @_;
4180 :     }
4181 :    
4182 :    
4183 :     #-----------------------------------------------------------------------------
4184 :     # Intersection of two or more sets:
4185 :     #
4186 :     # @intersection = intersection( \@set1, \@set2, ... )
4187 :     #-----------------------------------------------------------------------------
4188 :     sub intersection
4189 :     {
4190 :     my $set = shift;
4191 :     my @intersection = @$set;
4192 :    
4193 :     foreach $set ( @_ )
4194 :     {
4195 :     my %set = map { $_ => 1 } @$set;
4196 :     @intersection = grep { exists $set{ $_ } } @intersection;
4197 :     }
4198 :    
4199 :     @intersection;
4200 :     }
4201 :    
4202 :    
4203 :     #-----------------------------------------------------------------------------
4204 :     # Elements in set 1, but not set 2:
4205 :     #
4206 :     # @difference = set_difference( \@set1, \@set2 )
4207 :     #-----------------------------------------------------------------------------
4208 :     sub set_difference
4209 :     {
4210 :     my ($set1, $set2) = @_;
4211 :     my %set2 = map { $_ => 1 } @$set2;
4212 :     grep { ! ( exists $set2{$_} ) } @$set1;
4213 :     }
4214 :    
4215 :    
4216 : golsen 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3