R/root.R

   1 ## root.R (2007-12-21)
   2
   3 ##   Root of Phylogenetic Trees
   4
   5 ## Copyright 2004-2007 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 is.rooted <- function(phy)
  11 {
  12     if (!"phylo" %in% class(phy))
  13       stop('object "phy" is not of class "phylo"')
  14     if (!is.null(phy$root.edge)) return(TRUE)
  15     else
  16       if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2)
  17         return(FALSE)
  18       else return(TRUE)
  19 }
  20
  21 unroot <- function(phy)
  22 {
  23     if (class(phy) != "phylo")
  24       stop('object "phy" is not of class "phylo"')
  25     if (dim(phy$edge)[1] < 3)
  26       stop("cannot unroot a tree with two edges.")
  27     ## delete FIRST the root.edge (in case this is sufficient to
  28     ## unroot the tree, i.e. there is a multichotomy at the root)
  29     if (!is.null(phy$root.edge)) phy$root.edge <- NULL
  30     if (!is.rooted(phy)) return(phy)
  31     ## We remove one of the edges coming from the root, and
  32     ## eventually adding the branch length to the other one
  33     ## also coming from the root.
  34     ## In all cases, the node deleted is the 2nd one (numbered
  35     ## nb.tip+2 in `edge'), so we simply need to renumber the
  36     ## nodes by adding 1, except the root (this remains the
  37     ## origin of the tree).
  38     nb.tip <- length(phy$tip.label)
  39     ROOT <- nb.tip + 1
  40     EDGEROOT <- which(phy$edge[, 1] == ROOT)
  41     ## j: the target where to stick the edge
  42     ## i: the edge to delete
  43     if (phy$edge[EDGEROOT[1], 2] == ROOT + 1) {
  44         j <- EDGEROOT[2]
  45         i <- EDGEROOT[1]
  46     } else {
  47         j <- EDGEROOT[1]
  48         i <- EDGEROOT[2]
  49     }
  50     ## This should work whether the tree is in pruningwise or
  51     ## cladewise order.
  52     phy$edge <- phy$edge[-i, ]
  53     nodes <- phy$edge > ROOT # renumber all nodes except the root
  54     phy$edge[nodes] <- phy$edge[nodes] - 1
  55     if (!is.null(phy$edge.length)) {
  56         phy$edge.length[j] <- phy$edge.length[j] + phy$edge.length[i]
  57         phy$edge.length <- phy$edge.length[-i]
  58     }
  59     phy$Nnode <- phy$Nnode - 1
  60     if (!is.null(phy$node.label))
  61       phy$node.label <- phy$node.label[-2]
  62     phy
  63 }
  64
  65 root <- function(phy, outgroup)
  66 {
  67     if (class(phy) != "phylo")
  68       stop('object "phy" is not of class "phylo"')
  69     if (is.character(outgroup))
  70       outgroup <- which(phy$tip.label %in% outgroup)
  71     nb.tip <- length(phy$tip.label)
  72     if (length(outgroup) == nb.tip) return(phy)
  73
  74     ## First check that the outgroup is monophyletic--
  75     ## unless there's only one tip specified of course
  76     ROOT <- nb.tip + 1
  77     if (length(outgroup) > 1) {
  78         msg <- "the specified outgroup is not monophyletic!"
  79         ## If all tips in `tip' are not contiguous, then
  80         ## no need to go further:
  81         if (!all(diff(outgroup) == 1)) stop(msg)
  82         seq.nod <- .Call("seq_root2tip", phy$edge, nb.tip,
  83                          phy$Nnode, PACKAGE = "ape")
  84         sn <- seq.nod[outgroup]
  85         ## We go from the root to the tips: the sequence of nodes
  86         ## is identical until the MRCA:
  87         newroot <- ROOT
  88         i <- 2 # we start at the 2nd position since the root
  89                # of the tree is a common ancestor to all tips
  90         repeat {
  91             x <- unique(unlist(lapply(sn, "[", i)))
  92             if (length(x) != 1) break
  93             newroot <- x
  94             i <- i + 1
  95         }
  96         ## Check that all descendants of this node
  97         ## are included in the outgroup
  98         ## (1st solution... there may be something smarter)
  99         desc <- which(unlist(lapply(seq.nod,
 100                                     function(x) any(x %in% newroot))))
 101         if (length(outgroup) != length(desc)) stop(msg)
 102         if (!all(sort(outgroup) == sort(desc))) stop(msg)
 103
 104     } else newroot <- phy$edge[which(phy$edge[, 2] == outgroup), 1]
 105
 106     if (newroot == ROOT) return(phy)
 107
 108 ### <FIXME>
 109 ### The remaining part of the code has not been improved; this
 110 ### does not seem obvious. This is delayed...     (2006-09-23)
 111 ### </FIXME>
 112
 113     ## Invert all branches from the new root to the old one
 114     i <- which(phy$edge[, 2] == newroot)
 115     nod <- phy$edge[i, 1]
 116     while (nod != ROOT) {
 117         j <- which(phy$edge[, 2] == nod)
 118         phy$edge[i, 1] <- phy$edge[i, 2]
 119         phy$edge[i, 2] <- nod
 120         i <- j
 121         nod <- phy$edge[i, 1]
 122     }
 123
 124     i.oroot <- which(phy$edge[, 1] == ROOT)
 125     ## Unroot the tree if there's a basal dichotomy...
 126     if (length(i.oroot) == 2) {
 127         j <- i.oroot[which(i.oroot != i)]
 128         phy$edge[j, 1] <- phy$edge[i, 2]
 129         phy$edge <- phy$edge[-i, ]
 130         if (!is.null(phy$edge.length)) {
 131             phy$edge.length[j] <- phy$edge.length[j] + phy$edge.length[i]
 132             phy$edge.length <- phy$edge.length[-i]
 133         }
 134         phy$edge[which(phy$edge == newroot)] <- ROOT
 135     } else {
 136         ## ... otherwise just invert the root with the newroot
 137         phy$edge[which(phy$edge == newroot)] <- ROOT
 138         phy$edge[i.oroot] <- newroot
 139         ## ... and invert finally! (fixed 2005-11-07)
 140         phy$edge[i, ] <- rev(phy$edge[i, ])
 141     }
 142     if (!is.null(phy$node.label))
 143       ## It's important to not delete the label of the newroot
 144       ## to keep the positions of the other nodes
 145       phy$node.label[1] <- phy$node.label[newroot - nb.tip]
 146     ## Not needed: phy$Nnode <- phy$Nnode - 1
 147     read.tree(text = write.tree(phy))
 148 }