Monday, March 7, 2011

Creating a star phylogeny

I just posted up a function to create a star phylogeny (that is, a completely unresolved multifurcating tree) on my R phylogenetics page (code here). Very simple, and it seems likely that this can already be done with an existing R function that I am not aware of.

Basically, recalling back to the structure of a "phylo" object in memory, we just need to create an edge matrix containing n+1 (for n species) in the left column and 1:n in the right column:

edge<-matrix(NA,n,2)
edge[,1]<-n+1; edge[,2]<-1:n


Then we can optionally assign branch lengths, if they have been provided by the user:

if(!is.null(branch.lengths)) edge.length=branch.lengths

Finally we combine these elements into a list and assign the class "phylo":

phy<-list(edge=edge,edge.length=edge.length,Nnode=1,
   tip.label=as.vector(species))
class(phy)<-"phylo"


and we're done.

3 comments:

  1. Well, there is stree() in ape which already does this, but...

    > system.time(for(i in 1:100000){x<-stree(100)})
    user system elapsed
    15.13 0.00 15.46
    > system.time(for(i in 1:100000){x<-starTree(100)})
    user system elapsed
    2.91 0.00 2.98

    ...it appears your function runs faster. stree() is no slouch though, at ~15 seconds for 100000 star with 100 tips.
    -Dave

    ReplyDelete
  2. Thanks Dave. stree() looks like a handy function as it also can be used to create balanced and pectinate trees. Thanks for bringing it to my attention. My function also allows the inclusion of branch lengths. [Of course, this would also be possible with stree() and just one additional step.]

    ReplyDelete
    Replies
    1. Hi Liam,
      I wonder how can I add branch length to the starTree function?

      Thanks

      Delete