One of the most commonly used methods for identifying heterogeneity in diverisifcation rates across phylogenies is BAMM: Bayesian Analysis of Macroevolutionary Mixtures. One of the goals of BAMM is to identify where shifts in diversification dynamics occur. Within the package BAMMtools this can be done using the function getBestShiftConfiguration. This function returns a list of nodes associated with rate shifts, but these transitions can actually occur earlier than the node of interest. You can see how rate shifts can occur at different times above specific nodes on the BAMM gallery webpage. Therefore, it may be useful to know, for any given analysis, how long before each of these nodes the rate shifts occur. Fortunately, this is simple to achieve using the information output form getBestShiftConfiguration. Here, I present a simple function, time.to.rateshift.node (which can be downloaded here) to summarize this information.
source("time.to.rateshift.node.R")To demonstrate the function, let’s first load some example data generated from BAMM, and estimate where the rate shifts occur.
subset.edata <- readRDS("subset.BAMM.data.rds")
library(BAMMtools)## Loading required package: apelibrary(phytools)## Loading required package: mapsbss <- getBestShiftConfiguration(subset.edata, 30)## Processing event data from data.frame
##
## Discarded as burnin: GENERATIONS < 0
## Analyzing 1 samples from posterior
##
## Setting recursive sequence on tree...
##
## Done with recursive sequenceplot.bammdata(subset.edata, breaksmethod = 'jenks')
nodelabels(node = bss$eventData[[1]]$node, pch = 19, col = "deeppink")
As arguments, the function only needs the output from getBestShiftConfiguration and the phylogeny the data was generated from. Let’s run the function and take a look. The function returns a data.frame of three columns: index is the same index as identified by getBestShiftConfiguration, time.before.node is the amount of time (millions of years) that the rate shift occurs before the specific node, and proportion.parent.branch.length is the time.before.node value divided by the total length of the parent branch.
phy <- as.phylo(subset.edata)
time.to.rateshift.node(bss, phy)## index time.before.node proportion.parent.branch.length
## 1 1 0.00000000 NA
## 2 2 19.25917462 0.433450188
## 3 3 16.32883718 0.385986174
## 4 4 0.30395157 0.036436381
## 5 5 1.75604805 0.600443422
## 6 6 2.25943909 0.456942295
## 7 7 6.24197718 0.849351557
## 8 8 1.10119889 0.608641748
## 9 9 1.81448047 0.555943487
## 10 10 8.59514414 0.901724546
## 11 11 1.69652565 0.057007759
## 12 12 4.91560193 0.815643750
## 13 13 0.18611592 0.179125324
## 14 14 0.94562662 0.237660330
## 15 15 1.37846078 0.191423238
## 16 16 1.16483508 0.252709061
## 17 17 1.53319240 0.089619669
## 18 18 2.80743599 0.487194058
## 19 19 0.49922086 0.068233044
## 20 20 0.86577695 0.141016361
## 21 21 3.74034957 0.813210245
## 22 22 0.35609084 0.348244517
## 23 23 0.22057241 0.090408103
## 24 24 0.25886004 0.021171293
## 25 25 0.72339213 0.228584034
## 26 26 0.41529934 0.059645182
## 27 27 3.33322072 0.393277080
## 28 28 0.82591412 0.926579468
## 29 29 0.05292258 0.008028540
## 30 30 1.23864800 0.365707894
## 31 31 0.16991341 0.009379508
## 32 32 0.17814343 0.089927294
## 33 33 0.04167718 0.009537871
## 34 34 0.20156612 0.101461765These results show that the rateshifts are largely in close proximity to their associated nodes, but some are quite far away. Note, in this case the first node examined is the root node. Therefore, there is no time before this node nor is there a proportion of the parent branch length, as it has none!