In ms609/TreeDistData: Analysis of Phylogenetic Tree Distance Measures

library('Quartet', exclude = 'RobinsonFoulds')
library('TreeTools', quietly = TRUE, warn.conflicts = FALSE)
library('TreeDist', warn.conflicts = FALSE)
library('TreeDistData')

origPar <- par()

data('pectinateDistances11', package = 'TreeDistData')
pd11 <- pectinateDistances11
abbrevs <- TreeDistData::tdAbbrevs
pectinateTree <- Quartet::sq_trees$caterpillar

pm1l <- ape::read.tree(text = 
                       '(((((2, 3), 4), 5), 6), (7, (8, (9, (10, (11, 1))))));')
pm1l$edge.lengths <- c(rep(1, 4L), 2, rep(1, 12), .5, .5, 1.5)
threeMoves <- ape::read.tree(text = 
                       '(((((11, 2), 3), 4), 5), (7, (6, (8, (9, (10, 1))))));')
threeMoves$edge.lengths <- rep(1, 20)
threeMoves$edge.lengths[1] <- 2
threeMoves$edge.lengths[c(12, 14)] <- 0.5

topTail <- Quartet::sq_trees$top_and_tail

pd11['qd', ] <- pd11['qd', ] * 330L

BoxPlot <- function (method) {
  distances <- pd11[method, ]
  faintCol <- TreeDistCol(method, '44')
  opaqueCol <- TreeDistCol(method)
  if (method %in% c('rf', 'spr', 'mast', 'masti', 'mafi', 'tbr_u', 'nni_t',
                    'nni_u')) {
    histy <- hist(distances - 1e-7, plot = FALSE,
                  breaks = min(80L, length(unique(distances))))
    histMax <- max(histy$counts)
    xMax <- switch(method,
                   'mafi' = LnUnrooted(11) / log(2),
                   'mast' = 11,
                   'masti' = LnUnrooted(11) / log(2),
                   max(histy$breaks))
    xLim <- if(method %in% c('mast', 'masti')) c(xMax, 0) else c(0, xMax)
    plot(histy, axes = FALSE, main = NA, xlab = NA, ylab = NA, cex = 0.8,
         col = faintCol, border = NA, xlim = xLim)
  } else {
    # KDE not appropriate when not gamma distributed.
#    dists <- sample(distances, min(1000, length(distances)))
#    kde <- kdensity::kdensity(dists, tolerance = 5, kernel = 'gamma')
#    nBins <- 1024
#    binMiddle <- max(distances) / nBins / 2
#    kRange <- seq(0, max(distances), length.out = nBins)
#    histMax <- max(kde(kRange))
#    plot(kRange, kde(kRange), col = faintCol, type = 'n',
#         axes = FALSE, ylab = NA, xlab = NA,
#         xlim = c(0, max(distances, TDPair[[method]](pectinateTree, topTail))))
#    rect(kRange - binMiddle, 0, kRange + binMiddle, kde(kRange), col = faintCol,
#         border = NA)
    histy <- hist(distances - 1e-7, plot = FALSE,
                      breaks = min(512L, length(unique(distances))))
    histMax <- max(histy$counts)
    xLim <- c(0, max(histy$breaks, TDPair[[method]](pectinateTree, topTail)))
    plot(histy, axes = FALSE, main = NA, xlab = NA, ylab = NA, cex = 0.8,
         col = faintCol, border = NA, xlim = xLim)
  }
  yAt <- histMax / 2L
  boxplot(distances, border = opaqueCol, horizontal = TRUE,
          axes = FALSE, add = TRUE, at = yAt, outline = FALSE, 
          boxwex = histMax * 0.6, col = '#ffffff00', lwd = 2)
  legend('left', bty = 'n', legend = tdAbbrevs[method])
  points(TDPair[[method]](pectinateTree, pm1l), histMax * 0.8,
         pch = 25L, cex = 1.5, col = 'black', bg = opaqueCol)
  points(TDPair[[method]](pectinateTree, threeMoves), histMax * 0.5,
         pch = 3L, cex = 1.5, col = 'black', bg = opaqueCol)
  points(TDPair[[method]](pectinateTree, topTail), histMax * 0.2,
         pch = 24L, cex = 1.5, col = 'black', bg = opaqueCol)
  axis(1, col = opaqueCol)
  invisible()
}

TreePlotsPlotFunc <- function() {
  par(mar = c(1, 0.2, 1, 1), xpd = NA)
  nMethod <- length(tdPlotSequence)
  nRow <- nMethod / 2
  layout(matrix(c(rep(1:4, each = nRow),
                  4 + rep(seq_len(nMethod), each = 4)), ncol = 3),
         widths = c(1, 2, 2))

  par(mar = c(0.2, 0, 1, 1))
  TreeDistPlot(pectinateTree, '', leaveRoom = FALSE, 
               cex = 1, xpd = NA)
  legend('bottomleft', legend = "Pectinate reference\ntree", bty = 'n',
         inset = c(-0.1, -0.04))

  TreeDistPlot(pm1l, '', bold = 1L, leaveRoom = FALSE, cex = 1, graft = 20)
  legend('bottomleft', legend = "Move one leaf (A)", pch = 25, bty = 'n',
         inset = c(-0, -0.02), pt.cex = 1.5)

  TreeDistPlot(topTail, '', leaveRoom = FALSE, cex = 1,
               bold = c(1, 11), graft = c(6, 20))
  legend('bottomleft', legend = "Exchange two\nleaves (A, K)", pch = 24,
         bty = 'n', inset = c(-0.0, -0.005), pt.cex = 1.5)

  TreeDistPlot(threeMoves, '', leaveRoom = FALSE, cex = 1,
               bold = c(1, 6, 11), graft = c(5, 13, 20))
  legend('bottomleft', legend = "Three leaves\nmoved (A, K, F)", pch = 3, 
         bty = 'n', inset = c(-0.0, -0.005), pt.cex = 1.5)

  par(mar = c(1.9, 0.2, 0.2, 0.2))
  lapply(tdPlotSequence, BoxPlot) -> XX
}

PlotFunc <- function() {
  nMethod <- length(tdPlotSequence)
  nCol <- 3L
  nRow <- ceiling(nMethod / nCol)
  par(mfrow = c(nRow, nCol),
      mar = c(1, 0.2, 1, 1),
      xpd = NA,
      cex = 0.7,
      mar = c(1.9, 0.2, 0.2, 0.2)
  )

  lapply(tdPlotSequence, BoxPlot) -> XX
}

PlotFunc()

suppressWarnings(par(origPar))

Figure 1. Tree differences produced by moving one leaf (▽), exchanging two leaves (△), and moving a third leaf (+) in a pectinate tree, compared to the distances between the pectinate tree and 100 000 trees drawn at random from the uniform distribution of binary trees on the same leaves.

ms609/TreeDistData documentation built on May 21, 2021, 6:53 a.m.