Nothing
#' Find tip position differences
#'
#' Find the topologicial differences between two trees with the same tip labels. The function returns a data frame of the tips and the number of differences in their ancestry between the two trees.
#' Called by \code{\link{plotTreeDiff}}, which highlights the differing tips in a plot of the two trees.
#'
#' @author Michelle Kendall \email{michelle.louise.kendall@@gmail.com}
#'
#' @import ape
#' @importFrom combinat combn2
#'
#' @param tr1 an object of the class \code{phylo}: the first tree to compare.
#' @param tr2 an object of the class \code{phylo}: the second tree to compare.
#' @param vec1 an optional input, the result of \code{treeVec(tr1, lambda=0)}, to speed up the computation.
#' @param vec2 an optional input, the result of \code{treeVec(tr2, lambda=0)}, to speed up the computation.
#'
#' @return
#' A data frame of the tree tips and the number of ancestral differences between them in the two trees, in order of increasing difference.
#' A tip is said to have zero difference if each of its ancestral nodes admits the same tip partition in the two trees.
#'
#' @seealso \code{\link{medTree}} \code{\link{plotTreeDiff}}
#'
#' @examples
#' ## simple example on trees with five tips:
#' tr1 <- read.tree(text="((A:1,B:1):1,((C:1,D:1):1,E:1):1):1;")
#' tr2 <- read.tree(text="((A:1,B:1):1,(C:1,(D:1,E:1):1):1):1;")
#' tipDiff(tr1,tr2)
#'
#' ## example on larger woodmice trees
#' data(woodmiceTrees)
#' tipDiff(woodmiceTrees[[1]],woodmiceTrees[[2]])
#'
#' @export
tipDiff <- function(tr1,tr2,vec1=NULL,vec2=NULL) {
l <- length(tr1$tip.label)
lchoose2 <- l*(l-1)/2
if( l != length(tr2$tip.label)) stop("Trees must have the same number of tips")
if(setequal(tr1$tip.label,tr2$tip.label) == FALSE) stop("Trees must have the same tip label sets")
# get vec1
if (is.null(vec1)) {
vec1 <- treeVec(tr1) # emphasise.tips, emphasise.weight)
}
if (is.null(vec2)) {
vec2 <- treeVec(tr2) # emphasise.tips, emphasise.weight)
}
# trim pendant edge entries from vectors
vec1 <- vec1[1:lchoose2] # emphasise.tips, emphasise.weight)
vec2 <- vec2[1:lchoose2]
# find the positions where the vectors are different
vecDiff <- (vec1!=vec2)
# combine the tip labels (in order) and whether the vectors are different
treedf <- as.data.frame(cbind(combn2(tr1$tip.label[order(tr1$tip.label)]),vecDiff))
# list the tips which appear in "TRUE" rows
allTipDiffs <- c(as.character(treedf[,1][vecDiff]),as.character(treedf[,2][vecDiff]))
# find the number of times each tip appears
tipSignificance <- sapply(tr1$tip.label, function(x)
length(which(allTipDiffs==x)))
# prepare output as data frame of tips and their differences
tipDiff <- as.data.frame(cbind(names(tipSignificance[order(tipSignificance)]),tipSignificance[order(tipSignificance)]), stringsAsFactors = FALSE)
class(tipDiff[,2]) <- "numeric"
rownames(tipDiff) <- NULL
colnames(tipDiff) <- c("Tip","No. of differences")
return(tipDiff)
}
#' Plot tree differences
#'
#' Highlight the topologicial differences between two trees, plotted side by side.
#' This function is useful for comparing representative "median" trees - see \code{\link{medTree}}.
#' It relies on the function \code{\link{tipDiff}}.
#'
#' @author Michelle Kendall \email{michelle.louise.kendall@@gmail.com}
#'
#' @import ape
#' @importFrom adegenet lightseasun
#' @importFrom adegenet num2col
#' @importFrom grDevices colorRampPalette
#' @importFrom graphics layout
#'
#' @param tr1 an object of the class \code{phylo}: the first tree to plot.
#' @param tr2 an object of the class \code{phylo}: the second tree to plot.
#' @param tipDiff an optional input, the result of \code{\link{tipDiff}}. Supplying this will save time if calling \code{plotTreeDiff} repeatedly, for example with different aesthetics.
#' @param vec1 an optional input, the result of \code{treeVec(tr1, lambda=0)}. This argument is ignored if \code{tipDiff} is supplied; otherwise supplying this will save time if calling \code{plotTreeDiff} repeatedly, for example with different aesthetics.
#' @param vec2 an optional input, the result of \code{treeVec(tr2, lambda=0)}. This argument is ignored if \code{tipDiff} is supplied; otherwise supplying this will save time if calling \code{plotTreeDiff} repeatedly, for example with different aesthetics.
#' @param baseCol the colour used for tips with identical ancestry in the two trees. Defaults to "grey".
#' @param colourMethod the method to use for colouring. Default is "ramp", corresponding to the original implementation, where the function \code{colorRampPalette} is used to create a palette which ranges from \code{col1} to \code{col2}. For large trees this can be hard to interpret, and method \code{palette} may be preferred, which permits the selection of a palette to use in \code{adegenet}'s function \code{num2col}.
#' @param col1 the first colour used to define the colour spectrum for tips with differences. This colour will be used for tips with minor differences. Defaults to "peachpuff". Ignored if \code{colourMethod="palette"}.
#' @param col2 the second colour used to define the colour spectrum for tips with differences. This colour will be used for tips with major differences. Defaults to "red2". Ignored if \code{colourMethod="palette"}.
#' @param palette the colour palette to be used if \code{colourMethod="palette"}. For a list of available palettes see \code{?num2col}.
#' @param ... further arguments passed to \code{\link{plot.phylo}}
#'
#' @return
#' A plot of the two trees side by side. Tips are coloured in the following way:
#' \itemize{
#' \item if each ancestor of a tip in tree 1 occurs in tree 2 with the same partition of tip descendants, then the tip is coloured grey (or supplied "baseCol")
#' \item if not, the tip gets coloured pale orange to red on a scale according to how many differences there are amongst its most recent common ancestors with other tips. The colour spectrum can be changed according to preference.
#' }
#'
#' @seealso \code{\link{medTree}}, \code{\link{tipDiff}}
#'
#' @examples
#' ## simple example on trees with five tips:
#' tr1 <- read.tree(text="((A:1,B:1):1,((C:1,D:1):1,E:1):1):1;")
#' tr2 <- read.tree(text="((A:1,B:1):1,(C:1,(D:1,E:1):1):1):1;")
#' plotTreeDiff(tr1,tr2)
#'
#' ## example on larger woodmice trees
#' data(woodmiceTrees)
#' # find the tip differences in advance, to avoid recalculating with each plot
#' wmTipDiff <- tipDiff(woodmiceTrees[[1]],woodmiceTrees[[2]])
#' plotTreeDiff(woodmiceTrees[[1]],woodmiceTrees[[2]], tipDiff=wmTipDiff)
#' ## change aesthetics:
#' plotTreeDiff(woodmiceTrees[[1]],woodmiceTrees[[2]], tipDiff=wmTipDiff,
#' baseCol="grey2", col1="cyan", col2="navy",
#' edge.width=2, type="radial", cex=0.5, font=2)
#' ## use colour palette from adegenet:
#' plotTreeDiff(woodmiceTrees[[1]],woodmiceTrees[[2]], tipDiff=wmTipDiff,
#' baseCol="black", colourMethod="palette",
#' edge.width=2, type="cladogram", cex=0.5, font=2)
#'
#' @export
plotTreeDiff <- function(tr1,tr2,tipDiff=NULL,vec1=NULL,vec2=NULL,baseCol="grey",col1="peachpuff",col2="red2",colourMethod="ramp",palette=lightseasun,...) {
l <- length(tr1$tip.label)
if (is.null(tipDiff)){
# call tipDiff
tipDiff <- tipDiff(tr1,tr2,vec1,vec2)
}
# find the number of times each tip appears, in the order that the tip labels are read
tipSignificance1 <- sapply(tr1$tip.label, function(x)
tipDiff[which(tipDiff[,1]==x),2])
tipSignificance2 <- sapply(tr2$tip.label, function(x)
tipDiff[which(tipDiff[,1]==x),2])
if (colourMethod=="ramp") {
colfunc <- colorRampPalette(c(col1,col2))
if (min(tipDiff[,2])==0) { # make sure tips with no differences are coloured baseCol
tipSignificance1 <- tipSignificance1 + 1
tipSignificance2 <- tipSignificance2 + 1
if (max(tipDiff[,2])==0) {numCols <- 0}
else {numCols <- max(tipDiff[,2]) - min(tipDiff[,2][which(tipDiff[,2]!=0)]) + 1}
pal <- c(baseCol,colfunc(numCols))
}
else {
numCols <- max(tipDiff[,2]) - min(tipDiff[,2]) + 1
pal <- colfunc(numCols)
}
tipCols1 <- pal[as.factor(tipSignificance1)]
tipCols2 <- pal[as.factor(tipSignificance2)]
}
else {
tipCols1 <- num2col(tipSignificance1, col.pal=palette)
tipCols2 <- num2col(tipSignificance2, col.pal=palette)
if (min(tipDiff[,2])==0) { # make sure tips with no differences are coloured baseCol
tipCols1[which(tipSignificance1==0)] <- baseCol
tipCols2[which(tipSignificance2==0)] <- baseCol
}
}
# plot
layout(matrix(1:2, 1, 2))
plot.phylo(tr1, tip.color=tipCols1, no.margin=TRUE, ...)
plot.phylo(tr2, tip.color=tipCols2, no.margin=TRUE, ...)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.