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
-## DNA.R (2010-07-13)
+## DNA.R (2010-07-14)
## Manipulations and Comparisons of DNA Sequences
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:",
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
}
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]]
-## nodelabels.R (2010-03-12)
+## nodelabels.R (2010-07-17)
## Labelling Trees
}
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)
}
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)) {
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)) {
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", ...)
\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}
\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.
\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{
\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
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)
-/* dist_dna.c 2010-03-29 */
+/* dist_dna.c 2010-07-14 */
/* Copyright 2005-2010 Emmanuel Paradis
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)
}
}
-#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);\
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;
+
}
}
projects / ape.git / commitdiff