From a3ddfc06dd47c560b3ec5869ac104b0c68441eb1 Mon Sep 17 00:00:00 2001 From: paradis Date: Thu, 17 Apr 2008 12:46:51 +0000 Subject: [PATCH] new files by Damien + a few bug fixes git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@25 6e262413-ae40-0410-9e79-b911bd7a66b7 --- Changes | 11 +++ DESCRIPTION | 5 +- R/plot.cophylo.R | 173 ++++++++++++++++++++++++++++++++++++++++++++ R/plot.phylo.coor.R | 134 ++++++++++++++++++++++++++++++++++ R/subtreeplot.R | 47 ++++++++++++ R/subtrees.R | 55 ++++++++++++++ R/zoom.R | 8 +- Thanks | 22 +++--- man/plot.cophylo.Rd | 53 ++++++++++++++ man/subtreeplot.Rd | 39 ++++++++++ man/subtrees.Rd | 39 ++++++++++ 11 files changed, 569 insertions(+), 17 deletions(-) create mode 100644 R/plot.cophylo.R create mode 100644 R/plot.phylo.coor.R create mode 100644 R/subtreeplot.R create mode 100644 R/subtrees.R create mode 100644 man/plot.cophylo.Rd create mode 100644 man/subtreeplot.Rd create mode 100644 man/subtrees.Rd diff --git a/Changes b/Changes index 19e350a..5199d01 100644 --- a/Changes +++ b/Changes @@ -1,6 +1,14 @@ CHANGES IN APE VERSION 2.1-4 +NEW FEATURES + + o Four new functions have been written by Damien de Vienne for the + graphical exploration of large trees (plot.cophylo, subtrees, + subtreeplot), and to return the graphical coordinates of tree + (without plotting). + + BUG FIXES o read.dna() failed if "?" occurred in the first 10 sites of the @@ -12,6 +20,9 @@ BUG FIXES o Drawing the tip labels sometimes failed when plotting circular trees. + o zoom() failed when tip labels were used instead of their numbers + (thanks to Yan Wong for the fix). + OTHER CHANGES diff --git a/DESCRIPTION b/DESCRIPTION index d3be1eb..510eef4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,11 +1,12 @@ Package: ape Version: 2.1-4 -Date: 2008-03-28 +Date: 2008-04-17 Title: Analyses of Phylogenetics and Evolution Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Gangolf Jobb, Christoph Heibl, Vincent Lefort, Jim Lemon, - Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer + Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer, + Damien de Vienne Maintainer: Emmanuel Paradis Depends: R (>= 2.6.0) Suggests: gee, nlme, lattice diff --git a/R/plot.cophylo.R b/R/plot.cophylo.R new file mode 100644 index 0000000..7589f6a --- /dev/null +++ b/R/plot.cophylo.R @@ -0,0 +1,173 @@ +## plot.cophylo.R (2008-04-14) + +## Plots two phylogenetic trees face to +## face with the links between the tips + +## Copyright 2008 Damien de Vienne + +## This file is part of the R-package `ape'. +## See the file ../COPYING for licensing issues. + +plot.cophylo <- + function (x, y, assoc = NULL, use.edge.length = FALSE, space = 0, + length.line = 1, gap = 2, type = "phylogram", + rotate = FALSE, col = "red", show.tip.label = TRUE, + font = 3) +{ + if (is.null(assoc)) { + assoc <- matrix(ncol = 2) + print("No association matrix specified. Links will be omitted.") + } + if (rotate == TRUE) { + cat("\n Click on a node to rotate (right click to exit)\n\n") + repeat { + res <- plot.cophylo2(x, y, assoc = assoc, use.edge.length = use.edge.length, + space = space, length.line = length.line, gap = gap, + type = type, return = TRUE, col = col, show.tip.label = show.tip.label, + font = font) + click <- identify(res$c[, 1], res$c[, 2], n = 1) + if (click < length(res$a[, 1]) + 1) { + if (click > res$N.tip.x) + x <- rotate(x, click) + } + else if (click < length(res$c[, 1]) + 1) { + if (click > length(res$a[, 1]) + res$N.tip.y) + y <- rotate(y, click - length(res$a[, 1])) + } + plot.cophylo2(x, y, assoc = assoc, use.edge.length = use.edge.length, + space = space, length.line = length.line, gap = gap, + type = type, return = TRUE, col = col, show.tip.label = show.tip.label, + font = font) + } + on.exit(print("done")) + } + else plot.cophylo2(x, y, assoc = assoc, use.edge.length = use.edge.length, + space = space, length.line = length.line, gap = gap, + type = type, return = FALSE, col = col, show.tip.label = show.tip.label, + font = font) +} + +plot.cophylo2 <- + function (x, y, assoc = assoc, use.edge.length = use.edge.length, + space = space, length.line = length.line, gap = gap, + type = type, return = return, col = col, + show.tip.label = show.tip.label, font = font) +{ + res <- list() + +###choice of the minimum space between the trees### + left <- max(nchar(x$tip.label, type = "width")) + length.line + right <- max(nchar(y$tip.label, type = "width")) + length.line + space.min <- left + right + gap * 2 + if ((space <= 0) || (space < space.min)) + space <- space.min + + N.tip.x <- Ntip(x) + N.tip.y <- Ntip(y) + res$N.tip.x <- N.tip.x + res$N.tip.y <- N.tip.y + a <- plot.phylo.coor(x, use.edge.length = use.edge.length, + type = type) + res$a <- a + b <- plot.phylo.coor(y, use.edge.length = use.edge.length, + direction = "leftwards", type = type) + +###for the two trees to have the extreme leaves at the same ordinate. + a[, 2] <- a[, 2] - min(a[, 2]) + b[, 2] <- b[, 2] - min(b[, 2]) + + res$b <- b + + b2 <- b + b2[, 1] <- b[1:nrow(b), 1] * (max(a[, 1])/max(b[, 1])) + + space + max(a[, 1]) + b2[, 2] <- b[1:nrow(b), 2] * (max(a[, 2])/max(b[, 2])) + + res$b2 <- b2 + + c <- matrix(ncol = 2, nrow = nrow(a) + nrow(b)) + c[1:nrow(a), ] <- a[1:nrow(a), ] + c[nrow(a) + 1:nrow(b), 1] <- b2[, 1] + c[nrow(a) + 1:nrow(b), 2] <- b2[, 2] + res$c <- c + + plot(c, type = "n", xlim = NULL, ylim = NULL, log = "", main = NULL, + sub = NULL, xlab = NULL, ylab = NULL, ann = FALSE, axes = FALSE, + frame.plot = FALSE) + +###segments for cladograms + if (type == "cladogram") { + for (i in 1:(nrow(a) - 1)) + segments(a[x$edge[i, 1], 1], a[x$edge[i, 1], 2], + a[x$edge[i, 2], 1], a[x$edge[i, 2], 2]) + for (i in 1:(nrow(b) - 1)) + segments(b2[y$edge[i, 1], 1], b2[y$edge[i, 1], 2], + b2[y$edge[i, 2], 1], b2[y$edge[i, 2], 2]) + } + +###segments for phylograms + if (type == "phylogram") { + for (i in (N.tip.x + 1):nrow(a)) { + l <- length(x$edge[x$edge[, 1] == i, ][, 1]) + for (j in 1:l) { + segments(a[x$edge[x$edge[, 1] == i, ][1, 1], + 1], a[x$edge[x$edge[, 1] == i, 2], 2][1], a[x$edge[x$edge[, + 1] == i, ][1, 1], 1], a[x$edge[x$edge[, 1] == + i, 2], 2][j]) + segments(a[x$edge[x$edge[, 1] == i, ][1, 1], + 1], a[x$edge[x$edge[, 1] == i, 2], 2][j], a[x$edge[x$edge[, + 1] == i, 2], 1][j], a[x$edge[x$edge[, 1] == + i, 2], 2][j]) + } + } + for (i in (N.tip.y + 1):nrow(b)) { + l <- length(y$edge[y$edge[, 1] == i, ][, 1]) + for (j in 1:l) { + segments(b2[y$edge[y$edge[, 1] == i, ][1, 1], + 1], b2[y$edge[y$edge[, 1] == i, 2], 2][1], + b2[y$edge[y$edge[, 1] == i, ][1, 1], 1], b2[y$edge[y$edge[, + 1] == i, 2], 2][j]) + segments(b2[y$edge[y$edge[, 1] == i, ][1, 1], + 1], b2[y$edge[y$edge[, 1] == i, 2], 2][j], + b2[y$edge[y$edge[, 1] == i, 2], 1][j], b2[y$edge[y$edge[, + 1] == i, 2], 2][j]) + } + } + } + if (show.tip.label) { + text(a[1:N.tip.x, ], cex = 0, font = font, pos = 4, + labels = x$tip.label) + text(b2[1:N.tip.y, ], cex = 1, font = font, pos = 2, + labels = y$tip.label) + } + +###links between associated taxa. Takes into account the size of the character strings of the taxa names. + lsa <- 1:N.tip.x + lsb <- 1:N.tip.y + decx <- array(nrow(assoc)) + decy <- array(nrow(assoc)) + for (i in 1:nrow(assoc)) { + if (show.tip.label) { + decx[i] <- strwidth(x$tip.label[lsa[x$tip.label == + assoc[i, 1]]]) + decy[i] <- strwidth(y$tip.label[lsb[y$tip.label == + assoc[i, 1]]]) + } else { + decx[i] <- decy[i] <- 0 + } + segments(a[lsa[x$tip.label == assoc[i, 1]], 1] + decx[i] + + gap, a[lsa[x$tip.label == assoc[i, 1]], 2], a[lsa[x$tip.label == + assoc[i, 1]], 1] + gap + left, a[lsa[x$tip.label == + assoc[i, 1]], 2], col = col) + segments(b2[lsb[y$tip.label == assoc[i, 2]], 1] - (decy[i] + + gap), b2[lsb[y$tip.label == assoc[i, 2]], 2], b2[lsb[y$tip.label == + assoc[i, 2]], 1] - (gap + right), b2[lsb[y$tip.label == + assoc[i, 2]], 2], col = col) + segments(a[lsa[x$tip.label == assoc[i, 1]], 1] + gap + + left, a[lsa[x$tip.label == assoc[i, 1]], 2], b2[lsb[y$tip.label == + assoc[i, 2]], 1] - (gap + right), b2[lsb[y$tip.label == + assoc[i, 2]], 2], col = col) + } + if (return == TRUE) + return(res) +} diff --git a/R/plot.phylo.coor.R b/R/plot.phylo.coor.R new file mode 100644 index 0000000..b69fb95 --- /dev/null +++ b/R/plot.phylo.coor.R @@ -0,0 +1,134 @@ +## plot.phylo.coor.R (2008-04-14) + +## Coordinates of a Tree Plot + +## Copyright 2008 Damien de Vienne + +## This file is part of the R-package `ape'. +## See the file ../COPYING for licensing issues. + +plot.phylo.coor <- + function (x, type = "phylogram", use.edge.length = TRUE, node.pos = NULL, + direction = "rightwards") +{ + Ntip <- length(x$tip.label) + if (Ntip == 1) + stop("found only one tip in the tree!") + Nedge <- dim(x$edge)[1] + if (any(tabulate(x$edge[, 1]) == 1)) + stop("there are single (non-splitting) nodes in your tree; you may need to use collapse.singles().") + Nnode <- x$Nnode + if (is.null(x$edge.length)) use.edge.length <- FALSE + phyloORclado <- type %in% c("phylogram", "cladogram") + horizontal <- direction %in% c("rightwards", "leftwards") + if (phyloORclado) { + if (!is.null(attr(x, "order"))) + if (attr(x, "order") == "pruningwise") + x <- reorder(x) + yy <- numeric(Ntip + Nnode) + TIPS <- x$edge[x$edge[, 2] <= Ntip, 2] + yy[TIPS] <- 1:Ntip + + } + + xe <- x$edge + x <- reorder(x, order = "pruningwise") + ereorder <- match(x$edge[, 2], xe[, 2]) + + if (phyloORclado) { + if (is.null(node.pos)) { + node.pos <- 1 + if (type == "cladogram" && !use.edge.length) + node.pos <- 2 + } + if (node.pos == 1) + yy <- .C("node_height", as.integer(Ntip), as.integer(Nnode), + as.integer(x$edge[, 1]), as.integer(x$edge[, + 2]), as.integer(Nedge), as.double(yy), DUP = FALSE, + PACKAGE = "ape")[[6]] + else { + ans <- .C("node_height_clado", as.integer(Ntip), + as.integer(Nnode), as.integer(x$edge[, 1]), as.integer(x$edge[, + 2]), as.integer(Nedge), double(Ntip + Nnode), + as.double(yy), DUP = FALSE, PACKAGE = "ape") + xx <- ans[[6]] - 1 + yy <- ans[[7]] + } + if (!use.edge.length) { + if (node.pos != 2) + xx <- .C("node_depth", as.integer(Ntip), as.integer(Nnode), + as.integer(x$edge[, 1]), as.integer(x$edge[, + 2]), as.integer(Nedge), double(Ntip + Nnode), + DUP = FALSE, PACKAGE = "ape")[[6]] - 1 + xx <- max(xx) - xx + } + else { + xx <- .C("node_depth_edgelength", as.integer(Ntip), + as.integer(Nnode), as.integer(x$edge[, 1]), as.integer(x$edge[, + 2]), as.integer(Nedge), as.double(x$edge.length), + double(Ntip + Nnode), DUP = FALSE, PACKAGE = "ape")[[7]] + } + } + ##if (type == "fan") { + ## TIPS <- xe[which(xe[, 2] <= Ntip), 2] + ## xx <- seq(0, 2 * pi * (1 - 1/Ntip), 2 * pi/Ntip) + ## theta <- double(Ntip) + ## theta[TIPS] <- xx + ## theta <- c(theta, numeric(Nnode)) + ## theta <- .C("node_height", as.integer(Ntip), as.integer(Nnode), + ## as.integer(x$edge[, 1]), as.integer(x$edge[, 2]), + ## as.integer(Nedge), theta, DUP = FALSE, PACKAGE = "ape")[[6]] + ## if (use.edge.length) { + ## r <- .C("node_depth_edgelength", as.integer(Ntip), + ## as.integer(Nnode), as.integer(x$edge[, 1]), as.integer(x$edge[, + ## 2]), as.integer(Nedge), as.double(x$edge.length), + ## double(Ntip + Nnode), DUP = FALSE, PACKAGE = "ape")[[7]] + ## } + ## else { + ## r <- .C("node_depth", as.integer(Ntip), as.integer(Nnode), + ## as.integer(x$edge[, 1]), as.integer(x$edge[, + ## 2]), as.integer(Nedge), double(Ntip + Nnode), + ## DUP = FALSE, PACKAGE = "ape")[[6]] + ## r <- 1/r + ## } + ## xx <- r * cos(theta) + ## yy <- r * sin(theta) + ##} + ##if (type == "unrooted") { + ## XY <- if (use.edge.length) + ## unrooted.xy(Ntip, Nnode, x$edge, x$edge.length) + ## else unrooted.xy(Ntip, Nnode, x$edge, rep(1, Nedge)) + ## xx <- XY$M[, 1] - min(XY$M[, 1]) + ## yy <- XY$M[, 2] - min(XY$M[, 2]) + ##} + ##if (type == "radial") { + ## X <- .C("node_depth", as.integer(Ntip), as.integer(Nnode), + ## as.integer(x$edge[, 1]), as.integer(x$edge[, 2]), + ## as.integer(Nedge), double(Ntip + Nnode), DUP = FALSE, + ## PACKAGE = "ape")[[6]] + ## X[X == 1] <- 0 + ## X <- 1 - X/Ntip + ## yy <- c((1:Ntip) * 2 * pi/Ntip, rep(0, Nnode)) + ## Y <- .C("node_height", as.integer(Ntip), as.integer(Nnode), + ## as.integer(x$edge[, 1]), as.integer(x$edge[, 2]), + ## as.integer(Nedge), as.double(yy), DUP = FALSE, PACKAGE = "ape")[[6]] + ## xx <- X * cos(Y) + ## yy <- X * sin(Y) + ##} + if (phyloORclado && direction != "rightwards") { + if (direction == "leftwards") { + xx <- -xx + xx <- xx - min(xx) + } + if (!horizontal) { + tmp <- yy + yy <- xx + xx <- tmp - min(tmp) + 1 + if (direction == "downwards") { + yy <- -yy + yy <- yy - min(yy) + } + } + } + cbind(xx, yy) +} diff --git a/R/subtreeplot.R b/R/subtreeplot.R new file mode 100644 index 0000000..4e6f81f --- /dev/null +++ b/R/subtreeplot.R @@ -0,0 +1,47 @@ +## subtreeplot.R (2008-04-14) + +## Zoom on a Portion of a Phylogeny by Successive Clicks + +## Copyright 2008 Damien de Vienne + +## This file is part of the R-package `ape'. +## See the file ../COPYING for licensing issues. + +subtreeplot<-function(x, wait=FALSE, ...) { + + sub<-subtrees(x, wait=wait) + y<-NULL + plot.default(0, type="n",axes=FALSE, ann=FALSE) + repeat { + split.screen(c(1,2)) + screen(2) + if (is.null(y)) plot(x,...) + else plot(y,sub=paste("Node :", click),...) + screen(1) + plot(x,sub="Complete tree",main="Type ESC or right click to exit", cex.main=0.9, ...) + + N.tip<-Ntip(x) + N.node<-Nnode(x) + + coor<-plot.phylo.coor(x) + tips<-x$tip.label + nodes<-x$node.label + if (is.null(x$node.label)) nodes<-(N.tip+1):(N.tip+N.node) + labs<-c(rep("",N.tip), nodes) + + click<-identify(coor[,1], coor[,2], labels=labs, n=1) + if (length(click) == 0) {return(y)} + if (click > N.tip) { + close.screen(c(1,2),all.screens = TRUE) + split.screen(c(1,2)) + screen(1) #selects the screen to plot in + plot(x, sub="Complete tree", ...) # plots x in screen 1 (left) + screen(2) + for (i in 1:length(sub)) if (sub[[i]]$name==click) break + y<-sub[[i]] + } + else cat("this is a tip, you have to choose a node\n") + + } + on.exit(return(y)) +} diff --git a/R/subtrees.R b/R/subtrees.R new file mode 100644 index 0000000..bd5bab9 --- /dev/null +++ b/R/subtrees.R @@ -0,0 +1,55 @@ +## subtrees.R (2008-04-14) + +## All subtrees of a Phylogenetic Tree + +## Copyright 2008 Damien de Vienne + +## This file is part of the R-package `ape'. +## See the file ../COPYING for licensing issues. + +subtrees<-function(tree, wait = FALSE) +{ +N.tip<-Ntip (tree) +N.node<-Nnode(tree) +limit<-N.tip+N.node +sub<-list(N.node) +u<-0 + + for (k in (N.tip+1):limit) + { + u<-u+1 + if (wait==TRUE) cat("wait... Node",u,"out of", N.node, "treated\n") + fils<-NULL + pere<-res <- k + repeat + { + for (i in 1: length(pere)) fils<-c(fils, tree$edge[,2][tree$edge[,1]==pere[i]]) + res<-c(res, fils) + pere<-fils + fils<-NULL + if (length(pere)==0) break + } + + len<-res[res>N.tip] + if (u==1) { + tree2<-tree + len<-(N.tip+1):limit + } + else { + len.tip<-res[res