Joint inference of microsatellite mutation models, population history and genealogies using transdimensional Markov Chain Monte Carlo.

We provide a framework for Bayesian coalescent inference from microsatellite data that enables inference of population history parameters averaged over microsatellite mutation models. To achieve this we first implemented a rich family of microsatellite mutation models and related components in the software package BEAST. BEAST is a powerful tool that performs Bayesian MCMC analysis on molecular data to make coalescent and evolutionary inferences. Our implementation permits the application of existing nonparametric methods to microsatellite data. The implemented microsatellite models are based on the replication slippage mechanism and focus on three properties of microsatellite mutation: length dependency of mutation rate, mutational bias toward expansion or contraction, and number of repeat units changed in a single mutation event. We develop a new model that facilitates microsatellite model averaging and Bayesian model selection by transdimensional MCMC. With Bayesian model averaging, the posterior distributions of population history parameters are integrated across a set of microsatellite models and thus account for model uncertainty. Simulated data are used to evaluate our method in terms of accuracy and precision of estimation and also identification of the true mutation model. Finally we apply our method to a red colobus monkey data set as an example.

MICROSATELLITES, also called short tandem repeats (STRs), are repetitions of a DNA sequence motif with length between 1 and 6 bp. Because they are abundant, widely distributed in the genome, and highly polymorphic, microsatellites have become one of the most popular genetic markers for making inferences on molecular evolution and population genetics (S; S).

Unequal crossing over (S; R) and replication slippage (L) are the two main mechanisms proposed that potentially provide an explanation for the high mutation rate of microsatellites. The study by L using Escherichia coli showed that replication slippage is the predominant mutation mechanism of microsatellite DNA. Replication slippage occurs when the replicating strand and the template strand disassociate and then realign out of register, forming a loop in one of the strands. If the process of replication continues, a loop formed by the replicating strand gives rise to an insertion while that by the template strand results in a deletion.

The simplest microsatellite model is the stepwise mutation model (SMM) proposed by O, which states that the length of the microsatellite increases or decreases by 1 repeat unit at a rate independent of the microsatellite length. Although the SMM has been employed to devise commonly used statistics for estimating genetic divergence (S) and effective population size (W), the model has some drawbacks. Under the SMM, there is no stationary distribution and under this process the repeat length will eventually grow arbitrarily long, which is inconsistent with empirical microsatellite length distributions from genomic data (K). Moreover the SMM ignores various properties of microsatellite mutation that have been observed in empirical data. Many different models have been developed in attempts to capture some of these properties.

Observations from many studies support the fact that longer microsatellites have a higher mutation rate (G; W; S). A longer microsatellite allele has more locations for potential slippage errors and hence possesses a greater chance of experiencing a mutation event during replication, as demonstrated by S using bacteriophage T4. This is the motivation behind rate-dependent models such as the proportional slippage model (K) and others (C; S), which describe the mutation rate as a polynomial function of length in repeat units.

Another property is mutational bias, which exists when the probability of expansion and contraction is unequal for a mutation event. Evidence for this phenomenon has been found in genomes of several species including humans (R), which exhibit a preference for expansion, and the bacterium Mycoplasma gallisepticum (M), which tends to contract. Models proposed by C and W have accounted for this rate asymmetry (see original citations for the stationary distributions of these models).

In one-phase models, a mutation leads to expansion or contraction of the microsatellite by 1 repeat unit only. However, empirical data suggest that mutations can occasionally result in a change in the microsatellite length of >1 unit. According to the two-phase model (TPM), proposed by D, there is a probability of p that a mutation changes the microsatellite length by 1 unit and has a probability of 1 – p that a change in length is ≥1 repeat unit(s), where the length of change is given by a geometric distribution. The generalized mutational model (F) is a simplified version of this mixed model, which sets p to 0, and consequently the length of change is entirely governed by the geometric distribution.

Many population genetics inference methods for microsatellite data require the adoption of a mutation model such as those described above. These approaches can be divided into three categories. The first category involves moment estimators based on summary statistics, including sample homozygosity (K; X) and allele length variance (W; K), to estimate θ = 4Neμ (four times the product of effective population size and the mutation rate).

The second category consists of likelihood-based approaches to the estimation of θ. As it is not in general possible to evaluate the likelihood function analytically, it is approximated by computational methods including Markov chain Monte Carlo (MCMC) (B) and importance sampling (N). Significant progress has been made in the development of methods that employ importance sampling and composite likelihoods for microsatellite inference, allowing the maximum-likelihood estimate of demographic parameters to be computed efficiently (I; R). On the other hand, W, B, and others have applied MCMC to provide Bayesian inference of demographic history from microsatellite data, in which case population parameters are treated as random variables instead of unknown fixed parameters as in a maximum-likelihood approach. C investigates the underlying mutation process of microsatellites using reversible-jump MCMC (G) of microsatellite models.

The third category includes likelihood-free approaches such as approximate Bayesian computation (ABC) (W; B; B). The application of ABC to microsatellite data (B; C; T) aims to increase computation efficiency as the method uses summary statistics instead of the full data set and employs simulation to circumvent the likelihood computation step.

Many of the likelihood approaches mentioned above are based on the coalescent theory (K; G). Rather than assuming a parametric model for the population history, for example exponential growth or logistic growth models (P), advanced coalescent-based methods provide inference of the demographic history by estimating population as a function of time directly from the data (D; O; H; M), but most of them have not been accessible for microsatellite inference.

To extend previous work on Bayesian coalescent inference of microsatellite data, we develop a method that provides inference of the demographic history averaged over a nested set of microsatellite mutation models that incorporate length dependency, mutation bias, and step size. Our method can handle multiple loci and these are assumed to be unlinked or in independent blocks of linkage. The implementation of this method consists of two main parts. The first part is to introduce the implementation of a rich family of microsatellite mutational models and other necessary components to the BEAST software package (D), which provides microsatellite inference access to sophisticated coalescent models (D; H; M). The second part deals with model uncertainty by employing the product space formulation of transdimensional MCMC (S) as described in section 2.5 of G, which facilitates Bayesian model selection by producing posterior probabilities of the microsatellite mutation models and Bayesian model averaging for estimates of population history and genealogies over those models. The transdimensional MCMC technique chosen here uses techniques from Bayesian variable selection (BVS) (G; K) sensu G. The BVS-inspired scheme is preferred over other transdimensional MCMC techniques because it trades a small increase in MCMC state space for a high degree of simplification and flexibility in programming.

