Wright–Fisher exact solver (WFES): scalable analysis of population genetic models without simulation or diffusion theory

Abstract

MotivationThe simplifying assumptions that are used widely in theoretical population genetics may not always be appropriate for empirical population genetics. General computational approaches that do not require the assumptions of classical theory are therefore quite desirable. One such general approach is provided by the theory of absorbing Markov chains, which can be used to obtain exact results by directly analyzing population genetic Markov models, such as the classic bi-allelic Wright–Fisher model. Although these approaches are sometimes used, they are usually forgone in favor of simulation methods, due to the perception that they are too computationally burdensome. Here we show that, surprisingly, direct analysis of virtually any Markov chain model in population genetics can be made quite efficient by exploiting transition matrix sparsity and by solving restricted systems of linear equations, allowing a wide variety of exact calculations (within machine precision) to be easily and rapidly made on modern workstation computers.ResultsWe introduce Wright–Fisher Exact Solver (WFES), a fast and scalable method for direct analysis of Markov chain models in population genetics. WFES can rapidly solve for both long-term and transient behaviours including fixation and extinction probabilities, expected times to fixation or extinction, sojourn times, expected allele age and variance, and others. Our implementation requires only seconds to minutes of runtime on modern workstations and scales to biological population sizes ranging from humans to model organisms.Availability and ImplementationThe code is available at https://github.com/dekoning-lab/wfesSupplementary information Supplementary data are available at Bioinformatics online.

The Rate of Molecular Evolution When Mutation May Not Be Weak

“…Theoretical predictions from the full rate of evolution were consistent with uniform mutation simulations, regardless of θ. In contrast to the weak-mutation rate of evolution, which is not even consistent with the one-at-a-time model it assumes, when Wright-Fisher simulations were performed by allowing mutation to occur at any time ( Table 1, Simulation II), simulation averages were very close to their exact predictions computed from equation 4 ('General Theory') using WFES (Krukov et al, 2016). Notably, under both models, the cumulative time spent waiting for extinctions can be even larger than the time spent waiting for mutations when θ is moderately large (e.g., when θ = 0.25; Table 1, Simulations I and II).…”

“…S1). Although simple closed-form expressions for the component quantities are not available using diffusion theory, they can be easily calculated directly from the underlying Markov model using efficient computational techniques we recently described (De Sanctis et al, 2017;Krukov et al, 2016).…”

“…In all cases, a population size of N = 100 was used to make simulations tractable; even for such a small population, simulations took many CPU days to run. Exact predictions were made by directly computing the expected time between fixations and its variance with WFES (Krukov et al, 2016), which took less than one second. The number of mutational origins for IBS/IBD fixations are shown in Table S1.…”

“…Let x be the current number of mutants in a Wright-Fisher population of size N , and p the initial number of mutants to arise on a background of x = 0. To directly calculate the substitution rate without invoking the theory presented above, we modified our program WFES (Krukov et al, 2016) by making x = 0 a transient state so that computation of the mean time to absorption from a starting count of zero (p = 0) will represent the mean time it takes to go from being 100% wildtype to 100% mutant (i.e., the time between fixations). Similarly, we calculated the variance of the time between fixations as the variance of the time to absorption under these same conditions.…”

“…Mean times and probabilities in equation 4 were computed directly from a computationally efficient analysis of the appropriate Wright-Fisher Markov model, using absorbing Markov chain methods we previously described (Krukov et al, 2016). To measure properties of expected allele frequency trajectories demarcated by visits to either of the extinction or fixation boundaries, the boundary states were treated as absorbing (except when directly calculating the mean and variance of the time between fixations; see above), even though they may be escaped by mutation.…”