Hi Todd.
Challenge accepted ;)
Try the attached code. In my tests it returns the exact same result as
midpoint in phangorn, but it does not use any phangorn code internally.
(It does depend on phytools & dependencies, particularly ape.) It is
probably much less elegant & slower than phangorn (I haven't tested
that). If it works, please let me know & I will add to phytools.
All the best, Liam
Liam J. Revell, Assistant Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://blog.phytools.org
On 1/30/2014 1:00 PM, Todd Oakley wrote:
Colleagues,
Does anyone know of an existing midpoint rooting routine? I am
displaying trees and would like to show them as midpoint rooted.
I've been using the Phangorn package, which does the midpoint rooting
perfectly, but it has some dependencies that make it unstable on the
linux machines I've been using. When I looked last, Phangorn was the
only package I could find with midpoint routing.
Thanks for any ideas...
Todd Oakley, UCSB
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at
http://www.mail-archive.com/r-sig-phylo@r-project.org/
## function for midpoint rooting
## written by Liam J. Revell 2014
midpoint.root<-function(tree){
D<-cophenetic(tree)
dd<-max(D)
ii<-which(D==dd)[1]
ii<-c(ceiling(ii/nrow(D)),ii%%nrow(D))
if(ii[2]==0) ii[2]<-nrow(D)
spp<-rownames(D)[ii]
nn<-which(tree$tip.label==spp[2])
tree<-reroot(tree,nn,tree$edge.length[which(tree$edge[,2]==nn)])
aa<-getAncestors(tree,which(tree$tip.label==spp[1]))
D<-c(0,dist.nodes(tree)[which(tree$tip.label==spp[1]),aa])
names(D)[1]<-which(tree$tip.label==spp[1])
i<-0
while(D[i+1]<(dd/2)) i<-i+1
tree<-reroot(tree,as.numeric(names(D)[i]),D[i+1]-dd/2)
tree
}
## function gets ancestor node numbers
## written by Liam J. Revell
getAncestors<-function(tree,node,type=c("all","parent")){
type<-type[1]
if(type=="all"){
aa<-vector()
rt<-length(tree$tip.label)+1
currnode<-node
while(currnode!=rt){
currnode<-getAncestors(tree,currnode,"parent")
aa<-c(aa,currnode)
}
return(aa)
} else if(type=="parent"){
aa<-tree$edge[which(tree$edge[,2]==node),1]
return(aa)
} else stop("do not recognize type")
}_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
