Dear Paolo,

I modified the phylomorphospace3d code some time ago to be able to plot
colored branches (see below). You can use pointcol and edgecol to color the
plot. These can be simple color specifications (such as "black") or vectors
of colors that have the length of your number of tips or edges,
respectively.

Thank you Liam for this piece of code! It's really usefull for me! I hope I
didn't spoil your code too much :)

Jonas




ColoredPhylomorphospace3d <- function (tree, X, A = NULL, label = TRUE,
pointcol="black", edgecol="black", control = list(),
          method = c("dynamic", "static"), ...)
{
  method <- method[1]
  if (class(tree) != "phylo")
    stop("tree object must be of class 'phylo.'")
  # controls festlegen
  con = list(spin = TRUE, axes = TRUE, box = TRUE, simple.axes = FALSE,
             lwd = 1, ftype = "reg")
  con[(namc <- names(control))] <- control
  if (con$simple.axes)
    con$box <- con$axes <- FALSE
  con$ftype <- which(c("off", "reg", "b", "i", "bi") == con$ftype) -
    1
  if (is.null(A))
    A <- apply(X, 2, function(x, tree) fastAnc(tree, x),
               tree = tree)
  else A <- A[as.character(1:tree$Nnode + length(tree$tip)),
              ]
  x <- y <- z <- matrix(NA, nrow(tree$edge), 2)
  X <- X[tree$tip.label, ]
  for (i in 1:length(tree$tip)) {
    x[tree$edge[, 2] == i, 2] <- X[i, 1]
    y[tree$edge[, 2] == i, 2] <- X[i, 2]
    z[tree$edge[, 2] == i, 2] <- X[i, 3]
  }
  for (i in length(tree$tip) + 1:tree$Nnode) {
    x[tree$edge[, 1] == i, 1] <- x[tree$edge[, 2] == i,
                                   2] <- A[as.character(i), 1]
    y[tree$edge[, 1] == i, 1] <- y[tree$edge[, 2] == i,
                                   2] <- A[as.character(i), 2]
    z[tree$edge[, 1] == i, 1] <- z[tree$edge[, 2] == i,
                                   2] <- A[as.character(i), 3]
  }
  if (is.null(colnames(X)))
    colnames(X) <- c("x", "y", "z")
  if (method == "dynamic") {
    params <- get("r3dDefaults")
    plot3d(rbind(X, A), xlab = colnames(X)[1], ylab = colnames(X)[2],
           zlab = colnames(X)[3], axes = con$axes, box = con$box,
           params = params)
    spheres3d(X, radius = 0.02 * mean(apply(X, 2, max) -
                                        apply(X, 2, min)), col=pointcol)
    for (i in 1:nrow(tree$edge)) segments3d(x[i, ], y[i, ], z[i, ], lwd =
con$lwd, col=edgecol[i])
    ms <- colMeans(X)
    rs <- apply(rbind(X, A), 2, range)
    if (con$simple.axes) {
      lines3d(x = rs[, 1], y = c(rs[1, 2], rs[1, 2]),
              z = c(rs[1, 3], rs[1, 3]))
      lines3d(x = c(rs[1, 1], rs[1, 1]), y = rs[, 2],
              z = c(rs[1, 3], rs[1, 3]))
      lines3d(x = c(rs[1, 1], rs[1, 1]), y = c(rs[1, 2],
                                               rs[1, 2]), z = rs[, 3])
    }
    rs <- rs[2, ] - rs[1, ]
    for (i in 1:length(tree$tip)) {
      adj <- 0.03 * rs * (2 * (X[i, ] > ms) - 1)
      if (con$ftype)
        text3d(X[i, ] + adj, texts = tree$tip.label[i],
               font = con$ftype)
    }
    if (con$spin) {
      xx <- spin3d(axis = c(0, 0, 1), rpm = 10)
      play3d(xx, duration = 5)
      invisible(xx)
    }
    else invisible(NULL)
  }
  else if (method == "static") {
    if (hasArg(angle))
      angle <- list(...)$angle
    else angle <- 30
    if (hasArg(xlim))
      xlim <- list(...)$xlim
    else xlim <- NULL
    if (hasArg(ylim))
      ylim <- list(...)$ylim
    else ylim <- NULL
    if (hasArg(zlim))
      zlim <- list(...)$zlim
    else zlim = NULL
    xx <- scatterplot3d(X, xlab = colnames(X)[1], zlab = colnames(X)[3],
                        pch = 19, angle = angle, ylab = colnames(X)[2],
                        asp = 1, cex.symbols = 1.3, xlim = xlim, ylim =
ylim,
                        zlim = zlim)
    aa <- xx$xyz.convert(A)
    points(aa$x, aa$y, pch = 19, cex = 0.8)
    for (i in 1:nrow(tree$edge)) {
      aa <- xx$xyz.convert(x[i, ], y[i, ], z[i, ])
      lines(aa$x, aa$y, lwd = 2)
    }
    for (i in 1:length(tree$tip.label)) {
      aa <- xx$xyz.convert(x[which(tree$edge[, 2] == i),
                             2], y[which(tree$edge[, 2] == i), 2],
z[which(tree$edge[,

2] == i), 2])
      text(tree$tip.label[i], x = aa$x, y = aa$y, pos = 2)
    }
    invisible(xx)
  }
}

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Reply via email to