########################################################### ##Load required libraries ########################################################### library(ape) library(phytools) library(geiger) ########################################################### #Read in character data for Anolis, including geographic range, microhabitat (i.e., ecomorph), macrohabitat, svlCat, SVL, lnSVL ########################################################### read.csv("anolis_data_island.csv") -> anolis.dat #File includes categorical data on island (4 categories), microhabitat (6 categories), macrohabitat (3 categories), SVL (6 categories) and quantitative data on SVL data.frame(anolis.dat[,2:7]) -> anolisData rownames(anolisData) <- anolis.dat[,1] attach(anolisData) names(micro) <- rownames(anolisData) names(macro) <- rownames(anolisData) names(SVL) <- rownames(anolisData) ########################################################### #Read in trees for anoles and trim those taxa that are not found in the Greater Antilles ########################################################### read.tree("anolisFinalMCC.tre") -> anolisTree #Obtains MCC tree resulting from BEAST analyses name.check(anolisTree, micro) -> overlap overlap drop.tip(anolisTree, overlap$Tree.not.data) -> GAanolisTree ########################################################### #Visualize results ########################################################### help(read.tree) GAanolisTree$tip.label GAanolisTree$edge GAanolisTree$edge.length plot(GAanolisTree, show.tip.label=FALSE) nodelabels(cex=0.5) #Plot node labels text(rep(105, length(tree$tip.label)), 1:length(tree$tip.label), 1:length(tree$tip.label), cex=0.5) #Plot the numerical designations for the tip taxa #Now we're ready to plot our data plot(GAanolisTree, label.offset=8, cex=0.25, y.lim=c(-8,90), x.lim=c(0, 160)) #We've offset the text so we have room to plot character data and expanded the plotting area for the same reason #Add the microhabitat data microLabel <- character(length(GAanolisTree$tip.label)) #We're going to make a new matrix to store the colors will use to label our tip taxa names(microLabel) <- names(micro) microLabel[micro==0] <- "red" microLabel[micro==1] <- "yellow" microLabel[micro==2] <- "green" microLabel[micro==3] <- "light blue" microLabel[micro==4] <- "blue" microLabel[micro==5] <- "purple" points(rep(105, length(GAanolisTree$tip.label)), 1:length(GAanolisTree$tip.label), pch=21, bg=microLabel[GAanolisTree$tip.label], cex=0.5) text(105, length(GAanolisTree$tip.label)+2, "Microhabitat", cex=0.5) #Add the macrohabitat data macroLabel <- character(length(GAanolisTree$tip.label)) names(macroLabel) <- names(macro) macroLabel[macro==0] <- "red" macroLabel[macro==1] <- "green" macroLabel[macro==2] <- "blue" points(rep(109, length(GAanolisTree$tip.label)), 1:length(GAanolisTree$tip.label), pch=21, bg=macroLabel[GAanolisTree$tip.label], cex=0.5) text(109, length(GAanolisTree$tip.label)+5, "Macrohabitat", cex=0.5) #Add the SVL data abline(v=135, col="gray") segments(135, 1:length(SVL), 135 + (SVL[GAanolisTree$tip.label]/10), 1:length(SVL)) #Add scale bar segments(135, -1, 135+(max(SVL)/10), -1) text(c(135, 135+(max(SVL)/10)), c(-3,-3), c("0","191 mm"), cex=0.5) text(143, -6, "SVL", cex=0.5) #################################################################################### #Conduct basic ML reconstruction of microhabitat and macrohabitat in ape and visualize results #################################################################################### #Do preliminary analyses on microhabitat ace(micro, GAanolisTree, type="discrete", model="ER") -> GAanolisMicroAce GAanolisMicroAce #Check out the output you just saved GAanolisMicroAce$rates #You can call several elements of the output independently, such as the estimates rates. You can see a complete list of the output available in this manner by checking out the elements under "Value" in the help documentation for ace. GAanolisMicroAce$index.matrix #This is the index matrix describing our model GAanolisMicroAce$lik.anc #This matrix shows us the marginal likelihoods of each state at each internal node. You can confirm this list by asking whether the number of rows (nrow(GAanolisMicroAce$lik.anc)) is the same as the number of internal nodes in our tree (GAanolisTree$Nnode). plot(GAanolisTree, show.node.label=FALSE, label.offset=3, cex=0.5) points(rep(105, length(GAanolisTree$tip.label)), 1:length(GAanolisTree$tip.label), pch=21, bg=microLabel[GAanolisTree$tip.label]) title(main="ML Reconstructions for Microhabitat from 1 Rate Model in ACE") nodelabels(pie=GAanolisMicroAce$lik.anc, cex=1) #Use pie charts to illustrate the marginals on the tree #Do preliminary analyses on macrohabitat ace(macro, GAanolisTree, type="discrete", model="ER") -> GAanolisMacroAce GAanolisMacroAce #Check out the output you just saved plot(GAanolisTree, show.node.label=FALSE, label.offset=3, cex=0.5) title(main="ML Reconstructions for Macrohabitat from 1 Rate Model in ACE") points(rep(105, length(GAanolisTree$tip.label)), 1:length(GAanolisTree$tip.label), pch=21, bg=macroLabel[GAanolisTree$tip.label]) nodelabels(pie=GAanolisMacroAce$lik.anc, cex=1) #################################################################################### #Generate stochastic maps for microhabitat and macrohabitat #################################################################################### make.simmap(GAanolisTree, micro, nsim=100, model="ER") -> GAanolisMicroSimmap #Generate 100 stochastic maps for microhabitat make.simmap(GAanolisTree, macro, nsim=100, model="ER") -> GAanolisMacroSimmap #Generate 100 stochastic maps for macrohabitat macro_colors <-c("blue","red", "green"); names(macro_colors)<-c("0","1", "2") micro_colors <-c("red","yellow", "green", "light blue", "blue", "purple"); names(micro_colors)<-c("0","1", "2", "3", "4", "5") #Visualize one of your stochastic maps for microhabitat par(mfcol=c(2,2)) plotSimmap(GAanolisMicroSimmap[[1]], micro_colors, pts=FALSE, lwd=5) plotSimmap(GAanolisMicroSimmap[[2]], micro_colors, pts=FALSE, lwd=5) plotSimmap(GAanolisMicroSimmap[[3]], micro_colors, pts=FALSE, lwd=5) plotSimmap(GAanolisMicroSimmap[[4]], micro_colors, pts=FALSE, lwd=5) #Visualize one of your stochastic maps for macrohabitat plotSimmap(GAanolisMacroSimmap[[1]], macro_colors, pts=FALSE, lwd=5) plotSimmap(GAanolisMacroSimmap[[2]], macro_colors, pts=FALSE, lwd=5) plotSimmap(GAanolisMacroSimmap[[3]], macro_colors, pts=FALSE, lwd=5) plotSimmap(GAanolisMacroSimmap[[4]], macro_colors, pts=FALSE, lwd=5) #Let's extract some basic information from our maps. #Get numbers of changes for each trait # function counts transitions from a mapped history # written by Liam J. Revell 2013 countSimmap<-function(tree,states=NULL,message=TRUE){ n<-sum(sapply(tree$maps,length))-nrow(tree$edge) if(is.null(states)) states<-colnames(tree$mapped.edge) m<-length(states) TT<-matrix(NA,m,m,dimnames=list(states,states)) foo<-function(map,a,b){ if(length(map)==1) zz<-0 else { zz<-0; i<-2 while(i<=length(map)){ if(names(map)[i]==b&&names(map)[i-1]==a) zz<-zz+1 i<-i+1 } } return(zz) } for(i in 1:m) for(j in 1:m) if(i==j) TT[i,j]<-0 else TT[i,j]<-sum(sapply(tree$maps,foo,a=states[i],b=states[j])) return(list(N=n,Tr=TT,message=c( "N is the total number of character changes on the tree", "Tr gives the number of transitions from row state->column state"))) } countSimmap(GAanolisMicroSimmap[[1]]) phenogram(GAanolisMicroSimmap[[1]], micro, colors=micro_colors) map.overlap(GAanolisMicroSimmap[[1]],GAanolisMicroSimmap[[2]]) map.overlap(GAanolisMacroSimmap[[10]],GAanolisMacroSimmap[[21]]) #################################################################################### #Learning what SIMMAP does #Nothing in R should be a black box. Our goal here is to open the veil on stochastic mapping in R by going through a function that performs stochastic mapping step by step. #################################################################################### #make.simmap #> make.simmap #function (tree, x, model = "SYM", nsim = 1) #"model" here indicates the types and relative rates of changes we are permitting. See help(ace) for additional details on how to specify this matrix using precise commands or shortcuts. #{ #When we run the make.simmap function, we need to provide this function with four input variables, the tree (tree), the character data (x), the model for character reconstruction (model), and the number of simulations we'd like to do (nsim). We're going to run this function step by step and need to start by renaming our variables. x <- micro tree <- GAanolisTree model="ER" tree <- reorder.phylo(tree, "cladewise") #For more detail see help for reorder.phylo. Perhaps this is necessary because it is assumed by make.simmap that a edge leading from the root is the first edge? Not necessary if trees were read in using read.tree because this is the default #################################################################################### #1. Calculate conditional likelihood for each character state at each node of the tree, including the root. #################################################################################### XX <- ace(x, tree, type = "discrete", model = model) #Run maximum likelihood reconstruction of character evolution using the ace function of the ape package and save the output as a new object called XX. XX #Check out ancestral reconstruction par(mfcol=c(1,1)) #This just tells R to plot a single figure in the plotting window. plot(tree, show.node.label=FALSE, show.tip.label=TRUE, label.offset=3, cex=0.5) #Plot the tree title(main="ML Reconstructions for Microhabitat from 1 Rate Model in ACE") points(rep(105, length(GAanolisTree$tip.label)), 1:length(GAanolisTree$tip.label), pch=21, bg=microLabel[GAanolisTree$tip.label]) #Plot data at tips nodelabels(pie=GAanolisMicroAce$lik.anc, cex=1) #Plot marginals in the form of pie diagrams. #Get marginal values ready for analyses in make.simmap L <- XX$lik.anc #Store conditional likelihoods estimated via ML by ace L rownames(L) <- length(tree$tip) + 1:nrow(L) #Simply names the rows of the marginal likelihood data in L matrix by their node names L nodelabels(cex=0.5) #Confirm that the nodes are properly labeled text(rep(120, length(tree$tip.label)), 1:length(tree$tip.label), 1:length(tree$tip.label), cex=0.5) #Plot the numerical designations for the tip taxa #Prepare Q matrix for analyses in make.simmap XX$rates #Shows the stored output of ace and the rate estimate(s) element in particular. XX$index.matrix #Shows model of evolution Q <- matrix(XX$rates[XX$index.matrix], nrow(XX$index.matrix), ncol(XX$index.matrix), dimnames = list(colnames(XX$lik.anc), colnames(XX$lik.anc))) #Generate the Q matrix from rate parameters estimated in ace and associate model specification indices Q diag(Q) <- -colSums(Q, na.rm = TRUE) #Fills in diagonal of the Q matrix Q #Check out the final Q matrix #################################################################################### #2. Simulate ancestral states at each internal node by sampling from the posterior distribution of states. #################################################################################### mtrees <- list() #Generates an empty list that will eventually store or maps class(mtrees) <- "multiPhylo" #Tells R that the mtrees object is going to contain multiple phylogenetic trees #for (i in 1:nsim) #Initiates a loop that will take us through the generation of each stochastic map #{ mtree <- tree #Puts on input tree in the mtree multiPhylo object mtree$maps <- list() #Create a new element of the mtree multiPhylo that will store the data from our stochastic mapping node.states <- matrix(NA, nrow(tree$edge), ncol(tree$edge)) #Creates matrix to store nodal states at beginning and end of each edge node.states #Notice that this is just an empty matrix that will be used to store character states at the beginning [,1] and end [,2] of each edge. node.states[which(tree$edge[, 1] == (length(tree$tip) + 1)), 1] <- rstate(L[1, ]) #This line begins the process of assigning character states to nodes. Let's unpack exactly what is going on here. tree$edge[,1] #This is just the first column of the edge elements and indicates the starting point for each edge. length(tree$tip) + 1 #This should give us the number indicating the root of the tree. node.states[which(tree$edge[, 1] == (length(tree$tip) + 1)), 1] #In combination, this tells us to find the edges in our tree that start at the root. L[1, ] #Recall that L is our matrix of marginal values (i.e., the values that we've plotted in the form of the pie charts on our internal nodes) rstate(L[1, ]) #The internal phytools function rstate, which "picks a random element in a vector according to the probability assigned that element. It returns the name." node.states[which(tree$edge[, 1] == (length(tree$tip) + 1)), 1] <- rstate(L[1, ]) #All together then, this step identifies root and determines state based on marginals for each state at the root node.states #See which states have been filled in our map tree$edge #Confirm that node.states matrix is filled in for the root node #for (j in 1:nrow(tree$edge)) #This is the start of a loop that will determine the state for each other (non-root) node in our tree. #{ j <- 1 #Instead of running the whole loop, we're going to practice with just a few edges. By setting j to 1, we start with the first edge in our list of edges. #if (tree$edge[j, 2] > length(tree$tip)) #Asks if the edge we're considering - tree$edge[j,] - ends at an internal node rather than an tip. Edges end with the tip or terminal taxon in the second column of tree$edge. We therefore want to ask if the value in the second column of the edge row we're interested in tree$edge[j,2] is an internal node. We can do this by asking if the number at the end of the edge is greater than the number of terminals in our tree (internal nodes are numbered following tip taxa and therefore the names for internal nodes are always going to be greater than length(tree$tip)). #{ tree$edge[j, 2] #See which node we're dealing with at the end of the first edge. p <- expm(Q * tree$edge.length[j])[node.states[j, 1], ] * L[as.character(tree$edge[j, 2]), ] #We're going to unpack each element of this line. expm(Q * tree$edge.length[j]) #Exponentiate the Q matrix given the length of the edge in question (tree$edge.length[j]) [node.states[j, 1], ] #This is the state at the beginning of the edge in question expm(Q * tree$edge.length[j])[node.states[j, 1], ] #Extract the row of this exponentiated matrix that is relevant given the reconstructed state at the beginning of thh edge in question (node.states[j, 1]) L[as.character(tree$edge[j, 2]), ] #Obtain marginal likelihoods for node at the end of the edge in question. expm(Q * tree$edge.length[j])[node.states[j, 1], ] * L[as.character(tree$edge[j, 2]), ] node.states[j, 2] <- rstate(p/max(p)) #This is where we assign the state of the node at the end of our edge (tree$edge[j, 2]). To do this, we use the phytools function rstate, which "picks a random element in a vector [the vector that we just stored above called p] according to the probability assigned that element. It returns the name." rstate #You can learn more about this function by simply typing the name. rstate(p/max(p)) #See what you get when you run rstate. rstate(p/max(p)) #Do this a bunch of times. Do you always get the same value? p #See if you can figure out why you're always getting the same value. node.states[which(tree$edge[, 1] == tree$edge[j, 2]), 1] <- node.states[j, 2] #Having obtained the mapping for the node at the edge of a branch we can now propogate this mapping to the beginning(s) or the edges that extend away from this node (i.e., are in the first column of tree$edge) node.states #Check out the fruits of our labors #} This is just the end of the section that identifies internal nodes using the "if" function. #else { #The line starting with "if" above idenifies internal nodes, here we deal with the "else" case where nodes are not internal. Here we will merely assign known values to the terminal taxa in our matrix containing reconstructed trait valeus. j <- 10 node.states[j, 2] <- x[tree$tip.label[tree$edge[j, 2]]] #If the edge in question, ends with a terminal taxon (i.e., it does not meet the if (tree$edge[j, 2] > length(tree$tip)) condition above), we can merely assign this terminal with its known state in our map. x #We're using the tip data stored in x to fill in the node state for the edges that end with a terminal. tree$edge #Find the first edge in list that involves a tip. node.states[j, 2] <- x[tree$tip.label[tree$edge[j, 2]]] node.states #Notice that we've filled in the state for the terminal in the matrix storing our mapping. #} #The end of the section diagnosting terminals and assigning them known values. #Now that we know how make.simmap is going to assign mapped states to each node, we can run the whole loop and apply ancestral states to each node in our tree. node.states <- matrix(NA, nrow(tree$edge), ncol(tree$edge)) node.states[which(tree$edge[, 1] == (length(tree$tip) + 1)), 1] <- rstate(L[1, ]) for (j in 1:nrow(tree$edge)) { if (tree$edge[j, 2] > length(tree$tip)) { p <- expm(Q * tree$edge.length[j])[node.states[j, 1], ] * L[as.character(tree$edge[j, 2]), ] node.states[j, 2] <- rstate(p/max(p)) node.states[which(tree$edge[, 1] == tree$edge[j, 2]), 1] <- node.states[j, 2] } else { node.states[j, 2] <- x[tree$tip.label[tree$edge[j, 2]]] } } node.states #Check out the list of reconstructed states in your map #Visualize your node mapping. nodeMappings <- vector() #rownames(nodeLabels) <- node.states[,1][tree$edge[,1][(length(tree$tip.label)+1):(length(tree$tip.label)+tree$Nnode)] count <- 86 for(i in 1:tree$Nnode) { nodeMappings[i] <- node.states[which(tree$edge[,1] == count),1][1] count <- count + 1 } micro_colors <-c("red","yellow", "green", "light blue", "blue", "purple"); names(micro_colors)<-c("0","1", "2", "3", "4", "5") nodeMappings[nodeMappings[]=="0"] <- "red" nodeMappings[nodeMappings[]=="1"] <- "yellow" nodeMappings[nodeMappings[]=="2"] <- "green" nodeMappings[nodeMappings[]=="3"] <- "light blue" nodeMappings[nodeMappings[]=="4"] <- "blue" nodeMappings[nodeMappings[]=="5"] <- "purple" #micro_colors <-c("red","yellow", "green", "light blue", "blue", "purple"); names(micro_colors)<-c("0","1", "2", "3", "4", "5") plot(tree, show.node.label=FALSE, label.offset=3, cex=0.5) title(main="ML Reconstructions for Microhabitat from 1 Rate Model in ACE") points(rep(105, length(GAanolisTree$tip.label)), 1:length(GAanolisTree$tip.label), pch=21, bg=microLabel[GAanolisTree$tip.label]) nodelabels(pie=GAanolisMicroAce$lik.anc, cex=1) #Use pie charts to illustrate the marginals on the tree nodelabels(pch=22, bg=nodeMappings) #If we just run the function up to this point, we would end up with the reconstructed states at each ancestor in our map. This is helpful, but we're not done yet. Even with this information, we don't even know the basic features of our reconstruction, like how many changes happened (recall that multiple changes are possible on a single edge in ML. Now we want to figure out exactly how many changes occur in our map and where they occur along each edge. Again, changes may occur along an edge even if the start and end states are the same. #################################################################################### #3 For each simulation, we're going to need to identify changes, or lack thereof on each edge #################################################################################### else { accept = FALSE while (!accept) { time = 0 #Set starting time. state <- node.states[j, 1] #Sets the state at the start of the edge in question. new.state <- state #Just giving state another name. dt <- vector() #Creates a new vector that we'll use in a minute. map <- vector() #Creates a new vector that we'll use in a minute. k <- 1 #A count variable while (time < tree$edge.length[j]) #Starts running the continuous time simulation and continues this process until we get to the end of the edge in question. { dt[1] <- time #Sets the first element of our the dt vector we just created to time, which will start at 0, but may be changed before we loop through this for a second time. dt[2] <- dt[1] + rexp(n = 1, rate = -Q[state, state]) #Obtains the exponentially distributed interval of time before change occurs. dt[2] #Check out the resulting value if (dt[2] < tree$edge.length[j]) #Asks if length of time we stay in the original state dt[2] is less than the length of the edge in question, if so, we're going to want to pick a new state to transition to along this edge dt[2] < tree$edge.length[j] #Gives TRUE or FALSE answer about whether a change happens on the edge in question new.state <- rstate(Q[, state][-match(state,rownames(Q))]/sum(Q[, state][-match(state,rownames(Q))])) #If change happens, this step decides which stat we're going to change to based on the probability of each type of change expressed in the Q matrix. dt[2] <- min(dt[2], tree$edge.length[j]) #Assigns to dt[2] the lower of the two values obtained for when the change happens and how long the edge is. map[k] <- dt[2] - dt[1] #Add timing of the change to our map vector. names(map)[k] <- state k <- k + 1 #Adds 1 to counting variable state <- new.state #Resets the state value for the next run of this loop. time <- dt[2] #Resets time to reflect new value } #End of while loop that runs until we reach the end of the edge in question if (names(map)[length(map)] == node.states[j, 2]) accept = TRUE } #End of while loop } #End of else loop mtree$maps[[j]] <- map } allstates <- vector() for (j in 1:nrow(mtree$edge)) allstates <- c(allstates, names(mtree$maps[[j]])) allstates <- unique(allstates) mtree$mapped.edge <- matrix(data = 0, length(mtree$edge.length), length(allstates), dimnames = list(apply(mtree$edge, 1, function(x) paste(x, collapse = ",")), state = allstates)) for (j in 1:length(mtree$maps)) for (k in 1:length(mtree$maps[[j]])) mtree$mapped.edge[j,names(mtree$maps[[j]])[k]] <- mtree$mapped.edge[j, names(mtree$maps[[j]])[k]] + mtree$maps[[j]][k] mtrees[[i]] <- mtree } if (nsim == 1) mtrees <- mtrees[[1]] return(mtrees) } #Let's try this in SIMMAP #################################################################################### #Create an xml formatted character data for SIMMAP #################################################################################### tree <- GAanolisTree write.table(c("", "", paste("", sep="")), file="AnolisMicro_SIMMAP1.5.xml", quote=FALSE, col.names=FALSE, row.names=FALSE) for(i in 1:length(tree$tip.label)) { write.table(paste("", micro[tree$tip.label[i]], macro[tree$tip.label[i]], "", sep=""), file="AnolisMicro_SIMMAP1.5.xml", quote=FALSE, col.names=FALSE, row.names=FALSE, append=TRUE) } write.table("", file="AnolisMicro_SIMMAP1.5.xml", quote=FALSE, col.names=FALSE, row.names=FALSE, append=TRUE) write.table("", file="AnolisMicro_SIMMAP1.5.xml", quote=FALSE, col.names=FALSE, row.names=FALSE, append=TRUE) count <- 1 for(i in 1:length(tree$tip.label)) { write.table(paste("", tree$tip.label[i],"", sep=""), file="AnolisMicro_SIMMAP1.5.xml", quote=FALSE, col.names=FALSE, row.names=FALSE, append=TRUE) count<-count+1 } for(i in 1:length(GAanolis_500_beast_trees)) { write.table(paste("", write.tree(GAanolis_500_beast_trees[[i]]),"",sep=""), file="AnolisMicro_SIMMAP1.5.xml", quote=FALSE, col.names=FALSE, row.names=FALSE, append=TRUE) } write.table("", file="AnolisMicro_SIMMAP1.5.xml", quote=FALSE, col.names=FALSE, row.names=FALSE, append=TRUE) write.table("", file="AnolisMicro_SIMMAP1.5.xml", quote=FALSE, col.names=FALSE, row.names=FALSE, append=TRUE) write.nexus(GAanolis_500_beast_trees, file="test.nex") \", subTreeDefinitions <- function(phy) { subtrees(phy) -> phySubTrees subtreeMatrix <- matrix(NA, length(phySubTrees[[1]]$node.label), length(phySubTrees[[1]]$tip.label)+1) for(i in 1:length(phySubTrees)) { subtreeMatrix[i, 1] <- phySubTrees[[i]]$node.label[1] for(j in 1:length(phySubTrees[[i]]$tip.label)) { subtreeMatrix[i, j+1] <- match(phySubTrees[[i]]$tip.label[j], phy$tip.label) } subtreeMatrix[i, 1] <- paste("AddMRCA mrca_", subtreeMatrix[i,1], sep="") } return(subtreeMatrix) } read.tree("anolis_500_beast_trees.tre") -> anolis_500_beast_trees #Obtains trees resulting from BEAST time calibration analysis #read.nexus("Anolis_Bayes_Trees") -> AnolisTrees #Obtains ultrametric trees from BEAST posterior distribution #Check for overlap between tree and data, drop tips from tree that aren't in data (i.e., non-Greater Antillean species) name.check(anolisTree, micro) -> overlap drop.tip(anolisTree, overlap$Tree.not.data) -> GAanolisTree list() -> GAanolis_500_beast_trees for(i in 1:length(anolis_500_beast_trees)) { drop.tip(anolis_500_beast_trees[[i]], c("Bplumifrons", "Pacutirostris", overlap$Tree.not.data)) -> GAanolis_500_beast_trees[[i]] #GAanolis_500_beast_trees[[i]] <- reorder.phylo(GAanolis_500_beast_trees[[i]], "cladewise") #Necessary? } #write.nexus(GAanolis_500_beast_trees, file="GAanolis_500_beast_trees")