This may be unnecessarily complicated, but for reasons that I will address in a future post, I wanted to write some code to resolve polytomies in all possible ways.
Now, there exists a function in ape, multi2di, that resolves
polytomies randomly and by grouping taxa in the order that they appear,
but (so far as I know), no prior R function finds all possible resolutions of
a polytomy.
The code chunk below does this in the following way:
1) Given a tree and a particular node, it first finds all daughters of that
node. If this number is greater than 2, it generates all possible subclades
descended from that node using allTrees.
2) It then cuts the input tree below the polytomy (using
splitTree) and above the polytomy (using
extract.clade).
3) Finally it reassembles the trees using the set of all possible resolutions
obtained using allTrees.
Note that the number of possible resolutions for a polytomy with N descendants will be equal to the number of possible rooted trees with N taxa: 3 for 3; 15 for 4; 105 for 5; and so on.
Here is an example.
Load packages:
library(phytools)
library(phangorn)
Our function:
resolveNode<-function(tree,node){
dd<-Children(tree,node)
if(length(dd)>2){
EL<-!is.null(tree$edge.length)
if(!EL) tree<-compute.brlen(tree)
n<-length(dd)
tt<-lapply(allTrees(n,TRUE,dd),untangle,"read.tree")
ROOT<-node==(Ntip(tree)+1)
SPL<-if(!ROOT) splitTree(tree,split=list(node=node,
bp=tree$edge.length[which(tree$edge[,2]==node)])) else
list(NULL,tree)
KIDS<-Children(SPL[[2]],SPL[[2]]$edge[1,1])
KIDS<-setNames(KIDS,dd)[KIDS>Ntip(SPL[[2]])]
SUBS<-list()
if(length(KIDS)>0)
for(i in 1:length(KIDS))
SUBS[[i]]<-extract.clade(SPL[[2]],KIDS[i])
obj<-list()
for(i in 1:length(tt)){
tt[[i]]$edge.length<-rep(0,nrow(tt[[i]]$edge))
for(j in 1:Ntip(tt[[i]]))
tt[[i]]$edge.length[which(tt[[i]]$edge[,2]==j)]<-
tree$edge.length[which(tree$edge[,2]==
as.numeric(tt[[i]]$tip.label[j]))]
ind<-as.numeric(tt[[i]]$tip.label)<=Ntip(tree)
tt[[i]]$tip.label[ind]<-
tree$tip.label[as.numeric(tt[[i]]$tip.label[ind])]
if(length(KIDS)>0)
for(j in 1:length(KIDS))
tt[[i]]<-bind.tree(tt[[i]],SUBS[[j]],
where=which(tt[[i]]$tip.label==
names(KIDS)[j]))
obj[[i]]<-if(!ROOT) bind.tree(SPL[[1]],tt[[i]],
where=which(SPL[[1]]$tip.label=="NA")) else
tt[[i]]
if(!EL) obj[[i]]$edge.length<-NULL
}
class(obj)<-"multiPhylo"
} else obj<-tree
obj
}
Now let's generate a tree with polytomies:
tree<-read.tree(text="(((A1,A2),(B1,B2),C,D),E,F);")
plot(tree,type="cladogram",edge.width=1,no.margin=TRUE)
nodelabels()
First, resolve node 9. This should result in three trees:
trees<-resolveNode(tree,9)
trees
## 3 phylogenetic trees
par(mfrow=c(3,1))
par(mar=c(1,8,1,8))
nulo<-lapply(trees,plot,type="cladogram",edge.width=1,cex=1.1)
Alternatively, we could resolve node 10. Since this is a quadrifurcation, this resolution should result in 15 trees:
trees<-resolveNode(tree,10)
trees
## 15 phylogenetic trees
par(mfrow=c(5,3))
par(mar=rep(0.5,4))
nulo<-lapply(trees,plot,type="cladogram",edge.width=1,cex=0.8)
It is straightforward to also see how this could be used to obtain all possible trees - in this case the product of 3 & 15 = 45.
For instance:
resolveAllNodes<-function(tree){
foo<-function(node,tree) length(Children(tree,node))
nodes<-1:tree$Nnode+Ntip(tree) ## all nodes
nodes<-nodes[sapply(1:tree$Nnode+Ntip(tree),foo,
tree=tree)>2]
for(i in 1:length(nodes)){
if(i==1) obj<-resolveNode(tree,nodes[i])
else {
for(j in 1:length(obj)){
MM<-matchNodes(tree,obj[[j]])
NODE<-MM[which(MM[,1]==nodes[i]),2]
if(j==1) tmp<-resolveNode(obj[[j]],NODE)
else tmp<-c(tmp,resolveNode(obj[[j]],NODE))
}
obj<-tmp
}
}
obj
}
trees<-resolveAllNodes(tree)
par(mfrow=c(9,5))
par(mar=rep(0.25,4))
nulo<-lapply(trees,plot,type="cladogram",edge.width=1,cex=0.8)
One nuance that we have to keep in mind as we iterate over the unresolved
nodes in our input tree, is that the node numbers will update as we resolve
the tree! This is why I have to match nodes back to the original tree using
the function matchNodes.
That's it.
