cophyloplot2colors.R

balbuena — May 29, 2014, 12:35 PM

# Demostration of cophyloplot where congruent and incongruent links
# are displayed in two different colors  
# Hosts(H) and parasites (P) share a common topology (com.tree)
# with 1 to 1 (congruent) H-P links.
# Other part of the tree is different between H and P with random
# (incongruent) H-P links. 
library(ape)
# set numbers of congruent and incongruent links
Nc = 6
Ni = 4
com.tree <- rtree(Nc) # topology shared by H and P trees
H.tree <- rcoal(n=Ni, tip.label=paste("I",1:Ni)) # H different topology
P.tree <- rcoal(n=Ni, tip.label=paste("I",1:Ni)) # P different topology
#
# bind trees
com.tree$root.edge <- H.tree$root.edge <- P.tree$root.edge <- runif(1)
H.tree <- com.tree + H.tree
P.tree <- com.tree + P.tree
#
# HPgamma creates HP matrix with coevl links in common part of tree and
# random links RESTRICTED TO random part of the tree
HPgamma <- function (NH, NP, Nrand) { # input # rows and # random links
  HP <- diag(1, NH, NP) #diagonal 1 matrix of NH rows and NP cols
  HP2 <- diag(1, Nrand)
  #The following randomizes the links of the non-coevolving part of the tree:
  i <- sapply(seq_len(Nrand), function(x) sample(seq_len(Nrand)[-x], 1))
  HP2 <- matrix(0, Nrand, Nrand)
  HP2[cbind(seq_len(Nrand), i)] <- 1
  HP3 <- which(HP2 ==1, arr.in=TRUE) #get coordinates of random links   
  HP[HP3] <- 1
  for (i in 1:Nrand) HP[i,i] <-0
  return(HP)
}
# proper script starts here:
# 1. Sort assoc= input of cophyloplot. Incongruent links go first 
HP <- HPgamma(Nc+Ni,Nc+Ni,Ni)
rownames(HP)<- rev(H.tree$tip.label)
colnames(HP)<- rev(P.tree$tip.label)
HP.LUT <- which(HP ==1, arr.in=TRUE)
links <- cbind(rownames(HP)[HP.LUT[,1]],colnames(HP)[HP.LUT[,2]])
# 2. create a binary vector where 1 = incongruent and 0 = congruent
# order follow 'links' above
x <- c(rep(1, Ni), rep(0, Nc))
# 3. color vector red = incongruent, blue = congruent
cols <- c("blue","red")[(x==1)+1] 
# 4. plot trees:
cophyloplot(H.tree, P.tree, assoc=links, use.edge.length=TRUE,
            gap=0.5, space=8, col= cols)

plot of chunk unnamed-chunk-1