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

Diff of /FigKernelPackages/gjonewicklib.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.11, Sat Mar 29 05:18:19 2008 UTC revision 1.12, Mon Jun 22 19:14:50 2009 UTC
# Line 221  Line 221 
221  #  $newtree = reroot_newick_to_node( $tree, @node )  #  $newtree = reroot_newick_to_node( $tree, @node )
222  #  $newtree = reroot_newick_to_node_ref( $tree, $noderef )  #  $newtree = reroot_newick_to_node_ref( $tree, $noderef )
223  #  $newtree = reroot_newick_between_nodes( $tree, $node1, $node2, $fraction )  #  $newtree = reroot_newick_between_nodes( $tree, $node1, $node2, $fraction )
224    #  $newtree = reroot_newick_to_midpoint( $tree )           # unweighted
225    #  $newtree = reroot_newick_to_midpoint_w( $tree )         # weight by tips
226  #  $newtree = reroot_newick_to_approx_midpoint( $tree )    # unweighted  #  $newtree = reroot_newick_to_approx_midpoint( $tree )    # unweighted
227  #  $newtree = reroot_newick_to_approx_midpoint_w( $tree )  # weight by tips  #  $newtree = reroot_newick_to_approx_midpoint_w( $tree )  # weight by tips
228  #  $newtree = uproot_tip_rooted_newick( $tree )  #  $newtree = uproot_tip_rooted_newick( $tree )
229  #  $newtree = uproot_newick( $tree )  #  $newtree = uproot_newick( $tree )
230  #  #
231  #  $newtree = prune_from_newick( $tree, $tip )  #  $newtree = prune_from_newick( $tree, $tip )
232    #  $newtree = rooted_newick_subtree( $tree,  @tips )
233    #  $newtree = rooted_newick_subtree( $tree, \@tips )
234  #  $newtree = newick_subtree( $tree,  @tips )  #  $newtree = newick_subtree( $tree,  @tips )
235  #  $newtree = newick_subtree( $tree, \@tips )  #  $newtree = newick_subtree( $tree, \@tips )
236    #  $newtree = newick_covering_subtree( $tree,  @tips )
237    #  $newtree = newick_covering_subtree( $tree, \@tips )
238  #  #
239  #  $newtree = collapse_zero_length_branches( $tree )  #  $newtree = collapse_zero_length_branches( $tree )
240  #  #
# Line 355  Line 361 
361          reroot_newick_to_node          reroot_newick_to_node
362          reroot_newick_to_node_ref          reroot_newick_to_node_ref
363          reroot_newick_between_nodes          reroot_newick_between_nodes
364            reroot_newick_to_midpoint
365            reroot_newick_to_midpoint_w
366          reroot_newick_to_approx_midpoint          reroot_newick_to_approx_midpoint
367          reroot_newick_to_approx_midpoint_w          reroot_newick_to_approx_midpoint_w
368          uproot_tip_rooted_newick          uproot_tip_rooted_newick
369          uproot_newick          uproot_newick
370    
371          prune_from_newick          prune_from_newick
372            rooted_newick_subtree
373          newick_subtree          newick_subtree
374            newick_covering_subtree
375          collapse_zero_length_branches          collapse_zero_length_branches
376    
377          newick_insert_at_node          newick_insert_at_node
# Line 473  Line 483 
483  #  #
484  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
485    
486  sub newick_desc_ref { $_[0]->[0] }  # = ${$_[0]}[0]  sub newick_desc_ref { ref($_[0]) ? $_[0]->[0] : Carp::confess() }  # = ${$_[0]}[0]
487  sub newick_lbl      { ref($_[0]) ? $_[0]->[1] : Carp::confess() }  sub newick_lbl      { ref($_[0]) ? $_[0]->[1] : Carp::confess() }
488  sub newick_x        { $_[0]->[2] }  sub newick_x        { ref($_[0]) ? $_[0]->[2] : Carp::confess() }
489  sub newick_c1       { ref($_[0]) ? $_[0]->[3] : Carp::confess() }  sub newick_c1       { ref($_[0]) ? $_[0]->[3] : Carp::confess() }
490  sub newick_c2       { $_[0]->[4] }  sub newick_c2       { ref($_[0]) ? $_[0]->[4] : Carp::confess() }
491  sub newick_c3       { $_[0]->[5] }  sub newick_c3       { ref($_[0]) ? $_[0]->[5] : Carp::confess() }
492  sub newick_c4       { $_[0]->[6] }  sub newick_c4       { ref($_[0]) ? $_[0]->[6] : Carp::confess() }
493  sub newick_c5       { $_[0]->[7] }  sub newick_c5       { ref($_[0]) ? $_[0]->[7] : Carp::confess() }
494    
495  sub newick_desc_list {  sub newick_desc_list {
496      my $node = $_[0];      my $node = $_[0];
# Line 1165  Line 1175 
1175      my $imax = newick_n_desc( $node );      my $imax = newick_n_desc( $node );
1176      for ( my $i = 1; $i <= $imax; $i++ ) {      for ( my $i = 1; $i <= $imax; $i++ ) {
1177         @path = path_to_node_ref( newick_desc_i( $node, $i ), $noderef, ( @path0, $i ) );         @path = path_to_node_ref( newick_desc_i( $node, $i ), $noderef, ( @path0, $i ) );
1178         if ( @path ) { return @path }          return @path if @path;
1179      }      }
1180    
1181      ();  #  Not found      ();  #  Not found
# Line 1921  Line 1931 
1931      my @path1 = path_to_node( $tree, $node1 ) or return undef;      my @path1 = path_to_node( $tree, $node1 ) or return undef;
1932      my @path2 = path_to_node( $tree, $node2 ) or return undef;      my @path2 = path_to_node( $tree, $node2 ) or return undef;
1933    
1934        reroot_newick_between_nodes_by_path( $tree, \@path1, \@path2, $fraction )
1935    }
1936    
1937    
1938    #-------------------------------------------------------------------------------
1939    #  Reroot a newick tree along the path between 2 nodes:
1940    #
1941    #  $tree = reroot_newick_between_node_refs( $tree, $node1, $node2, $fraction )
1942    #-------------------------------------------------------------------------------
1943    sub reroot_newick_between_node_refs
1944    {
1945        my ( $tree, $node1, $node2, $fraction ) = @_;
1946        array_ref( $tree ) or return undef;
1947    
1948        #  Find the paths to the nodes:
1949    
1950        my @path1 = path_to_node_ref( $tree, $node1 ) or return undef;
1951        my @path2 = path_to_node_ref( $tree, $node2 ) or return undef;;
1952    
1953        reroot_newick_between_nodes_by_path( $tree, \@path1, \@path2, $fraction )
1954    }
1955    
1956    
1957    #-------------------------------------------------------------------------------
1958    #  Reroot a newick tree along the path between 2 nodes defined by paths:
1959    #
1960    #  $tree = reroot_newick_between_nodes_by_path( $tree, $path1, $path2, $fraction )
1961    #-------------------------------------------------------------------------------
1962    sub reroot_newick_between_nodes_by_path
1963    {
1964        my ( $tree, $path1, $path2, $fraction ) = @_;
1965        array_ref( $tree ) and array_ref( $path1 ) and  array_ref( $path2 )
1966           or return undef;
1967        $fraction >= 0 && $fraction <= 1 or return undef;
1968    
1969        my @path1 = @$path1;
1970        my @path2 = @$path2;
1971    
1972      #  Trim the common prefix, saving it:      #  Trim the common prefix, saving it:
1973    
1974      my @prefix = ();      my @prefix = ();
# Line 1984  Line 2032 
2032    
2033      my $dists1 = average_to_tips_1( $tree );      my $dists1 = average_to_tips_1( $tree );
2034    
2035      #  Compile average tip to node distances descending, returning midpoint node      #  Compile average tip to node distances descending, returning midpoint
2036        #  cadidates as a list of [ $node1, $node2, $fraction ]
2037    
2038      my $node = average_to_tips_2( $dists1, undef, undef );      my @mids = average_to_tips_2( $dists1, undef, undef );
2039    
2040      #  Reroot      #  Reroot to first midpoint candidate
2041    
2042      $node ? reroot_newick_to_node_ref( $tree, $node ) : $tree      return $tree if ! @mids;
2043        my ( $node1, $node2, $fraction ) = @{ $mids[0] };
2044        reroot_newick_to_node_ref( $tree, $fraction >= 0.5 ? $node2 : $node1 );
2045  }  }
2046    
2047    
2048    #-------------------------------------------------------------------------------
2049    #  Move root of tree to a midpoint.
2050    #
2051    #  $newtree = reroot_newick_to_midpoint( $tree )
2052    #-------------------------------------------------------------------------------
2053    sub reroot_newick_to_midpoint {
2054        my ( $tree ) = @_;
2055    
2056        #  Compile average tip to node distances assending
2057    
2058        my $dists1 = average_to_tips_1( $tree );
2059    
2060        #  Compile average tip to node distances descending, returning midpoint
2061        #  [ $node1, $node2, $fraction ]
2062    
2063        my @mids = average_to_tips_2( $dists1, undef, undef );
2064    
2065        @mids ? reroot_newick_between_node_refs( $tree, @{ $mids[0] } ) : $tree;
2066    }
2067    
2068    
2069    #-------------------------------------------------------------------------------
2070    #  Compile average tip to node distances assending
2071    #-------------------------------------------------------------------------------
2072  sub average_to_tips_1 {  sub average_to_tips_1 {
2073      my ( $node ) = @_;      my ( $node ) = @_;
2074    
# Line 2004  Line 2079 
2079          foreach ( @desc_dists ) { $x_below += $_->[0] }          foreach ( @desc_dists ) { $x_below += $_->[0] }
2080          $x_below /= @desc_dists;          $x_below /= @desc_dists;
2081      }      }
2082    
2083      my $x = newick_x( $node ) || 0;      my $x = newick_x( $node ) || 0;
2084      my $x_net = $x_below + $x;      my $x_net = $x_below + $x;
2085    
# Line 2011  Line 2087 
2087  }  }
2088    
2089    
2090    #-------------------------------------------------------------------------------
2091    #  Compile average tip to node distances descending, returning midpoint as
2092    #  [ $node1, $node2, $fraction_of_dist_between ]
2093    #-------------------------------------------------------------------------------
2094  sub average_to_tips_2 {  sub average_to_tips_2 {
2095      my ( $dists1, $x_above, $anc_node ) = @_;      my ( $dists1, $x_above, $anc_node ) = @_;
2096      my ( undef, $x, $x_below, $desc_list, $node ) = @$dists1;      my ( undef, $x, $x_below, $desc_list, $node ) = @$dists1;
2097    
2098      #  Are we done?  Root is in this node's branch, or "above"?      #  Are we done?  Root is in this node's branch, or "above"?
2099    
2100      # defined( $x_above ) and print STDERR "x_above = $x_above\n";      my @mids = ();
     # print STDERR "x       = $x\n";  
     # print STDERR "x_below = $x_below\n";  
     # print STDERR "n_desc  = ", scalar @$desc_list, "\n\n";  
   
2101      if ( defined( $x_above ) && ( ( $x_above + $x ) >= $x_below ) )      if ( defined( $x_above ) && ( ( $x_above + $x ) >= $x_below ) )
2102      {      {
2103          #  At this point the root can only be in this node's branch,          #  At this point the root can only be in this node's branch,
# Line 2033  Line 2109 
2109    
2110          if ( ( $x_below + $x ) >= $x_above )          if ( ( $x_below + $x ) >= $x_above )
2111          {          {
2112              return ( $x_above >= $x_below ) ? $anc_node : $node;              #  We will need to make a new node for the root, $fract of
2113          }              #  the way from $node to $anc_node:
2114          else              my $fract = ( $x > 0 ) ? 0.5 * ( ( $x_above - $x_below ) / $x + 1 )
2115          {                                     : 0.5;
2116              return undef;              push @mids, [ $node, $anc_node, $fract ];
2117          }          }
2118      }      }
2119    
2120      #  The root must be somewhere below this node:      #  The root might be somewhere below this node:
2121    
2122      my $n_1      =   @$desc_list - ( $anc_node ? 0 : 1 );      my $n_1      =   @$desc_list - ( $anc_node ? 0 : 1 );
2123      my $ttl_dist = ( @$desc_list * $x_below ) + ( defined( $x_above ) ? ( $x_above + $x ) : 0 );      my $ttl_dist = ( @$desc_list * $x_below ) + ( defined( $x_above ) ? ( $x_above + $x ) : 0 );
# Line 2051  Line 2127 
2127          #  If input tree is tip_rooted, $n-1 can be 0, so:          #  If input tree is tip_rooted, $n-1 can be 0, so:
2128    
2129          my $above2 = $n_1 ? ( ( $ttl_dist - $_->[0] ) / $n_1 ) : 0;          my $above2 = $n_1 ? ( ( $ttl_dist - $_->[0] ) / $n_1 ) : 0;
2130          my $root = average_to_tips_2( $_, $above2, $node );          push @mids, average_to_tips_2( $_, $above2, $node );
         if ( $root ) { return $root }  
2131      }      }
2132    
2133      #  Was not anywhere below this node (oh-oh):      return @mids;
   
     return undef;  
