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.
Well, there is stree() in ape which already does this, but...
ReplyDelete> 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
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.]
ReplyDeleteHi Liam,
DeleteI wonder how can I add branch length to the starTree function?
Thanks