OBAMA: OBAMA for Bayesian amino-acid model averaging

Site model

We consider phylogenetic site models consisting of four parts:

• a symmetric rate matrix Q of size 20 × 20 specifying the rate of substitution among amino acids. During the MCMC, we average over the 15 empirical substitution models listed in Table 1, using a model indicator I_M that identifies the model by an integer number.

• a frequency model F specifying the base frequencies of the 20 amino acids. Empirical models come with empirical frequencies, but we can also estimate them using MCMC.

• a model of rate variation over sites. We consider no variation at all (homogenous rates) and gamma distributed heterogeneous rates using 4 categories (Yang, 1994) with gamma shape parameter α estimated.

• the choice whether to include a category of invariable sites and estimate the proportion of invariable sites p_inv or not.

The rate matrix Q and frequency model F together specify a rate matrix R = Q × F that determine a transition probability matrix P = e^Rt, which can be used to calculate the likelihood of a tree given the alignment efficiently using Felsenstein’s pruning algorithm (Felsenstein, 1981).

Prior

For the prior on frequencies, we collected all frequencies from the 15 empirical models, and Fig. 1A shows a histogram of these 15 × 20 = 300 frequencies. The range of frequencies is from 0.006 to 0.169 with mean 0.050001. A Dirichlet prior is appropriate for a set of variables that is constrained to sum to 1. Figure 1B shows a sample of 2000 cases from a Dirichlet(4,4, …,4), which has range for samples from 0.005 to 0.167 with mean 1/20=0.05, and this statistic resembles the empirical distribution for frequencies.

For bModelTest (Bouckaert & Drummond, 2017), we observed some lack of identifiability when both gamma rate heterogeneity and invariable sites were included in the true site model used to simluate data. This was earlier observed for maximum likelihood estimates of these parameters (p. 120, Yang (2014)). This phenomenon is particularly problematic when the shape parameter α governing gamma rate heterogeneity is small and a large number of invariable sites can be expected. To demonstrate this, Fig. 2 shows the probability of observing an invariable site for the set of trees (from our simulation study described in detail later) and substitution models with gamma rate heterogeneity using 4 categories, but no invariable sites. This is calculated using Felsenstein’s algorithm as the probability that each of the tips in a tree have exactly the same amino acid, summed over all 20 amino acids.

When we look at how rates behave with small α and 4 categories (Fig. 2), we notice that when α < 0.2 the slowest category has a rate less than 0.00015. This implies for a branch of length 1 (i.e., where the expected number of substitutions per site is 1) that for every 100 sites we expect more than 99 sites to be invariable under such rate. When α < 0.1 the second slowest rate decreases to 0.007 and now for the second rate we expect over 99 out of 100 sites to be invariable as well, while the remaining rate heterogeneity must be captured by the two fastest rates. At α < 0.05, the third rate decreases below 0.005, and almost all rate heterogeneity has to be captured by the remaining rate, which with decreasing α becomes 4 − ϵ with ϵ < 0.001. In this situation, the probability of invariable sites becomes 0.75 in the case of 4 categories.

This does not seem to be reasonable behaviour for this model, so our prior on α should exclude α below 0.1. Empirical results (Table 2.2 in Yang (2014)) appear to fall in a range above 0.1. Even at α = 0.1 there is still a high probability of invariable sites, so the identifiability problem is not completely solved, but at least somewhat diminished. Since we would like gamma rate heterogeneity only to be included in the model if there is a reasonable amount of heterogeneity in the data (which coincides with low values of α), by default we use an exponential prior with rate 1, truncated to exclude values below 0.1.

For the other items we chose uninformative priors. In summary, the prior consists of the product of the following independent priors:

• a uniform prior on empirical models through a uniform prior on the model indicator I_M.

• a uniform prior on using empirical frequencies (from the empirical model indicated by I_M) or estimated frequencies.

• a uniform prior on including gamma rate heterogeneity or not.

• a uniform prior on including a proportion of invariable sites or not.

• a Beta(1,4) prior on the proportion of invariable sites. This skews the proportion of invariable sites rather to lower than to higher proportions of invariable sites.

• a Dirichlet(4,4,…,4) prior on estimated frequencies.

• an exponential prior on the shape parameter for gamma rate heterogeneity with rate 1 and lower bound 0.1.

MCMC proposals

For switching between substitution models, we use a uniform proposal for the model indicator I_M, which uniformly selects an indicator value in the range of 1 to the number of empirical subsitution models included in the analysis.

