Module `dadi.Triallele.demographics`

Equilibrium, two epoch, and three epoch (bottleneck) model with selection (sig1, sig2) on the two derived alleles

Expand source code

"""
Equilibrium, two epoch, and three epoch (bottleneck) model with selection (sig1, sig2) on the two derived alleles
"""

from . import numerics, integration
import numpy as np
import dadi
from numpy import newaxis as nuax

from dadi.Triallele.TriSpectrum_mod import TriSpectrum

def equilibrium(params, ns, pts, sig1 = 0.0, sig2 = 0.0, theta1 = 1.0, theta2 = 1.0, misid = 0.0, dt = 0.005, folded = False):
    """
    Integrate the density function to equilibrium
    params = unused
    sig1, sig2 - population scaled selection coefficients for the two derived alleles
    theta1, theta2 - population scaled mutation rates
    misid - ancestral misidentification parameter
    dt - time step to use for integration
    folded = True - fold the frequency spectrum (if we assume we don't know the order that derived alleles appeared)
    """
    x = np.linspace(0,1,pts+1)
    sig1,sig2 = np.float(sig1),np.float(sig2)

    y1 = dadi.PhiManip.phi_1D(x,gamma=sig1)
    y2 = dadi.PhiManip.phi_1D(x,gamma=sig2)
    phi = np.zeros((len(x),len(x)))
    if sig1 == sig2 == 0.0 and theta1 == theta2 == 1:
        phi = integration.equilibrium_neutral_exact(x)
    else:
        phi = integration.equilibrium_neutral_exact(x)
        phi,y1,y2 = integration.advance(phi, x, 2, y1, y2, nu=1., sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)

    dx = numerics.grid_dx(x)

    try:
        ns = int(ns)
    except TypeError:
        ns = ns[0]

    fs = numerics.sample(phi, ns, x)
    fs.extrap_t = dt

    if folded == True:
        fs = fs.fold_major()

    if misid > 0.0:
        fs = numerics.misidentification(fs,misid)
    
    return fs
    
def two_epoch(params, ns, pts, sig1 = 0.0, sig2 = 0.0, theta1 = 1.0, theta2 = 1.0, misid = 0.0, dt = 0.005, folded = False):
    """
    Two epoch demography - a single population size change at some point in the past
    params = [nu,T,sig1,sig2,theta1,theta2,misid,dt]
    nu - relative poplulation size change to ancestral population size
    T - time in past that size change occured (scaled by 2N generations)
    sig1, sig2 - population scaled selection coefficients for the two derived alleles
    theta1, theta2 - population scaled mutation rates
    misid - ancestral misidentification parameter
    dt - time step to use for integration
    """
    nu,T = params

    x = np.linspace(0,1,pts+1)
    sig1,sig2 = np.float(sig1),np.float(sig2)

    y1 = dadi.PhiManip.phi_1D(x,gamma=sig1)
    y2 = dadi.PhiManip.phi_1D(x,gamma=sig2)
    phi = np.zeros((len(x),len(x)))

    # integrate to equilibrium first
    if sig1 == sig2 == 0.0 and theta1 == theta2 == 1:
        phi = integration.equilibrium_neutral_exact(x)
    else:
        phi = integration.equilibrium_neutral_exact(x)
        phi,y1,y2 = integration.advance(phi, x, 2, y1, y2, nu=1., sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)

    phi,y1,y2 = integration.advance(phi, x, T, y1, y2, nu=nu, sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)

    dx = numerics.grid_dx(x)

    try:
        ns = int(ns)
    except TypeError:
        ns = ns[0]

    fs = numerics.sample(phi, ns, x)
    fs.extrap_t = dt

    if folded == True:
        fs = fs.fold_major()
    
    if misid > 0.0:
        fs = numerics.misidentification(fs,misid)
    
    return fs


