3 ## Root of Phylogenetic Trees
5 ## Copyright 2004-2007 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 is.rooted <- function(phy)
12 if (!"phylo" %in% class(phy))
13 stop('object "phy" is not of class "phylo"')
14 if (!is.null(phy$root.edge)) return(TRUE)
16 if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2)
21 unroot <- function(phy)
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)
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) {
50 ## This should work whether the tree is in pruningwise or
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]
59 phy$Nnode <- phy$Nnode - 1
60 if (!is.null(phy$node.label))
61 phy$node.label <- phy$node.label[-2]
65 root <- function(phy, outgroup)
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)
74 ## First check that the outgroup is monophyletic--
75 ## unless there's only one tip specified of course
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:
88 i <- 2 # we start at the 2nd position since the root
89 # of the tree is a common ancestor to all tips
91 x <- unique(unlist(lapply(sn, "[", i)))
92 if (length(x) != 1) break
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)
104 } else newroot <- phy$edge[which(phy$edge[, 2] == outgroup), 1]
106 if (newroot == ROOT) return(phy)
109 ### The remaining part of the code has not been improved; this
110 ### does not seem obvious. This is delayed... (2006-09-23)
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
121 nod <- phy$edge[i, 1]
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]
134 phy$edge[which(phy$edge == newroot)] <- ROOT
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, ])
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))