To apply BVS a composite model space must be constructed that nests all submodels of interest over which inference of population history and genealogies should be averaged. In our case the submodels have a natural nesting by variable selection, because each model represents a special case of the most general microsatellite mutation model in our family. However, models that do not nest naturally can still be averaged over using BVS by introducing a simple index parameter alongside the union of all submodel parameters to construct a product space over models (C). However, an advantage of our nested model space is that it is able to indicate which of the three microsatellite mutation properties has a strong signal in the data.

MATERIALS AND METHODS

The basic global model:

Here we give an overview of the global model and the framework within which our implementation is developed. The microsatellite data, D, consist of L loci, D = {D1, … , D} and each locus is composed of a collection of microsatellite repeat lengths from the population of interest, where l = 1, … , L and n is the number of copies of locus l collected. In haploid data, n is the number of sampled individuals, while in diploid data n is twice the number of individuals from which the samples have been collected. We assume that the data have been generated by an underlying continuous time Markov chain (CTMC), along an unknown genealogy τ, which is a rooted bifurcating tree. In the simulations and analyses, the mutation rate is assumed constant along the tree within a locus, i.e., a strict molecular clock rate. The time intervals between successive coalescence events in the genealogy are modeled by the coalescent, which requires a demographic model component containing parameters Θ. The mutation process is defined by the microsatellite mutation model (more details in the Microsatellite models section) with parameters φ. Let τ = {τ, … , τ} and assuming the loci are independent and identically distributed given (τ, φ), the joint posterior distribution of τ, φ, and Θ is(D). The tree likelihood of locus l is Pr(D|τ, ) and can be evaluated using the peeling/pruning algorithm described by F, although we employ augmentation of internal nodes with repeat lengths. The coalescent models come into play by serving as priors for the tree topology and coalescent/divergence times. The form of the coalescent likelihood, f(τ | Θ), depends on a demographic model specified a priori and its parameters (Θ) are jointly estimated. The prior distributions for parameters of the demographic model and mutational model are selected from various standard univariate and multivariate distributions.

Microsatellite models:

The models of interest in this study were ones that could be approximated by finite state space continuous-time Markov chains to readily incorporate them into existing software for likelihood calculations on trees. We first need to decide on the coding of the data. Unlike nucleotides or amino acids that have a finite state space, the size of the microsatellite state space is ambiguous, because a universal upper bound for length of microsatellites probably cannot be defined (K). Yet, according to previous observations, the number of repeats in a microsatellite allele rarely exceeds a few tens (G). In addition, it seems sensible to impose a lower bound on repeat length, above which we can expect the characteristic behavior of microsatellite mutation to occur. In this article, an allele with i repeats is denoted as i. The imposed maximum and minimum lengths of a microsatellite are denoted as imax and imin, respectively. Therefore there are s = imax − imin + 1 possible states.

Once boundaries are set, it is easy to define the infinitesimal rate matrix of an ergodic Markov chain. The infinitesimal rate matrix is a square matrix wherein each element q, specifies the relative instantaneous rate of allele i mutating to allele j, and the shared lower and upper bounds of i and j are imin and imax, respectively. Given the mutation rate (μ), the Markov chain has the transition probability matrix (P),where p, (μt) is the probability of mutating from allele i to allele j, in time t.

In our implementation, the infinitesimal rate matrix of the most complex model is parameterized:orThe rate matrix is normalized so that the total mutational outflow is 1.0; i.e., let and find c so that −∑q′π = 1.0.

The function α(u0, u1, u2, d0, d1, d2, i) is the truncated version of the asymmetric quadratic model proposed by C and accounts for the length dependency of mutation rate and mutational bias by modeling the rate of expansion and contraction as separate quadratic equations. The rate can be symmetric if expansion and contraction share exactly the same quadratic equation; in other words u0 = d0, u1 = d1, and u2 = d2. Equal rate for all lengths is obtained by setting the coefficients of the linear terms and quadratic terms to zero in both equations. Similarly, the rates are modeled as linear functions of the length when u2 and d2 are set to 0.

In Equation 3 the α-function has symmetric rates. It is standardized so that parameters u0 and d0 are fixed to 1.0, because the rate of imin must be a positive real number and for any α(1, a1, a2, 1, a1, a2, i) it is equivalent to constant × α(1, a1, a2, 1, a1, a2, i), because of the normalization of the rate matrix.

The focal length is equal to if if imin ≤ if ≤ imax and is the repeated root of the equation (u0 − d0) + (u1 − d1)(if − imin) + (u2 – d2)(if − imin)2 = 0. At the focal length the rate of expansion and contraction is the same, so given a mutation event, there is equal probability of expansion and contraction.

Although the function α(u0, u1, u2, d0, d1, d2, i) has taken mutational bias into account, the parameterization may not necessarily provide answers to questions regarding the relationship between mutational bias and microsatellite length. As mentioned earlier, mutational bias can be quantified by the probability of expansion given a mutation event. The function β(b0, b1, i) models the probability of expansion by a simple logistic regression.

Both the bias constant parameter, b0, and the bias linear parameter, b1, take real values from the range (−∞, +∞). The probability of contraction is 1 − β(b0, b1, i). This is a modification of the parameterization adopted by S. They model expansion probability by simple linear regression. The probability of expansion then becomesIn Equation 4 the bias constant parameter is while the bias linear parameter is

It is worth noting that β(b0, b1, i) ∈ (0, 1) whereas The difference may not seem significant computationally, but can give rise to numerical instability, when there are several consecutive rows of zero expansion (or contraction) rates in the infinitesimal rate matrix. This situation is very rare since in reality it is unlikely that probability of expansion will be close to 0 or 1; it can occur, however, if the user does not specify the appropriate starting values for and . In addition, it is conventionally more appropriate to model a categorical variable with logistic regression (A).