def three_epoch(params, ns, pts, sig1 = 0.0, sig2 = 0.0, theta1 = 1.0, theta2 = 1.0, misid = 0.0, dt = 0.005, folded = False):
    """
    Three epoch demography - two instantaneous population size changes in the past
    params = [nu1,nu2,T1,T2]
    nu1,nu2 - relative poplulation size changes to ancestral population size (nu1 occurs before nu2, historically)
    T1,T2 - time for which population had relative sizes nu1, nu2 (scaled by 2N generations)
    sig1, sig2 - population scaled selection coefficients for the two derived alleles
    theta1, theta2 - population scaled mutation rates
    misid - ancestral misidentification parameter
    dt - time step to use for integration
    """
    nu1,nu2,T1,T2 = params
    x = np.linspace(0,1,pts+1)
    sig1,sig2 = np.float(sig1),np.float(sig2)
    
    y1 = dadi.PhiManip.phi_1D(x,gamma=sig1)
    y2 = dadi.PhiManip.phi_1D(x,gamma=sig2)
    phi = np.zeros((len(x),len(x)))
    
    # integrate to equilibrium first
    if sig1 == sig2 == 0.0 and theta1 == theta2 == 1:
        phi = integration.equilibrium_neutral_exact(x)
    else:
        phi = integration.equilibrium_neutral_exact(x)
        phi,y1,y2 = integration.advance(phi, x, 2, y1, y2, nu=1., sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)
    
    phi,y1,y2 = integration.advance(phi, x, T1, y1, y2, nu1, sig1, sig2, theta1, theta2, dt=dt)
    phi,y1,y2 = integration.advance(phi, x, T2, y1, y2, nu2, sig1, sig2, theta1, theta2, dt=dt)

    dx = numerics.grid_dx(x)
    
    try:
        ns = int(ns)
    except TypeError:
        ns = ns[0]
    
    fs = numerics.sample(phi, ns, x)
    fs.extrap_t = dt

    if folded == True:
        fs = fs.fold_major()
    
    if misid > 0.0:
        fs = numerics.misidentification(fs,misid)
    
    return fs

def bottlegrowth(params, ns, pts, sig1 = 0.0, sig2 = 0.0, theta1 = 1.0, theta2 = 1.0, misid = 0.0, dt = 0.005, folded = False):
    """
    Three epoch demography - two instantaneous population size changes in the past
    params = [nu1,nu2,T1,T2]
    nu1,nu2 - relative poplulation size changes to ancestral population size (nu1 occurs before nu2, historically)
    T1,T2 - time for which population had relative sizes nu1, nu2 (scaled by 2N generations)
    sig1, sig2 - population scaled selection coefficients for the two derived alleles
    theta1, theta2 - population scaled mutation rates
    misid - ancestral misidentification parameter
    dt - time step to use for integration
    """
    nuB,nuF,T = params
    if nuB == nuF:
        nu = nuB
    else:
        nu = lambda t: nuB*np.exp(np.log(nuF/nuB) * t/T)
    
    x = np.linspace(0,1,pts+1)
    sig1,sig2 = np.float(sig1),np.float(sig2)
    
    y1 = dadi.PhiManip.phi_1D(x,gamma=sig1)
    y2 = dadi.PhiManip.phi_1D(x,gamma=sig2)
    phi = np.zeros((len(x),len(x)))
    
    # integrate to equilibrium first
    if sig1 == sig2 == 0.0 and theta1 == theta2 == 1:
        phi = integration.equilibrium_neutral_exact(x)
    else:
        phi = integration.equilibrium_neutral_exact(x)
        phi,y1,y2 = integration.advance(phi, x, 2, y1, y2, nu=1., sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)
    
    phi,y1,y2 = integration.advance(phi, x, T, y1, y2, nu, sig1, sig2, theta1, theta2, dt=dt)

    dx = numerics.grid_dx(x)
    
    try:
        ns = int(ns)
    except TypeError:
        ns = ns[0]
    
    fs = numerics.sample(phi, ns, x)
    fs.extrap_t = dt

    if folded == True:
        fs = fs.fold_major()
    
    if misid > 0.0:
        fs = numerics.misidentification(fs,misid)
    
    return fs

Functions

def bottlegrowth(params, ns, pts, sig1=0.0, sig2=0.0, theta1=1.0, theta2=1.0, misid=0.0, dt=0.005, folded=False)

Three epoch demography - two instantaneous population size changes in the past params = [nu1,nu2,T1,T2] nu1,nu2 - relative poplulation size changes to ancestral population size (nu1 occurs before nu2, historically) T1,T2 - time for which population had relative sizes nu1, nu2 (scaled by 2N generations) sig1, sig2 - population scaled selection coefficients for the two derived alleles theta1, theta2 - population scaled mutation rates misid - ancestral misidentification parameter dt - time step to use for integration

Expand source code