For gamma rate heterogeneity and proportion of invariable sites we use the same birth and death operators as for bModelTest. For gamma rate heterogeneity, the birth and death proposal sets or unsets the category count flag and samples a new value for shape parameter α from the prior (by default an exponential with lower bound 0.1) when the flag is set. The proposal ratio is d(α′) for birth and 1∕d(α) for death where d(.) is the density used to sample from. Likewise, for setting the proportion of invariable sites flag and sampling p_inv from the prior (by default a Beta(1,4)). The Jacobian determinant is 1 for these proposals, so have no impact on the Hastings ratio. For moving α, we use the standard scale operator in BEAST 2 (Bouckaert et al., 2014), adapted so it only samples if the category count flag is set for α. In the same way, for p_inv such a scale operator is used only if the proportion of invariable sites flag is set.

The estimation of frequencies is done through Bayesian stochastic variable selection (Kuo & Mallick, 1998; Lemey et al., 2009), where we use a bit flip operator on a boolean parameter to switch between empirical model frequencies and estimated frequencies.

To estimate frequencies, efficiency of operators for sampling was conducted by comparing two operators: the delta exchange operator and the adaptive variance multivariate normal (AVMN) operator (Baele et al., 2017). A delta exchange operator randomly selects two frequencies f_a and f_b, then randomly selects a value δ and proposes new frequencies $f_a^{\prime }=f_a+\delta$, $f_b^{\prime }=f_b-\delta$ under the condition that $0\le f_a^{\prime }\le 1$ and $0\le f_b^{\prime }\le 1$. An AVMN operator parameterises the space spanned by frequencies as a log-constrained sum transformed multivariate parameter and approximates the distribution of that transformed parameter using a multivariate normal distribution N(M, Σ) for which the mean vector M and covariance matrix Σ are learned during the MCMC run (see Baele et al. (2017) for more details). We use the settings recommended in Baele et al. (2017) for the AVMN operators in our experiments.

Validation of model implementation

We performed a well calibrated simulation study to verify correctness of the implementation and investigated the behaviour of OBAMA on a published alignment.

Identifiability of gamma shape and proportion of invariable sites

When both gamma rate heterogeneity and invariable sites are in the true model under an exponential prior with lower bound zero (see Table S4), the estimates are slightly diminished due to lack of identifiability as explained before and observed by Bouckaert & Drummond (2017). However, when a lower bound of 0.1 is used as motivated by our discussion, identifiability increases as witnessed by the expected coverage of gamma shape and proportion of invariable sites estimates in Table 2.

Furthermore, as shown in Fig. 3, the estimates of these parameters are informed by the data and not only the prior: when values used to simulate the data increase, their estimates increase as well.

Frequency operator

The well calibrated study from Table 2 uses the AVMN operator. The simulation experiment was repeated with delta exchange and AVMN operator, both with a weight of 1, giving them equal probability of being selected, and all other operators unaltered. Run times for both experiments were very similar, even though the AVMN operator performs a relatively large amount of work to produce a proposal. The effective sample sizes (ESSs) for 100 runs for the simulation experiment are shown in Fig. 4A, showing the AVMN operator outperforming the delta exchange operator: in all cases the average ESS increases, though the variance of ESSs from MCMC is usually quite high.

Site models matter

To determine the effectiveness of the OBAMA model, we investigated an amino acid alignment, M200 from TreeBase¹ (Simmons et al., 2002). This alignment for 20 flowering plant taxa is 466 characters long. ProtTest (Darriba et al., 2011) (run with settings -all-matrices -all-distributions -F -threads 2) suggests JTT+G+I+F fits best (where +G stands for using gamma rate heterogeneity, +I using invariable sites, and +F using observed frequencies), followed by WAG+G+I+F and LG+G+I+F as third runner up, both based on the BIC and the Ln criteria.

For a Bayesian analysis, a relaxed log normal clock (Drummond et al., 2006) can be rejected based on a coefficient of variation being distributed with mean less than 0.1 skewed towards zero. So, we use a strict clock in the remainder, a Yule prior with birth rate estimated and vary the site model using the three models suggested by ProtTest. All runs were verified to have ESSs exceeding 200 and four independent runs converging to the same distribution. Data files used can be found at https://github.com/rbouckaert/obama/releases/tag/data. When we compare the posterior clade support for these models (Figs. 5A–5C), we see considerable differences. The largest clade support difference between WAG and JTT is 17.3%, while between JTT and LG as well as WAG and LG these exceed 20%. As comparison, two independent runs of OBAMA only show a difference of just 6.7% in clade support. This indicates the three different models sample from different distributions and the choice of substitution model matters for the analysis.

With the OBAMA model, about 84% of the time JTT is sampled and 16% WAG while the other substitution models have no recorded posterior support. Gamma rate heterogeneity, proportion of invariable sites and estimated frequencies were supported 100% of the time. Comparing the three individual substitution models with OBAMA (Figs. 5E and 5F) shows that JTT and WAG have differences in clade support of 11.0% and 11.8% respectively, which shows that model averaging has an effect on the posterior distribution, and matters as well.