Both parameterizations β and were implemented. Parameters of a logistic regression are not as easy to interpret as those of a linear regression, so for a more straightforward interpretation, can be chosen for the analysis.

To account for larger steps in state space by a single mutation, we employ the parameterization used by S, which is similar to the TPM proposed by D. Under this model, single-repeat mutations have a probability of p whereas multirepeat mutations (length of change ≥1 repeat) have a probability of 1 − p. For multirepeat mutations, the distribution of step size, |i − j|, is given by a truncated geometric distribution γ(g, i, j). The symbol g is the failure probability of the truncated geometric distribution.

Stationary distribution:

For an ergodic Markov chain, as time, t, approaches positive infinity, its transition probability matrix converges to a matrix in which every single row is the stationary distribution, lim→∞ P(μt) = 1π. As mentioned by S the stationary distributions of all one-phase models are special cases of the general birth–death chain.

Bayesian model uncertainty:

The output of a Bayesian analysis is the posterior distribution of the parameters given the data. However, the high-dimensional parameter space in a genealogy-based analysis dictates simulation of the posterior distribution by computationally intensive Monte Carlo methods such as MCMC or importance sampling. Here, the posterior distribution is produced by the Metropolis–Hastings MCMC algorithm (M).

In a Bayesian framework, the standard procedure to compare two models is by computing their Bayes factor (BF), which is the ratio of the marginal likelihoods of the two models (M1 and M2):If the space of potential models is large, then some techniques for model comparison are very time consuming. In addition and more importantly, the mutational model may not be of prime interest, i.e., nuisance, and therefore it is not ideal to perform a separate analysis for every mutation model. The solution to this problem is Bayesian model averaging (BMA). We employ transdimensional MCMC to provide joint inference via sampling the microsatellite model indicator, ν, to produce the posterior distributionHere, φ is a union of parameter vectors of all n models of interest, and where is the parameter vector of model M. The marginal posterior distribution of model indicator ν can be obtained from posterior samples of Pr(Θ, τ, φ, v|D), representing the posterior distribution of the microsatellite model. The joint posterior distribution of Θ and τ integrated over the models iswhich can also be expressed as the model-averaged posterior joint distribution of Θ and τ.

Our implementation of transdimensional MCMC combines the techniques of BVS (K) and pseudopriors or linking densities (C). Early BVS applications have aimed to solve the problem of variable selection encountered when building a linear regression model. Initially, there is a large number of potential predictors X, … , X, with values x, j = 1, … , p and the focus is on determining which of these predictors are linearly associated with the response variable Y. The full model describes the response y as a linear combination of the explanatory variables x:The term is the intercept, and the error term ε ∼ N (0, ). A coefficient φ that is (statistically) significantly different from 0 suggests predictor X may help in predicting the response. Conversely, a φ that does not significantly differ from 0 indicates X provides little additional information and can be excluded from the model. The variable selection method by K uses an auxiliary binary indicator variable, δ, of P dimension. δ = 1 indicates the presence and δ = 0 indicates the absence of the parameter φ. The full linear model becomesThe term ψ can be considered as the outcome of the function ψ = g(φ, δ) = φ. Setting δ to 0 forces ψ to 0, so that φ is effectively excluded from the model. However, even though in such a case the value of φ has no effect in the likelihood, it is still sampled by the MCMC machinery, but according to its prior distribution only. This means the dimensions of φ and δ, and hence the model parameter dimension, are not changed even though mutation model parameters are effectively included and excluded in the likelihood during the MCMC. Therefore it does not require the computation of the Jacobian ratio unlike transdimensional MCMC techniques such as reversible-jump (RJ)MCMC (G).

We augment the parameters in Equation 3 with a set of indicators, each associated with one of the parameters in the most general microsatellite mutation model, to produce a natural nesting of the described microsatellite mutation models.

Our most complex (full) microsatellite model (defined by Equation 3) contains the parameters, φ = {a1, a2, b0, b2, g, p}, and any submodel has only a subset of is the kth element in φ. We augment φ with a binary indicator variable δ = {δ1, … , δ6} to provide a set of toggle switches that can be used to define all nested models of φ. Letting q(φ) represent Equation 3, the equation that defines the instantaneous rate matrix of the new model is and our function for ψ is given byThus the instantaneous rate matrix depends on both φ and δ. Even when the value of has no effect on the likelihood as δ = 0, is still sampled by the MCMC machinery, but according to its prior distribution only. The joint posterior distribution when using this model iswhere φ = φ. If φ and δ are assumed independent, fΦ|(φ|δ)fΔ(δ) is replaced by fΦ(φ)fΔ(δ).

Prior on model space:

There are six free parameters in the full model, which theoretically give us 64 submodels. However, parameter p cannot be estimated for submodels in which g is not a free parameter (i.e., fixed to 0). This is because if g is fixed to 0, p does not have any effect on the likelihood. Furthermore, when modeling with regression, it is convention to estimate all polynomial terms in the model up to the largest degree considered in the model. We apply this convention to functions α(1, a1, a2, 1, a1, a2) and β(b0, b1, i) in Equation 3. The application of these restrictions to the model space results in a connected subspace of 27 models, and we apply a uniform prior so that the prior probability on each model is 1/27, while the remaining 37 models have a prior probability of 0.0. Figure A1 in shows the restricted model space.

F

Reduced model space produced by the restrictions described in the Prior on model space section in materials and methods.

Proposal distributions:

Model switching is performed by two proposal moves, the flip move and the pick move. The flip move uniformly picks an index of the bit vector δ at random and performs a flip, whereby the value at that index changes from 0 to 1 or vice versa. For a bit vector of n dimension, the probability is 1/n for both the flip move and its reverse. The Hastings ratio for M → M is q(M|M)/q(M|M) = (1/n)/(1/n) = 1.0. This proposal distribution over bit vector δ is symmetric and therefore no Hastings correction is required for this proposal. Since, in our case, the 27 models with nonzero prior probability form a connected component, the flip move will produce an ergodic and irreducible Markov chain. In effect, the nonhomogeneity of the restricted model space is corrected for by rejection of the neighboring models with zero prior probability, rather than defining a Hastings ratio specifically for the restricted model space (shown in Figure A1 of ).