def bottlegrowth(params, ns, pts, sig1 = 0.0, sig2 = 0.0, theta1 = 1.0, theta2 = 1.0, misid = 0.0, dt = 0.005, folded = False):
    """
    Three epoch demography - two instantaneous population size changes in the past
    params = [nu1,nu2,T1,T2]
    nu1,nu2 - relative poplulation size changes to ancestral population size (nu1 occurs before nu2, historically)
    T1,T2 - time for which population had relative sizes nu1, nu2 (scaled by 2N generations)
    sig1, sig2 - population scaled selection coefficients for the two derived alleles
    theta1, theta2 - population scaled mutation rates
    misid - ancestral misidentification parameter
    dt - time step to use for integration
    """
    nuB,nuF,T = params
    if nuB == nuF:
        nu = nuB
    else:
        nu = lambda t: nuB*np.exp(np.log(nuF/nuB) * t/T)
    
    x = np.linspace(0,1,pts+1)
    sig1,sig2 = np.float(sig1),np.float(sig2)
    
    y1 = dadi.PhiManip.phi_1D(x,gamma=sig1)
    y2 = dadi.PhiManip.phi_1D(x,gamma=sig2)
    phi = np.zeros((len(x),len(x)))
    
    # integrate to equilibrium first
    if sig1 == sig2 == 0.0 and theta1 == theta2 == 1:
        phi = integration.equilibrium_neutral_exact(x)
    else:
        phi = integration.equilibrium_neutral_exact(x)
        phi,y1,y2 = integration.advance(phi, x, 2, y1, y2, nu=1., sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)
    
    phi,y1,y2 = integration.advance(phi, x, T, y1, y2, nu, sig1, sig2, theta1, theta2, dt=dt)

    dx = numerics.grid_dx(x)
    
    try:
        ns = int(ns)
    except TypeError:
        ns = ns[0]
    
    fs = numerics.sample(phi, ns, x)
    fs.extrap_t = dt

    if folded == True:
        fs = fs.fold_major()
    
    if misid > 0.0:
        fs = numerics.misidentification(fs,misid)
    
    return fs

def equilibrium(params, ns, pts, sig1=0.0, sig2=0.0, theta1=1.0, theta2=1.0, misid=0.0, dt=0.005, folded=False)

