Approximate Likelihood Estimation of Divergence Time Range Using a Coalescent-based Model.

We present an estimation of divergence time-range based on a coalescent model. This model sum the probability of coalescent trees, taking into account the effect of incomplete lineage sorting. Maximum likelihood estimate based on this model has been computed previously; however, a formula for divergence time-range estimate or confidence interval has never been presented for this model as the expression of the likelihood makes this estimation difficult to compute. Our formula for the divergence time-range estimate can be readily coded into a program. We did not use a simulation or resampling-based approach and therefore our method is fast and less computationally intensive. We demonstrate that our method is much faster and as accurate a simulation-based approach.

Introduction

Estimation of divergence time is an integral part of population genetics and evolutionary biology. It is a common practice1,2,3 to estimate a range or a confidence interval of divergence time rather than a point-estimation (ie, a single value).

Despite being a common practice, maximum likelihood estimation of the time-range poses a challenge for coalescent-based models. This is partly due to the complicated expression of likelihoods of the coalescent-based models, and partly due to the fact that a likelihood-based formula for confidence interval requires the computation of estimated Fisher’s information matrix, which has an even more complicated expression than the likelihood.

In this article, we present a divergence time range or a confidence interval based on the maximum likelihood estimation (MLE) of the divergence time from a coalescent-based model. We focus on the coalescent model of population evolution by Nielsen et al,4 where the likelihood is computed by summing up over all possible coalescent trees between the present day and most recent common ancestor (MRCA). This one-step procedure takes into account the uncertainties of estimating the coalescent trees, as well as the effect of incomplete lineage sorting.

The coalescent framework is consistent with several models of population evolution, including a diffusion model, a Wright-Fisher model, a continuous-time, or a discrete-time Moran model (see5,6). Specifically, the coalescent framework in4 requires a model with finite population size mating randomly within each population. The framework works with both monoecious and dioecious organisms and for both haploids and diploids. A complete separation between the populations is assumed at the point of divergence, and consequently, between-population mating are not allowed after the divergence. This also means that the coalescent events cannot take place between two individuals of different populations (between the present and the point of divergence).

It is further assumed that each locus has two alleles, and the allele-frequency spectrum for these biallelic loci is modeled with a symmetric Beta distribution as in.5–9 This assumption of a symmetric Beta distribution results in a symmetric Beta-Binomial distribution for the allele-count. (Note that, the assumption of Beta distribution comes from a diffusion approximation; see, for example).10

The model is valuable for estimating the divergence time from related populations, as it produces the likelihood directly from the data, bypassing the gene trees. This is achieved by computing the exact probability of each coalescent tree as well as probabilities of allele frequencies given these trees and then summing over them with a closed-form mathematical formula. Thus, potential misestimation due to gene tree incongruence (eg, incomplete lineage sorting, see for example),11 is avoided.

While computing the probabilities of the allele counts, it is assumed in this model that the effect of new mutations in allele frequencies between the MRCA and the present are negligible. This is an appropriate assumption for divergence estimation from related populations. This is because the divergence times in such populations are typically short, making the number of mutations very small. The variation in allele type is assumed to come from the mutations at or before the MRCA. A detailed mathematical description of the model probabilities is provided below.

There have been significant developments made in this model since its inception by4 and a number of methods of inference have been proposed. Later12 introduced a Markov Chain Monte Carlo-based on this model. More recently,6 introduced a pruning algorithm to systematically compute the likelihood under this model for inference on a population tree. This approach builds on the approach of4 and makes it possible to compute the likelihood of a large tree under this model. An innovative two-stage pruning algorithm is introduced for simultaneously keeping track of probability of the number of lineages, and allele counts among them. Later7 introduced a composite likelihood method based on this model that can be used to analyze dependent data. This method treats the dependent allele counts from nearby loci as independently by multiplying the marginal likelihoods obtained from each locus. As a result, a composite likelihood is computed, which is then maximized to obtain a maximum composite likelihood estimator. This method makes it possible to use the information in physically close loci, whereas previous approaches could only use a set of independent loci to compute the likelihood.13 Introduced a variation of the model that takes into account the effect ascertainment correction in the likelihood and MLE. This method is useful for data from loci that were selected because of their observed allele-frequencies (ie, ascertained). This method modifies the pruning algorithm of6 so that the likelihood is corrected for ascertainment bias. In the same year13 introduced a computational method for incorporating the effect of mutation at each branch in the pruning algorithm of.6 This method makes it possible to apply the pruning algorithm to data from different species. This is because although the effect of new mutations can be ignored in closely related populations, their effect s are large when comparing different species. Finally,9 has shown this model to be identifiable. The identifiability is an important desirable property of a statistical model; unidentifiable model parameters are ill-defined and therefore inference on such models could produce erroneous, confusing, and self-contradictory results.

