Coloring the edges of the plotted trees in a cophylo plot

The other day I was emailed the following question:

“I would like to know if the function recently implemented in phytools (cophylo) allows for colorizing internal clades in the tree. The reason I wanted to do this is to show a more detailed view of the relationships between the two trees that are being compared with the function…. If this feature is not available in phytools, would you be able to add it?”

Indeed, this feature is not available in the current version of phytools. plot.cophylo uses a different function for plotting trees - the internal phytools function, phylogram, rather than the function plotSimmap that is used by many of phytools' other plotting methods.

However, it is indeed possible to add. I have just pushed this update to GitHub.

The most complicated thing is that in the existing tree plotter vertical lines are plotted as a single line segment; however, for the case of coloring different edges, we can anticipate a case in which the vertical lines of our right- or left-facing phylograms have descendant edges with different colors. For this circumstance, we want to be able to split the line into two differently colored segments. I accomplished this by substituting the function segments for lines in plotting the vertical lines of the tree.

Even more complicated is the situation in which multiple edges descend from a single node - potentially with different colors. In this case, I decided to pick the midpoint of the vertical line & plot line segments to the north or south of this midpoint the same color as their respective descendant edges.

First, let's try it using random color assignments:

library(phytools)
packageVersion("phytools")

## [1] '0.6.6'

t1

## 
## Phylogenetic tree with 26 tips and 19 internal nodes.
## 
## Tip labels:
##  A, S, I, X, Q, K, ...
## 
## Unrooted; includes branch lengths.

t2

## 
## Phylogenetic tree with 26 tips and 19 internal nodes.
## 
## Tip labels:
##  Q, P, Y, O, H, T, ...
## 
## Rooted; includes branch lengths.

obj<-cophylo(t1,t2)

## Rotating nodes to optimize matching...
## Done.

plot(obj,link.type="curved")

edge.col<-list(
    left=sample(c("red","blue"),nrow(t1$edge),replace=TRUE),
    right=sample(c("red","blue"),nrow(t2$edge),replace=TRUE))
edge.col

## $left
##  [1] "red"  "blue" "blue" "red"  "red"  "blue" "blue" "blue" "blue" "blue"
## [11] "red"  "red"  "blue" "blue" "red"  "blue" "red"  "red"  "red"  "red" 
## [21] "red"  "red"  "red"  "blue" "blue" "red"  "red"  "red"  "blue" "red" 
## [31] "red"  "blue" "blue" "red"  "red"  "red"  "red"  "blue" "red"  "blue"
## [41] "red"  "blue" "blue" "blue"
## 
## $right
##  [1] "red"  "red"  "blue" "blue" "red"  "red"  "red"  "blue" "red"  "blue"
## [11] "red"  "blue" "blue" "blue" "red"  "red"  "red"  "red"  "blue" "blue"
## [21] "red"  "blue" "blue" "blue" "red"  "red"  "blue" "blue" "blue" "red" 
## [31] "blue" "red"  "red"  "red"  "red"  "blue" "red"  "blue" "blue" "blue"
## [41] "red"  "blue" "red"  "blue"

plot(obj,link.type="curved",edge.col=edge.col)

We can also color clades or parts of the tree as desired. This requires identifying the set of edges in the rotated trees that correspond with the clades we would like to color. For instance:

plot(obj,link.type="curved")
nodelabels.cophylo(which="left")
nodelabels.cophylo(which="right")

left<-rep("black",nrow(obj$trees[[1]]$edge))
nodes<-getDescendants(obj$trees[[1]],36)
left[sapply(nodes,function(x,y) which(y==x),y=obj$trees[[1]]$edge[,2])]<-
    "red"
nodes<-getDescendants(obj$trees[[1]],43)
left[sapply(nodes,function(x,y) which(y==x),y=obj$trees[[1]]$edge[,2])]<-
    "blue"
right<-rep("black",nrow(obj$trees[[2]]$edge))
nodes<-getDescendants(obj$trees[[2]],30)
right[sapply(nodes,function(x,y) which(y==x),y=obj$trees[[2]]$edge[,2])]<-
    "red"
nodes<-getDescendants(obj$trees[[2]],42)
right[sapply(nodes,function(x,y) which(y==x),y=obj$trees[[2]]$edge[,2])]<-
    "blue"
edge.col<-list(left=left,right=right)
edge.col

## $left
##  [1] "black" "black" "black" "black" "black" "black" "black" "black"
##  [9] "black" "black" "black" "black" "black" "black" "black" "black"
## [17] "black" "black" "black" "red"   "red"   "red"   "red"   "red"  
## [25] "red"   "red"   "red"   "red"   "red"   "red"   "red"   "red"  
## [33] "black" "black" "black" "black" "blue"  "blue"  "blue"  "blue" 
## [41] "blue"  "blue"  "blue"  "black"
## 
## $right
##  [1] "black" "black" "black" "black" "black" "red"   "red"   "red"  
##  [9] "red"   "red"   "red"   "red"   "red"   "red"   "red"   "red"  
## [17] "red"   "red"   "red"   "red"   "red"   "red"   "red"   "red"  
## [25] "red"   "red"   "red"   "red"   "black" "black" "black" "black"
## [33] "black" "black" "black" "blue"  "blue"  "blue"  "blue"  "blue" 
## [41] "blue"  "blue"  "blue"  "blue"

plot(obj,link.type="curved",edge.col=edge.col,lwd=2)

Hopefully, you get the idea.

The trees for this exercise were simulated as follows:

t1<-rtree(n=26,tip.label=LETTERS)
t2<-rtree(n=26,tip.label=LETTERS)
internal<-which(t1$edge[,2]>Ntip(t1))
t1$edge.length[sample(internal,6)]<-0
t1<-di2multi(t1)
internal<-which(t2$edge[,2]>Ntip(t2))
t2$edge.length[sample(internal,6)]<-0
t2<-di2multi(t2)