2134  }  }
2135    
2136    
# Line 2069  Line 2142 
2142  sub reroot_newick_to_approx_midpoint_w {  sub reroot_newick_to_approx_midpoint_w {
2143      my ( $tree ) = @_;      my ( $tree ) = @_;
2144    
2145        #  Compile average tip to node distances assending from tips
2146    
2147        my $dists1 = average_to_tips_1_w( $tree );
2148    
2149        #  Compile average tip to node distances descending, returning midpoints
2150    
2151        my @mids = average_to_tips_2_w( $dists1, undef, undef, undef );
2152    
2153        #  Reroot to first midpoint candidate
2154    
2155        return $tree if ! @mids;
2156        my ( $node1, $node2, $fraction ) = @{ $mids[0] };
2157        reroot_newick_to_node_ref( $tree, $fraction >= 0.5 ? $node2 : $node1 );
2158    }
2159    
2160    
2161    #-------------------------------------------------------------------------------
2162    #  Move root of tree to an approximate midpoint.  Weight by tips.
2163    #
2164    #  $newtree = reroot_newick_to_midpoint_w( $tree )
2165    #-------------------------------------------------------------------------------
2166    sub reroot_newick_to_midpoint_w {
2167        my ( $tree ) = @_;
2168    
2169      #  Compile average tip to node distances assending      #  Compile average tip to node distances assending
2170    
2171      my $dists1 = average_to_tips_1_w( $tree );      my $dists1 = average_to_tips_1_w( $tree );
2172    
2173      #  Compile average tip to node distances descending, returning midpoint node      #  Compile average tip to node distances descending, returning midpoint node
2174    
2175      my $node = average_to_tips_2_w( $dists1, undef, undef, undef );      my @mids = average_to_tips_2_w( $dists1, undef, undef, undef );
2176    
2177      #  Reroot      #  Reroot at first candidate midpoint
2178    
2179      $node ? reroot_newick_to_node_ref( $tree, $node ) : $tree      @mids ? reroot_newick_between_node_refs( $tree, @{ $mids[0] } ) : $tree;
2180  }  }
2181    
2182    
# Line 2100  Line 2197 
2197          }          }
2198          $x_below /= $n_below;          $x_below /= $n_below;
2199      }      }
2200    
2201      my $x = newick_x( $node ) || 0;      my $x = newick_x( $node ) || 0;
2202      my $x_net = $x_below + $x;      my $x_net = $x_below + $x;
2203    
# Line 2113  Line 2211 
2211    
2212      #  Are we done?  Root is in this node's branch, or "above"?      #  Are we done?  Root is in this node's branch, or "above"?
2213    
2214      # defined( $x_above ) and print STDERR "x_above = $x_above\n";      my @mids = ();
     # print STDERR "x       = $x\n";  
     # print STDERR "x_below = $x_below\n";  
     # print STDERR "n_desc  = ", scalar @$desc_list, "\n\n";  
   