Although the MLE of the divergence time is computed based on4,6 model by the previous authors, a formula for the range of divergence time has not been described. The computation of the MLE requires computation of the likelihood and then maximization of the likelihood over a period of parameter values. Computation of a range or confidence interval for divergence time requires estimation of variance and MLE, which is a more complicated procedure.4 computed an estimate for the variance through simulations, which could be used to compute a confidence interval. To estimate variance in this manner, one needs to first estimate the divergence time; then one needs to simulate the whole dataset a large number of times (say M = 100 times) using the estimated value of the divergence time as the real divergence time. Then, the divergence time is estimated from each of those M = 100 simulated datasets. The sample variance among the M estimated divergence time is taken as the estimated variance. However, this method takes a long time to compute, as one needs to simulate and estimate the divergence time M = 100 times (or a pre-determined large number of times). A resampling approach (eg, bootstrap) would have a similar large time requirement.

Here, we present a formula for computing an asymptotic approximation for divergence time range. Our formula is based on asymptotic approximation of the variance of the MLE. In statistical literature, this approximation is known to be a first converging approximation.14 Typically, the number of independent data-points (independent loci) is very high (>10, 000) in genetic data and therefore this will be a close approximation. In addition, we also provide a formula for first and second order mixed derivatives of the likelihood and log-likelihood. This formula is useful as it can be used for maximizing the likelihood with Newton Raphson Method. To evaluate the performance of our method, we also estimated the coverage probability of the confidence interval through simulations, and established that the approximation is indeed quite good, and as accurate as a purely simulation-based approach. However, as demonstrated in the article, our method is much faster and could be computed in a fraction of time it takes to estimate the range using a simulation-based approach.

For demonstration purposes, we have used our method for estimating the range of the divergence time between two populations in HapMap data.15 Our estimates are found to be of the same range as a recently published estimate of the divergence time of the same two populations.

Model and Definitions

In this section, we will briefly describe the model of.4,6 This model assumes that the divergence time is sufficiently small so that we can ignore the effect of mutation in the site-frequency-spectrum after divergence. Consider populations A and B with MRCA O (Fig. 1). Let us assume that each locus exhibits two alleles, arbitrarily named ‘0’ and ‘1’, by ‘allele count’ we will mean the count of allele ‘1’. The model of4,6 gives the probability distribution of the allele counts (r, r) from (haploid) sample of sizes (n, n) respectively. (Note that by “haploid” we mean a sampling unit. The data can be from either haploid or diploid organisms. A haploid sampling unit is as a set of chromosomes containing one chromosome of each type. Each diploid individual has two sampling units).

Figure 1

Variables associated with our model.

The parameters of this model are τ (>0), the divergence time in generations in the unit of effective population size, or the population scaled divergence time (τ = 2 Nt), and θ(>0), a mutation parameter or a population scaled mutation rate (θ = 4 Nμ) where t, N and μ are, respectively, the number of generations since the divergence, the effective population size and (raw) mutation rate at the time of divergence. Note that we assume a molecular clock; ie, the scaled time between A and O is same as the scaled time between B and O. However, this assumption is not binding, and we will discuss relaxing this assumption later in this article.

Next, we will describe how the probability distribution of (r, r) is computed. First, we follow n and n lineages at populations A and B respectively to the MRCA O, and compute the probability distribution of the number of coalescent events k and k, respectively. Let n0 and n0 be the number of lineages at O that are ancestral to the sampled lineages at A and B, respectively, and let r0 and r0, respectively, be the allele count out of them (Fig. 1). Note that n0 = n − k and n0 = n − k. The distributions of n0A and n0 can be computed from the following formula (first proposed by).16

m = A, B, where λ = j(j − 1)/2 (1) and

Let n0 be the total number of lineages at O that are ancestral to the all lineages sampled at A and B. Using the fact n0 = n0+n0 along with Eq. (1), one can compute the probability distribution of n0. Let r0 be (random) allele count out of the n0 lineages at the MRCA (Fig. 1). (r0 = r0+r0.) The probability of r0 given n0 is given by a ‘root distribution’ that varies in different versions of the model. We use the root distribution used in6 (symmetric beta-binomial)

