Special reviews
On the use of kernel approximate Bayesian computation to infer population history
2015 Volume 90 Issue 3 Pages 153-162
Details
2015 Volume 90 Issue 3 Pages 153-162
Genetic variation is a product of mutation, recombination and natural selection along with a complex history involving population subdivision, gene flow and changes in population size. Elucidating the evolutionary forces that shape genetic differences among populations is a major objective of evolutionary genetics. Recent advances in high-throughput technology enable genomic data to be obtained from samples at a population-based scale. Further, the growth in computational power has facilitated extensive efforts to develop intensive simulation-based approaches with the aim of analyzing such large-scale data and making inferences about population history. Approximate Bayesian computation (ABC) provides a quantitative way to assess the goodness-of-fit of complex models that are based on previous knowledge and to estimate the parameters of interest that produce the observed data. The practical advantage of ABC is the application of Bayesian inference to any model without the need to derive a likelihood function. ABC has rapidly become popular in ecology and evolutionary studies due to the contribution it has made to improving computational efficiency over the past decade. This review provides a brief overview of the background of ABC, including potential biases in estimation due to the assumptions and approximation involved, followed by an in-depth review of one of the recently developed ABCs, “kernel ABC,” with an explanation of how to overcome these biases. Finally, the application of kernel ABC to the inference of demographic history is summarized.
One of the main goals in evolutionary biology is to make inferences about population history. All living individuals can be traced back to a common ancestor through their genomes. Genetic variation is the product of natural selection and stochastic processes involving demographic events, such as population subdivision, migration and changes in population size. Statistical inferences provide useful and quantitative methods for understanding such a complex historical process. Model-based approaches aim to estimate parameters of interest that may explain how the observed data were generated. Several methods have been proposed for demographic and phylogenetic inferences (Cavalli-Sforza and Edwards, 1967; Thompson, 1973; Felsenstein, 1981). In particular, population genetics has been strongly influenced by the development of the coalescent theory that describes a genealogy (i.e., the path from present individual sequences to a most recent common ancestor) backwards in time under a stochastic model (Fig. 1A) (Kingman, 1982; Tajima, 1983). Bayesian inference based on the coalescent theory allows for a powerful and flexible framework for model-based approaches (Beaumont and Rannala, 2004; Beaumont et al., 2010) and has been successfully applied to inferring demographic history (Wilson and Balding, 1998; Nielsen and Wakeley, 2001; Drummond et al., 2002; Rosenberg and Nordborg, 2002; Rannala and Yang, 2003; Wilson and Rannala, 2003; Hey and Nielsen, 2007).
Illustration of the coalescent process and genealogy. (A) A simplified schematic of DNA evolution. Each row corresponds to a single generation, and each circle denotes a DNA sequence. This is a constant size model, so every generation has the same number of samples. All five sequences colored with orange were sampled from the current generation. The parental sequences that lead to the present day samples are also indicated in orange. A pair of sequences coalesces backwards in time because they share an ancestor. The last point of coalescence is a common ancestor for all sequences, which is the most recent common ancestor (MRCA). (B) The genealogy for the five sequences sampled at the present day. The waiting time before which the first coalescence happens among the five sequences is shown as T5. The waiting times of the subsequent coalescence events correspond to T4, T3 and T2. The time to MRCA (TMRCA) is the sum of T5, T4, T3 and T2. (C) Since neutral mutations do not affect the genealogy, they are superimposed on the tree. The five DNA sequences are denoted as A–E. Six mutations (Seg), shown as stars, occur in the genealogy. The site frequency spectrum (SFS) is used to calculate how many mutations are present on one (ξ1), two (ξ2), three (ξ3) and four (ξ4) sequences in the genealogy. There are three mutations of ξ1, two mutations of ξ2 and one mutation of ξ3.
Bayesian inference consists of three main steps: (1) setting up a model, (2) conditioning the model on observed data, and (3) evaluating the fit of the model (Gelman et al., 2004). The first step is performed by constructing models based on previous knowledge and experience. Such information can be used to define a likelihood function (i.e., the probability of the observed data under a set of parameters) and prior distribution (i.e., information about a parameter before obtaining the observed data). The second step is to estimate a posterior distribution that is proportional to the product of the likelihood and prior. The posterior distribution is characterized by a point estimate (mode or mean) and an interval estimate (a set of true parameters with a certain probability), which can be used to interpret how the parameter changes by conditioning using the observed data. However, several hypotheses are often proposed of an evolutionary history that produces the observed data. Hence, comparing models and evaluating their fit with the observed data is the fundamental step in making inferences about population history. The third step is to evaluate the goodness-of-fit of the models and to identify the most parsimonious scenario. Models are a general approximation of a true history and can be improved by re-formulating (i.e., going back to the first step) and repeating these steps.
One of the difficulties in applying Bayesian inference to evolutionary studies is computing the likelihood function. This difficulty is due to the complexity of the models (e.g., history from multiple populations and changes in population size) and of the data (e.g., genome-wide data from multiple populations). In population genetics, a genealogy is treated as a nuisance parameter and the likelihood is averaged over the possible genealogy states under a set of parameters (e.g., population size and time of population subdivision) (Rosenberg and Nordborg, 2002). However, if the model of interest is complex, the structure of the genealogy becomes complicated and the parameters constituting the genealogy become broader, making it difficult to determine the likelihood function. Further, the explosive growth of high-throughput technology has allowed genomic data to be obtained at a population-based scale. Since previous techniques were based on simple models (e.g., a panmictic population of constant size) and applied to data of a limited size (e.g., single locus or multiple loci from a limited number of samples), these techniques are computationally expensive or completely infeasible for evaluating likelihood in genomic data. Therefore, simulation-based methods have been developed that can handle such large-scale data without the need to derive the likelihood function (Marjoram and Tavare, 2006).
The development of simulation-based Bayesian computation has been facilitated by the coalescent theory, which provides efficient methods to simulate samples of DNA fragments under many genetic and demographic scenarios (Hudson, 1990). Here, θ is a parameter of interest and D is the observed data. The goal is to sample θ from a posterior distribution that is defined as π(θ|D)
One application of rejection sampling is to infer the time to the most recent common ancestor (TMRCA) based on the number of segregating sites (Seg) (Tavare et al., 1997). TMRCA is the sum of the coalescent times in a genealogy consisting of n samples (Fig. 1B). Ti is the waiting time of coalescence between i sequences (i = n to 2), and it follows an exponential distribution with a parameter of i(i − 1)/2. The likelihood represents the probability of an observation, Seg, under a set of parameters, T = (Tn, ..., T2), and it follows a Poisson distribution. Thus, the algorithm of the rejection method is composed of two steps (Tavare et al., 1997). First, a set of parameters, t = (tn, tn−1, …, t2), is generated from the exponential distribution, and the total branch length (l) in the genealogy is calculated from l = 2t2 + … + (n − 1)tn−1 + ntn. Second, the simulated data is accepted with the likelihood computed from the ratio between the Poisson point probabilities of Po(
Demand for testing complicated models and recent improvements in computational power have facilitated the development of simulation-based and likelihood-free methods in Bayesian inference. Approximate Bayesian computation (ABC) is a flexible framework for estimating a posterior distribution without making explicit use of the likelihood function (Fu and Li, 1997; Tavare et al., 1997; Weiss and von Haeseler, 1998; Pritchard et al., 1999; Beaumont et al., 2002). As the scale of genetic data expands, approximation must be involved to summarize the data into a set of statistics. Summary statistics are values that capture complicated data in a simpler way by reducing the dimensionality of the data. In population genetics, several useful statistics have been developed to evaluate patterns of genetic variation, including the number of segregating sites (Seg), nucleotide diversity (π) (Nei and Li, 1979), Tajima’s D (TD) (Tajima, 1989), and Fu and Li’s D (FLD) and F (FLF) (Fu and Li, 1997). The aim of ABC is to approximate the posterior distribution of parameters and determine the evolutionary history behind the observed data.
The most basic form of ABC, the rejection algorithm, involves sampling parameters from the posterior distribution approximated using a set of summary statistics, s = (s1, s2, …, sm), where m is the number (i.e., dimension) of summary statistics.
A1. Generate θ′ from π(·).
A2. Simulate D′ using θ′.
A3. Accept θ′ if s(D) = s(D′), and go to A1.
The resulting posterior distribution is π(θ|s(D)), which is the approximation of π(θ|D). Although summary statistics are used to reduce the computational cost, this algorithm requires a lot of time to accept the exact same statistics as the observation. Another approximation involves using a tolerance (η) that defines the degree of similarity between the observed and simulated summary statistics and relaxes the conditions for accepting the parameters. Step A3 is replaced with “Accept θ′ if d(s(D), s(D′)) < η, and go to 1”, where d(·,·) is a function measuring the similarity between the two sets of summary statistics. The practical outcome of the ABC rejection algorithm is sampling from π(θ|d(s(D), s(D′)) < η).
To correct discrepancies between the simulated and observed summary statistics, local linear regression was incorporated into ABC (Beaumont et al., 2002). After accepting parameters by ABC rejection, this method assigns weights to the accepted parameters according to the similarity between observed and simulated summary statistics. Each of the weighted parameters can be taken as a random draw from π(θ|s(D)), and the posterior mean of a parameter is given as
ABC can be also applied to evaluate the goodness-of-fit of a model. In the case where there are two alternative models of M1 and M2 (e.g., constant population size vs. population expansion, or no migration vs. migration), a Bayes factor (BF) representing how much more likely one model is than the other is given as
Several technical improvements to standard ABCs have increased the accuracy of posterior estimates and computational efficiency (Marjoram et al., 2003; Sisson et al., 2007; Beaumont et al., 2009; Toni et al., 2009). However, most ABC-based methods are subject to biases from (1) a lack of sufficient statistics and (2) the use of non-zero tolerance. Summarizing data into a set of statistics implies disregarding useful information to some extent. Hence, the posterior distribution given summary statistics is not equal to the posterior distribution given the actual observed data (Marjoram and Tavare, 2006). In addition, another layer of approximation is introduced to reduce computational time by setting up positive values of tolerance. The flexibility and caveats of using ABC have been well summarized in several review articles (Marjoram and Tavare, 2006; Beaumont, 2010; Beaumont et al., 2010; Csillery et al., 2010; Sunnaker et al., 2013).
For the first source of bias, accepted parameters are the samples from an approximation of π(θ|D), unless a set of summary statistics is sufficient for the parameters of interest. Since the likelihood function is unknown in ABC, it is almost impossible to obtain sufficient statistics for genetic data. The availability of informative statistics may be limited to particular models. For instance, the number of segregating sites can be useful for estimating mutation rates but provides little information about demographic history, such as population bottlenecks and expansion (Hey and Machado, 2003). An intuitive way to overcome the lack of sufficient statistics is to increase the number of summary statistics and keep the loss of information to a minimum. The approximation of posterior distribution given data can be improved by adding summary statistics if the new set of summary statistics refines an existing set of statistics and provides more detailed information about the data (Nakagome et al., 2013a). For example, the site frequency spectrum (SFS) is a function of statistics commonly used in summarizing the pattern of genetic variation (Fig. 1C). SFS consists of the site frequency, ξi, which is defined as the number of segregating sites where a mutant nucleotide is present on i sequences in the sample (i = 1, …, n). Seg, π, TD, FLD, and FLF are calculated as a function of SFS;
Another concern when increasing the number of summary statistics is a reduction in the accuracy and stability of ABC. The acceptance rate sharply decreases using high-dimensional summary statistics, which may increase the variance of posteriors due to the small number of accepted parameters (Beaumont et al., 2002; Nakagome et al., 2013a). The space defined by a tolerance becomes the multi-dimensional sphere. Even though a similarity between simulated and observed summary statistics is small enough to be included within the tolerance, most of the simulated statistics are located on the boundary of the sphere. Consequently, it would take an enormous amount of time to find parameters that could produce summary statistics truly close to the observation (Hastie et al., 2009). Alternatively, various methods have been developed for choosing useful statistics from many kinds of summary statistics (Joyce and Marjoram, 2008; Wegmann et al., 2009; Blum and Francois, 2010; Nunes and Balding, 2010; Fearnhead and Prangle, 2012) depending on the model.
For the second source of bias, the choice of a positive tolerance involves a trade-off between the bias and the variance of a posterior estimate. With a sufficiently small tolerance, the resulting posterior distribution should be close to the target posterior distribution, which is π(θ|s(D)). As the acceptance rate decreases, the small number of accepted samples increases the variance of posterior distribution. In contrast, a large tolerance allows prior samples to be easily accepted, which reduces the variance but increases the bias from the prior. To this end, regression ABC has been developed to reduce the variance of the posterior estimates without decreasing the accepting conditions. Other studies have addressed the bias introduced by non-zero tolerance by proposing schemes to successively reduce tolerance using a sequential Monte Carlo algorithm (Sisson et al., 2007; Beaumont et al., 2009; Drovandi and Pettitt, 2011). Given a measuring error in data or a modeling error in statistical hypotheses, the optimal tolerance has been shown to be η ≠ 0 (Fearnhead and Prangle, 2012; Wilkinson, 2013). Noisy ABC has been developed to characterize and compensate for the bias from non-zero tolerance by introducing a specific form of errors to sampling parameters from a posterior distribution.
To minimize the biases involved in the standard ABCs, more summary statistics can be used that do not involve any arbitrary choice of simulated data. Kernel methods provide systematic ways of analyzing high-dimensional data by transforming the data into high or potentially infinite space, which is called reproducing kernel Hilbert space (RKHS) (Hofmann et al., 2008; Fukumizu, 2010). In this space, the similarity of summary statistics is evaluated as the inner product, and each inner product is computed by a kernel function such as the Gaussian RBF kernel, which is defined as k(si,sj) = exp
B1. Generate θ′ from π(·).
B2. Simulate D′ using θ′.
B3. Calculate s(D′), and go to B1.
B4. Compute
In this last step (B4), kernel ABC is used to compute a posterior mean from the sum of weighted parameters, and each weight is calculated based on the similarity between the observed and simulated summary statistics.
There are two main advantages of kernel ABC. First, a posterior mean given observed summary statistics is estimated without rejection. Hence, the number of simulations required is much smaller (e.g., 10,000 or 20,000) than in the standard ABCs (e.g., one million to accept 100 or 1,000 samples). In addition, kernel ABC is able to accurately estimate posterior means with a small number of simulations, so it requires less computational time than other ABCs (Nakagome et al., 2013a). Second, the estimator of the posterior mean in kernel ABC is consistent with the posterior mean given observed summary statistics,
This idea of kernel methods has been extended to model selection based on an approximated Bayes factor (aBF) (Osada et al., 2013). Applying model selection for ABC needs a long time until enough simulations are accepted if high-dimensional summary statistics are used. Therefore, Osada et al. (2013) proposed an aBF that takes a ratio of P(s(D)|M1) to that of P(s(D)|M2) computed from a kernel density estimate of P(s(D)|M) =
A large number of summary statistics can provide detailed information about data. The site frequency spectrum (SFS) is a refinement of the number of segregating sites (Seg) (Fig. 1C). For example, in Fig. 1, the genealogy has Seg = 6, while SFS = {3,2,1,0,0} represents how many mutations occur in the tips leading directly to each type of sequence or in the internal branches (i.e., shared branches among some of the sequences). To evaluate how the use of a large number of summary statistics improves estimating posterior means, two examples are given in this section under a simple and a complex demographic model.
The first model is the constant population size (Fig. 2A). The observed data are simulated under a population size of 10,000 with a mutation rate of 2.5 × 10−9 (per bp per generation), a sample of 100 sequences, and a 100-kb non-recombining region. The sequence data are summarized into Seg = 49 or SFS = {28,6,4,3,2,1,5}, in which SFS is divided into seven bins with frequencies of 0–8%, 8–16%, 16–24%, 24–32%, 32–40%, 40–48% and 48–100%. The parameter of interest is the population size (N), and the prior is given as the log-normal distribution (LN) of N~LN(10000, 100002). Under the constant size model, the likelihood function is analytically tractable, and the true posterior mean,
(A) A simple model of constant population size of 10,000 (Nakagome et al., 2013a). The posterior means in the table below (A) are estimated by the Bayesian computation with the full set of data (D) and by kernel ABC with SFS (sSFS), or Seg (sSeg), respectively. (B) A complex model of population divergence in which the ancestral population with a size of 100,000 (NANC) split into two populations with the same population size of NANC 25,000 generations ago (Osada et al., 2013). The total number of sequences is 30 and 18 in populations A and B. Sequence data are summarized into haplotype frequency spectrum (HFS) as well as SFS. The bin sizes of SFS and HFS are defined as 0, 1–4, 5–8, 9–12 and 13–30 for population A and 0, 1–3, 4–6, 7–9 and 10–18 for population B. The posterior means estimated by kernel ABC are shown in the table below (B) that includes the four parameters.
The second example is a complex model that assumes that the population A split from the population B (Fig. 2B). This model is constructed as the example of macaque demography (Osada et al., 2013) and consists of four parameters: population size (NA, NB, and NANC) and the time of population divergence (T). DNA sequence data include 26 unlinked loci with a length of 700 bp from two populations (nA = 30 and nB = 18). The observed data are simulated under the following conditions: constant population size (NA = NB = NANC = 100,000), T = 25,000 generations ago, a mutation rate of 7.2 × 10−9, and a recombination rate of 7.3 × 10−9. A one-dimensional SFS may reflect the pattern of mutation in each population, while the amount of genetic variation shared between two populations could be also useful for estimating when the population split occurred. Further, the effects of recombination are incorporated into this model, so haplotype-based information may also be necessary to make inferences about the complex history. Thus, the observed data is summarized into a two-dimensional SFS and haplotype frequency spectrum (HFS) that represents how many mutations and haplotypes are specific to each population and how many are shared between the populations. The posterior means are almost consistent with the true values, NANC = 102,064, NA = 104,553, NB = 98,028, and T = 25,602, even though the prior means are twice as large as the true values (Osada et al., 2013). The number of summary statistics may need to be increased in correspondence with the complexity of the model, and kernel ABC is powerful enough to handle high-dimensional statistics and untangle complex evolutionary history.
One of the main goals in evolutionary biology is to understand the demographic processes that natural populations have experienced over time. The application of kernel ABC has yielded new insights into the evolution of wild animals and human populations. Two such examples include the divergence of brown bears and polar bears (Nakagome et al., 2013c) and the migration of humans into the Ryukyu Islands, which are the southernmost islands of the Japanese archipelago (Sato et al., 2014).
The speciation of polar bears from brown bears is of considerable interest to evolutionary geneticists. Their phylogenetic relationship is discordant between mitochondrial DNA (mtDNA) and autosomal DNA (aDNA). A paraphyletic model in which the polar bear lineage is derived from brown bear diversity is supported by mtDNA sequences (Talbot and Shields, 1996), while a study based on multiple loci in aDNA sequences proposed a monophyletic model in which polar bears split from the ancestral lineage of all brown bears (Hailer et al., 2012). Further, the majority of the discordance is due to stochastic processes in which ancestral polymorphisms and subsequent incomplete lineage sorting resulted in a unique gene genealogy (Nakagome et al., 2013b). These studies highlight the general issues in inferring evolutionary history when different genes produce different trees. This genealogical discordance represents a classic conundrum in ecology and evolutionary biology of whether phylogenetics or population genetics should be applied to DNA sequence data. Nakagome et al. (2013b) employed population genetics approaches based on the coalescent model. The estimates from kernel ABC show a large effective population size of mtDNA in brown bears and highlight significant deviation in the effective population size ratio of mtDNA to aDNA from the theoretically expected ratio (Nakagome et al., 2013c). This suggests that ancestral brown bear populations were structured by behavioral traits including female philopatry and male-biased migration (Fig. 3A). Thus, ancestral polymorphisms and sex-biased migration are likely to contribute to conflicting branching patterns of brown bears and polar bears across aDNA genes and mtDNA.
Applications of kernel ABC for the inference of population history. (A) Evolutionary history of brown bears and polar bears (Nakagome et al., 2013c). The initial divergence of the brown bear population occurred 1.37 million years ago (MYA). Geographic lineages of brown bears gradually diverged from the ancestral lineage of other brown bears and polar bears and expanded into broad geographic regions. Afterwards, the polar bear lineage appeared from the ancestral lineage at 0.31 MYA and differentiated into the arctic region. During this process, male brown bears migrated among the sub-populations (male-biased migration), while female brown bears remained in each of the sub-populations (female philopatry). Thus, geographic differentiations of brown bear lineages have been shown in the mtDNA genealogy, although these patterns were likely homogenized in aDNA genealogies by male-biased migration between geographic regions (Nakagome et al., 2008, 2013c). (B) Demographic models of population divergence between people living on Okinawa and the Miyako Islands (Sato et al., 2014). The two populations are assumed to have split at different times (T) sampled from different priors: T~LN(50, 502) in the recent model, T~LN(100, 1002) in the intermediate model and T~LN(400, 4002) in the ancient model. Migrations are also incorporated into the models. The table represents log10-scaled aBFs calculated from the ratios of the models to the highest aML model (“Recent” and “0”). The models of “Recent” and “0” and “Intermediate” and “10−3” show significantly higher aMLs than the other models. Assuming one generation is approximately 25 years, the divergence of the Okinawa and Miyako populations is likely to have occurred 1,800 or 2,600 years ago, supporting the “replacement” model.
It is also of great interest to understand how human populations have expanded from continents to islands since the human migration out of Africa. Contemporary island populations may be descendants of ancient people who left their oldest bones on the island (“population continuity”) or they may be recent migrants who have come from continents and replaced the ancient people (“replacement”). The peopling of the Ryukyu Islands is one of the models used to address this issue. Sato et al. (2014) found that the people living on the Miyako Islands had genetically diverged from those on Okinawa island. Meanwhile, 26,000-year-old human remains have been found on the Miyako Islands (Sakura, 1985). If the contemporary people on the Miyako Islands are direct descendants of these ancient people, then the time of divergence between the Miyako and Okinawa populations should have occurred within the Pleistocene. Sato et al. (2014) proposed three models in which the divergence occurred in recent, intermediate or ancient times (Fig. 3B), and the recent divergence model was significantly more likely than the other models based on aBF. Incorporating migration between the two populations would be expected to substantially push back the time of divergence. Although marginal likelihoods were estimated for several models with weak and strong migration, recent divergence with no or weak migration was still strongly supported by model selection (Sato et al., 2014). Assuming that one generation is approximately 25 years, the divergence of the Okinawa and Miyako populations is likely to have happened 1,800 or 2,600 years ago. Therefore, the first inhabitants of the Ryukyu Islands during the Pleistocene were likely replaced with recent migrants.
In conclusion, ABC is a flexible and powerful method for Bayesian inference. There is no need to formulate a likelihood function, and ABC can be applied for any complex model. However, caution is required for the reliable use of ABC due to potential biases in summarizing and choosing data. Kernel ABC enables high-dimensional summary statistics to be used without involving any arbitrary choice of simulation data. The recent development of high-throughput technology dramatically reduces the cost of collecting large data sets. Further, population-scale genomic data have become available in both model and non-model organisms. The challenges remain in extracting useful information from the data and analyzing the data statistically and quantitatively. ABC is a promising statistical framework, and the application of kernel methods to ABC is a starting point to address these questions and to understand the complex history of demographic processes and of adaptation to local environments.
Thanks to Shuhei Mano and two anonymous reviewers for helpful comments on this manuscript. The author was supported by the Japan Society for the Promotion of Science Postdoctoral Fellowships for Research Abroad (26-541).
