R/phylo-build.R

#' Create phylo class tree object
#'
#' This function converts the output of any of the sims4 simulators and creates the full case tree
#' first before either returning this as a phylo object, or the collapsed tree, whereby any single
#' nodes in the network have been collapsed. 
#'
#'
#' @param cas.tree the case tree generated from a simulator function within sims4  
#' @param full boolean specifying ehwther the full case network is to be returned, or the collapsed (default)
#' @return Returns a phylo object 
#' 
#'
#' 


Case_To_Phylo <- function(case.tree, full=FALSE){
  
  ## Create tips conforming to phylo class convention, i.e. tips are labelled 1:n, and nodes from n+1:N ##
  
  extended <- list(tips=0)
  
  ## create needed variables
  n.tips <- length(which(case.tree$IsFather==FALSE))
  n.nodes <- as.numeric(sum(case.tree$IsFather==TRUE))
  n.cases <- length(case.tree$ID)
  
  ## find which cases have not caused new infections and label them 1:n.tips
  extended$tips[which(case.tree$IsFather==FALSE)]=1:n.tips 
  ## find which cases have  caused new infections and label them up to n.cases
  extended$tips[which(case.tree$IsFather==TRUE)]=(1+n.tips):n.cases 
  
  ## create edge matrix and fill with the end of branches in column 2 and ancestors in column 1
  edge <- matrix(nrow=n.cases,ncol=2)
  edge[,2] <- extended$tips
  edge[,1] <- extended$tips[case.tree$Parent] 
  ## need to remove the 1st row which is the imaginary branch causing the root
  edge <- edge[-1,]
  
  ## create branch/edge length and remove row 1 again
  edge.length = case.tree$Branch.Length[-1]
  
  ## create tip and node labels in form n/t#_sequenceID_Time
  t <- paste("t",seq(1:n.tips),sep="")
  tip.label <- paste(t,case.tree$Sequence[which(case.tree$IsFather==FALSE)],case.tree$Time[which(case.tree$IsFather==FALSE)],sep="_")
  n <- paste("n",seq(1:n.nodes),sep="")
  node.label <- paste(n,case.tree$Sequence[which(case.tree$IsFather==TRUE)],case.tree$Time[which(case.tree$IsFather==TRUE)],sep="_")
  
  ## can not find anywhere in ape to create a phylo object from scratch, so first create one and change it
  tree <- rtree(n = 2)
  tree$Nnode <- n.nodes
  tree$edge <- edge
  tree$edge.length <- edge.length
  tree$tip.label <- tip.label
  tree$node.label <- node.label
  
  if (full == FALSE){
    ## The tree contains single nodes that the inbuilt plot functon can't handle (FigTree can tho...)
    ## Remove all sngle nodes 
    single_tree <- collapse.singles(tree);
    return(single_tree)
  } else {
    return(tree)
  }
  
}
OJWatson/sims4 documentation built on May 7, 2019, 8:33 p.m.