where β(., .) is the beta function; θ is the aforementioned mutation parameter which is to be estimated as well. Note that as mentioned in,6 beta-binomial distribution model is due to the fact that the alleles at the root are binomial draws from the allele frequencies, and the allele frequency spectrum is modeled with a beta distribution (see, for example).10

Given n0, n0, and r0, the distribution of (r0, r0) can be computed as

4,6. Then, given n0, n0, r0 and r0 the distribution of (r, r) can be computed as

m = A, B.6

Next we combine the Eqs. (1, 2, 3, 4). With the observed allele counts are (r, r), the likelihood of (τ, θ) can be computed as

Note that setting τ = 0 does not maximize the above expression, as the coeffcients of the exponential terms could be positive or negative.

Thus, we have described the full model. Maximum likelihood estimate (MLE) of the divergence time is computed by numerically maximizing right side of Eq. (5) above. In the next section we will present an estimator of the divergence time range, rather than a point-estimator of divergence time.

Methods

Let the MLE of (τ, θ) be (τ̂MLE, θ̂MLE). Using the standard statistical results,

where → denotes convergence in distribution (see, for example);7L is the number of independent data-points (independent loci in our case), Normal2 denotes a bivariate normal distribution, and Inf denotes Fisher’s Information matrix:

where L(τ, θ) = L(τ, θ; (r, r) = (j, j)) is the likelihood. Using Eq. (6), one can estimate a (1 − α) confidence interval for τ as

and a (1 − α) confidence interval for θ as

where Inf() is the element (i, j) of Inf−1, z1−/2 is the (1−α/2)th quantile of standard normal distribution. Next, we obtain a simpler expression for Eqs. (8, 9).

Using Lemma 1 in the Appendix A (Eqs. (13, 14)) and Eq. (5) it follows that

for l, m ∈ {1, 2}. Also, note that

(l, m) ∈ {(1, 1), (2, 0), (0, 2)} and hence,

(l, m) ∈ {(1, 1), (2, 0), (0, 2)}.

Estimation

Using Eqs. (5, 7–12), one can numerically compute confidence intervals of τ and θ as follows.

First, τ̂MLE and θ̂MLE are computed maximizing the likelihood given in Eq. (5). The maximization may be done over a grid of values of τ and θ; alternatively, the Newton Raphson Method may be used; the derivatives of likelihood or log-likelihood for use in Newton Raphson Method can be obtained by using the expressions of derivatives in the likelihood in Eq. (10) in the expression of Eq. (11).

Once we have τ̂MLE and θ̂MLE, we can compute τ̂MLE, θ̂MLE using Eqs. (10,12) with (τ, θ) substituted by τ̂MLE, θ̂MLE. Next, is Inf(τ̂MLE, θ̂MLE)−1 numerically computed numerically and subsequently the confidence intervals of Eqs. (8,9) are computed.

Results: Simulation and Comparison with Direct Simulation Method

We have simulated N = 1, 000 datasets for each combination of sample sizes n = n = n, divergence time τ and number of independent loci L; n can take values 4 and 8; τ can take values 0.01, 0.02, 0.05, 0.075, 0.1, 0.15 and 0.2; L can take values 10,000 and 100,000. The mutation parameter was kept fixed at θ = 4 × 10,000 × (1.1 × 10−8) to reflect a (human) effective population size of 10, 000 and a mutation rate of 1.1 × 10−8.

Mechanism of model generation

Following the model of4,6 we started with a divergence time τ, mutation parameter θ, sample size (identical for the two populations) n = n = n, and a predetermined number of independent loci L. For each locus the process is identical and independent (given the parameters). Therefore, we will describe the process for a single locus only.

For a given locus, we simulate the model (described in Section 2) as follows. First n0 and n0 are simulated from n and n using Eq. (1). Next, n0 is simulated (or computed) as n0 = n0+n0. Then, r0 is simulated from n0 and θ from symmetric beta-binomial distribution (Eq. (2)) as in.6 Symmetric beta-binomial distribution is used in6 to characterize the allele-frequency spectrum. That is, it is assumed in6 that the allele-frequency spectrum over the L loci has a symmetric Beta distribution, and r0 given n0 is a randomly drawn allele-count from a sample of n haploids where the allele frequency is P (which is a random draw from the allele-frequency spectrum and has a symmetric beta distribution with parameter θ). Thus, r0 given n0 has a symmetric beta-binomial distribution with parameter θ. Next, we simulate r0 and r0 from n0, r0, n0 and n0 using Eq. (3). Then r is simulated from n, n0 and r0 using Eq. (4) for m = A and B. This process is repeated independently (given the parameters) L times to generate the allele counts r and r for each of the L loci.

