R/mrca.R

   1 ## mrca.R (2009-05-10)
   2
   3 ##   Find Most Recent Common Ancestors Between Pairs
   4
   5 ## Copyright 2005-2009 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 mrca <- function(phy, full = FALSE)
  11 {
  12     if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"')
  13     ##    if (!is.rooted(phy)) stop("the tree must be rooted.")
  14     ## Get all clades:
  15     nb.tip <- length(phy$tip.label)
  16     nb.node <- phy$Nnode
  17     BP <- .Call("bipartition", phy$edge, nb.tip,
  18                 nb.node, PACKAGE = "ape")
  19     N <- nb.tip + nb.node
  20     ROOT <- nb.tip + 1
  21     ## In the following matrix, numeric indexing will be used:
  22     M <- numeric(N * N)
  23     dim(M) <- c(N, N)
  24
  25     ## We start at the root:
  26     next.node <- ROOT
  27     while (length(next.node)) {
  28         tmp <- numeric(0)
  29         for (anc in next.node) {
  30             ## Find the branches which `anc' is the ancestor...:
  31             id <- which(phy$edge[, 1] == anc)
  32             ## ... and get their descendants:
  33             desc <- phy$edge[id, 2]
  34             ## `anc' is itself the MRCA of its direct descendants:
  35             M[anc, desc] <- M[desc, anc] <- anc
  36             ## Find all 2-by-2 combinations of `desc': `anc'
  37             ## is their MRCA:
  38             for (i in 1:length(desc))
  39               M[cbind(desc[i], desc[-i])] <- anc
  40             ## If one element of `desc' is a node, then the tips it
  41             ## leads to and the other elements of `desc' have also
  42             ## `anc' as MRCA!
  43             for (i in 1:length(desc)) {
  44                 if (desc[i] < ROOT) next
  45                 ## (get the tips:)
  46                 tips <- BP[[desc[i] - nb.tip]]
  47                 ## Same thing for the nodes...
  48                 node.desc <- numeric(0)
  49                 for (k in 1:nb.node) {
  50                     if (k == desc[i] - nb.tip) next
  51                     ## If the clade of the current node is a
  52                     ## subset of desc[i], then it is one of its
  53                     ## descendants:
  54                     if (all(BP[[k]] %in% tips))
  55                       node.desc <- c(node.desc, k)
  56                 }
  57                 ## all nodes and tips which are descendants of
  58                 ## `desc[i]':
  59                 ALLDESC <- c(tips, node.desc + nb.tip)
  60                 M[ALLDESC, desc[-i]] <- M[desc[-i], ALLDESC] <- anc
  61                 for (j in 1:length(desc)) {
  62                     if (j == i || desc[j] < ROOT) next
  63                     tips2 <- BP[[desc[j] - nb.tip]]
  64                     node.desc <- numeric(0)
  65                     for (k in 1:nb.node) {
  66                         if (k == desc[j] - nb.tip) next
  67                         if (all(BP[[k]] %in% tips2))
  68                           node.desc <- c(node.desc, k)
  69                     }
  70                     ALLDESC2 <- c(tips2, node.desc + nb.tip)
  71                     M[ALLDESC, ALLDESC2] <- M[ALLDESC2, ALLDESC] <- anc
  72                 }
  73                 ## `anc' is also the MRCA of itself and its descendants:
  74                 M[ALLDESC, anc] <- M[anc, ALLDESC] <- anc
  75             }
  76             ## When it is done, `desc' i stored to become
  77             ## the new `next.node', if they are nodes:
  78             tmp <- c(tmp, desc[desc > nb.tip])
  79         }
  80         next.node <- tmp
  81     }
  82     M[cbind(1:N, 1:N)] <- 1:N
  83     if (full)
  84       dimnames(M)[1:2] <- list(as.character(1:N))
  85     else {
  86         M <- M[1:nb.tip, 1:nb.tip]
  87         dimnames(M)[1:2] <- list(phy$tip.label)
  88     }
  89     M
  90 }