#' Visualise a matching
#'
#' Depict the splits that are matched between two trees using a specified
#' [Generalized Robinson–Foulds](
#' https://ms609.github.io/TreeDist/articles/Generalized-RF.html)
#' similarity measure.
#'
#' Note that when visualizing a Robinson–Foulds distance (using
#' `Func = RobinsonFouldsMatching`), matched splits are assigned a _similarity_
#' score of 1, which is deducted from the total number of splits to calculate
#' the Robinson–Foulds _distance_. Unmatched splits thus contribute one to
#' total tree distance.
#'
#' @param Func Function used to construct tree similarity.
#' @param tree1,tree2 Trees of class `phylo`, with identical leaf labels.
#' @param setPar Logical specifying whether graphical parameters should be
#' set to display trees side by side.
#' @param precision Integer specifying number of significant figures to display
#' when reporting matching scores.
#' @param Plot Function to use to plot trees.
#' @param matchZeros Logical specifying whether to pair splits with zero
#' similarity (`TRUE`), or leave them unpaired (`FALSE`).
#' @param plainEdges Logical specifying whether to plot edges with a uniform
#' width and colour (`TRUE`), or whether to draw edge widths according to the
#' similarity of the associated splits (`FALSE`).
#' @param edge.width,edge.color,\dots Additional parameters to send to `Plot()`.
#'
#' @importFrom ape nodelabels edgelabels plot.phylo
#' @importFrom colorspace qualitative_hcl sequential_hcl
#' @importFrom graphics par
#' @importFrom TreeTools as.Splits
#'
#' @examples
#' tree1 <- TreeTools::BalancedTree(6)
#' tree2 <- TreeTools::PectinateTree(6)
#'
#' VisualizeMatching(RobinsonFouldsMatching, tree1, tree2)
#' VisualizeMatching(SharedPhylogeneticInfo, tree1, tree2, matchZeros = FALSE)
#' @template MRS
#' @encoding UTF-8
#' @export
VisualizeMatching <- function(Func, tree1, tree2, setPar = TRUE,
precision = 3L, Plot = plot.phylo,
matchZeros = TRUE, plainEdges = FALSE,
edge.width = 1, edge.color = "black",
...) {
splits1 <- as.Splits(tree1)
edge1 <- tree1[["edge"]]
child1 <- edge1[, 2]
nTip <- attr(splits1, "nTip")
splitEdges1 <- vapply(as.integer(rownames(splits1)),
function(node) which(child1 == node), integer(1))
splits2 <- as.Splits(tree2, tipLabels = tree1)
edge2 <- tree2[["edge"]]
child2 <- edge2[, 2]
splitEdges2 <- vapply(as.integer(rownames(splits2)),
function(node) which(child2 == node), integer(1))
matching <- Func(tree1, tree2, reportMatching = TRUE)
pairings <- attr(matching, "matching")
scores <- attr(matching, "pairScores")
pairScores <- signif(mapply(function(i, j) scores[i, j],
seq_along(pairings), pairings), precision)
adjNo <- c(0.5, -0.2)
adjVal <- c(0.5, 1.1)
faint <- "#aaaaaa"
if (setPar) {
origPar <- par(mfrow = c(1, 2), mar = rep(0.5, 4))
on.exit(par(origPar))
}
LabelUnpaired <- function(splitEdges, unpaired) {
if (any(unpaired)) {
#edgelabels(text="\u2012", edge=splitEdges[unpaired],
edgelabels(text = expression("-"), edge = splitEdges[unpaired],
frame = "n", col = faint, adj = adjNo)
edgelabels(text = "0", edge = splitEdges[unpaired],
frame = "n", col = faint, cex = 0.8, adj = adjVal)
}
}
if (plainEdges) {
Plot(tree1, edge.width = edge.width, edge.color = edge.color, ...)
} else {
Normalize <- function(scores, na.rm = FALSE) {
if (length(scores) == 0) return(scores)
if (na.rm) {
scores <- scores[!is.na(scores)]
} else {
scores[is.na(scores)] <- 0
}
if (any(scores < 0)) {
stop("Negative scores not supported") # nocov
}
if (max(scores) == 0) return(scores)
if (min(scores) == max(scores)) return(rep(1L, length(scores)))
scores / max(scores)
}
OtherRootEdge <- function(splitNodes, edge) {
parent <- edge[, 1]
child <- edge[, 2]
rootEdges <- which(parent == min(parent))
rootChildren <- child[rootEdges]
treeIsRooted <- length(rootChildren) < 3
if (treeIsRooted) {
splitEdges <- vapply(splitNodes, match, table = child, 0)
got <- rootChildren %in% splitNodes
if (any(got)) {
if(sum(got) != 1) {
warning("Unexpected polytomy")
}
c(score = as.integer(which(splitNodes %in% rootChildren[got])),
edge = rootEdges[!got])
} else {
# `edge` is not a root edge
c(score = NA_integer_, edge = NA_integer_)
}
} else {
# Tree is unrooted => there is no root edge at all
c(score = NA_integer_, edge = NA_integer_)
}
}
edgeColPalette <- sequential_hcl(n = 256L, palette = "Viridis")
EdgyPlot <- function(tree, splits, edge, splitEdges,
normalizedScores, ...) {
splitNodes <- as.integer(names(splits))
ore <- OtherRootEdge(splitNodes, edge)
if (length(normalizedScores) && !is.na(ore[[1]])) {
ns <- c(normalizedScores, normalizedScores[ore["score"]])
se <- c(splitEdges, ore[[2]])
} else {
ns <- normalizedScores
se <- splitEdges
}
edge.width <- rep(1, nrow(edge))
edge.width[se] <- 1 + (10 * ns)
edge.color <- rep("black", nrow(edge))
edge.color[se] <- edgeColPalette[1 + ceiling(255 * ns)]
Plot(tree, edge.width = edge.width, edge.color = edge.color, ...)
}
EdgyPlot(tree1, splits1, edge1, splitEdges1, Normalize(pairScores), ...)
}
paired1 <- if (matchZeros) {
!is.na(pairings)
} else {
!is.na(pairings) & pairScores > 0
}
palette <- qualitative_hcl(sum(paired1), c = 42, l = 88)
pairedPairScores <- pairScores[paired1]
pairLabels <- seq_len(sum(paired1))
if (any(pairLabels)) {
edgelabels(text = pairLabels, edge = splitEdges1[paired1],
bg = palette, adj = adjNo)
edgelabels(text = pairedPairScores, edge = splitEdges1[paired1],
frame = "n", adj = adjVal, cex = 0.8,
col = ifelse(pairedPairScores, "black", faint))
}
LabelUnpaired(splitEdges1, !paired1)
paired2 <- seq_along(splitEdges2) %in% pairings[paired1]
pairNames2 <- pairings[paired1]
if (plainEdges) {
Plot(tree2, edge.width = edge.width, edge.color = edge.color, ...)
} else {
EdgyPlot(tree2, splits2[[pairNames2]], edge2, splitEdges2[pairNames2],
Normalize(pairedPairScores, na.rm = TRUE), ...)
}
if (any(pairLabels)) {
edgelabels(text = pairLabels, edge = splitEdges2[pairNames2],
bg = palette, adj=adjNo)
edgelabels(text = pairedPairScores, edge = splitEdges2[pairNames2],
frame = "n", adj = adjVal, cex = 0.8,
col = ifelse(pairedPairScores, "black", faint))
}
LabelUnpaired(splitEdges2, !paired2)
# Return:
invisible()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.