2215      if ( defined( $x_above ) && ( ( $x_above + $x ) >= $x_below ) )      if ( defined( $x_above ) && ( ( $x_above + $x ) >= $x_below ) )
2216      {      {
2217          #  At this point the root can only be in this node's branch,          #  At this point the root can only be in this node's branch,
# Line 2125  Line 2219 
2219          #  would mean that the midpoint is actually down a different          #  would mean that the midpoint is actually down a different
2220          #  path from the root of the current tree).          #  path from the root of the current tree).
2221          #          #
2222          #  Is the root in the current branch?          #  Is their a root in the current branch?
2223    
2224          if ( ( $x_below + $x ) >= $x_above )          if ( ( $x_below + $x ) >= $x_above )
2225          {          {
2226              return ( $x_above >= $x_below ) ? $anc_node : $node;              #  We will need to make a new node for the root, $fract of
2227          }              #  the way from $node to $anc_node:
2228          else              my $fract = ( $x > 0 ) ? 0.5 * ( ( $x_above - $x_below ) / $x + 1 )
2229          {                                     : 0.5;
2230              return undef;              push @mids, [ $node, $anc_node, $fract ];
2231          }          }
2232      }      }
2233    
# Line 2153  Line 2247 
2247    
2248          my $x_above2 = $n_above2 ? ( ( $ttl_w_dist - $n_2 * $_->[0] ) / $n_above2 )          my $x_above2 = $n_above2 ? ( ( $ttl_w_dist - $n_2 * $_->[0] ) / $n_above2 )
2249                                   : 0;                                   : 0;
2250          my $root = average_to_tips_2_w( $_, $x_above2, $n_above2 || 1, $node );          push @mids, average_to_tips_2_w( $_, $x_above2, $n_above2 || 1, $node );
         if ( $root ) { return $root }  
