R/cophyloplot.R

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