-## 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)
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) {
}
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)) {
}
}
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)
}
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{
-/* 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. */
#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,
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
}
projects / ape.git / commitdiff