Here is a fixed version of the function I just posted to drop leaves while retaining all ancestral nodes of remaining extant taxa as singletons:
drop.tip.singleton<-function(tree,tip){
N<-length(tree$tip.label)
m<-length(tip)
ii<-sapply(tip,function(x,y) which(y==x),y=tree$tip.label)
tree$tip.label<-tree$tip.label[-ii]
ii<-sapply(ii,function(x,y) which(y==x),y=tree$edge[,2])
tree$edge<-tree$edge[-ii,]
tree$edge.length<-tree$edge.length[-ii]
tree$edge[tree$edge<=N]<-
as.integer(rank(tree$edge[tree$edge<=N]))
tree$edge[tree$edge>N]<-tree$edge[tree$edge>N]-m
N<-N-m
if(any(sapply(tree$edge[tree$edge[,2]>N,2],"%in%",
tree$edge[,1])==FALSE)) internal<-TRUE
else internal<-FALSE
while(internal){
ii<-which(sapply(tree$edge[,2],"%in%",
c(1:N,tree$edge[,1]))==FALSE)[1]
nn<-tree$edge[ii,2]
tree$edge<-tree$edge[-ii,]
tree$edge.length<-tree$edge.length[-ii]
tree$edge[tree$edge>nn]<-tree$edge[tree$edge>nn]-1
tree$Nnode<-tree$Nnode-length(ii)
if(any(sapply(tree$edge[tree$edge[,2]>N,2],
"%in%",tree$edge[,1])==FALSE)) internal<-TRUE
else internal<-FALSE
}
tree
}
N<-length(tree$tip.label)
m<-length(tip)
ii<-sapply(tip,function(x,y) which(y==x),y=tree$tip.label)
tree$tip.label<-tree$tip.label[-ii]
ii<-sapply(ii,function(x,y) which(y==x),y=tree$edge[,2])
tree$edge<-tree$edge[-ii,]
tree$edge.length<-tree$edge.length[-ii]
tree$edge[tree$edge<=N]<-
as.integer(rank(tree$edge[tree$edge<=N]))
tree$edge[tree$edge>N]<-tree$edge[tree$edge>N]-m
N<-N-m
if(any(sapply(tree$edge[tree$edge[,2]>N,2],"%in%",
tree$edge[,1])==FALSE)) internal<-TRUE
else internal<-FALSE
while(internal){
ii<-which(sapply(tree$edge[,2],"%in%",
c(1:N,tree$edge[,1]))==FALSE)[1]
nn<-tree$edge[ii,2]
tree$edge<-tree$edge[-ii,]
tree$edge.length<-tree$edge.length[-ii]
tree$edge[tree$edge>nn]<-tree$edge[tree$edge>nn]-1
tree$Nnode<-tree$Nnode-length(ii)
if(any(sapply(tree$edge[tree$edge[,2]>N,2],
"%in%",tree$edge[,1])==FALSE)) internal<-TRUE
else internal<-FALSE
}
tree
}
The previous version ran into trouble when the number of remaining internal edges that the function tried to remove in a step is >1. This version is less ambitious & seems to work fine.
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.