## multi2di.R (2010-01-23)
## Collapse and Resolve Multichotomies
## Copyright 2005-2010 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
multi2di <- function(phy, random = TRUE)
{
degree <- tabulate(phy$edge[, 1])
target <- which(degree > 2)
if (!length(target)) return(phy)
nb.edge <- dim(phy$edge)[1]
n <- length(phy$tip.label)
nextnode <- n + phy$Nnode + 1
new.edge <- edge2delete <- NULL
wbl <- FALSE
if (!is.null(phy$edge.length)) {
wbl <- TRUE
new.edge.length <- NULL
}
for (node in target) {
ind <- which(phy$edge[, 1] == node)
N <- length(ind)
desc <- phy$edge[ind, 2]
if (random) {
## if we shuffle the descendants, we need to eventually
## reorder the corresponding branch lenghts (see below)
## so we store the result of sample()
tmp <- sample(length(desc))
desc <- desc[tmp]
res <- rtree(N)$edge
} else {
res <- matrix(0, 2*N - 2, 2)
res[, 1] <- N + rep(1:(N - 1), each = 2)
res[, 2] <- N + rep(2:N, each = 2)
res[seq(1, by = 2, length.out = N - 1), 2] <- 1:(N - 1)
res[length(res)] <- N
}
if (wbl) {
## keep the branch lengths coming from `node'
el <- numeric(dim(res)[1]) # initialized with 0's
el[res[, 2] <= N] <-
if (random) phy$edge.length[ind][tmp] else phy$edge.length[ind]
}
## now substitute the nodes in `res'
## `node' stays at the "root" of these new
## edges whereas their "tips" are `desc'
Nodes <- c(node, seq(from = nextnode, length.out = N - 2))
res[, 1] <- Nodes[res[, 1] - N]
tmp <- res[, 2] > N
res[tmp, 2] <- Nodes[res[tmp, 2] - N]
res[!tmp, 2] <- desc[res[!tmp, 2]]
new.edge <- rbind(new.edge, res)
edge2delete <- c(edge2delete, ind)
if (wbl) new.edge.length <- c(new.edge.length, el)
nextnode <- nextnode + N - 2
phy$Nnode <- phy$Nnode + N - 2
}
phy$edge <- rbind(phy$edge[-edge2delete, ], new.edge)
if (wbl)
phy$edge.length <- c(phy$edge.length[-edge2delete], new.edge.length)
if (!is.null(attr(phy, "order"))) attr(phy, "order") <- NULL
if (!is.null(phy$node.label))
phy$node.label <-
c(phy$node.label, rep("", phy$Nnode - length(phy$node.label)))
phy <- reorder(phy)
## the node numbers are not in increasing order in edge[, 2]: this
## will confuse drop.tip and other functions (root), so renumber them
newNb <- integer(phy$Nnode)
newNb[1] <- n + 1L
sndcol <- phy$edge[, 2] > n
## reorder node labels before changing edge:
if (!is.null(phy$node.label)) {
o <- 1 + rank(phy$edge[sndcol, 2])
## the root's label is not changed:
phy$node.label <- phy$node.label[c(1, o)]
}
## executed from right to left, so newNb is modified before phy$edge:
phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2] - n] <-
n + 2:phy$Nnode
phy$edge[, 1] <- newNb[phy$edge[, 1] - n]
phy
}
di2multi <- function(phy, tol = 1e-8)
{
if (is.null(phy$edge.length)) stop("the tree has no branch length")
## We select only the internal branches which are
## significantly small:
ind <- which(phy$edge.length < tol & phy$edge[, 2] > length(phy$tip.label))
n <- length(ind)
if (!n) return(phy)
## recursive function to `propagate' node #'s in case
## there is a series of consecutive edges to remove
foo <- function(ancestor, des2del) {
wh <- which(phy$edge[, 1] == des2del)
for (k in wh) {
if (phy$edge[k, 2] %in% node2del) foo(ancestor, phy$edge[k, 2])
else phy$edge[k, 1] <<- ancestor
}
}
node2del <- phy$edge[ind, 2]
anc <- phy$edge[ind, 1]
for (i in 1:n) {
if (anc[i] %in% node2del) next
foo(anc[i], node2del[i])
}
phy$edge <- phy$edge[-ind, ]
phy$edge.length <- phy$edge.length[-ind]
phy$Nnode <- phy$Nnode - n
## Now we renumber the nodes that need to be:
sel <- phy$edge > min(node2del)
for (i in which(sel))
phy$edge[i] <- phy$edge[i] - sum(node2del < phy$edge[i])
if (!is.null(phy$node.label))
phy$node.label <- phy$node.label[-(node2del - length(phy$tip.label))]
phy
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.