Nothing
### Compute the height of the rooted tree.
get.rooted.tree.height <- function(rooted.tree, tol = .Machine$double.eps^0.5){
if(!is.rooted(rooted.tree)){
stop("A rooted tree is required")
}
if(rooted.tree$Nnode > 1){
if(!is.ultrametric(rooted.tree, tol = tol, option = 2)){
stop("A rooted tree is not ultrametric")
}
}
edge.from <- rooted.tree$edge[, 1]
edge.to <- rooted.tree$edge[, 2]
edge.id <- which(edge.from == (length(rooted.tree$tip.label) + 1))
n.from <- edge.to[edge.id[1]]
ret <- rooted.tree$edge.length[edge.id[1]]
repeat{
edge.id <- which(edge.from == n.from)
if(length(edge.id) >= 2){
n.from <- edge.to[edge.id[1]]
ret <- ret + rooted.tree$edge.length[edge.id[1]]
} else{
break
}
}
ret
} # End of get.rooted.tree.height().
### Rescale the rooted tree branche lengths by a constant.
rescale.rooted.tree <- function(rooted.tree, scale.height = 1){
if(!is.rooted(rooted.tree)){
stop("A rooted tree is required")
}
if(scale.height <= 0){
stop("The scale.height is not correct.")
}
rooted.tree$edge.length <- rooted.tree$edge.length * scale.height
rooted.tree
} # End of rescale.rooted.tree().
### Reshap a rooted tree to a star shape.
as.star.tree <- function(rooted.tree, keep.bifurcation = TRUE){
if(!is.rooted(rooted.tree)){
stop("A rooted tree is required")
}
n.tip <- length(rooted.tree$tip.label)
tree.height <- get.rooted.tree.height(rooted.tree)
if(keep.bifurcation){
rooted.tree$edge.length[rooted.tree$edge[, 2] > n.tip] <- 0
rooted.tree$edge.length[rooted.tree$edge[, 2] <= n.tip] <- tree.height
} else{
rooted.tree$edge <- rooted.tree$edge[rooted.tree$edge[, 2] <= n.tip,]
rooted.tree$edge[, 1] <- max(rooted.tree$edge[, 2]) + 1
rooted.tree$Nnode <- 1
rooted.tree$edge.length <- rep(tree.height, n.tip)
}
rooted.tree
} # End of as.star.tree().
### Generate a star tree.
gen.star.tree <- function(N, total.height = 1){
if(N <= 1){
stop("N > 1.")
}
ms.star <- ms(N, 1, opts = paste("-T", sep = " "))
tree.star <- read.tree(text = ms.star[3])
tree.star$tip.label <- paste("s", 1:N, sep = "")
tree.star <- as.star.tree(tree.star)
th <- get.rooted.tree.height(tree.star)
ret <- rescale.rooted.tree(tree.star, total.height / th)
ret
} # End of gen.star.tree().
### Generate a tree with unit heights and K clusters.
gen.unit.K <- function(K, N.K, rate.anc = 10, rate.dec = 10){
if(K <= 1){
stop("K > 1.")
}
if(any(N.K <= 1)){
stop("All N.K > 1.")
}
if(length(N.K) != K){
stop("The length of N.K is not equal to K.")
}
if(rate.anc <= 0 || rate.dec <= 0){
stop("The rate.anc or rate.dec <= 0.")
}
ms.anc <- ms(K, 1, opts = paste("-T -G", rate.anc, sep = " "))
tree.anc <- read.tree(text = ms.anc[3])
tree.anc$tip.label <- paste("a", 1:K, sep = "")
tree.dec <- NULL
tree.org <- tree.anc
tree.star <- tree.anc
for(k in 1:K){
ms.dec <- ms(N.K[k], 1, opts = paste("-T -G", rate.dec, sep = " "))
tree.dec[[k]] <- read.tree(text = ms.dec[3])
tree.dec[[k]]$tip.label <- paste("d", k, ".", 1:N.K[k], sep = "")
tree.org <- bind.tree(tree.org, tree.dec[[k]],
where = which(tree.org$tip.label ==
paste("a", k, sep = "")))
tree.star <- bind.tree(tree.star, as.star.tree(tree.dec[[k]]),
where = which(tree.star$tip.label ==
paste("a", k, sep = "")))
}
tree.anc.height <- get.rooted.tree.height(tree.anc)
tree.dec.height <- unlist(lapply(tree.dec, get.rooted.tree.height))
tree.dec.mean.height <- mean(tree.dec.height)
total.height <- tree.anc.height + tree.dec.mean.height
max.height <- tree.anc.height + max(tree.dec.height)
tree.equal <- rescale.rooted.tree(tree.anc, 1 / total.height)
tmp.scale <- tree.dec.mean.height / total.height
for(k in 1:K){
tree.dec.equal <- rescale.rooted.tree(tree.dec[[k]],
tmp.scale / tree.dec.height[k])
tree.equal <- bind.tree(tree.equal, tree.dec.equal,
where = which(tree.equal$tip.label ==
paste("a", k, sep = "")))
}
tree.max <- rescale.rooted.tree(tree.org, 1 / max.height)
tree.star <- rescale.rooted.tree(tree.star, 1 / max.height)
ret <- list(K = K, N.K = N.K, rate.anc = rate.anc, rate.dec = rate.dec,
height.anc = tree.anc.height,
height.dec = tree.dec.height,
anc = tree.anc, dec = tree.dec,
org = tree.org,
equal = tree.equal,
max = tree.max,
star = tree.star)
ret
} # End of gen.unit.K().
### Generate star ancestor tree binded with equal star descendant trees.
gen.equal.star.anc.dec <- function(K, N.K, rate.f = 0.5){
if(length(N.K) != K){
stop("The length of N.K is not equal to K.")
}
if(rate.f <= 0 || rate.f >= 1){
stop("The rate.f <= 0 or >= 1.")
}
if(K <= 0 || length(N.K) != K || any(N.K <= 0)){
stop("The K or N.K is not consistant.")
}
ms.anc <- ms(K, 1, opts = paste("-T -G", 1, sep = " "))
tree.anc <- read.tree(text = ms.anc[3])
tree.anc$tip.label <- paste("a", 1:K, sep = "")
tree.dec <- NULL
for(k in 1:K){
ms.dec <- ms(N.K[k], 1, opts = paste("-T -G", 1, sep = " "))
tree.dec[[k]] <- read.tree(text = ms.dec[3])
tree.dec[[k]]$tip.label <- paste("d", k, ".", 1:N.K[k], sep = "")
}
tmp.tree <- as.star.tree(tree.anc)
tmp.height <- get.rooted.tree.height(tmp.tree)
tree.equal.star <- rescale.rooted.tree(tmp.tree, rate.f / tmp.height)
for(k in 1:K){
tmp.tree <- as.star.tree(tree.dec[[k]])
tmp.height <- get.rooted.tree.height(tmp.tree)
tmp.tree <- rescale.rooted.tree(tmp.tree, (1 - rate.f) / tmp.height)
tree.equal.star <- bind.tree(tree.equal.star, tmp.tree,
where = which(tree.equal.star$tip.label ==
paste("a", k, sep = "")))
}
ret <- list(K = K, N.K = N.K, rate.f = rate.f,
anc = tree.anc, dec = tree.dec,
equal.star = tree.equal.star)
ret
} # End of gen.equal.anc.dec.star().
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.