Dropping the same tip from many trees in a "multiPhylo" object

Following the same theme as an earlier post, I thought I would take a closer look a Google search strings that have led readers to my site. A recent search string that has multiple hits is the following: "drop.tip multiphylo", by which I presume they mean to ask how one drops the same tip (or set of tips) from a set of phylogenies in an object of class "multiPhylo". It turns out that we can do this very easily, in fact.

First, let's get a set of trees as an object of the class "multiPhylo". [In most cases, of course, these would probably be trees we have read in from file (say, as a posterior sample of trees in a Bayesian analysis). In this case, we will just use the ape function rmtree().]

> require(ape)
> trees<-rmtree(N=10,n=20)
> trees
10 phylogenetic trees
> require(phytools)
> plotTree(trees[[6]]) # for instance



This gives a set of 10 trees, each containing 20 taxa, with the tip labels "t1" to "t20".

Now, a "multiPhylo" object is really just a list of trees, so to apply the ape function drop.tip to all the trees in the list, we can just use the function lapply. Let's drop, say "t13" (for luck) from all the trees in the list:

> ntrees<-lapply(trees,drop.tip,tip="t13")
> class(trees)<-"multiPhylo"
> plotTree(ntrees[[6]]) # for comparison



and taxon "t13" is gone, just as we'd hope. (The same is true, of course, of all the trees in our set.)

Naturally, the same technique works for a set of tips, for instance:

> ntrees<-lapply(trees,drop.tip,tip=c("t11","t12","t13"))
> class(ntrees)<-"multiPhylo"

We can even combine this with the trick, described in an earlier post that allows us to drop all but the set of taxa we have in a list. In this case, we will make the special function keep.tip. For instance if we wanted to keep only "t1" through "t10" in each of our trees above, we could just do the following:

> species<-paste("t",1:10,sep="")
> species
[1] "t1" "t2" "t3" "t4" "t5" "t6" "t7" "t8" "t9" "t10"
> keep.tip<-function(tree,tip) drop.tip(tree,setdiff(tree$tip.label,tip))
> ntrees<-lapply(trees,keep.tip,tip=species)
> "multiPhylo"->class(ntrees)
> plotTree(ntrees[[6]])



[Thanks to Dan Rabosky for the setdiff() suggestion.] Cool.