# This function takes a set of trees and a dataframe with species names as row.names and the clade they are assigned to in the first column. Using this it calculates the crown and stem age of the clade across a set of trees summarizing the results as the maximum, minimum and median clade age across all trees. This function requires ape and phytools packages. cladeage<-function(phylo, cladelist){ library(ape) library(phytools) uniquefam<-unique(cladelist[,1]) # this generates a vector of unique clade names # This loop generates a set of vectors with the family names of Acanthomorphs in the tree and assigns species names to them which I will use for estimating MRCA of all families in all the trees. familynames<-list() for(i in 1:length(uniquefam)) { familynames[[i]]<-row.names(cladelist)[cladelist ==as.character(uniquefam[i])] } maxcrowncladeage<-vector(length=length(familynames)) mincrowncladeage <-vector(length=length(familynames)) mediancrowncladeage <-vector(length=length(familynames)) maxstemcladeage <-vector(length=length(familynames)) minstemcladeage <-vector(length=length(familynames)) medianstemcladeage <-vector(length=length(familynames)) percentparaphyly<-vector(length=length(familynames)) # this loops over the clade MRCA function in phytools and calculates the age of the crown and stem if you have more than 1 species in a clade or just the stem age if you only have 1 species allcrownages<-matrix(nrow=length(phylo), ncol=length(uniquefam)) colnames(allcrownages)<-uniquefam allstemages<-matrix(nrow=length(phylo), ncol=length(uniquefam)) colnames(allstemages)<-uniquefam paraphyly<-matrix(nrow=length(phylo), ncol=length(uniquefam)) colnames(paraphyly)<-uniquefam for(i in 1:length(familynames)){ print(i) if (length(familynames[[i]])>1){ for (j in 1:length(phylo)) { print(j) bt<-branching.times(phylo[[j]]) foocrown<-findMRCA(phylo[[j]], familynames[[i]]) allcrownages[j,i]<-bt[names(bt)==foocrown] stemnodemulti<-phylo[[j]]$edge[,1][phylo[[j]]$edge[,2]==foocrown] allstemages[j,i]<-bt[names(bt)==stemnodemulti] len<-length(phylo[[j]]$tip.label) # finds the length of the tip labels so can identify tip nodes when trying to identify paraphyly foodesc<-getDescendants(phylo[[j]], foocrown) cladetips<-foodesc[foodesc<=len] # this identifies just the descendants that are tips paraphyly[j,i]<-if(length(familynames[[i]])