The pick move, on the other hand, allows larger moves. It selects a model uniformly at random from the set of 27 permitted models. All permitted models have equal probability to be selected, and thus the pick move is symmetric.

Pseudopriors:

Pseudopriors or linking densities, a technique used in transdimensional MCMC, were first introduced by C. Their method considers the situation when there is no overlap among the individual parameter vectors of the n models of interest. The parameter “pool” φ is therefore simply a concatenation of all model parameter vectors. The joint posterior distribution for v and φ given data D can be written aswhere the and their values do not affect the likelihood when the current model is M. The expression represents the “pseudopriors” by C and can be considered the prior distributions of parameters in when (by definition) their values are not being used by the likelihood. If is assumed independent of the pseudopriors become Unlike “real priors,” these pseudopriors have no effect on the joint posterior f(v, |Y), but govern the mixing of MCMC as they play the role of jumping distributions in RJMCMC (G). Appropriate pseudopriors resemble efficient proposal distributions and achieve efficient sampling by preventing extremely unlikely parameter values of from being sampled.

Selecting suitable pseudopriors can overcome the problem of poor mixing encountered when the prior is very different from the posterior for parameters being model averaged. When a parameter is not in the likelihood, values sampled from the prior may have little agreement with the data. Consequently, a parameter may have difficulty reentering the likelihood, resulting in poor mixing.

G extended the method by C to allow arbitrary overlap among parameter vectors In addition to a pool of parameters, φ, indicators, I(ν), map ν to the elements of φ used by model M. The posterior is expressed as

C suggested that the pseudoprior of a parameter φ in model M should match closely the model-specific posterior distributions Pr(φ|Mv). It has been observed in some trial runs that even though two models and share the parameter φ, the marginal posterior distributions and are quite different. However, a parameter can have only a single pseudoprior. To achieve model specificity of pseudopriors, we can augment the parameter space, so that φ in is a separate parameter from that in , thus allowing different pseudopriors.

