Module dadi.DFE.Vourlaki2022

Functions

def Vourlaki_mixture(params, ns, s1, s2, theta, pts)

Inference model from Vourlaki et al.

params = alpha, beta, ppos_wild, gamma_pos, pchange, pchange_pos

Expand source code

def Vourlaki_mixture(params, ns, s1, s2, theta, pts):
    """
    Inference model from Vourlaki et al.

    params = alpha, beta, ppos_wild, gamma_pos, pchange, pchange_pos
    """
    #params = alpha, beta, gamma_pos, ppos_wild, pchange, pchange_pos
    alpha, beta, ppos_wild, gamma_pos, pchange, pchange_pos = params

    # We'll scale by theta at the end of the script, so set theta=1 here.
    # Case in which gamma is negative and equal in the two pops
    m5 = s1.integrate([alpha, beta], None, PDFs.gamma, 1, None)
    # Case in which gamma is negative and indepenent in the two pops
    m6 = s2.integrate([alpha, beta], None, PDFs.biv_ind_gamma, 1, None,
                      exterior_int=True)

    # Cases in which gamma is positive in both pops.
    # For simplicity, using a single gamma value, rather than exponential dist
    #  that was simulated.
    try:
        m2 = s2.spectra[s2.gammas == gamma_pos,
                        s2.gammas == gamma_pos][0]
        # Even if selection coefficient changed to a different positive value,
        #   we just model them as equal.
        m3 = m2
    except IndexError:
        raise IndexError('Failed to find requested gamma_pos={0:.4f} '
                         'in cached spectra. Was it included '
                         'in additional_gammas during '
                         'cache generation?'.format(gamma_pos))

    # Cases in which gamma is positive in one pop and negative in the other.
    weights = PDFs.gamma(-s2.neg_gammas, [alpha, beta])

    # Case in which pop1 is positive and pop2 is negative.
    pos_neg_spectra = np.squeeze(s2.spectra[s2.gammas==gamma_pos, :len(s2.neg_gammas)])
    m4 = np.trapz(weights[:,np.newaxis,np.newaxis]*pos_neg_spectra,
                       s2.neg_gammas, axis=0)

    # Case in which pop2 is positive and pop1 is negative.
    neg_pos_spectra = np.squeeze(s2.spectra[:len(s2.neg_gammas), s2.gammas==gamma_pos])
    m7 = np.trapz(weights[:,np.newaxis,np.newaxis]*neg_pos_spectra,
                       s2.neg_gammas, axis=0)

    # Contributions to m4 and m7 for gammas that aren't covered by our cache
    # Probability toward gamma=0 that is not covered by our cache
    weight_neu, err = scipy.integrate.quad(PDFs.gamma, 0, -s2.neg_gammas[-1], args=[alpha, beta])
    # Probability toward gamma=-inf that is not covered by our cache
    weight_del, err = scipy.integrate.quad(PDFs.gamma, -s2.neg_gammas[0], np.inf, args=[alpha, beta])

    # In both cases, use the most neutral or deleterious spectra simulated
    m4 += pos_neg_spectra[0]*weight_del
    m4 += pos_neg_spectra[-1]*weight_neu
    m7 += neg_pos_spectra[0]*weight_del
    m7 += neg_pos_spectra[-1]*weight_neu

    # Weights for various parts of distribution checked against Figure S0
    fs = m5*(1-ppos_wild)*(1-pchange) +\
        m6*(1-ppos_wild)*pchange*(1-pchange_pos) +\
        m7*(1-ppos_wild)*pchange*pchange_pos +\
        m2*ppos_wild*(1-pchange) +\
        m3*ppos_wild*pchange*pchange_pos +\
        m4*ppos_wild*pchange*(1-pchange_pos)

    return theta*fs

def visual_validation()

Simulations with a simple model and extreme parameter values, to check functionality and build intuition.

Expand source code

