Evolutionary Rates & Selection Mode (BEAST2)

This vignette explains how to extract evolutionary rate parameters estimated from relaxed clock Bayesian inference analyses produced by BEAST2. It also shows how to use evolutionary rate based inference of selection mode (strength) adapted to clock-based rates, as introduced by Simões and Pierce (2021). See the sister vignette “Evolutionary Rates & Selection Mode (BEAST2)” for an equivalent workflow using output data produced by Mr. Bayes.

Load the EvoPhylo package

library(EvoPhylo)

Evolutionary Rates Statistics and Plots

In this section, we will extract evolutionary rate parameters from each node from a Bayesian clock (time-calibrated) summary tree. The functions below will store them in a data frame, produce summary statistics tables, and create different plots showing how rates are distributed across morphological partitions and clades. Note that step 0 needs to be performed before running the inference when using BEAST2.

0. Set your BEAST2 analysis to log multiple clocks

This step is specific to BEAST2 analyses - if using Mr. Bayes, skip directly to step 1.

When using a default configuration generated through BEAUti, only one of the clocks will be logged along the tree. For the other clocks, only summary statistics (such as the mean rate) will be logged. In order to obtain the clock rates for different partitions, we need to modify the setup before running the inference.

In order to do this, open the BEAST2 XML configuration file and find the section which logs the tree, which should look similar to this:

<logger id="treelog.t:26" spec="Logger" fileName="multiclock.trees" logEvery="10" mode="tree">
  <log id="TreeWithMetaDataLogger.t:26" spec="beast.evolution.tree.TreeWithMetaDataLogger" branchratemodel="@RelaxedClock.c:26" tree="@Tree.t:26"/>
</logger>

This section needs to be duplicated for each additional clock, modifying the name of the clock (branchratemodel), the file name (fileName) and the name of the logger (id) like in this example:

<logger id="treelog.t:26.2" spec="Logger" fileName="multiclock.clock2.trees" logEvery="10" mode="tree">
  <log id="TreeWithMetaDataLogger.t:26" spec="beast.evolution.tree.TreeWithMetaDataLogger" branchratemodel="@ExponentialRelaxedClock.c:261" tree="@Tree.t:26"/>
</logger>

The clock name needs to reference the name (id) of one of the clock models set up in the analysis, which can be found earlier in the XML file.

1. Get rates from the clock tree and create a rate table

Here we assume that a consensus or summary tree has been generated for each of the tree log files generated by the inference configured in step 0, for instance by using the software TreeAnnotator. First, import these summary clock trees using treeio’s function read.beast().

## Import all clock summary trees (.tre, .t, or .tree) produced by BEAST2 from your local directory

tree_clock1 <- treeio::read.beast("tree_clock1.tre")
tree_clock2 <- treeio::read.beast("tree_clock2.tre")
#etc

Below, we use the example BEAST2 clock trees from 3 morphological partitions and one molecular clock partition that accompany EvoPhylo.

tree_clock1 <- system.file("extdata", "Penguins_MCC_morpho_part1.t", package = "EvoPhylo")
tree_clock2 <- system.file("extdata", "Penguins_MCC_morpho_part2.t", package = "EvoPhylo")
tree_clock3 <- system.file("extdata", "Penguins_MCC_morpho_part3.t", package = "EvoPhylo")
tree_clock4 <- system.file("extdata", "Penguins_MCC_dna.t", package = "EvoPhylo")
tree_clock1 <- treeio::read.beast(tree_clock1)
tree_clock2 <- treeio::read.beast(tree_clock2)
tree_clock3 <- treeio::read.beast(tree_clock3)
tree_clock4 <- treeio::read.beast(tree_clock4)

Subsequently, using get_clockrate_table_BEAST2(), users can extract mean or median rate values for each node in the summary tree. These mean or median rate values are calculated by TreeAnnotator taking into account all trees from the posterior sample. Please note that analyses must have reached the stationarity phase and independent runs converging for the summary statistics in each node to be meaningful summaries of the posterior sample.

The example shown here uses two different clocks, however the function supports any number of clocks which can all be passed as separate arguments.

## Get table of clock rates with summary stats for each node in 
## the tree for each relaxed clock partition (4 partitions in this dataset) 
RateTable_Medians_4p <- get_clockrate_table_BEAST2(tree_clock1, tree_clock2, tree_clock3, tree_clock4, summary = "median")
RateTable_Means_4p <- get_clockrate_table_BEAST2(tree_clock1, tree_clock2, tree_clock3, tree_clock4, summary = "mean")

2. Export the rate table and plot tree with node values