To accommodate augmented mutation model parameter space for a model-specific pseudoprior, the new function for ψ becomeswhere the mapping function I(k, δ) returns the index of the element in φ according to k, which indicates the type of parameter, and δ, which specifies the currently active model. For example, if and is the constant bias parameter for models of the form * *10* *, whereas is that for * *11* * models, then φ( maps to the parameter , where * is either 0 or 1. Again, when the rest of the parameters in the pool are not used in the likelihood, they are sampled from their respective pseudopriors and so the parameter dimension remains the same.

Tree likelihood computation:

F pruning algorithm of tree likelihood computation implicitly sums over all possible ancestral states. For a data type with s states and a rooted gene tree with n taxa (n − 1 ancestral nodes), the pruning algorithm is O(ns2). The speed of this algorithm is sufficiently fast for analysis to be completed on nucleotide data, which have state space size of 4 (A, T, C, and G). For microsatellite DNA, however, the number of states is many times larger than that of the nucleotide data type, and therefore likelihood calculation is much more time consuming.

One solution to this problem is to avoid summation across all possible states at ancestral nodes by treating unknown ancestral allelic states, , as auxiliary parameters (W). After augmentation of the tree with ancestral states and fixing to a particular microsatellite model with parameters the tree likelihood of loci l is the product of all likelihoods of nodes in a tree,where x is one of the 2n − 2 nodes in the tree excluding the root and i is the state of node x. The parent of node x is denoted as anc(x) and t is the length of the branch that connects x to anc(x). Following W, we replace Felsenstein's tree likelihood with the likelihood in Equation 14. The prior probability of the ancestral state in the root node, is computed from the stationary distribution of the mutational process, as is standard in Felsenstein's likelihood of an independent ergodic Markov process on a tree. Ancestral states in the remaining internal nodes have a uniform prior.

For a discussion on the numerical stability see .

Proposal moves for ancestral state sampling:

The candidate allelic state of an ancestral node is proposed by a random-walk integer move, which makes a step from the current allelic state. This move randomly picks direction and step size, which is an integer between 1 and a maximum step size specified by the user. The maximum step size permitted is less than the difference between the upper and lower boundaries of the allelic state. If after a random-walk step the value proposed exceeds the boundaries, then the exceeding proportion of the step is reflected back. Due to the condition on the maximum step size and the type of reflection chosen, the result of a reflection will not be on either boundary. Given maximum step size, w, the number of possible combinations of direction and step size is 2w. The Hastings ratio is thus the ratio for a move from state i to j and is given by H(j, i)/H(i, j), where H(i, j) = h1(i, j) + h2(i, j) + h3(i, j). The equations h1(i, j), h2(i, j), and h3(i, j) are given below:The proposal mechanism is independent of the currently indicated model in transdimensional MCMC.

Simulations:

After developing the implementation for Bayesian microsatellite analysis, it is of interest to obtain some indication of the accuracy and precision of the estimates. We consider only a subset of the 27 models in the restricted model space. This subset is obtained by setting a2 and p to 0; therefore the most complex model considered here has only four parameters and the bit vector δ has four dimensions. Because of the restriction that b1 is a free parameter only when b0 is a free parameter, this subset has 12 permitted models instead of 16.

These 12 models resemble the set in S, except we use simple logistic regression to model mutational bias. For convenience, we use their model naming system. This restricted model space of 12 models is illustrated in Figure 1.

F

Restricted model space of the 12 models considered in the simulation analyses.

Simulated data were generated under the 12 different microsatellite models from the procedure described below:

One hundred replicate data sets were generated under each microsatellite model.

For each replicate, 30 random coalescent trees were generated, each with 15 individuals assuming a constant population size with Neμ = 2.0 (where Ne is the effective population size of chromosomes and μ is the mutation rate representing the number of mutations per microsatellite locus per generation).

A microsatellite data type was created with minimum length of 1 repeat unit and a maximum of 30 repeat units.

For each coalescent tree, a microsatellite site pattern was simulated under the microsatellite model with mutation rate equal to 1.0. All site patterns in a trial were simulated under the same microsatellite model. There were 15 sampled haploid individuals in each site pattern.

Each of the 1200 simulated unlinked 30-locus data sets were analyzed with transdimensional MCMC to demonstrate how well our method identifies the true microsatellite mutational model (Mtrue). We also compared the accuracy and precision in the demographic estimates between model averaging and when the true model was known. Analyses have chain lengths of 70 million steps with transdimensional MCMC and 50 million steps with the true model. Sampled parameters were recorded every 50,000 steps.

It is also of interest to investigate the effect of the number of taxa and that of loci on the precision of θ-estimation using transdimensional MCMC. Data sets were simulated under the PL2 model with different combinations of number of loci and number of taxa presented in supporting information, Table S1, which also includes the MCMC chain length for each combination. One hundred simulations were carried out for each combination. For this set of simulations, we recorded every 10,000th step of the MCMC. The convergence of each analysis was examined by the trace analyses including the computation of the effective sample size (ESS) of each estimated parameter. Table S2 is a summary of the model parameter values that are chosen for simulating the data. All simulated data sets are provided in .csv format in File S1.

Measure of accuracy:

The accuracy was measured by computing the relative bias. Here we define relative bias aswhere θ is the true population size value and is the posterior median population size for trial k.

Accuracy of estimates may also be indicated from the percentage of trials with 95% highest probability intervals (95% HPD) containing the true answer. The HPD is the smallest interval containing of the posterior distribution.

Measure of precision:

The relative error was used to measure the precision of the estimates obtained. We define the relative error asAnother measure of precision is the 95% HPD relative bound and is given asThe credible interval coverage, relative bias, error, and HPD bound defined here are similar to the corresponding measures used by H.

Prior distribution for microsatellite model parameters:

We used a normal(0, 10) prior on both b0 and b1, an exponential(1) on a1, and uniform(0,1) on g.

From test runs it appears that only b0, b1, and g require pseudopriors to reach reasonable convergence. The pseudopriors for each variable are chosen to be tight distributions centered around the true parameter values since they were known. Pseudopriors for each parameter are shown in Table S3.

Tree prior:

For the simulations, the tree prior used was the coalescent with constant population (see K or G for details on the coalescent likelihood calculation). For inference, the constant population size model was chosen to match the simulation conditions.

The prior density for θ was set to one-on-x prior, . The one-on-x prior is an improper prior; however, in the case of constant population size, it can be shown to be Jeffrey's prior, and its application in this context leads to a proper posterior distribution (D, 2004).

Microsatellite model prior:

Because we did not have any a priori information regarding the microsatellite model, a uniform prior is applied to the set of 12 models considered.

Sampling tree topologies:

The tree proposal moves subtree-slide, narrow exchange, wide exchange, and Wilson–Balding are used for tree topology sampling in all the analyses undertaken in this article. A nice summary of these proposal moves is presented in H.

Red colobus monkey data:

Data:

The red colobus monkey (Pilocolobus tephrosceles) data set was kindly provided by J. Allen (University of Florida, Gainesville). This unpublished data set consists of 62 samples from each of 16 loci. Each locus was typed for both homologous copies from each of 31 red colobus monkeys from the Kibale National Park of Uganda. These loci are treated as unlinked as no clear signal of linkage had been found (J. Allen, personal communication). The allele lengths and the PCR primers are presented in Table S4 and Table S5.

Analyses:

An upper bound of 33 and lower bound of 6 repeats were imposed. We made boundaries wider than the observed length range of the data to account for the possibility that the observed sample range is not the true range in the population. We ran two separate analyses of 200,000,000 for each of the 12 microsatellite models. We also ran two replicate analyses using transdimensional MCMC. To compare mixing and performance between transdimensional MCMC and fixing the microsatellite model, we estimated θ from this data set with transdimensional MCMC and each of the 12 models used for simulation. Values for ESS per MCMC step were computed for θ, tree likelihoods, coalescent likelihoods, and mutation rates. We accommodated the potential mutation rate variation across loci by estimating the relative rates but fixing the average rate to 1.0.

Prior selection:

A uniform Dirichlet prior for 16 dimensions was applied to the relative mutation rates. The tree prior and mutation model parameter priors used for the real data analyses are the same as the ones for analyses of the simulation data. However, because the true values of the mutational parameters are unknown, pseudopriors could not be picked as easily as for simulated data. We took Carlin and Chib's suggestion and obtained preliminary posterior distributions of mutation model parameters by running a short MCMC (of 40,000,000 steps) with each microsatellite model. The posterior densities obtained from these preliminary runs are still quite broad, but they provide sufficient guidance for the selection of pseudopriors. The first 10% of each chain is discarded as burn-in and the remaining chain is used to fit a standard parametric distribution to the posterior sample of each parameter.

The marginal posterior distributions of the microsatellite mutation parameters were fitted using the maximum-likelihood–based fitdistr function in the MASS package (V) of R (R D), a software environment for statistical computing. The fitdistr function returns parameter estimates for a parametric distribution that best describes the posterior sample. The quality of the fit is then examined by the one-sample Kolmogorov–Smirnov test with the null hypothesis that the posterior sample has come from the fitted distribution. Several different parametric distributions were fitted to the sample and the one with the largest Kolmogorov–Smirnov test P-value (least evidence against a bad fit) was chosen.

Table S6 presents all the pseudoprior distributions. In this case, the parameter space is augmented so that there is no overlap in model parameter vectors and each parameter has its own model-specific pseudoprior.

RESULTS

For each analysis of simulated data, the first 10% of the chain was discarded as burn-in and analysis of traces confirm that all parameters have ESS > 100. The accuracy of Bayesian model selection by our method is indicated by how often the maximum a posteriori model corresponds with the true model, MTrue. Table 1 presents the frequency distribution of MBest (the model with maximum posterior probability using our transdimensional MCMC method) for data sets simulated under each MTrue. The highest percentage value in each row is in italics. It is shown that all diagonal values are in italics, which means MBest = MTrue has the highest frequency for all 12 models of MTrue.

TABLE 1

Percentage of true model recovery computed from transdimensional MCMC (tdMCMC) analysis of simulated data

	MBest
MTrue	EU1	EU2	EC1	EC2	EL1	EL2	PU1	PU2	PC1	PC2	PL1	PL2
EU1	95	4	0	0	0	0	1	0	0	0	0	0
EU2	13	87	0	0	0	0	0	0	0	0	0	0
EC1	0	0	59	0	1	0	0	0	37	3	0	0
EC2	2	5	13	63	0	1	0	0	2	15	0	0
EL1	18	0	0	0	59	6	1	0	0	0	17	0
EL2	0	0	0	0	10	76	0	0	0	0	6	14
PU1	0	0	0	0	0	0	97	3	0	0	0	0
PU2	0	0	0	0	0	0	17	83	0	0	0	0
PC1	0	0	21	0	0	0	0	0	73	6	0	0
PC2	0	1	1	5	0	0	0	0	18	75	0	0
PL1	4	0	0	0	4	1	3	0	0	0	78	10
PL2	0	0	0	0	0	1	8	5	0	0	21	65

Accuracy is also indicated by computing the percentage of trials that has the MTrue contained in the 95% HPD set of models. These values are presented in Table 2. A very high proportion (>0.9) of the trials captures MTrue within the 95% HPD set. The median 95% HPD set size is between two and four models, and the majority of the models have a median set size of two, suggesting good precision.

TABLE 2

Accuracy and precision of true model recovery

MTrue	Inside 95% HPD	95% HPD set size
EU1	1.00	2
EU2	1.00	2
EC1	1.00	3
EC2	0.98	3
EL1	0.95	4
EL2	0.99	2
PU1	1.00	2
PU2	1.00	2
PC1	1.00	3
PC2	0.99	2
PL1	0.98	3
PL2	0.93	2

Often, the user is more concerned with the accuracy of other evolutionary or demographic estimates rather than the mutational model per se. The values for relative bias, relative error, and relative 95% HPD bounds for the demographic parameter of constant population size are calculated for each trial. Table 3 is a summary of the median relative bias, the median relative error, the median 95% HPD relative bound, and the percentage of trials in which the 95% HPD interval captured the true value of Neμ = 2.0 for each model.

TABLE 3

Measure of accuracy and precision of model-averaged θ-estimates from transdimensional analyses and of θ-estimates from analyses that fixed the microsatellite mutational model to the true model

	Inside 95% HPD		Median relative error		Median relative bias		Median relative bound
MTrue	BMA	TM	BMA	TM	BMA	TM	BMA	TM
EU1	0.97	0.98	0.10	0.10	−0.06	−0.02	0.58	0.55
EU2	0.94	0.97	0.12	0.10	0.07	0.04	0.83	0.73
EC1	0.96	0.93	0.12	0.12	0.02	0.01	0.68	0.61
EC2	0.89	0.91	0.21	0.14	0.09	0.11	0.97	0.77
EL1	0.89	0.92	0.10	0.10	−0.02	0.01	0.59	0.55
EL2	0.92	0.92	0.16	0.12	0.14	0.08	0.95	0.70
PU1	0.98	0.99	0.10	0.09	−0.06	−0.03	0.65	0.62
PU2	0.94	0.99	0.22	0.18	0.15	0.10	1.01	0.87
PC1	0.93	0.93	0.11	0.10	−0.01	0.03	0.65	0.63
PC2	0.86	0.92	0.16	0.14	0.10	0.08	0.95	0.87
PL1	0.90	0.93	0.11	0.12	−0.04	0.03	0.32	0.31
PL2	0.89	0.95	0.15	0.13	0.07	0.00	0.80	0.75

Estimates with high precision have small values of median relative error or median 95% HPD relative bound. Accurate estimates have small values of median relative bias and a high percentage of 95% HPD intervals containing the true value. Within each row, the values of the four statistics of accuracy and precision are close between BMA and when the true model is known. However, when the true model (TM) is known, the results have greater coverage of the true population parameter value, smaller median relative error, smaller median relative 95% HPD bound, and smaller absolute median relative bias than model-averaged estimates of θ.

All model–method combinations had high frequentist coverage varying from 0.86 to 0.99. Given the small number of replicates, these coverage statistics are not significantly different from each other and are all consistent with an underlying proportion in the region of 0.95, although in the Bayesian setting there is no reason to expect coverage to be at the 0.95 level.

For either method, there is a spectrum of median relative error values across the models of MTrue, where the median relative error value is the smallest when the data are simulated under PU1 and largest when data are simulated under PU2.

Similarly, there is quite some variation in the size of the median relative 95% HPD bounds across different models of MTrue for both BMA and analysis with the true model. Median relative 95% HPD bounds are the smallest for data simulated under PL1 and the largest for PU2, which are approximately three times of that for PL1.

Precision vs. number of loci:

Values of median relative error and median relative HPD range are plotted against the number of loci (Figure 2), where the number of taxa is fixed to 15. As the number of loci increases, the relative error (Figure 2, dashed blue line) and 95% HPD range (Figure 2, solid red line) reduce linearly (on a log-log scale). It appears that increasing the number of loci to eight times larger will halve the relative error and reduces the 95% HPD range by a factor >2.

F

Measures of precision θ-estimation vs. the number of loci.

Precision vs. number of taxa:

Figure 3 shows median relative error and 95% HPD credible intervals as a function of the number of taxa per loci for 10 unlinked loci. Both the error (Figure 3, dashed blue line) and the HPD interval (Figure 3, solid red line) decrease as the number of taxa increases, but they seem to asymptote to some positive limits. Further reduction of error and HPD range can be achieved only by sampling more loci. The function form y = a0+a1/(a2+x) is chosen merely to illustrate the general trend.

F

Measures of precision θ-estimation vs. the number of taxa.

Real data example—colobus monkey data:

Tracer (R) was used to decide the length of the chain to be discarded from a log file as burn-in. The burn-in length was ∼10–20% of the original length of each log file. The two logs from analyses with the same microsatellite model were subsequently combined. The combined log files were then examined again by Tracer to investigate mixing and convergence to stationarity. All ESS values were >150.

According to the results from transdimensional MCMC, the model with the highest posterior probability is EU1. The 95% highest posterior probability set consists of EU1, EU2, PL1, EL1, and EL2, with the respective probabilities 0.483, 0.324, 0.053, 0.049, and 0.048. The posterior probability of including a parameter a1 = 0.089, b0 = 0.193, b1 = 0.186, and g = 0.412. These posterior probabilities suggest that multistep mutation is the most evident feature followed by mutation bias and rate dependency.

The posterior mean, median, and 95% HPD interval from analyses with each model are presented in Table 4.

TABLE 4

Estimates of θ from red colobus monkey data

Model	Mean	Median	95% HPD lower	95% HPD upper
BMA	4.08	4.01	2.56	5.76
EU1	4.37	4.31	3.02	5.83
EU2	3.54	3.47	2.25	4.98
EC1	4.33	4.28	3.06	5.75
EC2	3.46	3.40	2.27	4.83
EL1	4.62	4.54	3.26	6.15
EL2	3.87	3.80	2.56	5.38
PU1	4.33	4.28	2.93	5.92
PU2	3.50	3.45	2.11	4.94
PC1	4.60	4.48	2.85	6.68
PC2	3.59	3.49	2.06	5.28
PL1	5.05	4.97	3.33	6.88
PL2	4.05	3.98	2.45	5.72

The posterior median of θ for single-step models ranges from 4.28 to 4.97 and that for multistep models ranges from 3.40 to 3.98. Since a multistep model was sampled almost half the time, the model-averaged posterior median of θ is somewhere in between.

Mixing and performance:

We use ESS values per MCMC step as an indication of the sampling efficiency. Figure S1, Figure S2, Figure S3, Figure S4, and Figure S5 are dot plots of ESS value per MCMC step for θ, tree likelihoods, coalescent likelihoods, root heights, and mutation rates. These plots suggest that the sampling efficiency of transdimensional MCMC is only slightly less than the average sampling efficiency of single-model analyses. The ratio of ESS per MCMC step for transdimensional MCMC vs. single-model analyses on average is 0.78 for θ. Averaging across loci, the ratio is 0.90 for tree likelihood, 0.74 for coalescent likelihood, 0.71 for root height, and 0.98 for relative mutation rate.

DISCUSSION

The focus of this research was on the implementation of a nested family of microsatellite mutation models in the BEAST software package (D). There are many analysis tools unique to BEAST, including nonparametric coalescent-based inference methods such as the extended Bayesian skyline plot (H) and the newly developed multispecies coalescent method *BEAST (H). Equipping BEAST with microsatellite models and other related software components permits the application of these methods to microsatellite data.

The microsatellite models implemented have all been previously described in the literature (D; F; C; S). It was not intended to introduce new models of microsatellite evolution in this work, except to make slight modifications where it made the models more suitable for Bayesian inference in an MCMC setting. Aside from their implementation, simulations were used to investigate their statistical properties.

Simulation results show moderate variation in performance across data sets simulated under different microsatellite models. However, when a smaller number of loci were simulated (10), more biased results were observed (results not shown).

As mentioned in parameter prior specification (see section on Prior distribution for microsatellite model parameters), free parameters shared by more than one microsatellite model have been given the same prior distribution for all models containing the parameter. For example, the constant bias parameter, b0, is in models EC1, EC2, PC1, and PC2, and the prior distribution for b0 is the same for all four models. The impact of the prior choices made for the mutational parameters has not been investigated in this study. It is quite possible that different prior choices would have altered the statistical properties of the estimators. For parameters such as g in the two-phase models, which are defined on [0, 1], the uniform prior is natural; however, for other scale parameters, a number of alternatives are feasible. Therefore, unsuitable prior distributions on the microsatellite model parameters may be partially responsible for the demographic estimation bias observed.

Red colobus monkey data example:

Our results suggest that the convergence speed of transdimensional MCMC is only slightly worse than for single microsatellite model analyses on average. In addition, only one analysis of transdimensional MCMC is required to perform model selection; however, for single-model analyses, we would need as many as the number of models of interest (12 in this case) and thus require a far longer time.

In this analysis we have selected pseudoprior densities by fitting univariate distributions to densities acquired from preliminary runs. The procedure can become much less time consuming if an empirical density function is used, coupled with the automation of the input file preparation for preliminary runs and the transdimensional MCMC. However, poorly mixed preliminary runs may still yield pseudopriors that are very different from the posterior and thus offer little improvement in mixing of the analysis with transdimensional MCMC. If human data are analyzed, then appropriate pseudopriors may also be constructed from the wealth of empirical data on the mutation parameters from sperm-typing (Z) and pedigree studies (W; X; W).

Model averaging and Bayes factors:

In addition to providing a model-averaged posterior distribution of population and genealogical parameters, our transdimensional MCMC method also facilitates robust estimates of the marginal likelihoods (and therefore Bayes factors) of the individual microsatellite models. These estimates are not subject to the large, even infinite variances (R; C) associated with the harmonic mean estimator of marginal likelihoods (N). Using transdimensional MCMC to estimate Bayes factors is also computationally efficient as it requires only a single MCMC run to determine the relative merits of all k mutational models, rather than the k (or more) independent runs required by other techniques, including thermodynamic integration (L).

Issues and improvement:

Low information content in a single microsatellite locus:

Low information content of a single microsatellite locus means inference results may be sensitive to poor prior choices. In comparison to microsatellite data, mtDNA sequence data possess much more information available for reconstruction of the genetic ancestry. As a result, an mtDNA tree has a higher level of resolution than a microsatellite tree. However, this does not mean inference on mtDNA sequence data is more accurate than that on microsatellite data. Given a population history, the coalescent admits wide variation in the topology and coalescent times of the gene trees. To make a more reliable inference, it is important to use multiple loci, each of which has an independent history (F). Even though mtDNA sequence data provide a clear view of the genealogy of the mtDNA sequences, the whole mtDNA genome is a completely linked locus. The mtDNA tree therefore provides only one of the many possible realizations of the coalescent process for a given population history. On the contrary, there are potentially thousands of independent microsatellite loci available to overcome the problem of stochastic variability of individual genealogies.

Speed and convergence:

A large microsatellite data set with hundreds of loci may give very accurate population size estimates, but is currently not practical in our implementation, due to slow convergence and computational inefficiencies. An analysis for a data set containing ∼60 unlinked loci and 100 taxa requires days to satisfy our heuristic diagnostic statistics for convergence when all the loci were used simultaneously in the same MCMC run. The slow convergence is due to the large joint parameter space when all loci are unlinked. The parameter space containing 60 independent 100-tip trees is much larger than that having a single 100-tip tree with 60 linked sites. One potential solution is sequential Monte Carlo methods (L), which take advantage of the independence structure of the likelihood to build up a full posterior distribution by sequential analysis of the loci (D). Besides the speed, there is also the need to improve the efficiency of sampling ancestral states. Our implementation samples ancestral states by a naive Metropolis–Hastings algorithm and therefore has low acceptance probability. Gibbs sampling (G) is an alternative MCMC algorithm that can sometimes produce more efficient sampling. It is a special case of the Metropolis–Hastings algorithm, whereby each proposed candidate is always accepted, since the components of the state that change in the proposal are drawn directly from the conditional posterior distribution. Gibbs sampling of internal nodes may improve the convergence.

Comparison with other software:

Well-known software programs such as BATWING (W) and Migrate (B) also provide Bayesian coalescent analysis on microsatellite data. However, these programs contain only a few simple microsatellite models. BATWING provides the SMM and the K-allele model; microsatellite model options in Migrate are the SMM (called the ladder model in Migrate) and Brownian motion (an approximation of the SMM). These models do not take into account many properties of microsatellite mutation and as mentioned in the Introduction there is much evidence for those properties. Therefore the simplifying assumptions of the SMM may not be adequate to perform inference on real data. Furthermore, we provide model averaging over a rich set of microsatellite models, which is absent in these programs.

In the case of BATWING, all microsatellites are assumed to be linked into a single locus. It was discussed earlier that incorporation of multiple loci is necessary for accurate inference. Using only a single locus overlooks the genome-wide distribution of microsatellites, a highly advantageous trait for coalescent inference.

Future directions:

In all analyses in this study, all loci shared the same (model-averaged) microsatellite model within an MCMC run. It is possible that the properties of mutation vary across loci. While our implementation allows for variation in both rates and microsatellite models across loci, we have not performed a systematic study of the properties of such models. A hierarchical prior structure can account for variation of the same component in different models. For example, every locus could have its own EC1 model, containing the parameter b0. During the MCMC, each b0 varies according to a given distribution, and the parameters of this distribution (hyperparameters) also have a prior.

Our framework for analysis of microsatellite data can be combined with the multispecies coalescent (H) to estimate the species tree using multiple microsatellite loci sampled from closely related species. Although microsatellite models can be used alongside various relaxed-clock models in BEAST (D) to estimate divergence times, we do not recommend this type of analysis, because each microsatellite locus does not have sufficient information to estimate rate heterogeneity among branches.

Model selection has always been an important problem in statistical inference. It is common to make inferences on the basis of the best model selected by a standard model comparison procedure. However, such a procedure may produce a subset of models that are not significantly different from one another in their goodness-of-fit and therefore create difficulty in deciding which model provides the most reliable inference. Our transdimensional method allows the data to speak for themselves and more importantly makes population inference on the basis of a set of microsatellite models, accounting for model uncertainty and avoiding model misspecification.

APPENDIX A: RESTRICTED MODEL SPACE

Figure A1 represents the restricted model space. The nodes represent the models, each labeled with its bit vector representation, and the arrow-edges are the Hastings ratios of a move from one model to it neighbor. Two models are neighbors if they have only a single difference. The nodes (models) are color coded according to the number of neighbors they have.

In this case the restricted model space is one connected component. However, some prior specification with zero probabilities on a subset of models may result in two or more disjointed components. In such cases the flip move alone cannot produce an ergodic Markov chain; therefore the pick move must be used so that all the models in the restricted space can be proposed.

APPENDIX B: NUMERICAL STABILITY

We used a few shortcuts in the tree likelihood calculation. If a proposal move results only in a change of likelihood on a few branches, then we subtract the initial logarithm of partial tree likelihood from the logarithm of full tree likelihood and add the new log partial tree likelihood. This is more efficient than computing the entire full log likelihood.

We have found that if only partial tree likelihood calculation is used, the difference in likelihood between partial and entire likelihood calculation at each step increases as the MCMC proceeds. However, a step that requires entire likelihood calculation will set this difference back to 0. Luckily, in a real analysis, there are many other parameters that will force the calculation of the entire full likelihood rather frequently. We have also provided the option so that the user can force the entire full likelihood computation every n number of likelihood computations.

Another component that is relevant to the numerical stability of the method is matrix exponentiation. Matrix exponentiation is achieved by using codes adapted from Cern Colt library 1.2 (http://acs.lbl.gov/software/colt/, more details in ). To ensure that matrix exponentiation of the CERN colt library is reliable, we compare the matrix exponentiation results computed from our codes with from the function MatrixExp in the msm package (J) of R.

Comparisons were made for five values of t (0.001, 0.01, 0.1, 1, and 10) and 12 different Q matrices, which are the instantaneous matrices under which our data were simulated. There was very little difference between the two exponential methods (all differences <10−13).

APPENDIX C: MATRIX EXPONENTIATION

The method of matrix exponentiation here requires the Q matrix to be diagonalizable. This requirement is checked by using the method described in G. If the method indicates that the requirement is not met, a 0 probability matrix is returned, which leads to the rejection of the proposal move. The code adapted from Cern COLT library 1.2 performs an eigen decomposition of the instantaneous rate matrix, so the matrix Q can be expressed as Q = VUV−1, where V is a matrix of eigenvectors and the U is a diagonal matrix of eigenvalues. Given some value of t, the exponential of Qt is then obtained by finding the matrix product e = Ve−1. The expression e is also a diagonal matrix, with diagonal values , where u represents the ith diagonal value.