Bayes Empirical Bayes Site Selection

Bayes Empirical Bayes (BEB) is a hierarchical Bayesian inference procedure introduced by Ziheng Yang, Wendy S. W. Wong, and Rasmus Nielsen in 2005 for identifying individual amino acid residues experiencing diversifying positive selection (nonsynonymous-to-synonymous rate ratio ω = d_N/d_S > 1) within a protein-coding sequence alignment. In random-sites codon substitution models such as M2a (PositiveSelection) and M8 (beta&ω), a global likelihood ratio test (LRT) first evaluates whether a gene contains a non-zero proportion of codons evolving with ω > 1. When the global LRT is statistically significant, investigators must still determine which specific codon columns in the alignment belong to the positively selected site class.

Whereas the earlier Naive Empirical Bayes (NEB) method (Nielsen & Yang, 1998) plugged single point estimates (maximum-likelihood estimates, θ̂_MLE) of the ω distribution parameters directly into Bayes' theorem as if they were known without error, Bayes Empirical Bayes assigns an explicit prior probability distribution over the parameters of the ω distribution and integrates numerically over their posterior uncertainty. By accounting for finite-sample estimation variance in the site-class proportions and ω values, BEB eliminates the high false-positive rates (>20–40%) that plagued NEB in small or low-divergence sequence alignments.

The Flaw of Point-Estimate Plug-In in Naive Empirical Bayes

Consider a codon alignment X = (x₁, x₂, ..., x_L) of L codon columns across N species, analyzed under a site model with K discrete selective classes having rates ω₀, ω₁, ..., ω_K−1 and prior mixing proportions p₀, p₁, ..., p_K−1 (denoted collectively by parameter vector η). Let θ = (η, b, κ, π) represent the complete set of model parameters, including branch lengths b, transition/transversion ratio κ, and equilibrium codon frequencies π.

In the Naive Empirical Bayes (NEB) approach, the parameter vector θ is first estimated by maximum-likelihood phylogenetics across all L sites to obtain θ̂. Treating θ̂ as fixed and exact, NEB calculates the posterior probability that codon column x_h belongs to the positive-selection class S_K−1 (where ω_K−1 > 1) via Bayes' theorem:

P_NEB(ω_h = ω_K−1 | x_h, θ̂) = [ p̂_K−1 f(x_h | ω̂_K−1) ] / [ ∑_c=0^K−1 p̂_c f(x_h | ω̂_c) ]

Here, the site-specific conditional probability f(x_h | ω_c) is computed by Felsenstein's tree-pruning algorithm. Although NEB is asymptotically valid when the alignment length L is very large and sequence divergence is high, it fails in small datasets or conserved alignments where θ̂ is subject to large sampling errors. For example, if the maximum-likelihood estimate accidentally places a boundary parameter on an extreme value (such as p̂₂ = 0.99 with ω̂₂ = 1.01, or conversely a tiny p̂₂ with huge variance), NEB ignores the flatness of the likelihood surface and reports artificially inflated posterior probabilities (P > 0.99) for neutral sites.

Hierarchical Parameter Priors in Bayes Empirical Bayes

To correct this overconfidence without incurring the heavy computational cost of full Markov chain Monte Carlo (MCMC) sampling over branch lengths and tree topologies, Yang, Wong, and Nielsen formulated Bayes Empirical Bayes (BEB) as a hybrid empirical-hierarchical approach:

1. Nuisance Parameters Fixed at MLEs: Parameters that are shared across all sites and well constrained by the entire alignment—specifically the branch lengths b̂, transition/transversion ratio κ̂, and equilibrium codon frequencies π̂—are fixed at their global maximum-likelihood estimates.
2. Bayesian Integration Over ω-Distribution Parameters (η): The parameters η that directly govern the site-to-site ω distribution (which suffer the largest sampling uncertainty) are assigned flat or symmetric Dirichlet prior distributions f(η) and integrated out over the full alignment data X:

• Under Model M2a (PositiveSelection): η = (p₀, p₁, p₂, ω₀, ω₂), where ω₁ = 1 is fixed. A flat Dirichlet(1, 1, 1) prior is placed on the mixing proportions (p₀, p₁, p₂) over the 2-simplex (p₀ + p₁ + p₂ = 1), while independent uniform priors are placed on the purifying rate ω₀ ~ U(0, 1) and the positive-selection rate ω₂ ~ U(1, 11).
• Under Model M8 (beta&ω): η = (p₀, p, q, ω_s), where a proportion p₀ of sites are drawn from a Beta(p, q) distribution on (0, 1) discretized into 10 equal-probability categories, and a proportion p₁ = 1 − p₀ evolve with rate ω_s ≥ 1. Independent uniform priors are assigned to p₀ ~ U(0, 1), p ~ U(0, 2), q ~ U(0, 2) (or over bounded logarithmic grids), and ω_s ~ U(1, 11).
• Under Branch-Site Model A: A Dirichlet(1, 1, 1, 1) prior is assigned to the four site-class proportions (p₀, p₁, p_2a, p_2b), along with uniform priors on the background purifying rate 0 < ω₀ < 1 and foreground adaptive rate ω₂ > 1.

Marginal Likelihood and Site-Class Posterior Integration

In Bayes Empirical Bayes, the posterior probability that codon site h belongs to site class c (out of K total classes, where class K − 1 corresponds to ω > 1) conditional on the entire alignment X is given by integrating over the parameter space Ω_η:

