Saturday, May 3, 2014

Fixed drop.tip.singleton

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
}

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