Integrate the density function to equilibrium params = unused sig1, sig2 - population scaled selection coefficients for the two derived alleles theta1, theta2 - population scaled mutation rates misid - ancestral misidentification parameter dt - time step to use for integration folded = True - fold the frequency spectrum (if we assume we don't know the order that derived alleles appeared)

Expand source code

def equilibrium(params, ns, pts, sig1 = 0.0, sig2 = 0.0, theta1 = 1.0, theta2 = 1.0, misid = 0.0, dt = 0.005, folded = False):
    """
    Integrate the density function to equilibrium
    params = unused
    sig1, sig2 - population scaled selection coefficients for the two derived alleles
    theta1, theta2 - population scaled mutation rates
    misid - ancestral misidentification parameter
    dt - time step to use for integration
    folded = True - fold the frequency spectrum (if we assume we don't know the order that derived alleles appeared)
    """
    x = np.linspace(0,1,pts+1)
    sig1,sig2 = np.float(sig1),np.float(sig2)

    y1 = dadi.PhiManip.phi_1D(x,gamma=sig1)
    y2 = dadi.PhiManip.phi_1D(x,gamma=sig2)
    phi = np.zeros((len(x),len(x)))
    if sig1 == sig2 == 0.0 and theta1 == theta2 == 1:
        phi = integration.equilibrium_neutral_exact(x)
    else:
        phi = integration.equilibrium_neutral_exact(x)
        phi,y1,y2 = integration.advance(phi, x, 2, y1, y2, nu=1., sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)

    dx = numerics.grid_dx(x)

    try:
        ns = int(ns)
    except TypeError:
        ns = ns[0]

    fs = numerics.sample(phi, ns, x)
    fs.extrap_t = dt

    if folded == True:
        fs = fs.fold_major()

    if misid > 0.0:
        fs = numerics.misidentification(fs,misid)
    
    return fs

def three_epoch(params, ns, pts, sig1=0.0, sig2=0.0, theta1=1.0, theta2=1.0, misid=0.0, dt=0.005, folded=False)

Expand source code

def three_epoch(params, ns, pts, sig1 = 0.0, sig2 = 0.0, theta1 = 1.0, theta2 = 1.0, misid = 0.0, dt = 0.005, folded = False):
    """
    Three epoch demography - two instantaneous population size changes in the past
    params = [nu1,nu2,T1,T2]
    nu1,nu2 - relative poplulation size changes to ancestral population size (nu1 occurs before nu2, historically)
    T1,T2 - time for which population had relative sizes nu1, nu2 (scaled by 2N generations)
    sig1, sig2 - population scaled selection coefficients for the two derived alleles
    theta1, theta2 - population scaled mutation rates
    misid - ancestral misidentification parameter
    dt - time step to use for integration
    """
    nu1,nu2,T1,T2 = params
    x = np.linspace(0,1,pts+1)
    sig1,sig2 = np.float(sig1),np.float(sig2)
    
    y1 = dadi.PhiManip.phi_1D(x,gamma=sig1)
    y2 = dadi.PhiManip.phi_1D(x,gamma=sig2)
    phi = np.zeros((len(x),len(x)))
    
    # integrate to equilibrium first
    if sig1 == sig2 == 0.0 and theta1 == theta2 == 1:
        phi = integration.equilibrium_neutral_exact(x)
    else:
        phi = integration.equilibrium_neutral_exact(x)
        phi,y1,y2 = integration.advance(phi, x, 2, y1, y2, nu=1., sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)
    
    phi,y1,y2 = integration.advance(phi, x, T1, y1, y2, nu1, sig1, sig2, theta1, theta2, dt=dt)
    phi,y1,y2 = integration.advance(phi, x, T2, y1, y2, nu2, sig1, sig2, theta1, theta2, dt=dt)

    dx = numerics.grid_dx(x)
    
    try:
        ns = int(ns)
    except TypeError:
        ns = ns[0]
    
    fs = numerics.sample(phi, ns, x)
    fs.extrap_t = dt

    if folded == True:
        fs = fs.fold_major()
    
    if misid > 0.0:
        fs = numerics.misidentification(fs,misid)
    
    return fs

def two_epoch(params, ns, pts, sig1=0.0, sig2=0.0, theta1=1.0, theta2=1.0, misid=0.0, dt=0.005, folded=False)

Two epoch demography - a single population size change at some point in the past params = [nu,T,sig1,sig2,theta1,theta2,misid,dt] nu - relative poplulation size change to ancestral population size T - time in past that size change occured (scaled by 2N generations) sig1, sig2 - population scaled selection coefficients for the two derived alleles theta1, theta2 - population scaled mutation rates misid - ancestral misidentification parameter dt - time step to use for integration

Expand source code

def two_epoch(params, ns, pts, sig1 = 0.0, sig2 = 0.0, theta1 = 1.0, theta2 = 1.0, misid = 0.0, dt = 0.005, folded = False):
    """
    Two epoch demography - a single population size change at some point in the past
    params = [nu,T,sig1,sig2,theta1,theta2,misid,dt]
    nu - relative poplulation size change to ancestral population size
    T - time in past that size change occured (scaled by 2N generations)
    sig1, sig2 - population scaled selection coefficients for the two derived alleles
    theta1, theta2 - population scaled mutation rates
    misid - ancestral misidentification parameter
    dt - time step to use for integration
    """
    nu,T = params

    x = np.linspace(0,1,pts+1)
    sig1,sig2 = np.float(sig1),np.float(sig2)

    y1 = dadi.PhiManip.phi_1D(x,gamma=sig1)
    y2 = dadi.PhiManip.phi_1D(x,gamma=sig2)
    phi = np.zeros((len(x),len(x)))

    # integrate to equilibrium first
    if sig1 == sig2 == 0.0 and theta1 == theta2 == 1:
        phi = integration.equilibrium_neutral_exact(x)
    else:
        phi = integration.equilibrium_neutral_exact(x)
        phi,y1,y2 = integration.advance(phi, x, 2, y1, y2, nu=1., sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)

    phi,y1,y2 = integration.advance(phi, x, T, y1, y2, nu=nu, sig1=sig1, sig2=sig2, theta1=theta1, theta2=theta2, dt=dt)

    dx = numerics.grid_dx(x)

    try:
        ns = int(ns)
    except TypeError:
        ns = ns[0]

    fs = numerics.sample(phi, ns, x)
    fs.extrap_t = dt

    if folded == True:
        fs = fs.fold_major()
    
    if misid > 0.0:
        fs = numerics.misidentification(fs,misid)
    
    return fs