[Phylobase-devl] [Fwd: Re: possible phylobase bug]

[Ben Bolker]

-------- Original Message --------
Subject: Re: possible phylobase bug
Date: Mon, 15 Dec 2008 13:13:19 -0800
From: Jim Regetz <regetz at nceas.ucsb.edu>
To: bolker at ufl.edu
References: <493D914F.8070709 at nceas.ucsb.edu> <49414945.6050200 at zoo.ufl.edu>

Ben Bolker wrote:
>   This is done.
>   We're planning a little "virtual hackathon" Dec 18-19, so if
> you have any other thoughts about design, usability, documentation,
> etc etc etc I'd love to hear them.
>

Thanks, Ben!

I'm actually collaborating with some folks on developing and
implementing a couple of new "community phylodiversity" metrics and
related ecological-relevant tree computations, all of which will be made
available in an R package. When we started this in the spring, I was
programming against the 'ape' tree implementation, but recently made the
switch to using phylobase as the base platform -- I gather that is the
raison d'etre of phylobase, analogous to what 'sp' provides for R
spatial? I've enjoyed working with it so far.

I've written a couple of new utility functions (code below) that make my
life easier when working with phylo4[d] objects. Not sure if you have
any interest in making either of these part of phylobase, but you're
welcome to do so. And for what it's worth, I'd be happy to help in any
way with the coding/development, now or in the future!

Quick synopsis: The 'nodes' function lets me easily extract a vector of
node numbers (with node names as vector names), for particular types of
nodes. The 'ancestralEdgeLength' (horrible name) function is basically
'edgeLength' with extra bells and whistles. Both include pirated pieces
of existing phylobase code, as you can see :) I haven't thoroughly
tested them or anything, and you might want to modify parameter names
and such for consistency, but they've been useful so far for me.

Cheers,
Jim

-----------

# Return vector of phylo4 nodes, with names if they exist; use
# parameters 'which' and 'include.root' to restrict results
nodes <- function(tree, which=c("all", "internal", "tip"),
  include.root=TRUE) {

  which <- match.arg(which)

  x <- as(tree, "data.frame")
  node.type <- x$node.type
  node <- x$node
  names(node) <- x$label

  tip.nodes <- node[node.type=="tip"]
  if (which=="tip") {
    return(tip.nodes)
  }

  int.nodes <- node[node.type=="internal"]
  if (include.root) {
    root <- node[node.type=="root"]
    int.nodes <- c(root, int.nodes)
  }

  if (which=="internal") {
    return(int.nodes)
  } else if (which=="all") {
    return(c(int.nodes, tip.nodes))
  }

  stop("internal error, contact package maintainer")

}

# Return length of the branch immediately ancestral to each specified
# node. Use 'node' parameter to restrict results to particular node(s).
ancestralEdgeLength <- function(tree, node=NULL) {
  E <- edges(tree)
  ancestor <- E[, 1]
  if (is.null(node)) {
    node <- E[,2]
  } else if (!all(node %in% E)) {
    stop("one or more nodes not found in tree")
  }
  idx <- match(node, E[, 2])
  ancestor <- E[idx, 1]
  branch.length <- edgeLength(tree)[idx]
  if (is.null(edgeLength(tree))) {
    branch.length <- rep(NA, length(node))
  }
  names(branch.length) <- node
  return(branch.length)
}

-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc