# Using R Studio. # Set working directory > setwd("/Users/michael.branstetter/Documents/aaProjects/Bees/orchard-bees/Phylogenetics/Results/bamm") # Get BAMM priors library(BAMMtools) setBAMMpriors(read.tree("my_tree.tre")) ############################################### # Prior block chosen by BAMMtools::setBAMMpriors expectedNumberOfShifts = 1.0 lambdaInitPrior = 6.22928773238883 lambdaShiftPrior = 0.00976797198764379 muInitPrior = 6.22928773238883 #### End Prior block ###################### # Run BAMM using Terminal > ~/Programs/bamm -c control.txt # Load tree and event data into BAMMtools tree <- read.tree("my_tree.tre") edata <- getEventData(tree, eventdata = "event_data.txt", burnin=0.1) # Assess convergence in R using BAMMtools and Coda > mcmcout <- read.csv("mcmc_out.txt", header=T) > plot(mcmcout$logLik ~ mcmcout$generation) > burnstart <- floor(0.1 * nrow(mcmcout)) > postburn <- mcmcout[burnstart:nrow(mcmcout), ] > library(coda) > effectiveSize(postburn$N_shifts) var1 730.7812 > effectiveSize(postburn$logLik) var1 605.0266 > plot(postburn$logLik ~ postburn$generation) # Analyze rate shifts postfile <- "mcmc_out.txt" bfmat <- computeBayesFactors(postfile, expectedNumberOfShifts=1, burnin=0.1) bfmat 1 2 3 4 5 6 7 1 1.000000000 1.011859583 1.600390156 3.83081897 16.25762195 33.32813 166.6406 2 0.988279419 1.000000000 1.581632653 3.78591954 16.06707317 32.93750 164.6875 3 0.624847632 0.632258065 1.000000000 2.39367816 10.15853659 20.82500 104.1250 4 0.261040788 0.264136622 0.417767107 1.00000000 4.24390244 8.70000 43.5000 5 0.061509611 0.062239089 0.098439376 0.23563218 1.00000000 2.05000 10.2500 6 0.030004688 0.030360531 0.048019208 0.11494253 0.48780488 1.00000 5.0000 7 0.006000938 0.006072106 0.009603842 0.02298851 0.09756098 0.20000 1.0000 plotPrior(postfile, expectedNumberOfShifts=1) edata <- getEventData(tree, eventdata = "event_data.txt", burnin=0.1) shift_probs <- summary(edata) cst=cumulativeShiftProbsTree(edata) plot.phylo(cst, cex=0.5) branch_priors = getBranchShiftPriors(tree,expectedNumberOfShifts = 1) plot.phylo(branch_priors, cex=0.5) mo=marginalOddsRatioBranches(edata,branch_priors) Analyzed 1801 posterior samples Shift posterior distribution: 1 0.58000 2 0.30000 3 0.09700 4 0.01600 5 0.00220 6 0.00110 7 0.00056 plot.bammdata(edata, lwd=2, labels = T, cex = 0.4, legend = T) css <- credibleShiftSet(edata, expectedNumberOfShifts=1, threshold=5, set.limit = 0.95) plot.credibleshiftset(css) css$number.distinct [1] 15 > summary(css) 95 % credible set of rate shift configurations sampled with BAMM Distinct shift configurations in credible set: 15 Frequency of 9 shift configurations with highest posterior probability: rank probability cumulative Core_shifts 1 0.51415880 0.5141588 1 2 0.16823987 0.6823987 2 3 0.05163798 0.7340366 1 4 0.04997224 0.7840089 1 5 0.02776235 0.8117712 2 6 0.02443087 0.8362021 1 7 0.02387562 0.8600777 1 8 0.01998890 0.8800666 3 9 0.01332593 0.8933926 2 ...omitted 6 additional distinct shift configurations from the credible set. You can access the full set from your credibleshiftset object # Get best rate shift configuration > best <- getBestShiftConfiguration(edata, expectedNumberOfShifts=1) Processing event data from data.frame Discarded as burnin: GENERATIONS < 0 Analyzing 1 samples from posterior Setting recursive sequence on tree... Done with recursive sequence > plot.bammdata(best, lwd = 2, legend = T, labels = T, cex = 0.4) > addBAMMshifts(best, cex=2.5) # Get pdf figure for supplemental methods > pdf(file = "orchard-bees-bamm-best-shift-config.pdf",useDingbats = F) > plot.bammdata(best, lwd = 2, legend = T, labels = T, cex = 0.45) > addBAMMshifts(best, cex=2.5) > dev.off()