2251      }      }
2252    
2253      #  Was not anywhere below this node (oh-oh):      return @mids;
   
     return undef;  
2254  }  }
2255    
2256    
# Line 2474  Line 2565 
2565    
2566    
2567  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
2568    #  Produce a potentially rooted subtree with the desired tips:
2569    #
2570    #     Except for (some) tip nodes, the tree produced is a copy.
2571    #     There is no check that requested tips exist.
2572    #
2573    #  $newtree = rooted_newick_subtree( $tree,  @tips )
2574    #  $newtree = rooted_newick_subtree( $tree, \@tips )
2575    #-------------------------------------------------------------------------------
2576    sub rooted_newick_subtree {
2577        my ( $tr, @tips ) = @_;
2578        if ( @tips == 1 && ref( $tips[0] ) eq "ARRAY" ) { @tips = @{ $tips[0] } }
2579    
2580        if ( @tips < 2 ) { return undef }
2581        my $keephash = { map { ( $_, 1 ) } @tips };
2582        my $tr2 = subtree1( $tr, $keephash );
2583        $tr2->[2] = undef if $tr2;                   # undef root branch length
2584        $tr2;
2585    }
2586    
2587    
2588    #-------------------------------------------------------------------------------
2589  #  Produce a subtree with the desired tips:  #  Produce a subtree with the desired tips:
2590  #  #
2591  #     Except for (some) tip nodes, the tree produced is a copy.  #     Except for (some) tip nodes, the tree produced is a copy.
# Line 2547  Line 2659 
2659  }  }
2660    
2661    
2662    #-------------------------------------------------------------------------------
2663    #  The smallest subtree of rooted tree that includes @tips:
2664    #
2665    #    $node = newick_covering_subtree( $tree,  @tips )
2666    #    $node = newick_covering_subtree( $tree, \@tips )
2667    #-------------------------------------------------------------------------------
2668    
2669    sub newick_covering_subtree {
2670        my $tree = shift;
2671        my %tips = map { $_ => 1 } ( ( ref( $_[0] ) eq 'ARRAY' ) ? @{ $_[0] } : @_ );
2672    
2673        #  Return smallest covering node, if any:
2674    
2675        ( newick_covering_subtree( $tree, \%tips ) )[ 0 ];
2676    }
2677    
2678    
2679    sub newick_covering_subtree_1 {
2680        my ( $node, $tips ) = @_;
2681        my $n_cover = 0;
2682        my @desc = newick_desc_list( $node );
2683        if ( @desc )
2684        {
2685            foreach ( @desc )
2686            {
2687                my ( $subtree, $n ) = newick_covering_subtree_1( $_, $tips );
2688                return ( $subtree, $n ) if $subtree;
2689                $n_cover += $n;
2690            }
2691        }
2692        elsif ( $tips->{ newick_lbl( $node ) } )
2693        {
2694            $n_cover++;
2695        }
2696    
2697        #  If all tips are covered, return node
2698    
2699        ( $n_cover == keys %$tips ) ? ( $node, $n_cover ) : ( undef, $n_cover );
2700    }
2701    
2702    
2703  #===============================================================================  #===============================================================================
2704  #  #
2705  #  Representative subtrees  #  Representative subtrees

Legend:
Removed from v.1.11  
changed lines
  Added in v.1.12

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3