### This function computes the height (distance to the root) of each node of a given tree
### If weight = TRUE, branch lengths are used; if a root edge is present, it is ignored!
computeHeights = function(tree, weight = FALSE) {
tree = checkPhylogeneticTree(tree)
n = ape::Ntip(tree)
N = 2 * n - 1
Tab = rep(0, N)
edges = tree$edge
for (ind in 1:(N - 1)) {
curRow = edges[ind,]
Tab[curRow[2]] = ifelse(weight, tree$edge.length[ind], 1) + Tab[curRow[1]]
}
Tab
}
### This function computes the sizes of the left and right subtrees rooted at internal nodes
computeLRSizes = function(tree) {
return(computeLRValues(tree, sum))
}
### This function computes the depths of the left and right subtrees rooted at internal nodes
### If weight = TRUE, the depth is the total weight of the path from the root to the node!
computeLRDepths = function(tree, weight = FALSE) {
return(computeLRValues(tree, max, weight = weight))
}
### This function factory recursively computes values for subtrees rooted at internal nodes
computeLRValues = function(tree, FUN, weight = FALSE) {
tree = checkPhylogeneticTree(tree)
n = ape::Ntip(tree)
N = 2 * n - 1
Tab = matrix(NA, n - 1, 2)
edges = tree$edge
for (ind in (N - 1):1) {
curRow = edges[ind,] - n
pos = Tab[curRow[1], 1]
W = ifelse(weight, tree$edge.length[ind], 1)
Tab[curRow[1], 2 - is.na(pos)] = W + ifelse(curRow[2] <= 0, 0, FUN(Tab[curRow[2],]))
}
Tab
}
### This function augments the tree with any subset of sizes and depths for each internal node
### and heights for each node; weight determines if the depths and heights are weighted.
augmentTree = function(tree, sizes, depths, heights, weight = FALSE) {
if (sizes && is.null(tree$subtreeSizes)) {
tree = addSubtreeSizes(tree)
}
if (depths && is.null(tree[[paste0("depths", rep("Weighted", weight))]])) {
tree = addDepths(tree, weight)
}
if (heights && is.null(tree[[paste0("heights", rep("Weighted", weight))]])) {
tree = addHeights(tree, weight)
}
tree
}
### This function computes the subtree sizes and adds them to a given tree
addSubtreeSizes = function(tree) {
treeSizes = computeLRSizes(tree)
tree$subtreeSizes = treeSizes
tree$subtreeSizesTips = (treeSizes + 1)/2
tree
}
### This function computes the weighted or unweighted depths and adds them to a given tree
addDepths = function(tree, weight) {
if (weight) {
tree$depthsWeighted = computeLRDepths(tree, TRUE)
} else {
tree$depths = computeLRDepths(tree, FALSE)
}
tree
}
### This function computes the weighted or unweighted heights and adds them to a given tree
addHeights = function(tree, weight) {
if (weight) {
tree$heightsWeighted = computeHeights(tree, TRUE)
} else {
tree$heights = computeHeights(tree, FALSE)
}
tree
}
### This function checks that the input tree is phylogenetic (i.e. rooted and binary)
### If not, arbitrarily roots it at the 1st internal node and/or binarizes it by multi2di
checkPhylogeneticTree = function(tree) {
if (!(ape::is.rooted(tree))) {
print("Warning: unrooted input tree; it will be rooted for further processing")
tree = ape::root(tree, node = ape::Ntip(tree) + 2, resolve.root = TRUE)
}
if (!(ape::is.binary.tree(tree))) {
print("Warning: non-binary input tree; it will be binarized for further processing")
tree = ape::multi2di(tree)
}
return(tree)
}
### This function is used to create a verbal description of a matrix based on parameters.
createName = function(Weight, Dist, Full, Lap, Norm) {
Name = ifelse(Lap, "Laplacian", "")
if (Dist) {
Name = paste0("Distance", ifelse(Name == "", "", " "), Name)
if (Full) {
Name = paste0(Name, ", Full")
}
}
if (Lap && Norm) {
Name = paste0(Name, ", Normalized")
}
if (Weight) {
Name = paste0(Name, ifelse(Name == "", "", ", "), "Weighted")
}
Name = paste(Name, "spectrum")
Name
}
### Utility function from the NeighborJoin.R file
createEmptyTree = function(numNodes) {
Tree = rep(list(c(NA, Inf, NA, NA)), numNodes)
names(Tree) = 1:numNodes
Tree
}
### Utility function from the NeighborJoin.R file
createNewNode = function(Tree, children, branchLengths, position) {
N = position
for (ind in 1:length(children)) {
curChild = children[ind]
Tree[[curChild]][1] = N
Tree[[curChild]][2] = branchLengths[ind]
}
Tree[[N]] = c(NA, Inf, children)
Tree
}
### Conversion function from the NeighborJoin.R file
convertPhylo = function(Tree) {
N = length(Tree)
n = (N + 1)/2
Labels = names(Tree)[1:n]
Edges = matrix(NA, N - 1, 3)
for (ind in 1:(N - 1)) {
curElement = Tree[[ind]]
Edges[ind,] = c(curElement[1], ind, curElement[2])
}
Edges[Edges[,1] > n, 1] = 3 * n - Edges[Edges[,1] > n, 1]
Edges[Edges[,2] > n, 2] = 3 * n - Edges[Edges[,2] > n, 2]
iTree = igraph::graph(t(Edges[, 1:2]), directed = TRUE)
orderT = igraph::graph.dfs(iTree, root = 3 * n - N)$order
Edges = Edges[order(Edges[,2]),]
newEdges = matrix(NA, nrow(Edges), ncol(Edges))
newEdges[order(orderT[-1]),] = Edges
newLengths = newEdges[,3]
newEdges = newEdges[,-3]
randTree = ape::rtree(n, rooted = TRUE)
randTree$edge = newEdges
randTree$edge.length = newLengths
randTree$tip.label = Labels
randTree
}
### This function creates the tree used to cross-check the statistics
createTestTree = function() {
Tr = createEmptyTree(numNodes = 13)
Tr = createNewNode(Tr, c(4,5), c(1,1), 8)
Tr = createNewNode(Tr, c(6,7), c(1,1), 9)
Tr = createNewNode(Tr, c(2,3), c(2,2), 10)
Tr = createNewNode(Tr, c(8,9), c(1,1), 11)
Tr = createNewNode(Tr, c(10,11), c(2,2), 12)
Tr = createNewNode(Tr, c(1,12), c(8,4), 13)
Tree = convertPhylo(Tr)
Tree$tip.label = c("a","f","g","d","e","b","c")
Tree$node.label = LETTERS[1:6]
Tree
}
### This function creates an unweighted version of a weighted tree
makeUnweighted = function(Tree) {
uTree = Tree
uTree$edge.length = rep(1, nrow(Tree$edge))
uTree
}
### This function scales the tree to have branch lengths that have a mean of 1; warning for any lengths <= 0
scaleTreeBranches = function(Tree) {
if (any(Tree$edge.length < 0)) {
warning("Some edges have negative branch lengths!")
}
if (any(Tree$edge.length == 0)) {
warning("Some edges have zero branch lengths!")
}
sTree = Tree
sTree$edge.length = sTree$edge.length/mean(sTree$edge.length)
sTree
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.