Once a rate table has been obtained from BEAST2 files, it is necessary to export it. This is a necessary step to subsequently open the rate table spreadsheet locally (e.g., using Microsoft Office Excel) and customize the table with clade names associated with with each node in the tree for downstream analyses. Note that the root node may include “NA” for rate value, so it must be removed from the rate table.

## Export the rate tables 
write.csv(RateTable_Means_4p, file = "RateTable_Means_4p.csv")
write.csv(RateTable_Medians_4p, file = "RateTable_Medians_4p.csv")

Plot tree node labels to customize clade names

To visualize the node values in the tree, you can use ggtree(). This can be done on any of the imported trees, as all BEAST2 Summary trees should have the same topology and divergence times, differing only on rate parameters.

## Plot tree node labels
library(ggtree)
tree_nodes <- ggtree(tree_clock1, branch.length = "none", size = 0.05) +
  geom_tiplab(size = 2, linesize = 0.01, color = "black", offset = 0.5) +
  geom_label(aes(label = node), size = 2, color="purple")
tree_nodes

## Save your plot to your working directory as a PDF
ggplot2::ggsave("Tree_nodes.pdf", width = 10, height = 10)

3. Import rate table with custom clade memberships

A new “clade” column has been added to the rates table. Below, we use the rate table with clade membership RateTable_Means_4p_Clades that accompanies EvoPhylo (clade names in this example are for demonstration purposes only).

## Import rate table with clade membership (new "clade" column added) 
## from your local directory
RateTable_Medians_Clades <- system.file("extdata", "RateTable_Medians_Clades.csv", package = "EvoPhylo")
RateTable_Medians_Clades <- read.csv(RateTable_Medians_Clades, header = TRUE)

head(RateTable_Medians_Clades, 5)
##               clade nodes     rates1     rates2     rates3     rates4
## 1 EarlyStemPenguins     1 0.01059049 0.01142373 0.04521063 0.03214791
## 2 EarlyStemPenguins     2 0.01127067 0.01117188 0.04558499 0.03297338
## 3 EarlyStemPenguins    58 0.01252320 0.01143236 0.04487146 0.03119615
## 4     CrownPenguins     3 0.01493279 0.01009941 0.04626338 0.04521892
## 5     CrownPenguins     4 0.01050220 0.01008585 0.04537996 0.03908424

4. Get summary stats for each clade/clock partition

Obtain summary statistics table and plots for each clade by clock partition using clockrate_summary(). Supplying a file path to file save the output to that file.

## Get summary statistics table for each clade by clock 
clockrate_summary(RateTable_Medians_Clades,
                  file = "Sum_RateTable_Medians_4p.csv")

Rate table summary statistics
clade	clock	n	mean	sd	min	Q1	median	Q3	max
CrownPenguins	1	47	0.012	0.004	0.006	0.009	0.011	0.013	0.024
EarlyStemPenguins	1	39	0.011	0.004	0.008	0.010	0.011	0.011	0.035
LateStemPenguins	1	22	0.015	0.015	0.007	0.010	0.011	0.012	0.081
CrownPenguins	2	47	0.011	0.001	0.010	0.011	0.011	0.012	0.013
EarlyStemPenguins	2	39	0.012	0.001	0.010	0.011	0.011	0.012	0.013
LateStemPenguins	2	22	0.012	0.001	0.011	0.011	0.011	0.012	0.014
CrownPenguins	3	47	0.045	0.001	0.043	0.045	0.045	0.046	0.051
EarlyStemPenguins	3	39	0.046	0.003	0.038	0.045	0.045	0.046	0.062
LateStemPenguins	3	22	0.045	0.002	0.040	0.045	0.045	0.046	0.052
CrownPenguins	4	47	0.033	0.006	0.019	0.029	0.032	0.034	0.058
EarlyStemPenguins	4	39	0.034	0.004	0.026	0.032	0.033	0.033	0.047
LateStemPenguins	4	22	0.032	0.002	0.025	0.032	0.032	0.033	0.036

5. Plot rates by clock partition and clade

Plot distributions of rates by clock partition and clade with clockrate_dens_plot().

## Overlapping plots
clockrate_dens_plot(RateTable_Medians_Clades, stack = FALSE,
                    nrow = 1, scales = "fixed")

Sometimes using stacked plots provides a better visualization as it avoids overlapping distributions.

## Stacked plots
clockrate_dens_plot(RateTable_Medians_Clades, stack = TRUE,
                    nrow = 1, scales = "fixed")

It is also possible to append extra layers using ggplot2 function, such as for changing the color scale. Below, we change the color scale to be the Viridis scale.

