Construct Phylogeny form Edgelist

Question

I am working on large scale biodiversity simulations and have collected a lot of data on simulated phylogenies. I have a list of child-parent edgelist species in the form:

list = [1, 1, 2, 2]

      /-- 3
1 -- 2
      \-- 4 

children = [1, 2, 3, 4]
parents = [1, 1, 2, 2]   

Where the index is the child and the value is the parent. 1 is it's own parent (the root).

I have such a tree for a large dataset upwards of 1000 individuals. I can build the tree programmatically in igraph or networkx easily and visualise it. However, I cannot find a way to convert the tree into a format for phylogenetic analysis. There are function to convert dendropy and ete trees from igraph/networkx into their own formats but no way to re-import the igraph/networkx trees. There also appear to be no way to write the igraph/networkx trees into phylo formats (newick).

There also seems to be no support in dendropy or ete to programmatically build trees and it would be infeasible to manually build all the trees I have.

Is there a way to build the trees so I can analyse and visualise them as phylogenies?

Tried approach

library(ape)
p <- c(NA, 1, 1, 3, 3, 4, 4)
x<-data.frame(name=as.character(c(1:7)) , id=1:7,  parent = c(NA, 1,   1, 3, 3, 4, 4))
x1<-na.omit(x)[,3:2]
n <- c(c(1, 3, 4), c(2, 5, 6, 7))
edges <- matrix(c( x1$parent, x1$id), ncol=2)

tree<-list(edge=edges, Nnode = 3, tip.label = as.character(x$name[c(2, 5, 6, 7)]), node.label = as.character(x$name[c(1, 3, 4)]))

class(tree)<-"phylo"

#plot.phylo(tree, show.node.label=TRUE)

The plot command just causes R to crash and the tree object has an invalid structure.

Answer 1

This will be a complete answer, but it's got code, so it won't work as a comment. Working backwards from rtree(3), here's a way to build a phylo object from component parts, which might be helpful:

children <- c(5, 1, 2, 3)
parents <- c(4, 5, 5, 4)
x <- matrix(c(parents, children), ncol = 2)
tr <- list(edge = x, tip.label = 1:3, Nnode = 2)
class(tr) <- "phylo"
plot(tr)

Doing the same with your children and parents, doesn't work, however.

children <- c(1, 2, 3, 4)
parents <- c(1, 1, 2, 2)
x <- matrix(c(parents, children), ncol = 2)
tr <- list(edge = x, tip.label = 1:3, Nnode = 2)
class(tr) <- "phylo"
tr

Gives

Error in if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2) FALSE else TRUE : missing value where TRUE/FALSE needed

I don't know hot to make the edge list conform to the expected order for a phylo.

Answer 2

@kmm's code works with a few tweaks. One problem your edgelist has is that it contains a singleton node, which ape does not like to deal with. Turning it into a proper phylo object means removing the singleton 'root' node and instead making the child of this singleton node the root. Here is an example using a slightly larger tree, which demonstrates how to make a tree with both leaf nodes and internal, non-root nodes.

The trickiest part is translating the parents and children list into ape's unusual internal numbering system.

phylo_from_el <- function(el){
  i <- 0
  for(n in el[2:length(el)]){
    i <- i+1
    parents[i] <- n 
    children[i] <- i+1
  }
  nnode <- sum(children %in% parents) + 1
  nleaf <- length(el)-nnode
  # Translate preorder sequence into ape sequence
  ttable <- matrix(c(1:length(el), rep(0, length(el))), ncol=2)
  ttable[unique(parents),2] <- 1:nnode + nleaf
  ttable[-unique(parents),2] <- 1:nleaf

  parents_t <- sapply(parents, FUN = function(x)ttable[x,2])
  children_t <- sapply(children, FUN = function(x)ttable[x,2])

  edge <- matrix(c(parents_t, children_t), ncol = 2)
  tr <- list(edge = edge, tip.label = 1:nleaf, Nnode = nnode)
  class(tr) <- "phylo"
  return(tr)
}
el <- c(1,1,1,3,3)

tr <- phylo_from_el(el)

plot(tr)

Answer 3

This code did not work for me so I wrote an R Script to generate a Newick code.

I have a edge list tree_data

 head(tree_data)
  parent node branch.length label
1    117  112      5.093319   124
2    117  116      9.090761   128
3    117  136      9.219770   153
4    117  109      9.292008   121
5    117  108      9.370598   120
6    117  118      9.552086   130

Here is the R code:

tree_data$newick = paste(tree_data$label,':',signif(tree_data$branch.length,3),sep="")
nparents = 2
while(nparents > 1){
leaf = setdiff(tree_data$node,tree_data$parent)
parent_leaf = unique(tree_data$parent[!is.na(match(tree_data$node,leaf))])
index_remove = integer(0)
for(i in 1:length(parent_leaf)){
  parent_children = tree_data$node[which(tree_data$parent == parent_leaf[i])]
  all_children_in_leaf = all(parent_children %in% leaf)
  if(all_children_in_leaf==T){
    index_remove = c(index_remove,match(parent_children,tree_data$node))
    tree_data$newick[which(tree_data$node==parent_leaf[i])] = paste("(",paste(tree_data$newick[match(parent_children,tree_data$node)],collapse=" , "),")",tree_data$newick[which(tree_data$node==parent_leaf[i])] ,sep = " ")
  }
}
tree_data = tree_data[-index_remove,]
parent_next = unique(tree_data$parent[!is.na(match(tree_data$node,leaf))])
nparents = length(parent_next)
if (nparents==1) newick_code = paste("(",paste(tree_data$newick,collapse=","),")",colnames(Adj)[parent_next],";",sep = " ")
}

tree <- ape::read.tree(text=newick_code)
p = ggtree(tree, layout="circular")