Nothing
## evonet.R (2017-07-28)
## Evolutionary Networks
## Copyright 2011-2012 Emmanuel Paradis, 2017 Klaus Schliep
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
evonet <- function(phy, from, to = NULL)
{
if (!inherits(phy, "phylo"))
stop('object "phy" is not of class "phylo".')
if (!is.rooted(phy))
warning("the tree is unrooted")
x <- phy
if (is.null(to)) {
if (is.data.frame(from))
from <- as.matrix(from)
if (!is.matrix(from))
stop("'from' must be a matrix or a data frame if 'to' is not given")
if (ncol(from) > 2) {
warning("'from' has more than two columns: only the first two will be used.")
ret <- from[, 1:2]
} else if (ncol(from) < 2) {
stop("'from' must have at least two columns")
} else ret <- from
} else {
from <- as.vector(from)
to <- as.vector(to)
if (length(from) != length(to))
stop("'from' and 'to' not of the same length after coercing as vectors")
ret <- cbind(from, to)
}
## check that values are not out of range:
storage.mode(ret) <- "integer"
if (any(is.na(ret)))
stop("some values are NA's after coercing as integers")
if (any(ret < 0) || any(ret > Ntip(phy) + phy$Nnode))
stop("some values are out of range")
x$reticulation <- ret
class(x) <- c("evonet", "phylo")
x
}
as.phylo.evonet <- function(x, ...)
{
x$reticulation <- NULL
class(x) <- "phylo"
x
}
plot.evonet <- function(x, col = "blue", lty = 1, lwd = 1, alpha = 0.5,
arrows = 0, arrow.type = "classical", ...)
{
## changed 5/24/17 by Klaus
plot.phylo(x, ...)
edges(x$reticulation[, 1], x$reticulation[, 2],
col = rgb(t(col2rgb(col)), alpha = 255 * alpha,
maxColorValue = 255),
lty = lty, lwd = lwd, arrows = arrows, type = arrow.type)
}
as.networx.evonet <- function(x, weight = NA, ...)
{
if (any(x$reticulation <= Ntip(x)))
stop("some tips are involved in reticulations: cannot convert to \"networx\"")
x <- reorder(x, "postorder")
ned <- Nedge(x)
nrt <- nrow(x$reticulation)
x$edge <- rbind(x$edge, x$reticulation)
colnames(x$edge) <- c("oldNodes", "newNodes")
x$reticulation <- NULL
x$edge.length <- c(x$edge.length, rep(weight, length.out = nrt))
x$split <- c(1:ned, 1:nrt)
class(x) <- c("networx", "phylo")
x
}
as.network.evonet <- function(x, directed = TRUE, ...)
{
class(x) <- NULL
x$edge <- rbind(x$edge, x$reticulation)
as.network.phylo(x, directed = directed, ...)
}
as.igraph.evonet <- function(x, directed = TRUE, use.labels = TRUE, ...)
{
class(x) <- NULL
x$edge <- rbind(x$edge, x$reticulation)
## added check by Klaus (2017-05-26)
if (use.labels) {
if (!is.null(x$node.label)){
tmp <- nchar(x$node.label)
if (any(tmp == 0)){
newLabel <- paste0("number", 1:x$Nnode)
x$node.label[tmp == 0] <- newLabel[tmp == 0]
}
}
if (any(duplicated(c(x$tip.label, x$node.label))))
stop("Duplicated labels!")
}
as.igraph.phylo(x, directed = directed, use.labels = use.labels, ...)
}
print.evonet <- function(x, ...)
{
nr <- nrow(x$reticulation)
cat("\n Evolutionary network with", nr, "reticulation")
if (nr > 1) cat("s")
cat("\n\n --- Base tree ---")
print.phylo(as.phylo(x))
}
## new stuff by Klaus (2017-05-26)
reorder.evonet <- function(x, order = "cladewise", index.only = FALSE, ...)
{
reticulation <- x$reticulation
y <- reorder(as.phylo(x), order = order, index.only = index.only, ...)
if (index.only) return(y)
y$reticulation <- reticulation
class(y) <- c("evonet", "phylo")
y
}
as.evonet <- function(x, ...)
{
if (inherits(x, "evonet")) return(x)
UseMethod("as.evonet")
}
as.evonet.phylo <- function(x, ...)
{
pos <- grep("#", x$tip.label)
ind <- match(pos, x$edge[, 2])
reticulation <- x$edge[ind, , drop = FALSE]
edge <- x$edge[-ind, , drop = FALSE]
nTips <- as.integer(length(x$tip.label))
reticulation[, 2] <- as.integer(match(x$tip.label[pos], x$node.label) + nTips)
for (i in sort(pos, TRUE)) {
edge[edge > i ] <- edge[edge > i] - 1L
reticulation[reticulation > i] <- reticulation[reticulation > i] - 1L
}
x$tip.label <- x$tip.label[-pos]
nTips <- as.integer(length(x$tip.label))
nn <- length(unique(edge[,1]))
if(nn < x$Nnode){
ne <- as.integer( x$Nnode - nn )
edge[edge > nTips] <- edge[edge > nTips] + ne
reticulation[reticulation > nTips] <- reticulation[reticulation > nTips] +ne
z <- logical(max(edge))
z[edge[, 2]] <- TRUE
z[seq_len(nTips)] <- FALSE
z[edge[, 1]] <- FALSE
pos2 <- which(z)
k <- 1
for (i in sort(pos2, TRUE)) {
nTips <- as.integer( nTips + 1L )
edge[edge==i] <- nTips
reticulation[reticulation == i] <- nTips
edge[edge > i] <- edge[edge > i] - 1L
reticulation[reticulation > i] <- reticulation[reticulation > i] - 1L
}
x$Nnode <- nn
x$node.label <- NULL
x$tip.label <- c(x$tip.label , rep("", ne))
}
x$edge <- edge
x$reticulation <- reticulation
if (!is.null(x$edge.length)) x$edge.length <- x$edge.length[-ind]
class(x) <- c("evonet", "phylo")
x
}
## requires new version of clado.build and tree.build
read.evonet <- function(file = "", text = NULL, comment.char = "", ...)
{
x <- read.tree(file = file, text = text, comment.char = comment.char, ...)
as.evonet.phylo(x)
}
.evonet2phylo <- function(x)
{
nTips <- as.integer(length(x$tip.label))
if (!is.null(x$edge.length)) {
nd <- node.depth.edgelength(x)
x$edge.length <- c(x$edge.length, nd[x$reticulation[, 2]] - nd[x$reticulation[, 1]])
}
if (!is.null(x$node.label))
x$tip.label <- c(x$tip.label, x$node.label[x$reticulation[, 2] - nTips])
else {
newLabels <- paste0("#H", x$reticulation[, 2])
x$tip.label <- c(x$tip.label, newLabels)
x$node.label <- rep("", x$Nnode)
ind <- which((x$reticulation[, 2] > nTips) & !duplicated(x$reticulation[, 2]))
x$node.label[x$reticulation[ind, 2] - nTips] <- newLabels[ind]
}
nrets <- as.integer(nrow(x$reticulation))
x$edge[x$edge > nTips] <- x$edge[x$edge > nTips] + nrets
x$reticulation[, 1] <- x$reticulation[, 1] + nrets
x$reticulation[, 2] <- nTips + (1L:nrets)
x$edge <- rbind(x$edge, x$reticulation)
x$reticulation <- NULL
attr(x, "order") <- NULL
class(x) <- "phylo"
x
}
write.evonet <- function(x, file = "", ...)
{
x <- .evonet2phylo(x)
write.tree(x, file = file, ...)
}
Nedge.evonet <- function(phy) dim(phy$edge)[1] + dim(phy$reticulation)[1]
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.