P_BEB(S_h = c | X) = [1 / f(X)] ∫_Ωη P(S_h = c | x_h, η) · f(X | η) · f(η) dη

where the marginal likelihood of the alignment X is the normalizing constant:

f(X) = ∫_Ωη f(X | η) f(η) dη = ∫_Ωη [ ∏_k=1^L ∑_j=0^K−1 p_j f(x_k | ω_j) ] f(η) dη

Because η contains four parameters under M8 (p₀, p, q, ω_s) and four independent parameters under M2a (p₀, p₁, ω₀, ω₂), the four-dimensional integral cannot be solved analytically. Instead, PAML's codeml evaluates the integral deterministically by discretizing each parameter into a 10-point grid:

• In M2a, the ternary simplex (p₀, p₁, p₂) is partitioned into 100 equilateral triangles (evaluating the centroid of each triangle), while ω₀ ∈ (0, 1) and ω₂ ∈ (1, 11) are each discretized into 10 equal-interval points, producing 100 × 10 × 10 = 10,000 grid points.
• Because the tree-pruning likelihoods f(x_h | ω_j) depend only on the single scalar value ω_j (with b̂, κ̂, and π̂ fixed), codeml pre-calculates f(x_h | ω) at the 21 distinct grid values of ω (10 values for ω₀, 1 value for ω₁ = 1, and 10 values for ω₂) in just 21 passes of Felsenstein's tree-pruning algorithm. Combining those pre-computed site likelihoods across all 10,000 grid points takes a fraction of a second.

Decision Thresholds, Posterior Mean ω, and Alignment Pitfalls

Once the posterior probability distribution P_BEB(S_h = c | X) is computed across all K site classes for codon h, two summary statistics are reported:

1. Probability of Positive Selection, P(ω_h > 1 | X): In standard practice (PAML rst output), codon sites with P(ω > 1) > 0.95 are flagged with a single asterisk (*) as statistically significant candidates for positive selection, and sites with P(ω > 1) > 0.99 are flagged with a double asterisk (**). Extensive Monte Carlo simulations confirm that when P_BEB(ω > 1) ≥ 0.95, the empirical false-positive rate remains below 5% even in alignments of only 6–10 taxa.
2. Posterior Mean and Standard Error of ω_h: The site-specific expected selective ratio E[ω_h | X] and its posterior variance Var(ω_h | X) are calculated by averaging the class rates ω_c weighted by their BEB posterior probabilities:

E[ω_h | X] = ∑_c=0^K−1 ω_c P_BEB(S_h = c | X)

Investigators applying BEB must ensure high upstream data quality: frameshift sequencing errors, non-homologous exon alignment around indel boundaries, or intragenic recombination (such as horizontal gene transfer mosaicism that violates the single-tree assumption) can generate clustered nonsynonymous mismatches that mimic ω > 1 unless filtered prior to codon analysis.

When should you report Bayes Empirical Bayes (BEB) results instead of Naive Empirical Bayes (NEB)?

You should always report Bayes Empirical Bayes (BEB) results rather than NEB when identifying positively selected sites under models M2a, M8, or Branch-Site Model A. When the alignment is large and informative, BEB and NEB converge to identical posterior probabilities; when the alignment is small or moderately conserved, NEB ignores parameter estimation uncertainty and suffers from high false-positive rates, whereas BEB integrates over the parameter prior to maintain reliable false-positive control.

Can you use BEB posterior probabilities if the global Likelihood Ratio Test (M1a vs. M2a or M7 vs. M8) is not significant?

Standard statistical protocol requires a two-stage hierarchical test: first verify that the global Likelihood Ratio Test (such as M1a vs. M2a or M7 vs. M8) rejects the neutral/purifying null model (p < 0.05), and then use BEB to identify individual sites with P(ω > 1) > 0.95. However, even if the null hypothesis holds true across the alignment, simulations by Yang et al. (2005) showed that BEB rarely produces P > 0.95 (false-positive rate < 1–3%), unlike NEB.

Why is the BEB 4D grid integration so fast if it evaluates 10,000 parameter combinations?

Evaluating a phylogenetic likelihood via Felsenstein's pruning algorithm is only needed when the rate matrix Q(ω) changes. Because the mixing proportions (p₀, p₁, p₂) do not appear inside Q(ω), codeml only runs the tree-pruning algorithm across the 21 discrete grid values of ω (under M2a) or 110 values (under M8), and then takes fast algebraic weighted sums over the 10,000 prior grid points.

1. Yang, Z., Wong, W. S. W., & Nielsen, R. (2005). "Bayes Empirical Bayes inference of amino acid sites under positive selection." Molecular Biology and Evolution, 22(4), 1107–1118. https://doi.org/10.1093/molbev/msi097
2. Nielsen, R., & Yang, Z. (1998). "Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene." Genetics, 148(3), 929–936. https://doi.org/10.1093/genetics/148.3.929
3. Zhang, J., Nielsen, R., & Yang, Z. (2005). "Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level." Molecular Biology and Evolution, 22(12), 2472–2479. https://doi.org/10.1093/molbev/msi237
4. Yang, Z. (2014). Molecular Evolution: A Statistical Approach. Oxford University Press. https://doi.org/10.1093/acprof:oso/9780199602605.001.0001