median.tree: Get the tree with median taxa ages among a list of trees

In jacobmaugoust/ULT: Useful Little Things

Get the tree with median taxa ages among a list of trees

Description

For a given list of trees of same structure (same tips, same branching pattern) but of various ages of taxa, this function computes (and returns) a single tree with median node ages. It also names the nodes by either taking the already provided node labels or by creating new ones.

Usage

median.tree(trees,node.labelling=FALSE,node.label.opt=NULL)

Arguments

trees	The list of phylogenetic trees to consider
node.labelling	Logical. Whether to set node labels or not.
node.label.opt	A list of arguments to be passed to the makeNodeLabel function

Examples

require(ape)
require(phytools)
# Get a random tree of 20 tips
set.seed(1)
tree<-rtree(20)
tree<-ULT::add.ages.phylo(tree,tip.ages=rep(0,Ntip(tree)))
# Create 50 trees with some variation in its taxa ages (except the root), between the actual age and parent node age
nH<-max(nodeHeights(tree))-nodeHeights(tree) # Get ages (from the root) for each taxon
ages_bounds<-t(mapply(function(x){
  if(x==(Ntip(tree)+1)){
    rep(max(nH),2)
  }
  else{
    c(nH[tree$edge[,2]==x,2],nH[tree$edge[,2]==x,1])
  }
},c(1:(Ntip(tree)+Nnode(tree))))) # Get bounds of variation for each taxon
trees<-list()
for(i in 1:50){
  set.seed(i+1)
  trees[[i]]<-tree
  new_nH<-nH
  for(j in 1:(Ntip(tree)+Nnode(tree))){
    new_age<-runif(1,ages_bounds[j,1],ages_bounds[j,2])
    for(k in 1:2){
      new_nH[which(tree$edge[,k]==j),k]<-new_age
    }
  }
  new_nH<-max(new_nH)-new_nH
  trees[[i]]$edge.length<-apply(new_nH,1,diff)
}
# Plot four of them to just check the discrepancies in taxa ages
par(mfrow=c(2,2),mar=c(2.1,0,0,0))
for(i in 1:4){
  plot(trees[[i]])
  axisPhylo()
}
# Get tree with median taxa ages
medtree<-median.tree(trees)

# Plot median tree
par(mfrow=c(1,1))
plot(medtree,show.tip.label = FALSE,x.lim=c(0,4.4))
for(i in 1:Ntip(medtree)){
  text(4.4,i,medtree$tip.label[i],adj=c(0,0.5),font=3)
}
axisPhylo()
# Add bars for the age variation of each taxon (light gray for tips, darker gray for nodes)
phylo.coords<-get("last_plot.phylo", envir = ape::.PlotPhyloEnv)[c("root.time","yy")]
for(i in 1:nrow(ages_bounds)){
  if(i==(Ntip(medtree)+1)){next}
  segments(phylo.coords$root.time-ages_bounds[i,1],phylo.coords$yy[i],phylo.coords$root.time-ages_bounds[i,2],phylo.coords$yy[i],col=rgb(t(rep(ifelse(i<Ntip(medtree),0.5,0),3)),alpha=0.25),lwd=10)
}
# Add original tree as shaded edges, to check whether bounds have been respected well
par(new=TRUE)
plot(tree,edge.col=rgb(1,0,0,0.25),show.tip.label=FALSE,x.lim=c(0,4.4))
par(mar=c(5,4,4,2)+0.1)

jacobmaugoust/ULT documentation built on May 16, 2023, 1:29 p.m.