R/cophyloplot.R

   1 ## cophyloplot.R (2012-02-14)
   2
   3 ##   Plots two phylogenetic trees face to
   4 ##   face with the links between the tips
   5
   6 ## Copyright 2008-2010 Damien de Vienne
   7
   8 ## This file is part of the R-package `ape'.
   9 ## See the file ../COPYING for licensing issues.
  10
  11 cophyloplot <-
  12     function(x, y, assoc = NULL, use.edge.length = FALSE, space = 0,
  13              length.line = 1, gap = 2, type = "phylogram", rotate = FALSE,
  14              col = par("fg"), lwd = par("lwd"), lty = par("lty"),
  15              show.tip.label = TRUE, font = 3, ...)
  16 {
  17     if (is.null(assoc)) {
  18         assoc <- matrix(ncol = 2)
  19         print("No association matrix specified. Links will be omitted.")
  20     }
  21     if (rotate == TRUE) {
  22         cat("\n    Click on a node to rotate (right click to exit)\n\n")
  23         repeat {
  24             res <- plotCophylo2(x, y, assoc = assoc, use.edge.length = use.edge.length,
  25                 space = space, length.line = length.line, gap = gap,
  26                 type = type, return = TRUE, col = col, lwd=lwd, lty=lty, show.tip.label = show.tip.label,
  27                 font = font)
  28             click <- identify(res$c[, 1], res$c[, 2], n = 1)
  29             if (click < length(res$a[, 1]) + 1) {
  30                 if (click > res$N.tip.x)
  31                   x <- rotate(x, click)
  32             } else if (click < length(res$c[, 1]) + 1) {
  33                 if (click > length(res$a[, 1]) + res$N.tip.y)
  34                   y <- rotate(y, click - length(res$a[, 1]))
  35             }
  36         }
  37         on.exit(print("done"))
  38     }
  39     else plotCophylo2(x, y, assoc = assoc, use.edge.length = use.edge.length,
  40         space = space, length.line = length.line, gap = gap,
  41         type = type, return = FALSE, col = col, lwd=lwd, lty=lty, show.tip.label = show.tip.label, font = font)
  42 }
  43
  44 plotCophylo2 <-
  45     function (x, y, assoc = assoc, use.edge.length = use.edge.length,
  46               space = space, length.line = length.line, gap = gap,
  47               type = type, return = return, col = col, lwd=lwd, lty=lty,
  48               show.tip.label = show.tip.label,
  49               font = font, ...)
  50 {
  51     res <- list()
  52 ###choice of the minimum space between the trees
  53     left <- max(nchar(x$tip.label, type = "width")) + length.line
  54     right <- max(nchar(y$tip.label, type = "width")) + length.line
  55     space.min <- left + right + gap * 2
  56     if ((space <= 0) || (space < space.min)) space <- space.min
  57     N.tip.x <- Ntip(x)
  58     N.tip.y <- Ntip(y)
  59     res$N.tip.x <- N.tip.x
  60     res$N.tip.y <- N.tip.y
  61     a <- plotPhyloCoor(x, use.edge.length = use.edge.length, type = type)
  62     res$a <- a
  63     b <- plotPhyloCoor(y, use.edge.length = use.edge.length,
  64                        direction = "leftwards", type = type)
  65 ###for the two trees to have the extreme leaves at the same ordinate.
  66     a[, 2] <- a[, 2] - min(a[, 2])
  67     b[, 2] <- b[, 2] - min(b[, 2])
  68     res$b <- b
  69     b2 <- b
  70     b2[, 1] <- b[1:nrow(b), 1] * (max(a[, 1])/max(b[, 1])) +
  71         space + max(a[, 1])
  72     b2[, 2] <- b[1:nrow(b), 2] * (max(a[, 2])/max(b[, 2]))
  73     res$b2 <- b2
  74     c <- matrix(ncol = 2, nrow = nrow(a) + nrow(b))
  75     c[1:nrow(a), ] <- a[1:nrow(a), ]
  76     c[nrow(a) + 1:nrow(b), 1] <- b2[, 1]
  77     c[nrow(a) + 1:nrow(b), 2] <- b2[, 2]
  78     res$c <- c
  79     plot(c, type = "n", xlim = NULL, ylim = NULL, log = "", main = NULL,
  80         sub = NULL, xlab = NULL, ylab = NULL, ann = FALSE, axes = FALSE,
  81         frame.plot = FALSE)
  82  ###segments for cladograms
  83    if (type == "cladogram") {
  84         for (i in 1:(nrow(a) - 1)) segments(a[x$edge[i, 1], 1],
  85             a[x$edge[i, 1], 2], a[x$edge[i, 2], 1], a[x$edge[i,
  86                 2], 2], col="red")
  87         for (i in 1:(nrow(b) - 1))
  88             segments(b2[y$edge[i, 1], 1], b2[y$edge[i, 1], 2],
  89                      b2[y$edge[i, 2], 1], b2[y$edge[i, 2], 2])
  90     }
  91 ###segments for phylograms
  92     if (type == "phylogram") {
  93         for (i in (N.tip.x + 1):nrow(a)) {
  94             l <- length(x$edge[x$edge[, 1] == i, ][, 1])
  95             for (j in 1:l) {
  96                 segments(a[x$edge[x$edge[, 1] == i, ][1, 1],
  97                   1], a[x$edge[x$edge[, 1] == i, 2], 2][1], a[x$edge[x$edge[,
  98                   1] == i, ][1, 1], 1], a[x$edge[x$edge[, 1] ==
  99                   i, 2], 2][j])
 100                 segments(a[x$edge[x$edge[, 1] == i, ][1, 1],
 101                   1], a[x$edge[x$edge[, 1] == i, 2], 2][j], a[x$edge[x$edge[,
 102                   1] == i, 2], 1][j], a[x$edge[x$edge[, 1] ==
 103                   i, 2], 2][j])
 104             }
 105         }
 106         for (i in (N.tip.y + 1):nrow(b)) {
 107             l <- length(y$edge[y$edge[, 1] == i, ][, 1])
 108             for (j in 1:l) {
 109                 segments(b2[y$edge[y$edge[, 1] == i, ][1, 1],
 110                   1], b2[y$edge[y$edge[, 1] == i, 2], 2][1],
 111                   b2[y$edge[y$edge[, 1] == i, ][1, 1], 1], b2[y$edge[y$edge[,
 112                     1] == i, 2], 2][j])
 113                 segments(b2[y$edge[y$edge[, 1] == i, ][1, 1],
 114                   1], b2[y$edge[y$edge[, 1] == i, 2], 2][j],
 115                   b2[y$edge[y$edge[, 1] == i, 2], 1][j], b2[y$edge[y$edge[,
 116                     1] == i, 2], 2][j])
 117             }
 118         }
 119     }
 120     if (show.tip.label) {
 121         text(a[1:N.tip.x, ], cex = 0, font = font, pos = 4,
 122              labels = x$tip.label)
 123         text(b2[1:N.tip.y, ], cex = 1, font = font, pos = 2,
 124              labels = y$tip.label)
 125     }
 126 ###links between associated taxa. Takes into account the size of the character strings of the taxa names.
 127     lsa <- 1:N.tip.x
 128     lsb <- 1:N.tip.y
 129     decx <- array(nrow(assoc))
 130     decy <- array(nrow(assoc))
 131
 132     #colors
 133     if (length(col)==1) colors<-c(rep(col, nrow(assoc)))
 134     else if (length(col)>=nrow(assoc)) colors<-col
 135     else  colors<-c(rep(col, as.integer(nrow(assoc)/length(col))+1))
 136
 137     #lwd
 138     if (length(lwd)==1) lwidths<-c(rep(lwd, nrow(assoc)))
 139     else if (length(lwd)>=nrow(assoc)) lwidths<-lwd
 140     else  lwidths<-c(rep(lwd, as.integer(nrow(assoc)/length(lwd))+1))
 141
 142     #lty
 143     if (length(lty) == 1) ltype <- c(rep(lty, nrow(assoc)))
 144     else if (length(lty) >= nrow(assoc)) ltype <- lty
 145     else  ltype <- c(rep(lty, as.integer(nrow(assoc)/length(lty))+1))
 146
 147     for (i in 1:nrow(assoc)) {
 148         if (show.tip.label) {
 149             decx[i] <- strwidth(x$tip.label[lsa[x$tip.label == assoc[i, 1]]])
 150             decy[i] <- strwidth(y$tip.label[lsb[y$tip.label == assoc[i, 2]]])
 151         } else {
 152             decx[i] <- decy[i] <- 0
 153         }
 154
 155         segments(a[lsa[x$tip.label == assoc[i, 1]], 1] + decx[i] + gap,
 156                  a[lsa[x$tip.label == assoc[i, 1]], 2],
 157                  a[lsa[x$tip.label == assoc[i, 1]], 1] + gap + left,
 158                  a[lsa[x$tip.label ==  assoc[i, 1]], 2],
 159                  col = colors[i], lwd = lwidths[i], lty = ltype[i])
 160
 161         segments(b2[lsb[y$tip.label == assoc[i, 2]], 1] - (decy[i] + gap),
 162                  b2[lsb[y$tip.label == assoc[i, 2]], 2],
 163                  b2[lsb[y$tip.label == assoc[i, 2]], 1] - (gap + right),
 164                  b2[lsb[y$tip.label ==  assoc[i, 2]], 2],
 165                  col = colors[i], lwd = lwidths[i], lty = ltype[i])
 166
 167         segments(a[lsa[x$tip.label == assoc[i, 1]], 1] + gap + left,
 168                  a[lsa[x$tip.label == assoc[i, 1]], 2],
 169                  b2[lsb[y$tip.label == assoc[i, 2]], 1] - (gap + right),
 170                  b2[lsb[y$tip.label == assoc[i, 2]], 2],
 171                  col = colors[i], lwd = lwidths[i], lty = ltype[i])
 172     }
 173     if (return == TRUE)  return(res)
 174 }