lucasnell · April 24, 2020 20:56 · Apr 24, 2020
diff --git a/add_outgroup.R b/add_outgroup.R
@@ -0,0 +1,48 @@
+#'
+#' 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)
+}