## Stacked plots with viridis color scale
clockrate_dens_plot(RateTable_Medians_Clades, stack = TRUE,
                    nrow = 1, scales = "fixed") +
  ggplot2::scale_color_viridis_d() +
  ggplot2::scale_fill_viridis_d()

6. Rate linear models

We can also plot linear model regressions between rates from two or more clocks with clockrate_reg_plot().

## Plot regressions of rates from multiple clocks
#Morpho-morpho rates
p1<- clockrate_reg_plot(RateTable_Medians_Clades, clock_x = 1, clock_y = 2)
p2<- clockrate_reg_plot(RateTable_Medians_Clades, clock_x = 1, clock_y = 3)
p3<- clockrate_reg_plot(RateTable_Medians_Clades, clock_x = 2, clock_y = 3)

#Morpho-Mol rates
p4<- clockrate_reg_plot(RateTable_Medians_Clades, clock_x = 1, clock_y = 4)
p5<- clockrate_reg_plot(RateTable_Medians_Clades, clock_x = 2, clock_y = 4)
p6<- clockrate_reg_plot(RateTable_Medians_Clades, clock_x = 3, clock_y = 4)


library(patchwork) #for combining plots
p1 + p2 + p3 + p4 + p5 + p6 + plot_layout(nrow = 2)

## Save your plot to your working directory as a PDF
ggplot2::ggsave("Plot_regs.pdf", width = 8, height = 8)

Selection Mode

In this section, we will use evolutionary rate based inference of selection mode, as first introduced by Baker et al. (2016) for continuous traits, and later adapted to clock-based rates by Simões and Pierce (2021).

Although the dataset used in the example here includes both morphological and molecular partitions, we will focus the selection mode inference using morphological partitions only. There are multiple available methods to infer the strength of selection using molecular data (e.g., dn/ds ratio) that are better suited to the intrinsic properties of molecular data, or which take into account additional variables available for extant taxa (generation times), and we refer users to those approaches.

1. Import posterior log file

The combined posterior log file from BEAST2 is outputted by LogCombiner from their software package. Our own function to combined log files combine_log is intended to work with Mr. Bayes posterior files only.

When using more two or more clock partitions, posterior log files from BEAST2 will have the mean evolutionary rates for each partition labeled as “rate.mean”. To be read and analyzed by Evophylo, users must ensure that all clock partitions that wish to be analyzed should include the word “part” just before the partition number in this label, as outputted by the Evophylo’s function write_partitioned_alignments. So, this label should read “ratepart.mean”. In our example provided here, with filename “KairukuEA2015”, the first morphological partition was exported with write_partitioned_alignments as a Nexus file labeled “KairukuEA2015_Morpho_3p_part1”, and so the mean evolutionary rates for the first partition is labeled “rate.KairukuEA2015_Morpho_3p_part1.mean”.

Below, we use the posterior dataset “Penguins_log.log” that accompanies EvoPhylo.

posterior <- system.file("extdata", "Penguins_log.log", package = "EvoPhylo")
posterior <- read.table(posterior, header = TRUE)

## Show first 10 lines of combined log file
head(posterior, 10)

###2. Check background rates distribution and if they need transformation

The output includes histograms showing the data distribution before and after data transformation for comparisons.

3. Plot selection gradient on the summary tree

Using different thresholds, Identify the selection mode across branches in the tree for each clock partition with plot_treerates_sgn().

Users must indicate the type of output file (between Mr. Bayes and BEAST2) and whether they would like to log transform the background rate to meet assumptions of normally distributed data, based on the results obtained from plot_back_rates. Users should also indicate in “clock” the number of the clock partition they would like to plot rates from and the desired significance threshold to interpret branch rates (we recommend number of standard deviations around the mean of background rates).Finally, a series of arguments enable users to customize the geological timescale to add to the tree.

## Plot tree using various thresholds for clock partition 1
A1<-plot_treerates_sgn(
  type = "BEAST2", trans = "log10",                         #Indicates software name output and type of transformation
  tree_clock1, posterior,                                   #Calls tree for clock partition 1 and posterior log file
  clock = 1,                                                #Use background rates for clock partition 1
  threshold = c("1 SD", "3 SD"),                            #sets threshold for selection mode
  branch_size = 1.5, tip_size = 3,                          #sets size for tree elements
  xlim = c(-70, 30), nbreaks = 8, geo_size = list(3, 3))    #sets limits and breaks for geoscale
A1

Plot tree using various thresholds for the other clock partitions and combine them.

