Phylogenetic Model.

The s species are related to one another through an unknown phylogeny, represented as a bifurcating tree. Possible phylogenies are denoted τ₁, τ₂, …, τ_B(s), where B(s) is the number of possible phylogenetic trees for s species. We assume a time-reversible model of nucleotide substitution (see below), which means that the phylogenetic trees can be arbitrarily rooted without changing the likelihood. Therefore, we assume unrooted phylogenetic trees throughout this study. The number of unrooted phylogenetic trees for s species is B(s) = (2s − 5)!/[2^s−3(s − 3)!] (28). Each unrooted tree has 2s − 3 branches. Each branch on the tree has a length, ν. The set of branch lengths for the jth tree is denoted ν_j = (ν₁, ν₂, …, ν_2s−3).

We assume that substitutions occur on the phylogenetic tree according to a continuous-time Markov process. The instantaneous rate of change from codon i to codon j is contained in a matrix, Q, the entries of which are given by

where ω is the nonsynonymous/synonymous rate ratio, κ is the transition/transversion rate ratio, and π_j is the stationary frequency of codon j (2, 14, 15). The rate matrix, Q, is 61 × 61 in dimension; although there are 64 possible codon triplets, the three termination codons of the universal code are excluded from the state space, leaving the 61 sense codons as the states of the substitution process. The diagonal elements of the rate matrix Q are specified such that each row sums to 0. We rescale the matrix such that the average rate of synonymous substitution is 1; this means that the branch lengths are in terms of expected number of synonymous substitutions per codon site. This scaling is different from the usual for phylogenies, where the rate matrix is rescaled such that the branch lengths are in terms of expected number of substitutions per site. Finally, note that the substitution process is time-reversible because π_iq_ij = π_jq_ji for all i ≠ j. Practically speaking, this means that we can use unrooted phylogenies instead of rooted trees when calculating the likelihood on the tree (29). The transition probabilities of the process are the probability of finding the process in codon j conditioned on the process starting in codon i and run over a branch of length ν. The transition probabilities can be calculated as P(ν) = e^Qν.

Each of the c sites can be in one of k selection categories, where each category differs in its d_N/d_S value. The information on the category membership is contained in an association vector, z. For example, for the alignment of four sites given above, one possible assignment of sites to selection categories is z = (1, 1, 2, 3). Here, k = 3, and the first two sites are in the same category. The allocation vector describes a partition of the sites into different classes, each class with its own d_N/d_S value. The number of ways to partition a set of c sites into k classes is described by the Stirling numbers of the second kind (30):

The sites for the example alignment of c = 4 sites can be assigned to k = 3 categories in a total of S(4, 3) = 6 ways: z = (1, 1, 2, 3); z = (1, 2, 1, 3); z = (1, 2, 3, 1); z = (1, 2, 2, 3); z = (1, 2, 3, 2); and z = (1, 2, 3, 3). We label possible partitions of the sites into classes by using the restricted growth function notation of Stanton and White (31), where elements are sequentially numbered with the constraint that the index numbers for two sites are the same if they are found in the same category. The number of sites in the ith category is denoted η_i. For the example given above, η₁ = 2, η₂ = 1, and η₃ = 1.

Likelihood.

Assuming that substitutions are independent across sites, the likelihood for an alignment can be calculated as the product of the site likelihoods:

We used the Felsenstein (29) pruning algorithm to calculate likelihoods for specific combinations of parameter values.

Statistical Model.

We estimate parameters in a Bayesian framework. All parameters of the model are considered random values with their own prior and posterior probability distributions. The joint posterior probability of all of the parameters of the model is

The posterior probability [f(· |X)] is equal to the product of the likelihood [(·) or, equivalently, f(X| ·)] and the prior probability distribution [f(·)] of the parameters divided by the marginal likelihood [f(X)].

Bayesian analysis requires that a prior probability distribution be specified for the parameters. Here we assume that the parameters take the prior probability distributions shown in Table 8.

