Extract edges and plot on a phylogenetic tree

get_edges.R

# PLOTTING WITH HIGHLIGHTED TIPS

# LIB
library(ape)

# FUNCTIONS FROM MoreTreeTools
getEdges <- function (tree, node = NULL, tips = NULL, type = 1) 
{
  if (!is.null(node)) {
    edges <- c()
    while (TRUE) {
      bool <- tree$edge[, 1] %in% node
      if (sum(bool) > 1) {
        node <- tree$edge[bool, 2]
        edges <- c(edges, which(bool))
      }
      else {
        break
      }
    }
    return(edges)
  }
  if (!type %in% c(1, 2, 3)) {
    stop("Type must be an integer: 1, 2 or 3.")
  }
  if (length(tips) == length(tree$tip.label)) {
    return(tree$edge)
  }
  if (type == 1 & length(tips) == 1) {
    stop("length (tips) == 1 :\n         Cannot return a single edge for type 1.")
  }
  if (is.character(tips)) {
    edges <- match(match(tips, tree$tip.label), tree$edge[, 
                                                          2])
  }
  else {
    edges <- match(tips, tree$edge[, 2])
  }
  end.nodes <- tree$edge[edges, 1]
  term.node <- length(tree$tip.label) + 1
  if (all(end.nodes %in% term.node)) {
    return(edges)
  }
  else {
    if (type == 3) {
      while (any(duplicated(end.nodes))) {
        start.node <- end.nodes[duplicated(end.nodes)][1]
        if (sum(tree$edge[, 1] %in% start.node) == sum(end.nodes %in% 
                                                       start.node)) {
          edge <- match(start.node, tree$edge[, 2])
          end.node <- tree$edge[edge, 1]
          edges <- c(edges, edge)
          end.nodes <- c(end.nodes[!end.nodes %in% start.node], 
                         end.node)
        }
        else {
          end.nodes <- end.nodes[end.nodes != start.node]
        }
      }
    }
    else {
      while (TRUE) {
        if (type == 2) {
          if (sum(term.node == end.nodes) == length(end.nodes)) {
            break
          }
        }
        else {
          if (sum(end.nodes[1] == end.nodes) == length(end.nodes)) {
            break
          }
        }
        end.nodes <- sort(end.nodes, TRUE)
        start.node <- end.nodes[1]
        edge <- match(start.node, tree$edge[, 2])
        end.node <- tree$edge[edge, 1]
        edges <- c(edges, edge)
        end.nodes <- c(end.nodes[!end.nodes %in% start.node], 
                       end.node)
      }
    }
    return(edges)
  }
}

getParent <- function (tree, node = NULL, tips = NULL, edges = NULL) 
{
  if (!is.null(node) & length(node) == 1) {
    if (!is.numeric(node)) {
      stop("Node must be numeric")
    }
    if (node > getSize(tree) + tree$Nnode) {
      stop("Node not in tree")
    }
    if ((node == getSize(tree) + 1) & is.rooted(tree)) {
      return(node)
    }
    return(tree$edge[tree$edge[, 2] == node, 1])
  }
  else if (!is.null(tips)) {
    if (is.character(tips)) {
      edges <- match(match(tips, tree$tip.label), tree$edge[, 
                                                            2])
    }
    else {
      edges <- match(tips, tree$edge[, 2])
    }
  }
  else if (!is.null(node)) {
    edges <- which(tree$edge[, 2] %in% node)
  }
  else if (!is.null(edges)) {
    if (is.character(edges) & !is.null(tree$edge.label)) {
      edges <- match(edges, tree$edge.label)
    }
  }
  else {
    stop("Must provide either edges, tips or nodes argument")
  }
  end.nodes <- tree$edge[edges, 1]
  term.node <- length(tree$tip.label) + 1
  while (TRUE) {
    if (sum(end.nodes[1] == end.nodes) == length(end.nodes)) {
      break
    }
    end.nodes <- sort(end.nodes, TRUE)
    start.node <- end.nodes[1]
    edge <- match(start.node, tree$edge[, 2])
    end.node <- tree$edge[edge, 1]
    edges <- c(edges, edge)
    end.nodes <- c(end.nodes[!end.nodes %in% start.node], 
                   end.node)
  }
  return(end.nodes[1])
}