## Plot tree using various thresholds for other clock partition and combine them
A2<-plot_treerates_sgn(
  type = "BEAST2", trans = "log10",                         #Indicates software name output and type of transformation
  tree_clock2, posterior,                                   #Calls tree for clock partition 2 and posterior log file
  clock = 2,                                                #Use background rates for clock partition 3
  threshold = c("1 SD", "3 SD"),                            #sets threshold for selection mode
  branch_size = 1.5, tip_size = 3,                          #sets size for tree elements
  xlim = c(-70, 30), nbreaks = 8, geo_size = list(3, 3))    #sets limits and breaks for geoscale


A3<-plot_treerates_sgn(
  type = "BEAST2", trans = "log10",                         #Indicates software name output and type of transformation
  tree_clock3, posterior,                                   #Calls tree for clock partition 3 and posterior log file
  clock = 3,                                                #Use background rates for clock partition 3
  threshold = c("1 SD", "3 SD"),                            #sets threshold for selection mode
  branch_size = 1.5, tip_size = 3,                          #sets size for tree elements
  xlim = c(-70, 30), nbreaks = 8, geo_size = list(3, 3))    #sets limits and breaks for geoscale


library(patchwork)
A1 + A2 + A3 + plot_layout(nrow = 1)

## Save your plot to your working directory as a PDF
ggplot2::ggsave("Tree_Sel_Morpho_3p.pdf", width = 18, height = 8)

4. Pairwise t-tests of rate values among clades

The function get_pwt_rates_BEAST2() complements the functionality of plot_treerates_sgn by producing a table of pairwise t-tests for differences between the mean background rate in the posterior and the absolute rate for each summary tree branches Should be used only for normally distributed data in which a CI=0.95 is considered a good threshold. In many cases, however, using multiple standard deviations as outputted using plot_treerates_sgn provides a more robust test of whether branch rates are significantly different from background rates.

4.1. Import rate table with clade membership (new “clade” column added) from your local directory with “mean” values

Below, we use the rate table with clade membership “RateTable_Means_Clades.csv” that accompanies EvoPhylo.

## Import rate table with clade membership (new "clade" column added) 
## from your local directory with "mean" values
RateTable_Means_Clades <- system.file("extdata", "RateTable_Means_Clades.csv", package = "EvoPhylo")
RateTable_Means_Clades <- read.csv(RateTable_Means_Clades, header = TRUE)

head(RateTable_Means_Clades, 5)
##               clade nodes     rates1     rates2     rates3     rates4
## 1 EarlyStemPenguins     1 0.01680336 0.01217294 0.04894430 0.03657966
## 2 EarlyStemPenguins     2 0.01782966 0.01181350 0.04985747 0.03746188
## 3 EarlyStemPenguins    58 0.02404947 0.01263853 0.04995735 0.03741211
## 4     CrownPenguins     3 0.01913400 0.01017825 0.05145375 0.04892464
## 5     CrownPenguins     4 0.01319288 0.01016964 0.04996688 0.04325669

4.2. Get and export table of pairwise t-tests

## Get table of pairwise t-tests for difference between the posterior 
## mean and the rate for each tree node
RateSign_Tests <- get_pwt_rates_BEAST2(RateTable_Means_Clades, posterior)

Combined log file
clade	node	clock	background.rate.mean	absolute.rate.mean	relative.rate.mean	p.value
EarlyStemPenguins	43	1	0.0260902	0.0156937	0.6015154	0e+00
LateStemPenguins	65	1	0.0260902	0.0215020	0.8241390	0e+00
LateStemPenguins	37	1	0.0260902	0.0153367	0.5878329	0e+00
LateStemPenguins	64	1	0.0260902	0.0450708	1.7274962	0e+00
EarlyStemPenguins	38	1	0.0260902	0.0175102	0.6711403	0e+00
EarlyStemPenguins	39	1	0.0260902	0.0161433	0.6187486	0e+00
EarlyStemPenguins	99	1	0.0260902	0.0189790	0.7274382	0e+00
EarlyStemPenguins	40	1	0.0260902	0.0178498	0.6841549	0e+00
EarlyStemPenguins	98	1	0.0260902	0.0130211	0.4990813	0e+00
EarlyStemPenguins	63	1	0.0260902	0.0223161	0.8553429	5e-07

## Export the table
write.csv(RateSign_Tests, file = "RateSign_Tests.csv")

References

Baker, Joanna, Andrew Meade, Mark Pagel, and Chris Venditti. 2016. “Positive Phenotypic Selection Inferred from Phylogenies.” Biological Journal of the Linnean Society 118 (1): 95–115.

Simões, Tiago R., and Stephanie E. Pierce. 2021. “Sustained High Rates of Morphological Evolution During the Rise of Tetrapods.” Nature Ecology & Evolution 5 (10): 1403–14. https://doi.org/10.1038/s41559-021-01532-x.