Source code for moments.LD.Demographics3D

Three-population demographic models.
import numpy as np
import moments.LD

[docs]def out_of_Africa(params, rho=None, theta=0.001, pop_ids=["YRI", "CEU", "CHB"]): """ The Gutenkunst et al (2009) out-of-Africa that has been reinferred a number of times. :param params: List of parameters, in the order (nuA, TA, nuB, TB, nuEu0, nuEuF, nuAs0, nuAsF, TF, mAfB, mAfEu, mAfAs, mEuAs). :param rho: Recombination rate or list of recombination rates (population-size scaled). :param theta: Population-size scaled mutation rate. :param pop_ids: List of population IDs. """ if pop_ids is not None and len(pop_ids) != 3: raise ValueError("pop_ids must be a list of three population IDs") ( nuA, TA, nuB, TB, nuEu0, nuEuF, nuAs0, nuAsF, TF, mAfB, mAfEu, mAfAs, mEuAs, ) = params y = moments.LD.Demographics1D.snm(rho=rho, theta=theta) # integrate modern human branch with expansion y.integrate([nuA], TA, rho=rho, theta=theta) # split into African and Eurasian branches and integrate y = y.split(0) mig_mat = [[0, mAfB], [mAfB, 0]] y.integrate([nuA, nuB], TB, m=mig_mat, rho=rho, theta=theta) # split Eurasian into CEU and CHB y = y.split(1) nu_func = lambda t: [ nuA, nuEu0 * np.exp(np.log(nuEuF / nuEu0) * t / TF), nuAs0 * np.exp(np.log(nuAsF / nuAs0) * t / TF), ] mig_mat = [[0, mAfEu, mAfAs], [mAfEu, 0, mEuAs], [mAfAs, mEuAs, 0]] y.integrate(nu_func, TF, m=mig_mat, rho=rho, theta=theta) y.pop_ids = pop_ids return y