# SUGGESTED FUNCTIONS
plotUniqueEdges <- function(tips)
{
  edgs <- getEdges(tree, tips=tips, type=1)
  edgs_cols <- ifelse(1:nrow(tree$edge) %in% edgs, "black", "grey")
  plot(tree, show.tip.label=FALSE, edge.color=edgs_cols, no.margin=TRUE,
       edge.width=1.5)
}

plotUniqueEdgeSets <- function(tip_set_1, tip_set_2)
{
  prnt_1 <- getParent(tree, tips=tip_set_1)
  edgs_1 <- getEdges(tree, node=prnt_1)
  prnt_2 <- getParent(tree, tips=tip_set_2)
  edgs_2 <- getEdges(tree, node=prnt_2)
  shrd_edgs <- getEdges(tree, tips=c(tip_set_1, tip_set_2),
                        type=1)
  edgs <- unique(c(edgs_1, edgs_2, shrd_edgs))
  edgs_cols <- ifelse(1:nrow(tree$edge) %in% edgs, "black", "grey")
  plot(tree, show.tip.label=FALSE, edge.color=edgs_cols, no.margin=TRUE,
       edge.width=1.5)
}

# TREE (25 TIPS) + PLOT
tree <- read.tree(text="((((t19:0.4166666667,(t4:0.375,((t12:0.125,(t7:0.08333333333,(t21:0.04166666667,t2:0.04166666667):0.04166666667):0.04166666667):0.2083333333,((t15:0.125,(t11:0.08333333333,(t20:0.04166666667,t6:0.04166666667):0.04166666667):0.04166666667):0.04166666667,t3:0.1666666667):0.1666666667):0.04166666667):0.04166666667):0.08333333333,(t24:0.04166666667,t16:0.04166666667):0.4583333333):0.04166666667,t25:0.5416666667):0.4583333333,((t9:0.08333333333,(t18:0.04166666667,t5:0.04166666667):0.04166666667):0.3333333333,(t1:0.2916666667,(((t13:0.04166666667,t8:0.04166666667):0.08333333333,(t10:0.04166666667,t17:0.04166666667):0.08333333333):0.125,(t22:0.08333333333,(t14:0.04166666667,t23:0.04166666667):0.04166666667):0.1666666667):0.04166666667):0.125):0.5833333333);")
# unique
tips <- sample(1:length(tree$tip.label), 2)
plotUniqueEdges(tree$tip.label[tips])
tiplabels(tip=tips, frame="none", col="red1", pch=19)
# unique sets
tip_set_1 <- sample(1:(length(tree$tip.label)-1), 1)
tip_set_1 <- c(tip_set_1, tip_set_1 + 1)  # keep random sets near each other for visual reasons
tip_set_2 <- sample(1:(length(tree$tip.label)-1), 1)
tip_set_2 <- c(tip_set_2, tip_set_2 + 1)
plotUniqueEdgeSets(tip_set_1, tip_set_2)
tiplabels(tip=tip_set_1, frame="none", col="red1", pch=19)
tiplabels(tip=tip_set_2, frame="none", col="cornflowerblue", pch=19)

