Friday, August 3, 2012

Plotting phylogenies with lineages leading to extinct taxa highlighted

I just created a new function, fancyTree, in which I intend to build in some custom wrappers to ape plot.phylo as well as phytools functions plotTree, plotSimmap, and phenogram to help phytools users draw special types of "fancy" phylogenetic trees.

I'm open to suggestions, of course, but my first idea was to create a function that would allow users to plot phylogenies containing extinct taxa (lineages terminating before the end of the tree) in a matter so as to highlight these branches. In particular, let's say we want to color branches of the tree leading either only to extinct taxa or proceeding in time the MRCA of all extant taxa in the tree a different color than the remaining branches (i.e, branches descended from the MRCA of all extant species and leading to extant species).

This should be pretty easy so long as we can identify all the branches that we want to color in this way. So let's think about how to do this. The first thing we need to do is identify the extant species and their common ancestor. We can do this using two functions in phytools, as follows:

extant<-getExtant(tree)
ca<-findMRCA(tree,extant)


Now, let's create a vector of all the edges in the tree (identified by the tipward node, i.e., the right column of tree$edge). We'll use this to keep score of which branches we want to flag:

edges<-rep(0,nrow(tree$edge))
names(edges)<-tree$edge[,2]


Now, after testing to see if the MRCA of extant species is the root node, and assuming that it is not, we can find all the edges preceding the MRCA of extant taxa using the phangorn function Descendants and the base function setdiff, as follows:

z<-setdiff(Descendants(tree,root.node,type="all"), Descendants(tree,ca,type="all"))
# now mark them in vector edges
edges[as.character(z)]<-1


Finally, we can go through the remaining nodes of the tree and determine if the descendants of each node above the MRCA of all extant taxa include extant taxa, both extant and extinct species, or only extinct species. In the final case, we should flag the branch in our vector edges. Here, we can just get a vector of all the nodes tipward of the MRCA; and then go through that vector and try to match the tips descended from each node in that vector to our vector of extant species, extant. This can be done as follows:

# get descendants of the MRCA
z<-Descendants(tree,ca,type="all")
# compute list of all descendants for each node in z
y<-Descendants(tree,z)
# finally, loop over y and flag edges not leaving
# extant descendants
for(i in 1:length(z))
   if(!any(tree$tip.label[y[[i]]]%in%extant))
      edges[as.character(z[i])]<-1


Note that there is some admittedly lazy programming above, because as soon as we visit an edge with no extant descendants, we should not visit any of its descendants, but we do anyway! This should be a simple thing to fix.

That's pretty much it. We can then use plot.phylo to create our tree, and voila!

plot.phylo(tree,edge.color=edges+1,edge.lty=edges+1,edge.width=2, no.margin=TRUE)


Code for this function can be downloaded from my phytools page and should be in the next version of the package.

1 comment: