To respond to a post by Matthew Vavrek (and for fun, besides), I just created a function, exhaustiveMP(), that performs exhaustive as well as branch & bound MP tree searches. To do the parsimony calculations and the tree search it calls functions from {ape} and (especially) {phangorn}. I have posted this to my R-phylogenetics page (direct link to code, here).
The main challenge in doing this was adding a tip in all possible places on the tree (the subject of my last post); once this was solved, the rest was easily done using the {phangorn} functions parsimony() and allTrees(); and the {ape} function bind.tree().
This function could easily be modified to conduct ML exhaustive and branch & bound tree searches, using, for instance, pml() from the {phangorn} package. I use parsimony here only because Matthew's inquiry was about MP; and because computing the parsimony score is much faster.
Of course, nothing about this function is fast. It will in fact be very slow for more than 7 or 8 species. Note that method="exhaustive" allows no more than 10 species; and method="branch.and.bound" no more than 15. I have not tried the function with that many species.
I believe this function works properly, but I welcome any experiences or feedback.
Hi Liam,
ReplyDeletenice code. Here are two small additions. The first is a heuristic of the input order of the labels. We want to add always the taxa which is furthest apart from the taxa in the tree, to push the parsimony score over the bound as quick as possible for "bad" trees.
The second function I added computes parsimonious uninformative sites. There are more similar tricks which should be incorporated.
In the example below summaryRprof() showed that there is lots of potential of improving, especially bind.tree(). I takes about 4 sec. (+1 sec to reorder the tree), from 14 sec. in total. Computing the parsimony score takes only ~2.5 sec.
Here is the code:
parsinfo = function(x){
nr <- attr(x, "nr")
labels = attr(x, "allLevels")
result = matrix(NA, nrow = length(x), ncol = nr)
for (i in 1:length(x)) result[i, ] <- x[[i]]
contrast <- attr(x, "contrast")
D = tcrossprod(contrast)
l = dim(contrast)[1]
tab = apply(result, 2, tabulate, l)
tab2 = tab
tab2[] = FALSE
tab2[tab>1] = TRUE
tab3 = tab
tab3[] = FALSE
tab3[tab>0] = TRUE ## ==1
tab3 = t(tab3)
ps = rowSums(tab3) - 1
diag(D) = 0
tab4 = tab3 %*% D
tab4 = tab3 * tab4 #
ind = which(colSums(tab2)==1 & rowSums(tab4)==0)
ps = rowSums(tab3) - 1
cbind(ind, ps[ind])
}
getorder <- function(x){
label = names(x)
dm = as.matrix(dist.hamming(x, FALSE))
ind = which(dm==max(dm), arr.ind=TRUE)[1,]
added = label[ind]
remaining<-label[-ind]
for(i in 3:length(x)){
tmp = which.max(rowSums(dm[-ind, ind, drop=FALSE]))
added = c(added, remaining[tmp])
ind = c(ind, match(remaining[tmp], label))
remaining<-remaining[-tmp]
}
added
}
branch.and.bound.2 <- function(data, tree=NULL, trace=1, ...){
# first, compute the parsimony score on the input tree
tree = pratchet(data, start=tree, trace=trace-1, ...)
# use only parsimony informative sites
nr <- attr(data, "nr")
pis <- parsinfo(data)
p0 <- sum(attr(data, 'weight')[pis[,1]] * pis[,2])
data.orig <- data
data <- phangorn:::getRows(data, c(1:nr)[-pis[,1]], TRUE)
#compute input order
label = getorder(data)
data <- phangorn:::prepareDataFitch(data)
#compute upper bound
bound<- phangorn:::fit.fitch(tree,data)
# now pick three species to start
new<-list(stree(n=3,tip.label=label[1:3]))
class(new)<-"multiPhylo"
added<-new[[1]]$tip.label; remaining<-setdiff(label,added)
# branch & bound
while(length(remaining)>0){
old<-new; new<-list()
new.tip<-remaining[1]
pscores<-vector()
for(i in 1:length(old)){
temp<-add.everywhere(old[[i]],new.tip)
temp <- unclass(temp)
score = sapply(temp, phangorn:::fit.fitch, data)
# score<-fitch(temp,data)
new<-unlist(list(new,temp[score<=bound]),recursive=FALSE); class(new)<-"multiPhylo"
pscores<-c(pscores,score[score<=bound])
}
added<-c(added,new.tip)
if(trace) print(paste(length(added),"species added;",length(new),"trees retained",collapse=""))
remaining<-setdiff(tree$tip.label,added)
}
# ok, done, now sort what needs to be returned
trees<-new[pscores==min(pscores)]
for(i in 1:length(trees)) attr(trees[[i]],"pscore")<-min(pscores)+p0
return(trees) # return all mp trees
}
library(phangorn)
data(yeast)
tree = pratchet(yeast)
system.time(tree1 <- branch.and.bound(yeast, tree))
system.time(tree2 <- branch.and.bound.2(yeast, tree))
Great post Klaus. I will try it out when I get a chance.
ReplyDelete