Skip to content

Instantly share code, notes, and snippets.

@lucasnell
Created April 24, 2020 20:56
Show Gist options
  • Save lucasnell/394ebabf698cec912f45d8679848892b to your computer and use it in GitHub Desktop.
Save lucasnell/394ebabf698cec912f45d8679848892b to your computer and use it in GitHub Desktop.
#'
#' Add outgroup to phylogenetic tree.
#' The tree must be ultrametric.
#'
#' `tree`: a `phylo` object
#' `taxon_name`: name of the taxon to add; should be a single string
#' `divergence` total divergence time between a species inside `tree` and the outgroup,
#' in whatever units branchlengths in `tree` are in
#'
add_outgroup <- function(tree, taxon_name, divergence) {
stopifnot(inherits(tree, "phylo"))
stopifnot(is.ultrametric(tree))
stopifnot(inherits(taxon_name, "character") && length(taxon_name) == 1)
stopifnot(is.numeric(divergence) && length(divergence) == 1)
tip <- list(edge = matrix(c(2,1),1,2),
tip.label = taxon_name,
edge.length = divergence,
root.edge = 0,
Nnode = 1)
class(tip) <- "phylo"
tree$root.edge <- divergence - max(node.depth.edgelength(tree))
new_tr <- bind.tree(tree, tip, position = divergence - max(node.depth.edgelength(tree)))
# Remove unnecessary node and its relevant info:
rm_node <- new_tr$edge[new_tr$edge[,2] == length(new_tr$tip.label), 1]
rm_to <- new_tr$edge[,2] == rm_node
rm_from <- new_tr$edge[,1] == rm_node
new_edgelength <- sum(new_tr$edge.length[rm_from | rm_to])
new_tr$edge.length[rm_to] <- new_edgelength
new_tr$edge[rm_to, 2] <- length(new_tr$tip.label)
new_tr$edge <- new_tr$edge[!rm_from,]
new_tr$edge.length <- new_tr$edge.length[!rm_from]
new_tr$Nnode <- new_tr$Nnode - 1
if (!is.null(new_tr$node.label)) {
new_tr$node.label <- new_tr$node.label[-(rm_node - length(new_tr$tip.label))]
}
return (new_tr)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment