From a03a8c554a6fde0dc4313688e3248bfae2e521e4 Mon Sep 17 00:00:00 2001 From: paradis Date: Fri, 19 Mar 2010 16:51:48 +0000 Subject: [PATCH] fixing various bugs + news in cophyloplot git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@115 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 14 ++++++ DESCRIPTION | 2 +- R/DNA.R | 14 +++++- R/cophyloplot.R | 106 +++++++++++++++++++++------------------------ R/plot.phylo.R | 10 +++-- man/bind.tree.Rd | 3 +- man/cophyloplot.Rd | 43 +++++++++++------- src/dist_dna.c | 69 +++++++++++++++++++++++++++-- 8 files changed, 178 insertions(+), 83 deletions(-) diff --git a/ChangeLog b/ChangeLog index c551ded..1088c13 100644 --- a/ChangeLog +++ b/ChangeLog @@ -11,6 +11,9 @@ NEW FEATURES o drop.tip(), extract.clade(), root(), and bind.tree() now have an 'interactive' option to make the operation on a plotted tree. + o cophyloplot() gains two new arguments 'lwd' and 'lty' for the + association links; they are recycled like 'col' (which wasn't before). + BUG FIXES @@ -29,6 +32,17 @@ BUG FIXES o plot.phylo(type = "p") sometimes failed to colour correctly the vertical lines representing the nodes. + o plot.phylo(direction = "l", x.lim = 30) failed to plot the branches + in the correct direction though the tip labels were displayed + correctly. + + +OTHER CHANGES + + o The c, cbind, and rbind methods for "DNAbin" objetcs now check that + the sequences are correctly stored (in a list for c, in a matrix + for the two other functions). + CHANGES IN APE VERSION 2.5 diff --git a/DESCRIPTION b/DESCRIPTION index c23ac10..981ea9f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape Version: 2.5-1 -Date: 2010-03-15 +Date: 2010-03-19 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, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer, Damien de Vienne Maintainer: Emmanuel Paradis diff --git a/R/DNA.R b/R/DNA.R index c751c50..e7c0e3a 100644 --- a/R/DNA.R +++ b/R/DNA.R @@ -1,8 +1,8 @@ -## DNA.R (2009-10-02) +## DNA.R (2010-03-16) ## Manipulations and Comparisons of DNA Sequences -## Copyright 2002-2009 Emmanuel Paradis +## Copyright 2002-2010 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -90,6 +90,9 @@ rbind.DNAbin <- function(...) obj <- list(...) n <- length(obj) if (n == 1) return(obj[[1]]) + for (i in 1:n) + if (!is.matrix(obj[[1]])) + stop("the 'rbind' method for \"DNAbin\" accepts only matrices") NC <- unlist(lapply(obj, ncol)) if (length(unique(NC)) > 1) stop("matrices do not have the same number of columns.") @@ -106,6 +109,9 @@ cbind.DNAbin <- obj <- list(...) n <- length(obj) if (n == 1) return(obj[[1]]) + for (i in 1:n) + if (!is.matrix(obj[[1]])) + stop("the 'cbind' method for \"DNAbin\" accepts only matrices") NR <- unlist(lapply(obj, nrow)) for (i in 1:n) class(obj[[i]]) <- NULL if (check.names) { @@ -144,7 +150,11 @@ cbind.DNAbin <- } c.DNAbin <- function(..., recursive = FALSE) +{ + if (!all(unlist(lapply(list(...), is.list)))) + stop("the 'c' method for \"DNAbin\" accepts only lists") structure(NextMethod("c"), class = "DNAbin") +} print.DNAbin <- function(x, ...) { diff --git a/R/cophyloplot.R b/R/cophyloplot.R index 3bb3702..04425c6 100644 --- a/R/cophyloplot.R +++ b/R/cophyloplot.R @@ -1,18 +1,18 @@ -## cophyloplot.R (2008-07-08) +## cophyloplot.R (2010-03-18) ## Plots two phylogenetic trees face to ## face with the links between the tips -## Copyright 2008 Damien de Vienne +## Copyright 2008-2010 Damien de Vienne ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. cophyloplot <- - 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, ...) + function(x, y, assoc = NULL, use.edge.length = FALSE, space = 0, + length.line = 1, gap = 2, type = "phylogram", rotate = FALSE, + col = par("fg"), lwd = par("lwd"), lty = par("lty"), + show.tip.label = TRUE, font = 3, ...) { if (is.null(assoc)) { assoc <- matrix(ncol = 2) @@ -23,7 +23,7 @@ cophyloplot <- repeat { res <- plotCophylo2(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, + type = type, return = TRUE, col = col, lwd=lwd, lty=lty, 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) { @@ -36,75 +36,62 @@ cophyloplot <- } plotCophylo2(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) + type = type, return = TRUE, col = col, lwd=lwd, lty=lty, show.tip.label = show.tip.label, font = font) } on.exit(print("done")) } else plotCophylo2(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) + type = type, return = FALSE, col = col, lwd=lwd, lty=lty, show.tip.label = show.tip.label, font = font) } plotCophylo2 <- 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, ...) + type = type, return = return, col = col, lwd=lwd, lty=lty, + 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 - + 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 <- plotPhyloCoor(x, use.edge.length = use.edge.length, - type = type) + a <- plotPhyloCoor(x, use.edge.length = use.edge.length, type = type) res$a <- a b <- plotPhyloCoor(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]) + 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]) + 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], col="red") 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)) { @@ -135,39 +122,46 @@ plotCophylo2 <- } } if (show.tip.label) { - text(a[1:N.tip.x, ], cex = 0, font = font, pos = 4, - labels = x$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)) + + #colors + if (length(col)==1) colors<-c(rep(col, nrow(assoc))) + else if (length(col)>=nrow(assoc)) colors<-col + else colors<-c(rep(col, as.integer(nrow(assoc)/length(col))+1)) + + #lwd + if (length(lwd)==1) lwidths<-c(rep(lwd, nrow(assoc))) + else if (length(lwd)>=nrow(assoc)) lwidths<-lwd + else lwidths<-c(rep(lwd, as.integer(nrow(assoc)/length(lwd))+1)) + + #lty + if (length(lty) == 1) ltype <- c(rep(lty, nrow(assoc))) + else if (length(lty) >= nrow(assoc)) ltype <- lty + else ltype <- c(rep(lty, as.integer(nrow(assoc)/length(lty))+1)) + for (i in 1:nrow(assoc)) { - if (show.tip.label) { + + 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, 2]]]) - } else { + } + 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) + + 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 = colors[i], lwd=lwidths[i], lty=ltype[i]) + 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 = colors[i], lwd=lwidths[i], lty=ltype[i]) + 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 = colors[i], lwd=lwidths[i], lty=ltype[i]) } - if (return == TRUE) - return(res) + if (return == TRUE) return(res) } diff --git a/R/plot.phylo.R b/R/plot.phylo.R index de7278e..631fb4d 100644 --- a/R/plot.phylo.R +++ b/R/plot.phylo.R @@ -1,4 +1,4 @@ -## plot.phylo.R (2010-03-15) +## plot.phylo.R (2010-03-19) ## Plot Phylogenies @@ -198,8 +198,6 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, tmp <- if (show.tip.label) max(xx.tips + strWi/alp) else max(xx.tips) } x.lim[2] <- tmp - if (direction == "leftwards") xx <- x.lim[2] - xx #max(xx[ROOT] + tmp) -# else max(xx[1:Ntip] + tmp) } else x.lim <- c(1, Ntip) } else switch(type, "fan" = { if (show.tip.label) { @@ -227,6 +225,8 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, if (show.tip.label) -1 - max(nchar(x$tip.label) * 0.03 * cex) else -1 } + ## mirror the xx: + if (phyloORclado && direction == "leftwards") xx <- x.lim[2] - xx if (is.null(y.lim)) { if (phyloORclado) { if (horizontal) y.lim <- c(1, Ntip) else { @@ -245,7 +245,6 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, tmp <- if (show.tip.label) max(yy.tips + strWi/alp) else max(yy.tips) } y.lim[2] <- tmp - if (direction == "downwards") yy <- y.lim[2] - yy } } else switch(type, "fan" = { if (show.tip.label) { @@ -271,6 +270,8 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, if (type == "radial") y.lim[1] <- if (show.tip.label) -1 - max(nchar(x$tip.label) * 0.018 * max(yy) * cex) else -1 } + ## mirror the yy: + if (phyloORclado && direction == "downwards") yy <- y.lim[2] - yy if (phyloORclado && root.edge) { if (direction == "leftwards") x.lim[2] <- x.lim[2] + x$root.edge if (direction == "downwards") y.lim[2] <- y.lim[2] + x$root.edge @@ -332,6 +333,7 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, "upwards" = segments(xx[ROOT], 0, xx[ROOT], x$root.edge), "downwards" = segments(xx[ROOT], yy[ROOT], xx[ROOT], yy[ROOT] + x$root.edge)) if (show.tip.label) { + if (is.expression(x$tip.label)) underscore <- TRUE if (!underscore) x$tip.label <- gsub("_", " ", x$tip.label) if (phyloORclado) diff --git a/man/bind.tree.Rd b/man/bind.tree.Rd index ff8cafe..58464e0 100644 --- a/man/bind.tree.Rd +++ b/man/bind.tree.Rd @@ -1,9 +1,10 @@ \name{bind.tree} \alias{bind.tree} +\alias{+.phylo} \title{Binds Trees} \usage{ bind.tree(x, y, where = "root", position = 0, interactive = FALSE) -\method{+}{phylo}(x, y) +\special{x + y} } \arguments{ \item{x}{an object of class \code{"phylo"}.} diff --git a/man/cophyloplot.Rd b/man/cophyloplot.Rd index 431cbb5..4f125dc 100644 --- a/man/cophyloplot.Rd +++ b/man/cophyloplot.Rd @@ -5,25 +5,38 @@ This function plots two trees face to face with the links if specified. It is possible to rotate the branches of each tree around the nodes by clicking. } \usage{ -cophyloplot(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, \dots) +cophyloplot(x, y, assoc = NULL, use.edge.length = FALSE, space = 0, + length.line = 1, gap = 2, type = "phylogram", rotate = FALSE, + col = par("fg"), lwd = par("lwd"), lty = par("lty"), + show.tip.label = TRUE, font = 3, \dots) } \arguments{ \item{x, y}{two objects of class \code{"phylo"}.} - \item{assoc}{a matrix with 2 columns specifying the associations between the tips. If NULL, no links will be drawn.} - \item{use.edge.length}{a logical indicating whether the branch lengths should be used to plot the trees; default is FALSE.} - \item{space}{a positive value that specifies the distance between the two trees.} - \item{length.line}{a positive value that specifies the length of the horizontal line associated to each taxa. Default is 1.} - \item{gap}{a value specifying the distance between the tips of the phylogeny and the lines.} - \item{type}{a character string specifying the type of phylogeny to be drawn; it must be one of "phylogram" (the default) or "cladogram".} - \item{rotate}{a logical indicating whether the nodes of the phylogeny can be rotated by clicking. Default is FALSE.} - \item{col}{a character string indicating the color to be used for the links. Default is red.} - \item{show.tip.label}{a logical indicating whether to show the tip labels on the phylogeny (defaults to 'TRUE', i.e. the labels are shown).} - \item{font}{an integer specifying the type of font for the labels: 1 - (plain text), 2 (bold), 3 (italic, the default), or 4 (bold - italic).} + \item{assoc}{a matrix with 2 columns specifying the associations + between the tips. If NULL, no links will be drawn.} + \item{use.edge.length}{a logical indicating whether the branch lengths + should be used to plot the trees; default is FALSE.} + \item{space}{a positive value that specifies the distance between the + two trees.} + \item{length.line}{a positive value that specifies the length of the + horizontal line associated to each taxa. Default is 1.} + \item{gap}{a value specifying the distance between the tips of the + phylogeny and the lines.} + \item{type}{a character string specifying the type of phylogeny to be + drawn; it must be one of "phylogram" (the default) or "cladogram".} + \item{rotate}{a logical indicating whether the nodes of the phylogeny + can be rotated by clicking. Default is FALSE.} + \item{col}{a character vector indicating the color to be used for the + links; recycled as necessary.} + \item{lwd}{id. for the width.} + \item{lty}{id. for the line type.} + \item{show.tip.label}{a logical indicating whether to show the tip + labels on the phylogeny (defaults to 'TRUE', i.e. the labels are + shown).} + \item{font}{an integer specifying the type of font for the + labels: 1 (plain text), 2 (bold), 3 (italic, the default), or 4 + (bold italic).} \item{\dots}{(unused)} } \details{ diff --git a/src/dist_dna.c b/src/dist_dna.c index 2aa6da6..9057b9b 100644 --- a/src/dist_dna.c +++ b/src/dist_dna.c @@ -1,6 +1,6 @@ -/* dist_dna.c 2009-09-16 */ +/* dist_dna.c 2010-03-16 */ -/* Copyright 2005-2008 Emmanuel Paradis +/* Copyright 2005-2010 Emmanuel Paradis /* This file is part of the R-package `ape'. */ /* See the file ../COPYING for licensing issues. */ @@ -350,14 +350,14 @@ void distDNA_K81_pairdel(unsigned char *x, int *n, int *s, double *d, #define COMPUTE_DIST_F84\ P = ((double) Ns/L);\ Q = ((double) (Nd - Ns)/L);\ - d[target] = -2*A*log(1 - (P/(2*A) - (A - B)*Q/(2*A*C))) + 2*(A - B - C)*log(1 - Q/(2*C));\ + d[target] = -2*A*log(1 - P/(2*A) - (A - B)*Q/(2*A*C)) + 2*(A - B - C)*log(1 - Q/(2*C));\ if (*variance) {\ t1 = A*C;\ t2 = C*P/2;\ t3 = (A - B)*Q/2;\ a = t1/(t1 - t2 - t3);\ b = A*(A - B)/(t1 - t2 - t3) - (A - B - C)/(C - Q/2);\ - var[target] = (a*a*P + b*b*Q - pow(a*P + b*Q, 2))/2;\ + var[target] = (a*a*P + b*b*Q - pow(a*P + b*Q, 2))/L;\ } void distDNA_F84(unsigned char *x, int *n, int *s, double *d, @@ -800,6 +800,67 @@ void distDNA_BH87(unsigned char *x, int *n, int *s, double *d, for (i2 = i1 + 1; i2 <= *n; i2++) { for (k = 0; k < 16; k++) Ntab[k] = 0; for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) { +#define PREPARE_BF_F84\ + A = (BF[0]*BF[2])/(BF[0] + BF[2]) + (BF[1]*BF[3])/(BF[1] + BF[3]);\ + B = BF[0]*BF[2] + BF[1]*BF[3];\ + C = (BF[0] + BF[2])*(BF[1] + BF[3]); + +#define COMPUTE_DIST_F84\ + P = ((double) Ns/L);\ + Q = ((double) (Nd - Ns)/L);\ + d[target] = -2*A*log(1 - P/(2*A) - (A - B)*Q/(2*A*C)) + 2*(A - B - C)*log(1 - Q/(2*C));\ + if (*variance) {\ + t1 = A*C;\ + t2 = C*P/2;\ + t3 = (A - B)*Q/2;\ + a = t1/(t1 - t2 - t3);\ + b = A*(A - B)/(t1 - t2 - t3) - (A - B - C)/(C - Q/2);\ + var[target] = (a*a*P + b*b*Q - pow(a*P + b*Q, 2))/2;\ + } + +void distDNA_F84(unsigned char *x, int *n, int *s, double *d, + double *BF, int *variance, double *var) +{ + int i1, i2, Nd, Ns, L, target, s1, s2; + double P, Q, A, B, C, a, b, t1, t2, t3; + + PREPARE_BF_F84 + L = *s; + + target = 0; + for (i1 = 1; i1 < *n; i1++) { + for (i2 = i1 + 1; i2 <= *n; i2++) { + Nd = Ns = 0; + for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) { + COUNT_TS_TV + } + COMPUTE_DIST_F84 + target++; + } + } +} + +void distDNA_F84_pairdel(unsigned char *x, int *n, int *s, double *d, + double *BF, int *variance, double *var) +{ + int i1, i2, Nd, Ns, L, target, s1, s2; + double P, Q, A, B, C, a, b, t1, t2, t3; + + PREPARE_BF_F84 + + target = 0; + for (i1 = 1; i1 < *n; i1++) { + for (i2 = i1 + 1; i2 <= *n; i2++) { + Nd = Ns = L = 0; + for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) { + CHECK_PAIRWISE_DELETION + COUNT_TS_TV + } + COMPUTE_DIST_F84 + target++; + } + } +} DO_CONTINGENCY_NUCLEOTIDES } -- 2.39.2