#' 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)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.