5

I would like to plot two phylogenies opposite each other in R using the ape package. One tree has 40 nodes and one has 26 nodes:

library(ape)
tree1 <- rtree(40)
tree2 <- rtree(26)

The cophyloplot function plots these face to face with specified links.

I'm having trouble specifying the links.

Note that in my actual nexus tree files, tip labels are text (and I'm unsure how to change these to numbers, if necessary...).

The links should be as follows:

If, in the tree1 nexus file, the tip labels of the sequences are 1-40. In the tree2 nexus file, tip labels are 1-26. Then links should be:

a <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40)
b <- c(14,1,4,1,9,12,2,10,6,3,13,5,14,15,18,19,19,7,14,9,10,11,25,22,21,16,23,24,26,17,1,12,12,21,15,16,21,8,20,21) 
association <- cbind(a, b)

(i.e. sequence 1 in tree1 is linked with sequence 14 in tree2)

So, I use something like this to plot the trees:

cophyloplot(tree1, tree2, assoc=association,length.line=4, space=28, gap=10, rotate=TRUE)

And calculate the distance matrix:

dist.topo(tree1, tree2, method = "PH85")

I'm not quite sure where I'm going wrong here. Any help would be appreciated!

SlowLoris
  • 995
  • 6
  • 28
user2861089
  • 1,205
  • 4
  • 22
  • 44

2 Answers2

5

To plot the trees, try this

library(ape)
set.seed(1)

# create trees
tree1 <- rtree(40)
tree2 <- rtree(26)

# modify tip labels
tree1$tip.label <- sub("t", "", tree1$tip.label, fixed = T)
tree2$tip.label <- sub("t", "", tree2$tip.label, fixed = T)

# create associations matrix 
a <- as.character(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40))
b <- as.character(c(14,1,4,1,9,12,2,10,6,3,13,5,14,15,18,19,19,7,14,9,10,11,25,22,21,16,23,24,26,17,1,12,12,21,15,16,21,8,20,21)) 
association <- cbind(a, b)

# plot
cophyloplot(tree1, tree2, assoc = association, length.line = 4, space = 28, gap = 3)

enter image description here

lukeA
  • 53,097
  • 5
  • 97
  • 100
  • thanks!! I'm getting the following error: `Error in decx[i] <- strwidth(x$tip.label[lsa[x$tip.label == assoc[i, 1]]]) : replacement has length zero` So I tried `tree1$tip.label<-(1:40) tree2$tip.label<-(1:26)` – user2861089 Jul 08 '15 at 23:57
  • When I run your code @lukeA I get the same trees but it is unreadable and needs to be expanded. How can I expand it to look like yours? – Jamie Nov 18 '20 at 00:38
1

The cophyloplot function does not require tip label indexing. One can refer to the taxa with their names. Note that lukeA's answer stores numbers from the associations as character. Transforming them to text that corresponds to tip labels and plotting the two trees shows the same result.

association <- apply(association, 2, function(x) sub("^","t", x))

head(association)
#      a    b    
# [1,] "t1" "t14"
# [2,] "t2" "t1" 
# [3,] "t3" "t4" 
# [4,] "t4" "t1" 
# [5,] "t5" "t9" 

cophyloplot(tree1, tree2, assoc=association, length.line=4, space=28, gap=3)

The order, in which the associations are listed in the matrix, is irrelevant. The best practice would be to import them from an external file with read.table().

nya
  • 2,138
  • 15
  • 29