For each combination of (n, τ, L), an MLE (τ̂MLE, θ̂MLE) was computed by maximizing the full likelihood over τ and θ for each N = 1, 000 repetitions. Thus, for each combination of (n, τ, L), we have N = 1, 000 estimates of (τ̂MLE, θ̂MLE). Then, applying our methods on these estimates, 95% confidence intervals for τ was estimated for each of N = 1, 000 repetitions for each combination of (n, τ, L). Then for each combination of (n, τ, L), we computed the average of the confidence intervals over N = 1, 000 estimated values (reported in Table 1 at the CIAsymp columns).

Table 1

Simulation: confidence interval for τ and estimated coverage probability P.

	L = 10,000CIASYMP	SIM. VAR. CI	L = 100,000CIASYMP	SIM. VAR. CI
τ = 0.01n = 4	CI: {0.01 ± 0.0010} P = 0.951	{0.01 ± 0.0009} P = 0.955	CI: {0.01 ± 0.0003} P = 0.949	{0.01 ± 0.0003} P = 0.955
τ = 0.01n = 8	CI: {0.01 ± 0.0004} P = 0.947	{0.01 ± 0.0005} P = 0.954	CI: {0.01 ± 0.0001} P = 0.954	{0.01 ± 0.0001} P = 0.945
τ = 0.02n = 4	CI: {0.02 ± 0.0015} P = 0.945	{0.02 ± 0.0014} P = 0.960	CI: {0.02 ± 0.0005} P = 0.950	{0.02 ± 0.0003} P = 0.953
τ = 0.02n = 8	CI: {0.02 ± 0.0006} P = 0.948	{0.02 ± 0.0005} P = 0.948	CI: {0.02 ± 0.0002} P = 0.942	{0.02 ± 0.0002} P = 0.944
τ = 0.05n = 4	CI: {0.05 ± 0.0025} P = 0.954	{0.05 ± 0.0029} P = 0.959	CI: {0.05 ± 0.0008} P = 0.955	{0.05 ± 0.0010} P = 0.953
τ = 0.05n = 8	CI: {0.05 ± 0.0012} P = 0.951	{0.05 ± 0.0010} P = 0.944	CI: {0.05 ± 0.0004} P = 0.951	{0.05 ± 0.0005} P = 0.945
τ = 0.075n = 4	CI: {0.075 ± 0.0032} P = 0.954	{0.075 ± 0.0036} P = 0.949	CI: {0.075 ± 0.0010} P = 0.942	{0.075 ± 0.0011} P = 0.950
τ = 0.075n = 8	CI: {0.075 ± 0.0016} P = 0.944	{0.075 ± 0.0015} P = 0.943	CI: {0.075 ± 0.0005} P = 0.952	{0.075 ± 0.0005} P = 0.961
τ = 0.1n = 4	CI: {0.10 ± 0.0040} P = 0.964	{0.10 ± 0.0041} P = 0.945	CI: {0.10 ± 0.0013} P = 0.956	{0.10 ± 0.0011} P = 0.961
τ = 0.1n = 8	CI: {0.10 ± 0.0022} P = 0.949	{0.10 ± 0.0018} P = 0.950	CI: {0.10 ± 0.0007} P = 0.949	{0.10 ± 0.0007} P = 0.948
τ = 0.15n = 4	CI: {0.15 ± 0.0055} P = 0.957	{0.15 ± 0.0049} P = 0.954	CI: {0.15 ± 0.0017} P = 0.948	{0.15 ± 0.0016} P = 0.955
τ = 0.15n = 8	CI: {0.15 ± 0.0034} P = 0.951	{0.15 ± 0.0038} P = 0.949	CI: {0.15 ± 0.0011} P = 0.957	{0.15 ± 0.0011} P = 0.956
τ = 0.2n = 4	CI: {0.2 ± 0.0071} P = 0.940	{0.2 ± 0.0076} P = 0.949	CI: {0.2 ± 0.0022} P = 0.954	{0.2 ± 0.0024} P = 0.957
τ = 0.2n = 8	CI: {0.2 ± 0.0050} P = 0.944	{0.2 ± 0.0054} P = 0.951	CI: {0.2 ± 0.0016} P = 0.957	{0.2 ± 0.0019} P = 0.948

