R/cophenetic.phylo.R

   1 ## cophenetic.phylo.R (2010-11-15)
   2
   3 ##   Pairwise Distances from a Phylogenetic Tree
   4
   5 ## Copyright 2006-2010 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 dist.nodes <- function(x)
  11 {
  12     if (is.null(x$edge.length))
  13         stop("your tree has no branch lengths")
  14
  15     trim <- FALSE
  16     if (!is.binary.tree(x) || !is.rooted(x)) {
  17         trim <- TRUE
  18         x <- makeNodeLabel(x)
  19         x <- multi2di(x, random = FALSE)
  20     }
  21
  22     n <- length(x$tip.label)
  23     n.node <- x$Nnode
  24     N <- n + n.node
  25     x <- reorder(x, order = "pruningwise")
  26
  27     res <- matrix(NA, N, N)
  28     res[cbind(1:N, 1:N)] <- 0 # implicit mode conversion
  29
  30     ## I like the simplicity of this one:
  31     res[x$edge] <- res[x$edge[, 2:1]] <- x$edge.length
  32
  33     ## compute the distances ...
  34     for (i in seq(from = 1, by = 2, length.out = n.node)) {
  35         j <- i + 1
  36         anc <- x$edge[i, 1]
  37         des1 <- x$edge[i, 2]
  38         des2 <- x$edge[j, 2]
  39
  40         ## If `des1' is a node, we look for the nodes and tips for
  41         ## which the distance up to `des1' has already been
  42         ## computed, including `des1' itself. For all these, we can
  43         ## compute the distance up to `anc' and all node(s) and
  44         ## tip(s) in `des2'.
  45         if (des1 > n) des1 <- which(!is.na(res[des1, ]))
  46
  47         ## id. for `des2'
  48         if (des2 > n) des2 <- which(!is.na(res[des2, ]))
  49
  50         ## The following expression is vectorized only on `des2' and
  51         ## not on `des1' because they may be of different lengths.
  52         for (y in des1)
  53           res[y, des2] <- res[des2, y] <- res[anc, y] + res[anc, des2]
  54         ## compute the distances between the tip(s) and node(s)
  55         ## in `des2' and the ancestor of `anc'; id. for `des2'
  56         ## (only if it is not the root)
  57         if (anc != n + 1) {
  58             ind <- which(x$edge[, 2] == anc)
  59             nod <- x$edge[ind, 1] # the ancestor of `anc'
  60             l <- x$edge.length[ind]
  61             res[des2, nod] <- res[nod, des2] <- res[anc, des2] + l
  62             res[des1, nod] <- res[nod, des1] <- res[anc, des1] + l
  63         }
  64     }
  65     if (trim) {
  66         i <- which(x$node.label == "") + n
  67         res <- res[-i, -i]
  68         N <- dim(res)[1]
  69     }
  70     dimnames(res)[1:2] <- list(1:N)
  71     res
  72 }
  73
  74 cophenetic.phylo <- function(x)
  75 {
  76     n <- length(x$tip.label)
  77     ans <- dist.nodes(x)[1:n, 1:n]
  78     dimnames(ans)[1:2] <- list(x$tip.label)
  79     ans
  80 }