Avoiding misleading estimates using mtDNA heteroplasmy statistics to study bottleneck size and selection

Abstract Mitochondrial DNA heteroplasmy samples can shed light on vital developmental and genetic processes shaping mitochondrial DNA populations. The sample means and sample variance of a set of heteroplasmy observations are typically used both to estimate bottleneck sizes and to perform fits to the theoretical “Kimura” distribution in seeking evidence for mitochondrial DNA selection. However, each of these applications raises problems. Sample statistics do not generally provide optimal fits to the Kimura distribution and so can give misleading results in hypothesis testing, including false positive signals of selection. Using sample variance can give misleading results for bottleneck size estimates, particularly for small samples. These issues can and do lead to false positive results for mitochondrial DNA mechanisms—all published experimental datasets we re-analyzed, reported as displaying departures from the Kimura model, do not in fact give evidence for such departures. Here we outline a maximum likelihood approach that is simple to implement computationally and addresses all of these issues. We advocate the use of maximum likelihood fits and explicit hypothesis tests, not fits and Kolmogorov–Smirnov tests via summary statistics, for ongoing work with mitochondrial DNA heteroplasmy.

An alternative, maximum likelihood fit gives us p = 0.5 (95% c.i.s 0.058-0.941), and a point estimate of b ≃ 0 (and corresponding n = 1) -all readily interpretable as well-behaved parameters of the Kimura distribution. Confidence intervals on nb can be estimated via taking likelihood derivatives or resampling, the latter of which is better behaved for this small, extreme dataset (see below), and gives nb between 1 and 13.6 with 95% confidence.

Numerical issues
An issue can arise when numerically maximising the likelihood associated with homoplasmic data. Because the parameter b is constrained to lie on the interval [0,1], we either bound the domain of optimisation or use a transformation to cast the real line onto that interval. In both cases, the gradients in the associated Hessian matrix behave poorly numerically when b approaches 0 -as is the case for homoplasmic data, where b = 0 is the point estimate. In these cases, we found that using bootstrap resampling to characterise confidence intervals for b was the better approach. Of course, bootstrapping could be used for confidence intervals in non-homoplasmic cases too, but for reasons of computational speed and accuracy we suggest Fisher information for estimating p in all cases and for estimating b in cases of heteroplasmy, and bootstrapping for estimating b in cases of homoplasmy.
Another numerical issue involves minimising the KS distance over parameterisations.
Because the KS distance surface as a function of parameter p and b can be non-convex (Fig. 1Aiii), it is hard to guarantee that the best solution has been found: the local optimum found may depend on initial conditions. In our implementation we perform two searches, using the method-of-moments and maximum likelihood parameter estimates as two different initial conditions, and reporting the best solution found from these two searches.
We found this to generally identify good solutions based on exhaustively visualising the full surface.

Kimura fourth central moment
The original article from Kimura (1955) gives an expression for the nth moment about zero of the Kimura distribution: where we have used the b = exp(-t / Neff) substitution from Wonnapinij et al. (2008).

Supplementary Fig. A. Influence of small observation perturbations on reported p-values from the KS test.
Ranges of estimated p-values from the WCS-K approach, with small amounts of noise added to given observations. Originally reported p-values are given by grey dots. MoM, method-of-moments fitting (as is typically used); minKS, minimum KS distance fitting. Even small amounts of noise in reported observations (much less than 1%) dramatically expand the range of p-values that can arise from the WCS-K approach. MinKS fits, providing better fits to the Kimura distribution than MoM fits, universally give higher p-values.
Supplementary Figure B. Example of synthetic data that significantly departs from the Kimura distribution. The heteroplasmy measurements consist of the pair (0,0.9) repeated 50 times. As in Fig 1, (i) probability distributions for each fit; (ii) comparison of cumulative distribution functions (CDFs) from the data and for each fit; (iii) the KS distance between empirical observations and theoretical distribution with parameters of the Kimura distribution. mKS, KS distance for each method of fitting; p-values are from the Monte Carlo KS test in WCS-K. As the sample size is higher in this case, the fit using the method of moments and that using maximum likelihood provide more similar parameter estimates.
Supplementary Figure C. Different ways to estimate the uncertainty of variance in samples from (A) synthetic normal distribution; (B) synthetic Kimura distribution; (C) example mouse data. For the synthetic distributions, samples are ordered by sample size; for the mouse data, 12 example LE type samples are depicted for which the ordering is (for just visual clarity) ascending according to the mean sample variance. Error bars correspond to 2x the standard error of variance as calculated by each method. Analytic (h-statistic) and bootstrapped results are typically comparable; Kimura fits often provide a more conservative estimate of V(h) variance.