Next, for each combination of (n, τ, L), we have also estimated the probability of the true τ falling into the estimated confidence interval (coverage probability) as the total number of time the true value was in the estimated interval (among N = 1, 000 repetitions) divided by N. The results are in Table 1 at the CIAsymp columns.

Next, we compared the performance of our method with the method of computing the confidence interval using simulation estimation of variance4 using the same simulated data. We have briefly described their method in Section 1. As we had done for our method, we estimated confidence interval for each of N = 1, 000 repetitions and for each combination of (n, τ, L). (Note that, this involved resimulating M = 100 datasets for each estimated τ̂MLE and estimating τ from the resimulated dataset, and thus a total of N × M = 105 simulations and estimations.) Then, we have computed the average estimated 95% confidence interval and coverage probability from N = 1, 000 simulated datasets for each combination of (n, τ, L). That is, after estimating the confidence intervals for each dataset, we computed the average of the confidence intervals over N = 1, 000 estimated values for each combination of (n, τ, L) (reported in Table 1). The coverage probabilities are estimated as before as the total number of time the true value was in the estimated interval divided by N. The results are shown in Table 1.

For both CIAsymp and simulated CI the estimated ranges have small lengths. As expected, the length becomes smaller with larger n and larger L. For τ = 0.01 the ranges have radii 0.0010 (for n = 4, L = 10, 000, CIAsymp) to 0.0001 (for n = 8, L = 100, 000, both methods). The length of the radii increases with τ. For τ = 0.02 the ranges have radii 0.0015 (for n = 4, L = 10, 000, CIAsymp) to 0.0002 (for n = 8, L = 100, 000, both methods). For τ = 0.05 the ranges have radii 0.0029 (for n = 4, L = 10, 000, simulated CI) to 0.0004 (for n = 8, L = 100, 000, CIAsymp). The radius increases up with an approximate proportionality with the true value of τ. For τ = 0.1 and 0.2 the ranges have radii 0.0041 and 0.0076, respectively, (for n = 4, L = 10, 000 simulated CI) to 0.0007 (for n = 8, L = 100, 000 both methods) and 0.0016 (for n = 8, L = 100, 000 CIASYMP), respectively.

A comparison of the two methods reveals no significant difference in the coverage probabilities of the lengths of the intervals. In Wilcoxon signed-rank tests we found no significant difference between coverage probabilities of the two methods (P > 0.3) as well as between lengths of the intervals (P > 0.9). However, computing time was an order of magnitude faster for CIAsymp than for the simulated CI method. Once we have τ̂MLE and θ̂MLE, it takes less than 30 minutes to compute the CI using CIAsymp using a R v3.0.2013-05-1217 code in a 2.54 GHz Dual Core processor. Using the same computer and same version of R, simulated CI (with at least M = 100 simulations for estimating the variance) takes more than a day to compute. This is because each simulation involves simulating L independent loci and maximization of the multi-locus likelihood over τ and θ using the allele-counts from these loci.

Results: Applications to HapMap Data

For the purpose of demonstrating our method and comparing its performance with known results, we applied our method to a subset of HapMap data.15 Specifically, we estimated a confidence interval for the divergence time between HCB (Han Chinese from Be-jing, China) and CEU (United States residents of northern and western European ancestry) populations using 112 independent SNP loci in Chromosome 19. To reduce the computational load, we only considered a random sub-sample of 8 unrelated haploids from each population. Our data consist of the allele-count in 8 individuals in each population for 112 SNP loci. The 112 SNPs were selected at least 0.5 MB apart from each other to ensure the independence of the coalescent trees in the lineages of 8 haploids.

The estimated scaled divergence time was τ̂ = 0.16, and the estimated 95% confidence intervals was. [0.125, 0.195]. These numbers, when transformed in years using an overall effective population size of 3,100 (see, for example)18 for HCB and CEU populations and generation time of 25 years, produce an estimated range of between 19,375–30,225 years ago. If an overall effective population size of 4,000 is assumed, then this produces as estimated range of between 25,000–39,000 years ago. Note that these ranges roughly match with other recent estimates of divergence time between HCB and CEU (see, for example).19

Discussion

We have presented a formula for computing divergence time range using the MLE on a coalescent model. As MLE is asymptotically effcient, our estimated confidence interval has a high degree of accuracy. We have also presented formulae for first and second order mixed derivatives of likelihood and log-likelihood, which is useful for computation of the MLE.