The model, as described, is fairly typical of codon models implemented to date (32). However, the Dirichlet process prior is new to the problem of detecting positive natural selection and requires more explanation. The Dirichlet process prior is widely used in clustering problems as a nonparametric alternative, where the data elements are drawn from a mixture distribution with a countably infinite number of components (33, 34). [Ewens and Tavaré (35) describe the relationship between the Dirichlet process prior and the Ewen’s sampling formula from population genetics.] The Dirichlet process prior has two parameters: a concentration parameter, which, loosely speaking, determines the degree to which data elements are grouped together, and a base measure, G₀(·), which describes the probability distribution of the parameter value associated with each cluster. The model can be described as follows: First, the number of classes, k, and an allocation variable, z, are randomly drawn from the probability distribution

Once the sites have been partitioned into k categories, a ω = d_N/d_S value is assigned to each category by drawing randomly from the following probability distribution k times, G₀(ω_i) = f(ω_i) = [1/(1 + ω_i)²] which is the probability density of the ratio of two identically distributed exponential random variables. Here, ω_i is the d_N/d_S rate ratio for the ith selection category. The probability of having k categories is obtained by summing over all possible partitions for k categories, and is

where _na_k is the absolute values of the Stirling numbers of the first kind. The expected number of categories is

Finally, the probability of finding two sites grouped together into the same cluster is f(z_i = z_j|α, c) = 1/(1 + α)].

MCMC.

We approximated posterior probabilities by using MCMC (36, 37). The general idea is to construct a Markov chain that has as its state space the parameters of the statistical model and a stationary distribution that is the posterior probability of the parameters. Samples drawn from the Markov chain when at stationarity are valid, albeit dependent, draws from the posterior probability distribution of the parameters (38). The fraction of the time the Markov chain visits any particular combination of parameter values is an approximation of the joint posterior probability of those parameters.

We wrote a computer program that implements the MCMC method for the phylogenetic model described above. For each iteration of the Markov chain, the program randomly picks a parameter to change, proposes a new value for the parameter, and then accepts or rejects the proposed state as the next state of the chain by using the Metropolis–Hastings formula. Most of the proposal mechanisms have been described elsewhere. We use the LOCAL mechanism of Larget and Simon (39) to simultaneously change the tree and branch lengths (τ, ν). We update the codon frequencies (π) by using a Dirichlet proposal mechanism and change the transition/transversion rate ratio (κ) and the category d_N/d_S values (ω) by using a rate multiplier [i.e., we multiply the old parameter value by e^λ(u−½) to obtain a new value for the parameter, where u is a uniform(0, 1) random number and λ is a tuning parameter]. Finally, we update the allocation vector by using a Gibbs sampling algorithm for nonconjugate models proposed by Neal (40). We implement the method of Neal (40) as follows: First, we randomly choose a site, i, and remove it from the set of k selection classes currently in computer memory. If the site is in a selection class by itself (i.e., η_zi for the class was 1), then the entire class is deleted, and k is decreased by 1. Otherwise, η_zi for the selection class that the site is associated with is decreased by 1. We then construct five new (temporary auxiliary) classes by drawing the d_N/d_S value for each from the prior probability distribution. These auxiliary classes represent components that have no other sites associated with them. Gibbs sampling is performed to update site i but with respect to the distribution that includes the temporary auxiliary parameters. The site is assigned to the jth previously existing class with probability Zη_−if(x_i|ω_j) and to the jth of five auxiliary classes with probability Z(α/5)f(x_i|ω_j), where Z is the appropriate normalizing constant. After this update, all d_N/d_S values not associated with a selection class are discarded.

All analyses were run for a total of 2 million update cycles. Chains were thinned, with samples taken every 100 cycles. Samples taken during the first 100,000 cycles were discarded as the burn-in phase. The MCMC procedure was repeated, resulting in two sets of samples for each gene. Convergence was assessed by using the program tracer (http://evolve.zoo.ox.ac.uk/software.html?id=tracer).

J.P.H. was supported by National Institutes of Health Grant GM-069801 and National Science Foundation Grants DEB-0445453 and CACTTOL-22035A.