# TREE (100 TIPS) + PLOT
tree <- read.tree(text="((((((t7:0.0303030303,(t18:0.0202020202,(t73:0.0101010101,t4:0.0101010101):0.0101010101):0.0101010101):0.1616161616,((((((t96:0.0101010101,t32:0.0101010101):0.0101010101,t10:0.0202020202):0.06060606061,((t64:0.0303030303,(t80:0.0202020202,(t56:0.0101010101,t44:0.0101010101):0.0101010101):0.0101010101):0.0202020202,(t85:0.0101010101,t50:0.0101010101):0.0404040404):0.0303030303):0.05050505051,((t41:0.0303030303,(t66:0.0202020202,(t60:0.0101010101,t89:0.0101010101):0.0101010101):0.0101010101):0.0101010101,t86:0.0404040404):0.09090909091):0.0101010101,t57:0.1414141414):0.0101010101,t30:0.1515151515):0.0404040404):0.1515151515,(((t70:0.05050505051,(t53:0.0404040404,(t99:0.0303030303,(t16:0.0202020202,(t19:0.0101010101,t59:0.0101010101):0.0101010101):0.0101010101):0.0101010101):0.0101010101):0.0303030303,(t29:0.0202020202,(t27:0.0101010101,t46:0.0101010101):0.0101010101):0.06060606061):0.06060606061,(((t1:0.0101010101,t52:0.0101010101):0.0202020202,(t2:0.0101010101,t45:0.0101010101):0.0202020202):0.0202020202,(t20:0.0101010101,t55:0.0101010101):0.0404040404):0.09090909091):0.202020202):0.4545454545,((((((t8:0.0101010101,t22:0.0101010101):0.0202020202,(t5:0.0101010101,t3:0.0101010101):0.0202020202):0.06060606061,((((t37:0.0101010101,t54:0.0101010101):0.0202020202,(t72:0.0101010101,t100:0.0101010101):0.0202020202):0.0101010101,t9:0.0404040404):0.0101010101,t21:0.05050505051):0.0404040404):0.1212121212,(t47:0.1111111111,((((t63:0.0101010101,t12:0.0101010101):0.0101010101,t58:0.0202020202):0.0101010101,t95:0.0303030303):0.07070707071,((t43:0.0101010101,t79:0.0101010101):0.05050505051,(((t92:0.0101010101,t36:0.0101010101):0.0101010101,t34:0.0202020202):0.0202020202,(t87:0.0101010101,t67:0.0101010101):0.0303030303):0.0202020202):0.0404040404):0.0101010101):0.101010101):0.05050505051,((((t88:0.0101010101,t84:0.0101010101):0.0101010101,t48:0.0202020202):0.0101010101,t31:0.0303030303):0.0101010101,t6:0.0404040404):0.2222222222):0.1818181818,((t28:0.0101010101,t97:0.0101010101):0.1616161616,((t69:0.101010101,(((t75:0.0101010101,t24:0.0101010101):0.0202020202,(t71:0.0101010101,t25:0.0101010101):0.0202020202):0.06060606061,(t49:0.05050505051,(((t26:0.0101010101,t11:0.0101010101):0.0101010101,t62:0.0202020202):0.0202020202,(t65:0.0101010101,t83:0.0101010101):0.0303030303):0.0101010101):0.0404040404):0.0101010101):0.05050505051,(t35:0.0404040404,((t40:0.0101010101,t39:0.0101010101):0.0202020202,(t17:0.0101010101,t13:0.0101010101):0.0202020202):0.0101010101):0.1111111111):0.0202020202):0.2727272727):0.3535353535):0.0404040404,((t38:0.0101010101,t78:0.0101010101):0.0202020202,(t61:0.0101010101,t23:0.0101010101):0.0202020202):0.8080808081):0.1616161616,((((t98:0.0202020202,(t81:0.0101010101,t14:0.0101010101):0.0101010101):0.0303030303,(t33:0.0202020202,(t76:0.0101010101,t74:0.0101010101):0.0101010101):0.0303030303):0.06060606061,(((t82:0.0101010101,t93:0.0101010101):0.0202020202,(t51:0.0101010101,t42:0.0101010101):0.0202020202):0.0202020202,(t15:0.0101010101,t91:0.0101010101):0.0404040404):0.06060606061):0.0404040404,(((t94:0.0101010101,t77:0.0101010101):0.0101010101,t68:0.0202020202):0.0101010101,t90:0.0303030303):0.1212121212):0.8484848485);")
# unique
tips <- sample(1:length(tree$tip.label), 2)
plotUniqueEdges(tree$tip.label[tips])
tiplabels(tip=tips, frame="none", col="red1", pch=19)
# unique sets
tip_set_1 <- sample(1:(length(tree$tip.label)-1), 1)
tip_set_1 <- c(tip_set_1, tip_set_1 + 1)  # keep random sets near each other for visual reasons
tip_set_2 <- sample(1:(length(tree$tip.label)-1), 1)
tip_set_2 <- c(tip_set_2, tip_set_2 + 1)
plotUniqueEdgeSets(tip_set_1, tip_set_2)
tiplabels(tip=tip_set_1, frame="none", col="red1", pch=19)
tiplabels(tip=tip_set_2, frame="none", col="cornflowerblue", pch=19)