Our simulation study shows that our method produces small confidence intervals with appropriate coverage of 95%. Thus, it reduces the amount of uncertainty regarding the actual divergence time. Although MLE of divergence time has been computed before using our model, an expression for the range of divergence time has never been derived because of the complexity of the expression. By deriving this expression we made it possible for the range of the divergence time to be estimated, and by evaluating its performance we established its usefulness.

The radii of the confidence intervals and the coverage probabilities produced by our methods have been shown to be statistically equivalent to that of the simulation-based estimator. However, we found that our method is an order of magnitude faster. This is expected because in a simulation or resampling-based method the data needs to be simulated a large number of times, and then confidence intervals needs to be reestimated for each resimulation or resample.

Certain assumptions are made about the underlying model. The accuracy of the esti-mated variance and range may be affected if the data do not fit the underlying model. For example, a molecular clock is assumed. That is, the scaled time τ between A and O is assumed to be same as that between B and O. However, if the effective population sizes in populations a and b are very different then this will induce a bias in the estimated variance, as well as the estimated range. However, it is straightforward to extend our method for models without molecular clock. If a molecular clock is not assumed then one needs to separately estimate scaled divergence times τ and τ for Populations a and b. It is straightforward, albeit tedious, to modify Eqs. (6–16)) for parameters (τ, τ, θ), rather than (τ,θ). Thus, analogous expressions for the variance and the range of (τ, τ, θ) can be computed if a molecular clock is not assumed.

Another departure from our model would be from the assumption of biallelic loci. Following the same principle of tracking coalescing lineages back to the MRCA, and then following the allele-types to the present time, one can modify Eqs. 1–5 for more than two alleles. However, for more than two alleles, one needs to replace the symmetric beta-binomial distribution at the MRCA (which arises from assuming a symmetric beta distribution for biallelic allele-frequency spectrum) with its “multiple-count” version: a symmetric Dirichlet-Multinomial distribution (which arises from assuming a symmetric Dirichlet distribution for multiallelic allele-frequency spectrum). Once the likelihood is computed using modified version of Eqs. 1–5, then Eqs. 6–16 can be modified for estimating the variance and the range from the multiallelic likelihood.

We assumed that the effect of mutation between the point of divergence and the present is negligible. This is an appropriate assumption if the two populations are closely related (and consequently the divergence time is small). For large divergence times this assumption is not appropriate, and the effect of mutation may create an amount of difference between the two population that is higher than expected from an ignore-mutation model”. As a result, an erroneously larger divergence time may be estimated, which will also induce an upward bias in the estimated variance and the estimated range.

previous studies have suggested that this model is appropriate for populations within the same species.5–8 This designation is intended to serve as an upper bound for the divergence time to be modeled. In the absence of a mathematically concrete upper bound for divergence time for this model, we also use this convention. Thus, although we do not have a more concrete upper bound, we suggest using our methods for doing inference on divergence time between populations within the same species. (A lower bound is not necessary as the model fits well for small divergence times as the amount of mutation will be very small in such cases.) Moreover, a recent article13 introduces a version of the4,6 model that takes into account the effect of mutation. An appropriate extension of our methods to the13 model needs to be derived for estimating the variance and the range of larger divergence times. With such an extension, our methods could be used to estimate variance and range of species divergence times.

A possible disadvantage of our method would be in a scenario where significant amount of migration has taken place between the two populations after the divergence. In the presence of migration, there will be less variation than expected in a “no migration” scenario. Consequently, divergence time will be underestimated, which will also induce a downward bias in the estimated variance and the estimated range. Another limitation is that our method uses asymptotic approximation. Thus, our method is only applicable when we have a large number of independent loci. To quantify a large number of loci, a rule of thumb used by statisticians is that the asymptotic methods are to be used if there are at least 30 independent data-points. Thus, as long as we have allele counts for at least 30 independent loci, our method can be used. As most modern datasets have more than 30 independent loci, this is not a real limitation. Note that, as the effect of mutation after divergence is assumed to be negligible and is ignored, the difference in mutation rate between the two populations does not play a part in our method.

We applied our method to estimate a range for the divergence time between CEU and HCB populations from the HapMap data. Our estimates are found to be in the same range as a recently estimated divergence time between these two populations.19

A possible extension of this method could be computation of the confidence interval of the branch-length of a phylogenetic tree given the tree-topology. This would require creating an algorithm that could effciently and systematically compute the derivatives of the likelihood of the tree. Another possible extension could be using dependent loci to estimate the range of divergence time. This will have potential applications in high-resolution NextGen Sequencing data.