def visual_validation():
    """
    Simulations with a simple model and extreme parameter values, to check
    functionality and build intuition.
    """
    demo_params = [2,2,0.2,0] #nu1, nu2, T, m
    ns, pts_l = [10,10], [30,35,40]

    s1 = Cache1D(demo_params, ns, DemogSelModels.split_mig_sel_single_gamma, pts=pts_l, 
                 gamma_pts=30, gamma_bounds=(1e-4, 200), mp=True, additional_gammas=[10])
    s2 = Cache2D(demo_params, ns, DemogSelModels.split_mig_sel, pts=pts_l, 
                 gamma_pts=30, gamma_bounds=(1e-4, 200), mp=True, additional_gammas=[10])

    alpha, beta, gamma_pos = 0.2, 10, 10

    import matplotlib.pyplot as plt
    fig = plt.figure(10, figsize=(10,8))
    fig.clear()

    axw = fig.add_subplot(3,2,5)
    axd = fig.add_subplot(3,2,6)

    # alpha, beta, ppos_wild, gamma_pos, pchange, pchange_pos
    fs = Vourlaki_mixture([alpha, beta, 0, gamma_pos, 0, 0], None, s1, s2, 1.0, None)
    ax = fig.add_subplot(3,2,1); dadi.Plotting.plot_single_2d_sfs(fs, ax=ax)
    ax.set_title('ppos_wild=0.0, pchange=0.0, pchange_pos=0.0')
    ax.text(0.02,0.98,'A', fontsize='x-large', va='top', transform=ax.transAxes)
    fsw, fsd = fs.filter_pops(tokeep=[1]), fs.filter_pops(tokeep=[2])
    axw.semilogy(fsw, '-o', label='A'); axd.semilogy(fsd, '-o', label='A')

    # alpha, beta, ppos_wild, gamma_pos, pchange, pchange_pos
    fs = Vourlaki_mixture([alpha, beta, 1.0, gamma_pos, 0, 0], None, s1, s2, 1.0, None)
    ax = fig.add_subplot(3,2,2); dadi.Plotting.plot_single_2d_sfs(fs, ax=ax)
    ax.set_title('ppos_wild=1.0, pchange=0.0, pchange_pos=0.0')
    ax.text(0.02,0.98,'B', fontsize='x-large', va='top', transform=ax.transAxes)
    fsw, fsd = fs.filter_pops(tokeep=[1]), fs.filter_pops(tokeep=[2])
    axw.semilogy(fsw, '-o', label='B'); axd.semilogy(fsd, '-o', label='B')

    # alpha, beta, ppos_wild, gamma_pos, pchange, pchange_pos
    fs = Vourlaki_mixture([alpha, beta, 0.0, gamma_pos, 1.0, 0], None, s1, s2, 1.0, None)
    ax = fig.add_subplot(3,2,3); dadi.Plotting.plot_single_2d_sfs(fs, ax=ax)
    ax.set_title('ppos_wild=0.0, pchange=1.0, pchange_pos=0.0')
    ax.text(0.02,0.98,'C', fontsize='x-large', va='top', transform=ax.transAxes)
    fsw, fsd = fs.filter_pops(tokeep=[1]), fs.filter_pops(tokeep=[2])
    axw.semilogy(fsw, '-o', label='C'); axd.semilogy(fsd, '-o', label='C')

    # alpha, beta, ppos_wild, gamma_pos, pchange, pchange_pos
    # This makes me think I got m4, m7 backwards
    fs = Vourlaki_mixture([alpha, beta, 0.0, gamma_pos, 1.0, 0.9], None, s1, s2, 1.0, None)
    ax = fig.add_subplot(3,2,4); dadi.Plotting.plot_single_2d_sfs(fs, ax=ax)
    ax.set_title('ppos_wild=0.0, pchange=1.0, pchange_pos=0.9')
    ax.text(0.02,0.98,'D', fontsize='x-large', va='top', transform=ax.transAxes)
    fsw, fsd = fs.filter_pops(tokeep=[1]), fs.filter_pops(tokeep=[2])
    axw.semilogy(fsw, '-o', label='D'); axd.semilogy(fsd, '-o', label='D')

    axw.set_title('wild'); axd.set_title('domesticate')
    axw.legend(); axd.legend()

    fig.tight_layout()
    fig.savefig('validation.pdf')
    plt.show()