From: paradis Date: Mon, 19 Jul 2010 11:50:40 +0000 (+0000) Subject: BOTHlabels(... hozir = TRUE) X-Git-Url: https://git.donarmstrong.com/?p=ape.git;a=commitdiff_plain;h=82dd3702485179ba5408f1e3e57eb856d025e16c BOTHlabels(... hozir = TRUE) git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@126 6e262413-ae40-0410-9e79-b911bd7a66b7 --- diff --git a/ChangeLog b/ChangeLog index 5e2a0fa..df3b32b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -5,6 +5,13 @@ NEW FEATURES o There is now an 'as.list' method for the class "DNAbin". + o dist.dna() can compute the number of transitions or transversions + with the option model = "Ts" or model = "Tv", respectively. + + o [node|tip|edge]labels() gain three options with default values to + control the aspect of thermometers: horiz = TRUE, width = NULL, + and height = NULL. + CHANGES IN APE VERSION 2.5-3 diff --git a/R/DNA.R b/R/DNA.R index 8b016ca..1f4f5d6 100644 --- a/R/DNA.R +++ b/R/DNA.R @@ -1,4 +1,4 @@ -## DNA.R (2010-07-13) +## DNA.R (2010-07-14) ## Manipulations and Comparisons of DNA Sequences @@ -313,7 +313,7 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE, as.matrix = FALSE) { MODELS <- c("RAW", "JC69", "K80", "F81", "K81", "F84", "T92", "TN93", - "GG95", "LOGDET", "BH87", "PARALIN", "N") + "GG95", "LOGDET", "BH87", "PARALIN", "N", "TS", "TV") imod <- pmatch(toupper(model), MODELS) if (is.na(imod)) stop(paste("'model' must be one of:", @@ -322,7 +322,7 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE, warning("computing variance temporarily not available for model BH87.") variance <- FALSE } - if (gamma && imod %in% c(1, 5:7, 9:12)) { + if (gamma && imod %in% c(1, 5:7, 9:15)) { warning(paste("gamma-correction not available for model", model)) gamma <- FALSE } @@ -331,7 +331,9 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE, n <- dim(x) s <- n[2] n <- n[1] - BF <- if (is.null(base.freq)) base.freq(x) else base.freq + if (imod %in% c(4, 6:8)) { + BF <- if (is.null(base.freq)) base.freq(x) else base.freq + } else BF <- 0 if (!pairwise.deletion) { keep <- .C("GlobalDeletionDNA", x, n, s, rep(1L, s), PACKAGE = "ape")[[4]] diff --git a/R/nodelabels.R b/R/nodelabels.R index 9d5b39c..40bfab1 100644 --- a/R/nodelabels.R +++ b/R/nodelabels.R @@ -1,4 +1,4 @@ -## nodelabels.R (2010-03-12) +## nodelabels.R (2010-07-17) ## Labelling Trees @@ -46,7 +46,7 @@ floating.pie.asp <- function(xpos, ypos, x, edges = 200, radius = 1, } BOTHlabels <- function(text, sel, XX, YY, adj, frame, pch, thermo, - pie, piecol, col, bg, ...) + pie, piecol, col, bg, horiz, width, height, ...) { if (missing(text)) text <- NULL if (length(adj) == 1) adj <- c(adj, 0.5) @@ -84,26 +84,50 @@ BOTHlabels <- function(text, sel, XX, YY, adj, frame, pch, thermo, } if (!is.null(thermo)) { parusr <- par("usr") - width <- CEX * (parusr[2] - parusr[1]) / 40 - height <- CEX * (parusr[4] - parusr[3]) / 15 + + if (is.null(width)) { + width <- CEX * (parusr[2] - parusr[1]) + width <- if (horiz) width/15 else width/40 + } + + if (is.null(height)) { + height <- CEX * (parusr[4] - parusr[3]) + height <- if (horiz) height/40 else height/15 + } + if (is.vector(thermo)) thermo <- cbind(thermo, 1 - thermo) - thermo <- height * thermo + thermo <- if (horiz) width * thermo else height * thermo + if (is.null(piecol)) piecol <- rainbow(ncol(thermo)) + xl <- XX - width/2 + adj[1] - 0.5 # added 'adj' from Janet Young (2009-09-30) xr <- xl + width yb <- YY - height/2 + adj[2] - 0.5 - if (is.null(piecol)) piecol <- rainbow(ncol(thermo)) - ## draw the first rectangle: - rect(xl, yb, xr, yb + thermo[, 1], border = NA, col = piecol[1]) - for (i in 2:ncol(thermo)) - rect(xl, yb + rowSums(thermo[, 1:(i - 1), drop = FALSE]), - xr, yb + rowSums(thermo[, 1:i]), - border = NA, col = piecol[i]) + yt <- yb + height + + if (horiz) { + ## draw the first rectangle: + rect(xl, yb, xl + thermo[, 1], yt, border = NA, col = piecol[1]) + for (i in 2:ncol(thermo)) + rect(xl + rowSums(thermo[, 1:(i - 1), drop = FALSE]), yb, + xl + rowSums(thermo[, 1:i]), yt, border = NA, col = piecol[i]) + } else { + ## draw the first rectangle: + rect(xl, yb, xr, yb + thermo[, 1], border = NA, col = piecol[1]) + for (i in 2:ncol(thermo)) + rect(xl, yb + rowSums(thermo[, 1:(i - 1), drop = FALSE]), + xr, yb + rowSums(thermo[, 1:i]), + border = NA, col = piecol[i]) + } + ## check for NA's before drawing the borders s <- apply(thermo, 1, function(xx) any(is.na(xx))) xl[s] <- xr[s] <- NA - rect(xl, yb, xr, yb + height, border = "black") - segments(xl, YY, xl - width/5, YY) - segments(xr, YY, xr + width/5, YY) + rect(xl, yb, xr, yt, border = "black") + + if (!horiz) { + segments(xl, YY, xl - width/5, YY) + segments(xr, YY, xr + width/5, YY) + } } ## from BB: if (!is.null(pie)) { @@ -122,33 +146,39 @@ BOTHlabels <- function(text, sel, XX, YY, adj, frame, pch, thermo, pch = pch, col = col, bg = bg, ...) } -nodelabels <- function(text, node, adj = c(0.5, 0.5), frame = "rect", - pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, - col = "black", bg = "lightblue", ...) +nodelabels <- + function(text, node, adj = c(0.5, 0.5), frame = "rect", + pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, + col = "black", bg = "lightblue", horiz = FALSE, + width = NULL, height = NULL, ...) { lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv) if (missing(node)) node <- (lastPP$Ntip + 1):length(lastPP$xx) XX <- lastPP$xx[node] YY <- lastPP$yy[node] BOTHlabels(text, node, XX, YY, adj, frame, pch, thermo, - pie, piecol, col, bg, ...) + pie, piecol, col, bg, horiz, width, height, ...) } -tiplabels <- function(text, tip, adj = c(0.5, 0.5), frame = "rect", - pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, - col = "black", bg = "yellow", ...) +tiplabels <- + function(text, tip, adj = c(0.5, 0.5), frame = "rect", + pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, + col = "black", bg = "yellow", horiz = FALSE, + width = NULL, height = NULL, ...) { lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv) if (missing(tip)) tip <- 1:lastPP$Ntip XX <- lastPP$xx[tip] YY <- lastPP$yy[tip] BOTHlabels(text, tip, XX, YY, adj, frame, pch, thermo, - pie, piecol, col, bg, ...) + pie, piecol, col, bg, horiz, width, height, ...) } -edgelabels <- function(text, edge, adj = c(0.5, 0.5), frame = "rect", - pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, - col = "black", bg = "lightgreen", ...) +edgelabels <- + function(text, edge, adj = c(0.5, 0.5), frame = "rect", + pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, + col = "black", bg = "lightgreen", horiz = FALSE, + width = NULL, height = NULL, ...) { lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv) if (missing(edge)) { @@ -171,7 +201,7 @@ edgelabels <- function(text, edge, adj = c(0.5, 0.5), frame = "rect", YY <- (lastPP$yy[subedge[, 1]] + lastPP$yy[subedge[, 2]]) / 2 } BOTHlabels(text, sel, XX, YY, adj, frame, pch, thermo, - pie, piecol, col, bg, ...) + pie, piecol, col, bg, horiz, width, height, ...) } edges <- function(nodes0, nodes1, arrows = 0, type = "classical", ...) diff --git a/inst/doc/MoranI.Rnw b/inst/doc/MoranI.Rnw index fce8ab8..04de862 100644 --- a/inst/doc/MoranI.Rnw +++ b/inst/doc/MoranI.Rnw @@ -81,7 +81,7 @@ calculate the weights $w$: \begin{itemize} \item With phylogenetic distances among species, e.g., $w_{ij} = - 1/d_{ij}$, where $1/d_{ij}$ are distances measured on a tree. + 1/d_{ij}$, where $d_{ij}$ are distances measured on a tree. \item With taxonomic levels where $w_{ij} = 1$ if species $i$ and $j$ belong to the same group, 0 otherwise. \end{itemize} diff --git a/man/corBrownian.Rd b/man/corBrownian.Rd index 46c5753..14e082c 100644 --- a/man/corBrownian.Rd +++ b/man/corBrownian.Rd @@ -4,41 +4,44 @@ \alias{corMatrix.corBrownian} \title{Brownian Correlation Structure} \usage{ - corBrownian(value=1, phy, form=~1) - \method{coef}{corBrownian}(object, unconstrained = TRUE, ...) - \method{corMatrix}{corBrownian}(object, - covariate = getCovariate(object), corr = TRUE, ...) +corBrownian(value=1, phy, form=~1) +\method{coef}{corBrownian}(object, unconstrained = TRUE, ...) +\method{corMatrix}{corBrownian}(object, covariate = getCovariate(object), corr = TRUE, ...) } \arguments{ - \item{value}{The \eqn{\gamma}{gamma} parameter (default to 1)} - \item{phy}{An object of class \code{phylo} representing the phylogeny - (with branch lengths) to consider} - \item{object}{An (initialized) object of class \code{corBrownian}} - \item{corr}{a logical value. If 'TRUE' the function returns the correlation matrix, otherwise it returns - the variance/covariance matrix.} - \item{form}{ignored for now.} - \item{covariate}{ignored for now.} - \item{unconstrained}{a logical value. If 'TRUE' the coefficients are returned + \item{value}{The \eqn{\gamma}{gamma} parameter (default to 1)} + \item{phy}{An object of class \code{phylo} representing the phylogeny + (with branch lengths) to consider} + \item{object}{An (initialized) object of class \code{corBrownian}} + \item{corr}{a logical value. If 'TRUE' the function returns the + correlation matrix, otherwise it returns the variance/covariance matrix.} + \item{form}{ignored for now.} + \item{covariate}{ignored for now.} + \item{unconstrained}{a logical value. If 'TRUE' the coefficients are returned in unconstrained form (the same used in the optimization algorithm). If 'FALSE' the coefficients are returned in "natural", possibly constrained, form. Defaults to 'TRUE'} - \item{\dots}{some methods for these generics require additional arguments. - None are used in these methods.} + \item{\dots}{some methods for these generics require additional arguments. + None are used in these methods.} } \description{ - Expected covariance under a Brownian model (Felsenstein 1985, Martins and Hansen 1997): - \deqn{V_{ij} = \gamma \times t_a}{% - Vij = gamma . ta} - where \eqn{t_a}{ta} is the distance on the phylogeny between the root and the most recent common ancestor - of taxa \eqn{i}{i} and \eqn{j}{j} and \eqn{\gamma}{gamma} is a constant. + Expected covariance under a Brownian model (Felsenstein 1985, Martins + and Hansen 1997) + + \deqn{V_{ij} = \gamma \times t_a}{Vij = gamma . ta} + + where \eqn{t_a}{ta} is the distance on the phylogeny between the root + and the most recent common ancestor of taxa \eqn{i}{i} and \eqn{j}{j} + and \eqn{\gamma}{gamma} is a constant. } \value{ - An object of class \code{corBrownian} or coefficient from an object of this class (actually sends \code{numeric(0)}!) - or the correlation matrix of an initialized object of this class. + An object of class \code{corBrownian}, or the coefficient from an + object of this class (actually sends \code{numeric(0)}), or the + correlation matrix of an initialized object of this class. } \author{Julien Dutheil \email{julien.dutheil@univ-montp2.fr}} \seealso{ - \code{\link{corClasses}}. + \code{\link{corClasses}} } \references{ Felsenstein, J. (1985) Phylogenies and the comparative method. diff --git a/man/nodelabels.Rd b/man/nodelabels.Rd index fa81274..9a7f51a 100644 --- a/man/nodelabels.Rd +++ b/man/nodelabels.Rd @@ -11,13 +11,16 @@ \usage{ nodelabels(text, node, adj = c(0.5, 0.5), frame = "rect", pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, - col = "black", bg = "lightblue", ...) + col = "black", bg = "lightblue", horiz = FALSE, + width = NULL, height = NULL, ...) tiplabels(text, tip, adj = c(0.5, 0.5), frame = "rect", pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, - col = "black", bg = "yellow", ...) + col = "black", bg = "yellow", horiz = FALSE, + width = NULL, height = NULL, ...) edgelabels(text, edge, adj = c(0.5, 0.5), frame = "rect", pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, - col = "black", bg = "lightgreen", ...) + col = "black", bg = "lightgreen", horiz = FALSE, + width = NULL, height = NULL, ...) } \arguments{ @@ -57,6 +60,9 @@ edgelabels(text, edge, adj = c(0.5, 0.5), frame = "rect", \code{points} functions (e.g. \code{cex} to alter the size of the text or the symbols, or \code{font} for the text; see the examples below).} + \item{horiz, width, height}{parameters controlling the aspect of + thermometers; by default, their width and height are determined + automatically.} } \details{ These three functions have the same optional arguments and the same @@ -141,6 +147,10 @@ text(4.5, 15, "Are you \"pie\"...", font = 4, cex = 1.5) plot(tr, "c", FALSE, no.margin = TRUE) nodelabels(thermo = x, col = rainbow(4), cex = 1.3) text(4.5, 15, "... or \"thermo\"?", font = 4, cex = 1.5) +plot(tr, "c", FALSE, no.margin = TRUE) +nodelabels(thermo = x, col = rainbow(4), cex = 1.3) +plot(tr, "c", FALSE, no.margin = TRUE) +nodelabels(thermo = x, col = rainbow(4), width = 3, horiz = TRUE) layout(matrix(1)) plot(tr, main = "Showing Edge Lengths") edgelabels(round(tr$edge.length, 3), srt = 90) diff --git a/src/dist_dna.c b/src/dist_dna.c index 2cabb7d..a289f15 100644 --- a/src/dist_dna.c +++ b/src/dist_dna.c @@ -1,4 +1,4 @@ -/* dist_dna.c 2010-03-29 */ +/* dist_dna.c 2010-07-14 */ /* Copyright 2005-2010 Emmanuel Paradis @@ -64,21 +64,49 @@ double detFourByFour(double *x) if (KnownBase(x[s1]) && KnownBase(x[s2])) L++;\ else continue; -void distDNA_raw(unsigned char *x, int *n, int *s, double *d, int scaled) +#define COUNT_TS_TV\ + if (SameBase(x[s1], x[s2])) continue;\ + Nd++;\ + if (IsPurine(x[s1]) && IsPurine(x[s2])) {\ + Ns++;\ + continue;\ + }\ + if (IsPyrimidine(x[s1]) && IsPyrimidine(x[s2])) Ns++; + +void distDNA_TsTv(unsigned char *x, int *n, int *s, double *d, int Ts, int pairdel) { - int i1, i2, s1, s2, target, Nd; + int i1, i2, s1, s2, target, Nd, Ns; + + 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) { + if (pairdel && !(KnownBase(x[s1]) && KnownBase(x[s2]))) continue; + COUNT_TS_TV + } + if (Ts) d[target] = ((double) Ns); /* output number of transitions */ + else d[target] = ((double) Nd - Ns); /* output number of transversions */ + target++; + } + } +} - target = 0; - for (i1 = 1; i1 < *n; i1++) { - for (i2 = i1 + 1; i2 <= *n; i2++) { - Nd = 0; - for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) - if (DifferentBase(x[s1], x[s2])) Nd++; - if (scaled) d[target] = ((double) Nd / *s); - else d[target] = ((double) Nd); - target++; +void distDNA_raw(unsigned char *x, int *n, int *s, double *d, int scaled) +{ + int i1, i2, s1, s2, target, Nd; + + target = 0; + for (i1 = 1; i1 < *n; i1++) { + for (i2 = i1 + 1; i2 <= *n; i2++) { + Nd = 0; + for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) + if (DifferentBase(x[s1], x[s2])) Nd++; + if (scaled) d[target] = ((double) Nd / *s); + else d[target] = ((double) Nd); + target++; + } } - } } void distDNA_raw_pairdel(unsigned char *x, int *n, int *s, double *d, int scaled) @@ -150,15 +178,6 @@ void distDNA_JC69_pairdel(unsigned char *x, int *n, int *s, double *d, } } -#define COUNT_TS_TV\ - if (SameBase(x[s1], x[s2])) continue;\ - Nd++;\ - if (IsPurine(x[s1]) && IsPurine(x[s2])) {\ - Ns++;\ - continue;\ - }\ - if (IsPyrimidine(x[s1]) && IsPyrimidine(x[s2])) Ns++; - #define COMPUTE_DIST_K80\ P = ((double) Ns/L);\ Q = ((double) (Nd - Ns)/L);\ @@ -1013,39 +1032,33 @@ void dist_dna(unsigned char *x, int *n, int *s, int *model, double *d, switch (*model) { case 1 : if (pairdel) distDNA_raw_pairdel(x, n, s, d, 1); else distDNA_raw(x, n, s, d, 1); break; - case 2 : if (pairdel) distDNA_JC69_pairdel(x, n, s, d, variance, var, gamma, alpha); else distDNA_JC69(x, n, s, d, variance, var, gamma, alpha); break; - case 3 : if (pairdel) distDNA_K80_pairdel(x, n, s, d, variance, var, gamma, alpha); else distDNA_K80(x, n, s, d, variance, var, gamma, alpha); break; - case 4 : if (pairdel) distDNA_F81_pairdel(x, n, s, d, BF, variance, var, gamma, alpha); else distDNA_F81(x, n, s, d, BF, variance, var, gamma, alpha); break; - case 5 : if (pairdel) distDNA_K81_pairdel(x, n, s, d, variance, var); else distDNA_K81(x, n, s, d, variance, var); break; - case 6 : if (pairdel) distDNA_F84_pairdel(x, n, s, d, BF, variance, var); else distDNA_F84(x, n, s, d, BF, variance, var); break; - case 7 : if (pairdel) distDNA_T92_pairdel(x, n, s, d, BF, variance, var); else distDNA_T92(x, n, s, d, BF, variance, var); break; - case 8 : if (pairdel) distDNA_TN93_pairdel(x, n, s, d, BF, variance, var, gamma, alpha); else distDNA_TN93(x, n, s, d, BF, variance, var, gamma, alpha); break; - case 9 : if (pairdel) distDNA_GG95_pairdel(x, n, s, d, variance, var); else distDNA_GG95(x, n, s, d, variance, var); break; - case 10 : if (pairdel) distDNA_LogDet_pairdel(x, n, s, d, variance, var); else distDNA_LogDet(x, n, s, d, variance, var); break; - case 11 : distDNA_BH87(x, n, s, d, variance, var); break; - case 12 : if (pairdel) distDNA_ParaLin_pairdel(x, n, s, d, variance, var); else distDNA_ParaLin(x, n, s, d, variance, var); break; case 13 : if (pairdel) distDNA_raw_pairdel(x, n, s, d, 0); else distDNA_raw(x, n, s, d, 0); break; + case 14 : if (pairdel) distDNA_TsTv(x, n, s, d, 1, 1); + else distDNA_TsTv(x, n, s, d, 1, 0); break; + case 15 : if (pairdel) distDNA_TsTv(x, n, s, d, 0, 1); + else distDNA_TsTv(x, n, s, d, 0, 0); break; + } }