#' Resolve Polytomies by Order of First Appearance
#'
#' Resolves polytomies in trees with lineages arranged in a pectinate pattern
#' (i.e. a ladder-like subtree), ordered by the time of first appearance (FAD)
#' for each lineage.
#'
#' @details
#' This method of resolving polytomies assumes that the order of stratigraphic
#' appearance perfectly depicts the order of branching. This may not be a good
#' assumption for poorly sampled fossil records.
#'
#' This function is for resolving trees when a continuous time-scale is known.
#' For discrete time-scales, see the function \code{bin_timePaleoPhy}.
#'
#' Taxa with the same identical first appearance date will be ordered randomly.
#' Thus, the output is slightly stochastic, but only when ties exist. This is
#' probably uncommon with real data on continuous time-scales.
#'
#' Taxa not shared between the input tree and the \code{timeData} matrix, or listed as
#' having a FAD or LAD of \code{NA} in \code{timeData} will be dropped and will not be
#' included in the output tree.
#'
#' See this blog post for more information:
#'
#' \url{https://nemagraptus.blogspot.com/2012/07/resolving-polytomies-according-to.html}
#'
#' @param tree A phylogeny, as an object of class \code{phylo}.
#' @param timeData Two-column \code{matrix} of per-taxon first and last occurrences in
#' absolute continuous time.
#' @return Returns the modified tree as an object of class \code{phylo}, with no edge
#' lengths.
#' @author David W. Bapst
#' @seealso \code{\link[ape]{di2multi}}
#' @examples
#' \donttest{
#' set.seed(444)
#' record <- simFossilRecord(p = 0.1, q = 0.1, nruns = 1,
#' nTotalTaxa = c(100,200))
#' taxa <- fossilRecord2fossilTaxa(record)
#' tree <- taxa2cladogram(taxa)
#' ranges <- sampleRanges(taxa,r = 0.5)
#' tree1 <- timeLadderTree(tree,ranges)
#' layout(1:2)
#' plot(ladderize(tree),show.tip.label = FALSE)
#' plot(ladderize(tree1),show.tip.label = FALSE)
#'
#' #an example with applying timeLadderTree to discrete time data
#' rangeData <- binTimeData(ranges,int.len = 5) #sim discrete range data
#' tree2 <- bin_timePaleoPhy(tree,timeList = rangeData,timeres = TRUE)
#' plot(ladderize(tree),show.tip.label = FALSE)
#' plot(ladderize(tree2),show.tip.label = FALSE)
#' axisPhylo()
#'
#' layout(1)
#' }
#' @export timeLadderTree
timeLadderTree <- function(tree,timeData){
#resolves all polytomies in a tree as ladders to match FADs in timeData
#only applicable to continuous time data
#require(ape)
#first sanitize data
if(!inherits(tree,"phylo")){
stop("tree must be of class 'phylo'")
}
if(!inherits(timeData,"matrix")){
if(inherits(timeData,"data.frame")){
timeData <- as.matrix(timeData)
}else{
stop("timeData not of matrix or data.frame format")
}
}
if(ncol(timeData) == 6){ #also allow it to accept taxad objects
timeData <- timeData[,3:4,drop = FALSE]
}
#first clean out all taxa which are NA or missing in timeData
#remove taxa that are NA or missing in timeData
tree <- drop.tip(tree,tree$tip.label[is.na(match(tree$tip.label,names(which(!is.na(timeData[,1])))))])
if(Ntip(tree)<2){
stop("Less than two valid taxa shared between the tree and temporal data")
}
timeData <- timeData[!is.na(timeData[,1]),]
if(any(is.na(timeData))){
stop("Weird NAs in Data??")
}
if(any(timeData[,1]<timeData[,2])){
stop("timeData is not in time relative to modern (decreasing to present)")
}
#node IDs of polytomies
nodes <- Ntip(tree)+1:Nnode(tree)
polys <- nodes[sapply(nodes,function(x) sum(tree$edge[,1] == x)>2)]
#error if no polytomies
if(length(polys) == 0){stop("No Polytomies in Cropped Tree?")}
while(length(polys)>0){
node <- polys[1]
#make a list of the first FAD of each descendant
desc <- tree$edge[tree$edge[,1] == node,2]
dTips <- lapply(desc,function(x) if(x>Ntip(tree)){
tree$tip.label[prop.part(tree)[[x-Ntip(tree)]]]
}else{tree$tip.label[x]})
#get max FAD
dFADs <- sapply(dTips,function(x) max(timeData[x,1]))
#build a ladderized subtree of the node
subtree <- stree(length(dFADs),type = "left")
#figure out the order of tips, use a second vector of random numbers for ties
subtree$tip.label <- desc[order(-dFADs,sample(1:length(dFADs)))]
#stick on tip labels
for(i in desc){
dtip <- which(subtree$tip.label == i)
if(i>Ntip(tree)){ #if its a clade
subclade <- extract.clade(tree,i)
subtree <- bind.tree(subtree,subclade,where = dtip)
subtree <- collapse.singles(subtree)
}else{subtree$tip.label[dtip] <- tree$tip.label[i]} #if its a tip
}
#replace original node with new, resolved, scaled node
if(node != (Ntip(tree)+1)){ #if it isn't the root node
drtips <- prop.part(tree)[[node-Ntip(tree)]]
tip_lab <- tree$tip.label[drtips[1]] #cut out all but one tip, just to put it back together later
droptree <- collapse.singles(drop.tip(tree,drtips[-1]))
droptree <- bind.tree(droptree,subtree,where = which(droptree$tip.label == tip_lab)) #put in subtree at tip
tree <- droptree
}else{ #if it is the root node
tree <- subtree
}
nodes <- Ntip(tree)+1:Nnode(tree)
polys <- nodes[sapply(nodes,function(x) sum(tree$edge[,1] == x)>2)]
}
tree$edge.length <- NULL
return(tree)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.