Here, I am just going to redo the same 2D model fitting as I did before for the original 2D spectrum, but with the spectrum to which I applied Ludovic's correction.
In [1]:
from ipyparallel import Client
cl = Client()
cl.ids
Out[1]:
In [2]:
%%px --local
# run whole cell on all engines a well as in the local IPython session
import numpy as np
import sys
sys.path.insert(0, '/home/claudius/Downloads/dadi')
import dadi
In [3]:
def run_dadi(p_init): # for the function to be called with map, it needs to have one input variable
"""
p_init: initial parameter values to run optimisation from
"""
if perturb == True:
p_init = dadi.Misc.perturb_params(p_init, fold=fold,
upper_bound=upper_bound, lower_bound=lower_bound)
# note upper_bound and lower_bound variables are expected to be in the namespace of each engine
# run optimisation of paramters
popt = dadi_opt_func(p0=p_init, data=sfs, model_func=func_ex, pts=pts_l, \
lower_bound=lower_bound, upper_bound=upper_bound, \
verbose=verbose, maxiter=maxiter, full_output=full_output, \
fixed_params=fixed_params)
# pickle to file
import dill
name = outname[:] # make copy of file name stub!
for p in p_init:
name += "_%.4f" % (p)
with open(name + ".dill", "w") as fh:
dill.dump((p_init, popt), fh)
return p_init, popt
In [4]:
from glob import glob
import dill
from utility_functions import *
import pandas as pd
# turn on floating point division by default, old behaviour via '//'
from __future__ import division
In [5]:
lbview = cl.load_balanced_view()
In [6]:
from itertools import repeat
In [7]:
%matplotlib inline
import pylab
pylab.rcParams['figure.figsize'] = [12, 10]
pylab.rcParams['font.size'] = 14
In [8]:
%%px --local
# load spectrum modified with Ludovic's correction, p=35
sfs2d = dadi.Spectrum.from_file("EryPar_modified.2dsfs")
In [9]:
dadi.Plotting.plot_single_2d_sfs(sfs2d, vmin=1, cmap=pylab.cm.jet)
Out[9]:
In [10]:
# number of SNP's in the spectrum
sfs2d.S()
Out[10]:
In [11]:
def split_nomig(params, ns, pts):
"""
params = (nu1,nu2,T)
ns = (n1,n2)
Split into two populations of specifed size, no migration.
nu1: Size ratio of population 1 after split (with respect to ancestral population size Na)
nu2: Size ratio of population 2 after split (with respect to ancestral population size Na)
T: Time in the past of split (in units of 2*Na generations)
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,T = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
phi = dadi.Integration.two_pops(phi, xx, T, nu1, nu2, m12=0, m21=0)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [14]:
cl[:].push(dict(split_nomig=split_nomig))
Out[14]:
In [15]:
%%px --local
# create extrapolating version of the model function
func_ex = dadi.Numerics.make_extrap_log_func(split_nomig)
In [16]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 10 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_nomig" # set file name stub for opt. result files
fixed_params = None
In [17]:
%%px --local
# set lower and upper bounds to nu1, nu2 and T
upper_bound = [1e4, 1e4, 5]
lower_bound = [1e-4, 1e-4, 0]
In [ ]:
# perturb parameter neutral values
p0 = [0.5, 0.5, 0.1] # split into equal proportions at >200,000 generations ago
#ar_split_nomig = lbview.map(run_dadi, repeat(p0, 10))
In [28]:
% ll OUT_2D_models/split_nomig_[!p]*dill
In [47]:
ar_split_nomig = []
for filename in glob("OUT_2D_models/split_nomig_[!p]*dill"):
ar_split_nomig.append(dill.load(open(filename)))
In [48]:
get_flag_count(ar_split_nomig, NM=True)
In [93]:
returned = [flatten(out)[:7] for out in ar_split_nomig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'T_0', 'nu1_opt', 'nu2_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Out[93]:
Beautiful convergence.
In [94]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 3:6])
popt
Out[94]:
In [95]:
# get unscaled, best-fit model spectrum
model = func_ex(popt, ns, pts_l)
# get logL of best-fit model
ll_model = dadi.Inference.ll_multinom(model, sfs2d)
ll_model
Out[95]:
In [96]:
# get theta
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
# get total sequence length in spectrum
L = sfs2d.data.sum()
# get theta per site
print "The optimal value of theta per site for the ancestral population is {0:.4f}.".format(theta/L)
In [58]:
mu = 3e-9 # assumed per generation per site mutation rate
print "The total sequence length for the 2D spectrum is {0:,}.".format(int(L))
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ery and par (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
I am assuming that $\nu_1$ refers to ery and $\nu_2$ refers to par. The split_nomig model with its optimal parameter values suggests the following:
In [62]:
print "The ancestral population of ery and par split apart {0:,} generations ago.".format(int(popt[2]*2*N_ref)),
print "Since then until present the ERY population had a size of {0:,} and the PAR population of {1:,}.".format(int(popt[0]*N_ref), int(popt[1]*N_ref))
With the uncorrected spectrum, I inferred a split time of 400,188 generations and effective population sizes for ery and par of 438,776 and 847,045, respectively.
The following table compares the best-fit parameter values from the uncorrected and corrected spectrum:
| parameter | uncorrected | corrected |
|---|---|---|
| $N_a$ | 688,875 | 631,166 |
| $\nu_1$ | 438,776 | 481,209 |
| $\nu_2$ | 847,045 | 963,542 |
| T | 400,188 | 381,147 |
| -logL | 20387 | 15791 |
The values are given translated to absolute units, i. e. $\nu_x$ in individuals and T in generations.
In [63]:
dadi.Plotting.plot_2d_comp_multinom(model, sfs2d, vmin=1)
The best-fit model spectrum predicts far fewer SNP's with frequency [1,1] than the observed (corrected) spectrum.
The model fitting above was done with a 2D SFS corrected with a p of 0.35. I have created two more spectra corrected with different values of p (0.3 and 0.39). They were derived by using a different null model for optimising p (for details see Ludovics_correction.ipynb).
In [9]:
%%px --local
# load spectrum modified with Ludovic's correction
sfs2d_a = dadi.Spectrum.from_file("EryPar_modified_a.2dsfs") # p=0.30
sfs2d_b = dadi.Spectrum.from_file("EryPar_modified_b.2dsfs") # p=0.39
In [ ]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d_a.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d_a # use sfs corrected with p = 0.3
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 10 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_nomig_p0.30" # set file name stub for opt. result files
fixed_params = None
In [ ]:
%%px --local
# set lower and upper bounds to nu1, nu2 and T
upper_bound = [1e4, 1e4, 5]
lower_bound = [1e-4, 1e-4, 0]
In [ ]:
# perturb parameter neutral values
p0 = [0.5, 0.5, 0.1] # split into equal proportions at >200,000 generations ago
#ar_split_nomig = lbview.map(run_dadi, repeat(p0, 10))
In [34]:
% ll OUT_2D_models/split_nomig_p0.30_*dill
In [35]:
ar_split_nomig = []
for filename in glob("OUT_2D_models/split_nomig_p0.30_*dill"):
ar_split_nomig.append(dill.load(open(filename)))
In [36]:
get_flag_count(ar_split_nomig, NM=True)
In [37]:
returned = [flatten(out)[:7] for out in ar_split_nomig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'T_0', 'nu1_opt', 'nu2_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Out[37]:
The same $T$ is inferred as with the spectrum corrected with a $p$ of 0.35.
In [ ]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d_a.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d_b # use sfs corrected with p=0.39
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 10 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_nomig_p0.39" # set file name stub for opt. result files
fixed_params = None
In [ ]:
%%px --local
# set lower and upper bounds to nu1, nu2 and T
upper_bound = [1e4, 1e4, 5]
lower_bound = [1e-4, 1e-4, 0]
In [ ]:
# perturb parameter neutral values
p0 = [0.5, 0.5, 0.1] # split into equal proportions at >200,000 generations ago
#ar_split_nomig = lbview.map(run_dadi, repeat(p0, 10))
In [38]:
ar_split_nomig = []
for filename in glob("OUT_2D_models/split_nomig_p0.39_*dill"):
ar_split_nomig.append(dill.load(open(filename)))
In [39]:
get_flag_count(ar_split_nomig, NM=True)
In [40]:
returned = [flatten(out)[:7] for out in ar_split_nomig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'T_0', 'nu1_opt', 'nu2_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Out[40]:
The inferred $T$ is the same as before. So the degree of correction has no influence on the inferred divergence time. The divergence time was the only parameter that was varied in the null models used for optimising p.
In [18]:
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(dadi.Demographics2D.split_mig)
In [19]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 10 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_mig" # set file name stub for opt. result files
fixed_params = None
In [20]:
%%px --local
# set lower and upper bounds to nu1, nu2, T, m
upper_bound = [1e4, 1e4, 2, 10]
lower_bound = [1e-4, 1e-4, 0, 0]
In [21]:
# nu1, nu2, T, m:
p0 = [0.5, 0.5, 0.1, 0.1]
#ar_split_mig = lbview.map(run_dadi, repeat(p0, 10))
In [22]:
% ll OUT_2D_models/split_mig_[!p]*dill
In [25]:
% ll OUT_2D_models/split_mig_[0-9]*dill
In [26]:
ar_split_mig = []
for filename in glob("OUT_2D_models/split_mig_[0-9]*dill"):
ar_split_mig.append(dill.load(open(filename)))
In [27]:
get_flag_count(ar_split_mig, NM=True)
In [28]:
# get optimisation results
returned = [flatten(out)[:9] for out in ar_split_mig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'T_0', 'm_0', 'nu1_opt', 'nu2_opt', 'T_opt', 'm_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[28]:
In [ ]:
%%px
maxiter = 300
In [ ]:
# nu1, nu2, T, m:
p0 = [1.5, 2.8, 1.7, 0.2]
#ar_split_mig = lbview.map(run_dadi, repeat(p0, 10))
In [74]:
get_flag_count(ar_split_mig, NM=True)
In [75]:
ar_split_mig = []
for filename in glob("OUT_2D_models/split_mig_[!p]*dill"):
ar_split_mig.append(dill.load(open(filename)))
In [29]:
# get "unsuccessfull" optimisations
returned = [flatten(out)[:9] for out in ar_split_mig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'T_0', 'm_0', 'nu1_opt', 'nu2_opt', 'T_opt', 'm_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[29]:
The divergence time $T$ inferred with the uncorrected spectrum was only 0.92 ($2N_a$ generations).
| parameter | uncorrected | corrected |
|---|---|---|
| $\nu_1$ | 0.994185 | 1.584756 |
| $\nu_2$ | 1.766127 | 2.862212 |
| T | 0.922632 | 1.726329 |
| m | 0.250688 | 0.210858 |
| -logL | 18574 | 12809 |
The fact that with the corrected spectrum higher population sizes are inferred makes sense, since the correction shifts counts from higher frequency to lower frequency variants and a proportional increase in low frequency variants is the effect of population increase on the SFS. The almost doubling of the divergence time and reduction in migration is less obvious to me, but note that these parameters are relative to the $N_a$ and if inferred $N_a$ changed proportionally, then parameter values in absolute units may not be that divergent (see below).
In [30]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 4:8])
popt
Out[30]:
In [31]:
# get unscaled, best-fit model spectrum
model = func_ex(popt, ns, pts_l)
# get logL of best-fit model
ll_model = dadi.Inference.ll_multinom(model, sfs2d)
ll_model
Out[31]:
In [32]:
# get theta
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
# get total sequence length in spectrum
L = sfs2d.data.sum()
# get theta per site
print "The optimal value of theta per site for the ancestral population is {0:.4f}.".format(theta/L)
In [33]:
mu = 3e-9 # assumed per generation per site mutation rate
print "The total sequence length for the 2D spectrum is {0:,}.".format(int(L))
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ery and par (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
The uncorrected spectrum inferred an ancestral population size of 468,295.
In [34]:
print "The ancestral population split apart {0:,} generation ago.".format(int(popt[2]*2*N_ref)),
print "Since then, ERY and PAR had a constant population size of {0:,} and {1:,}, respectively.".format(int(popt[0]*N_ref), int(popt[1]*N_ref)),
print "Since the split, a fraction of {0:.2e} of the population size of ERY were made up each generation of new immigrant individuals from PAR".format(popt[3]/2/N_ref/popt[0]),
print "and a fraction of {0:.2e} of the population size of PAR were made up each generation of new immigtant individuals of ERY.".format(popt[3]/2/N_ref/popt[1])
print "Put another way:",
print "Since the split ERY received a constant number of {0:.2f} new immigrant alleles per generation, while PAR received a constant number of {1:.2f} per generation.".format(popt[3]*popt[0], popt[3]*popt[1])
The time of split inferred with the uncorrected spectrum was: 864,128 generation ago. So an almost doubling of the parameters in relative units translates to an increase by about 1/4 in absolute units (generations).
The following table compares the parameter estimates for the split-migration model derived with the uncorrected and the corrected spectrum in their absolute units:
| parameter | uncorrected | corrected |
|---|---|---|
| $N_a$ | 468,295 | 308,029 |
| $N_{ERY}$ | 465,572 | 488,151 |
| $N_{PAR}$ | 827,069 | 881,645 |
| T (gen.) | 864,128 | 1,063,520 |
| $p_{par->ery}$ | 2.69e-07 | 2.16e-07 |
| $p_{ery->par}$ | 1.52e-07 | 1.20e-07 |
| -logL | 18574 | 12809 |
$N_x$ have unit individuals, T has unit generations, $p_x$ are proportions of new immigrant alleles per generation.
There is a marked difference in the inferred ancestral population size ($N_a$ eq. to $N_{ref}$), which affects all other parameters. The estimated contemporary population sizes for ERY and PAR are slightly higher for the corrected spectrum, but not much so. The inferred divergence time is 1/4 higher with the corrected spectrum. The inferred migration rates (as proportion of new immigrant individuals per generation) are both slightly smaller with the corrected spectrum, but not dramatically.
In [109]:
dadi.Plotting.plot_2d_comp_multinom(data=sfs2d, model=model, vmin=1)
This model with gene flow has greatly reduced residuals. It still predicts too few SNP's with frequency [1,1], but this residual is much reduced when compared with the divergence in isolation model above.
In [110]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d_a.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d_a # use the sfs corrected with p=0.3
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 10 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_mig_p0.30" # set file name stub for opt. result files
fixed_params = None
In [111]:
%%px --local
# set lower and upper bounds to nu1, nu2, T, m
upper_bound = [1e4, 1e4, 2, 10]
lower_bound = [1e-4, 1e-4, 0, 0]
In [ ]:
# nu1, nu2, T, m:
p0 = [1.5, 2.8, 1.7, 0.2]
#ar_split_mig = lbview.map(run_dadi, repeat(p0, 10))
In [113]:
% ll OUT_2D_models/split_mig_p0.30*dill
In [114]:
ar_split_mig = []
for filename in glob("OUT_2D_models/split_mig_p0.30*dill"):
ar_split_mig.append(dill.load(open(filename)))
In [115]:
# get optimisations
returned = [flatten(out)[:9] for out in ar_split_mig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'T_0', 'm_0', 'nu1_opt', 'nu2_opt', 'T_opt', 'm_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[115]:
The inferred $T$ and $m$ are not very different. The values inferred with the spectrum corrected with a p of 0.35 (see above) are quite similar ($T=1.726$, $m=0.211$).
In [ ]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d_a.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d_b # use the sfs corrected with p=0.39
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 10 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_mig_p0.39" # set file name stub for opt. result files
fixed_params = None
In [ ]:
# nu1, nu2, T, m:
p0 = [1.5, 2.8, 1.7, 0.2]
#ar_split_mig = lbview.map(run_dadi, repeat(p0, 10))
In [116]:
ar_split_mig = []
for filename in glob("OUT_2D_models/split_mig_p0.39*dill"):
ar_split_mig.append(dill.load(open(filename)))
In [117]:
# get optimisations
returned = [flatten(out)[:9] for out in ar_split_mig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'T_0', 'm_0', 'nu1_opt', 'nu2_opt', 'T_opt', 'm_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[117]:
Again, the inferred $T$ is not very different, although it varies quite a bit. This model fitting requires a few more runs to achieve convergence (and maybe also a finer grid).
Since the degree of correction does not seem to have a strong influence on parameter estimation, I am going to use only the spectrum corrected with $p=0.35$ in the following.
In [17]:
def split_asym_mig(params, ns, pts):
"""
params = (nu1,nu2,T,m1,m2)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration.
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
T: Time in the past of split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,T,m1,m2 = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, T, nu1, nu2, m12=m2, m21=m1)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [22]:
cl[:].push(dict(split_asym_mig=split_asym_mig))
Out[22]:
In [23]:
%%px --local
func = split_asym_mig
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [24]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 300 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_asym_mig" # set file name stub for opt. result files
fixed_params = None
In [122]:
%%px --local
# set lower and upper bounds to nu1, nu2, T, m1, m2
upper_bound = [1e4, 1e4, 4, 10, 10] # note, I have increased the upper bound for T
lower_bound = [1e-4, 1e-4, 0, 0, 0]
In [17]:
# using optimal parameters from simple divergence with migration model (split_mig)
p0 = [1.572018, 2.838641, 1.703611, 0.211194, 0.211194]
In [ ]:
#ar_split_asym_mig = lbview.map(run_dadi, repeat(p0, 10))
In [13]:
% ll OUT_2D_models/split_asym_mig_[0-9]*dill
In [14]:
ar_split_asym_mig = []
for filename in glob("OUT_2D_models/split_asym_mig_[0-9]*dill"):
ar_split_asym_mig.append(dill.load(open(filename)))
In [15]:
get_flag_count(ar_split_asym_mig, NM=True)
In [18]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'T_0', 'm1_0', 'm2_0', 'nu1_opt', 'nu2_opt', 'T_opt', 'm1_opt', 'm2_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Out[18]:
This looks like convergence. Allowing for asymmetric migration rates improves the model fit by 450 logL units and is therefore highly significant.
In [19]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 5:10])
popt
Out[19]:
In [25]:
# calculate best-fit model spectrum
model = func_ex(popt, ns, pts_l)
ll_model = dadi.Inference.ll_multinom(model, sfs2d)
ll_model
Out[25]:
In [26]:
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
L = sfs2d.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.5f}.".format(theta/L)
In [27]:
mu = 3e-9
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ery and par (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
In [32]:
print "The ancestral population split apart {0:,} generations ago.".format(int(popt[2]*2*N_ref)),
print "Immediately after the split the ERY population changed to a size of {0:,} and the PAR population to {1:,}.".format(int(popt[0]*N_ref), int(popt[1]*N_ref)),
print "Since the split of the ancestral population, PAR received 1 individual from ERY every {0:.2f} generations,".format(1.0/(popt[3]/2*popt[1])),
print "while ERY received 1 PAR individual every {0:.2f} generations.".format(1.0/(popt[4]/2*popt[0])),
print "Put another way: The PAR population contained a constant proportion of {0:.2e} of new immigrant alleles each generation".format(popt[3]/2/N_ref/popt[1]),
print "and the ERY population contained a constant proportion of {0:.2e} of new immigrant alleles each generation.".format(popt[4]/2/N_ref/popt[0])
The following table compares the inferred parameters for the asymmetric migration model from the uncorrected and corrected spectrum in their absolute units:
| parameter | uncorrected | corrected |
|---|---|---|
| $N_{a}$ | 410,678 | 255,715 |
| $N_{ery}$ | 378,198 | 397,740 |
| $N_{par}$ | 919,352 | 975,012 |
| T | 1,022,329 | 1,198,692 |
| $p_{ery->par}$ | 4.58e-08 | 3.91e-08 |
| $p_{par->ery}$ | 6.54e-07 | 4.37e-07 |
| -logL | 18104 | 12359 |
$N_x$ have unit individuals, T has unit generations, $p_x$ are proportions of new immigrant alleles per generation.
There is marked difference in the inferred size of the ancestral population between corrected and uncorrected spectra. The inferred time of the split in number of generations for the unmodified spectrum was 1,022,329. Not a dramatic difference from the one inferred with this corrected spectrum, but still 176,000 years difference.
In [136]:
dadi.Plotting.plot_2d_comp_multinom(data=sfs2d, model=model, vmin=1)
Note the asymmetry in the model spectrum for frequency classes [1, x] and [x, 1].
At the end of the last Ice Age, about 10 Ky ago, both ERY and PAR expanded from refugial areas across Iberia and across northern, central and western Europe, respectively. This expansion likely happened by long-distance migration, where only very few individuals colonised a new habitat. This should lead to a series of founder events that successively reduced genetic diversity the further away from the refugium.
In [13]:
def split_asym_mig_2epoch(params, ns, pts):
"""
params = (nu1_1,nu2_1,T1,nu1_2,nu2_2,T2,m1,m2)
ns = (n1,n2)
Split into two populations of specified size, with potentially asymmetric migration.
The split coincides with a stepwise size change in the daughter populations. Then,
have a second stepwise size change at some point in time after the split. This is
enforced to happen at the same time for both populations. Migration is assumed to
be the same during both epochs.
nu1_1: pop size ratio of pop 1 after split (with respect to Na)
nu2_1: pop size ratio of pop 2 after split (with respect to Na)
T1: Time from split to second size change (in units of 2*Na generations)
nu1_2: pop size ratio of pop 1 after second size change (with respect to Na)
nu2_2: pop size ratio of pop 2 after second size change (with respect to Na)
T2: time in past of second size change (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1_1,nu2_1,T1,nu1_2,nu2_2,T2,m1,m2 = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration for time T1
phi = dadi.Integration.two_pops(phi, xx, T1, nu1_1, nu2_1, m12=m2, m21=m1)
# divergence with potentially asymmetric migration and different pop size for time T2
phi = dadi.Integration.two_pops(phi, xx, T2, nu1_2, nu2_2, m12=m2, m21=m1)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [14]:
cl[:].push(dict(split_asym_mig_2epoch=split_asym_mig_2epoch))
Out[14]:
In [15]:
%%px --local
func = split_asym_mig_2epoch
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [13]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_asym_mig_2epoch" # set file name stub for opt. result files
fixed_params = None
In [16]:
%%px --local
# set lower and upper bounds to nu1, nu2, T, m1, m2
upper_bound = [1e4, 1e4, 6, 1e4, 1e4, 6, 10, 10]
lower_bound = [1e-4, 1e-4, 0, 1e-4, 1e-4, 0, 0, 0]
In [10]:
p0 = [1.0, 3.0, 2.0, 1.5, 3.8, 0.343805, 0.076152, 0.347306]
In [23]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [24]:
ar_split_asym_mig_2epoch.done()
Out[24]:
In [25]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [26]:
get_flag_count(ar_split_asym_mig_2epoch, NM=True)
In [28]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True)
Out[28]:
This clearly needs refinement.
In [29]:
%%px --local
# reduce perturbation of starting values
fold = 1 # perturb randomly up to `fold` times 2-fold
# increase maximum number of iterations
maxiter = 300 # run a maximum of 300 iterations
In [30]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 8:16])
popt
Out[30]:
In [ ]:
p0 = popt
In [31]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [32]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [33]:
get_flag_count(ar_split_asym_mig_2epoch, NM=True)
In [34]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True)
Out[34]:
This is interesting and warrants further refinement.
In [35]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 8:16])
popt
Out[35]:
In [17]:
%%px --local
pts_l = [50, 60, 70]
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
In [36]:
dadi.Inference.optimize_log?
In [38]:
p0 = popt
In [39]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [11]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [12]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[12]:
This looks like convergence. Allowing for a second stepwise size change by adding three more parameters ($\nu_{ery_2}$, $\nu_{par_2}$ and $T2$) improves the likelihood by 325 logL units as compared to the asymmetric migration model (best logL -12359).
This recent bottleneck model can be reduced to the simpler asymmetric migration model by setting T2 to 0. So, in order to compare these two models, I just have to ask whether T2 is significantly different from 0.
I am going to fit the recent bottleneck model to the data with T2 fixed at 0. I expect that it will yield similar best-fit parameters as the asymmetric migration model from above.
I can then do a LRT of the full model vs. the nested model where T2 is set to 0.
In [43]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [50, 60, 70]
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
sfs = sfs2d
perturb = True
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 300 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_asym_mig_2epoch_fixedT2_at0" # set file name stub for opt. result files
fixed_params = [None, None, None, None, None, 0.0, None, None] # fix T2
In [42]:
popt = df.sort_values(by='-logL', ascending=True).iloc[0,8:16]
popt
Out[42]:
In [49]:
popt = np.array(popt)
In [39]:
asym_mig_popt = [1.5554055 , 3.81288413, 2.34380465, 0.07615248, 0.34730636]
p_names = ['ery_opt', 'par_opt', 'T_opt', 'm1_opt', 'm2_opt']
for i in range(len(p_names)):
print p_names[i].rjust(10) + str(asym_mig_popt[i])[:5].rjust(10)
The best-fit parameter values from the asymmetric migration model are fairly similar to the best-fit parameters of the recent bottleneck model.
In [50]:
# set p0, T2_0 will be ignored and change of ery_2 or par_2 should not change the likelihood
p0 = popt
In [51]:
#ar_split_asym_mig_2epoch_fixedT2_at0 = lbview.map(run_dadi, repeat(p0, 10))
In [52]:
ar_split_asym_mig_2epoch_fixedT2_at0.done()
Out[52]:
In [53]:
ar_split_asym_mig_2epoch_fixedT2_at0 = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_fixedT2_at0*dill"):
ar_split_asym_mig_2epoch_fixedT2_at0.append(dill.load(open(filename)))
In [54]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch_fixedT2_at0]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[54]:
This seems to be well converged and recovers the best-fit model spectrum from the previous asymmetric migration model fitting very well. I am therefore confident that I can do a simple LRT with weights (0.5, 0.5) for $\chi^2_0$ and $\chi^2_1$.
In [58]:
ll_c = -12034
ll_s = -12359
In [61]:
dadi.Godambe.sum_chi2_ppf?
In [59]:
D = 2*(ll_c - ll_s)
D
Out[59]:
In [60]:
dadi.Godambe.sum_chi2_ppf(D, weights=(0.5, 0.5))
Out[60]:
Allowing for a second epoch with different population sizes after the split significantly improves the fit to the observed data. I very much doubt that there could be enough linkage in the data to compromise this result.
In [13]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 8:16])
popt
Out[13]:
In [14]:
# calculate best-fit model spectrum
model = func_ex(popt, ns, pts_l)
ll_model = dadi.Inference.ll_multinom(model, sfs2d)
ll_model
Out[14]:
In [45]:
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
L = sfs2d.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.5f}.".format(theta/L)
In [47]:
mu = 3e-9
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ERY and PAR (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
In [49]:
ery_1, par_1, T1, ery_2, par_2, T2, m1, m2 = popt
In [58]:
print "The ancestral population split apart {0:,} generations ago.".format(int((T1+T2)*2*N_ref)),
print "Immediately after the split the ERY population changed to a size of {0:,} and the PAR population to {1:,}.".format(int(ery_1*N_ref), int(par_1*N_ref)),
print "Since the split of the ancestral population, PAR received 1 individual from ERY every {0:.2f} generations,".format(1.0/(m1*par_1/2)),
print "while ERY received 1 PAR individual every {0:.2f} generations.".format(1.0/(m2*ery_1/2)),
#print "Put another way: The PAR population contained a constant proportion of {0:.2e} of new immigrant alleles each generation".format(popt[3]/2/N_ref/popt[1]),
#print "and the ERY population contained a constant proportion of {0:.2e} of new immigrant alleles each generation.".format(popt[4]/2/N_ref/popt[0])
print "ERY and PAR underwent a second stepwise simultaneous population size change {0:,} generations in the past.".format(int(T2*2*N_ref)),
print "ERY changed to a size of {0:,} individuals and PAR to a size of {1:,} individuals.".format(int(ery_2*N_ref), int(par_2*N_ref))
The best fit time and strength of the second population size change seem unreasonably recent and severe. Fitting this model to the uncorrected spectrum yields very similar absolute parameter estimates.
Could there be more biologically reasonable parameter combinations with almost equally high likelihood, i. e. could there be a ridge in the parameter space as is common for different combinations of time and strength of population size changes?
In [70]:
6000/2/N_ref
Out[70]:
Let's fix $T2$ at 0.0104, which should correspond to roughly 6,000 generations, i. e. roughly the time after which the expansion of PAR and ERY out of their refugia should have completed, following the major reversal in climate (cooling) of the Younger Dryas (Hewitt1996).
In [75]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [50, 60, 70]
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
sfs = sfs2d
perturb = True
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 300 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_asym_mig_2epoch_fixedT2" # set file name stub for opt. result files
fixed_params = [None, None, None, None, None, 0.0104, None, None] # fix T2
In [76]:
dadi.Inference.optimize_log?
In [77]:
# set p0, T2_0 will be ignored
p0 = popt * [1, 1, 1, 1e2, 1e2, 1, 1, 1]
In [79]:
#ar_split_asym_mig_2epoch_fixedT2 = lbview.map(run_dadi, repeat(p0, 10))
In [80]:
ar_split_asym_mig_2epoch_fixedT2 = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_fixedT2*dill"):
ar_split_asym_mig_2epoch_fixedT2.append(dill.load(open(filename)))
In [81]:
get_flag_count(ar_split_asym_mig_2epoch_fixedT2, NM=False)
In [82]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch_fixedT2]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True)
Out[82]:
This seems to be very well converged, but the log likelihood of these parameters (-12052) is 18 logL units below the logL of the previous model, where I allowed T2 to vary freely (-12034).
Are these parameter combinations significantly worse than the previous ones?
In [84]:
ll_s = -12052
ll_c = -12034
In [87]:
D = 2*(ll_c - ll_s)
# using Chi^2 dist. with 1 d.f., since there is one parameter fewer to vary in the
# simple model (T2 is fixed) and the fixed parameter is not fixed to the boundary of
# the parameter space
p = dadi.Godambe.sum_chi2_ppf(D, weights = (0, 1))
p
Out[87]:
Yes, the 2 orders of magnitude lower T2 (corresponding to less than 100 generations) inferred above allows for a significantly better fit than a T2 fixed at 0.0104 (corresponding to about 6,000 generations). This should hold even though potential linkage between SNP's has not been taken into account here.
In [24]:
dadi.Plotting.plot_2d_comp_multinom(data=sfs2d, model=model, vmin=1)
There is a marked decrease in residuals for low frequency SNP's in PAR that are absent from ERY as compared to the asymmetric migration model. If they are absent from ERY they are likely to be either of recent origin or selected against in ERY. Bottlenecks lead to the extinction of low frequency variants and the fixation of high frequency variants. This can be seen in the SNP category ery:0, par:1, for which the recent bottleneck model predicts far fewer SNP's than the asymmetric migration model.
I would like to compare the best-fit model spectra of the asymmetric migration and the recent bottleneck model.
In [15]:
recent_bottleneck = dadi.Inference.optimally_scaled_sfs(model, sfs2d)
In [28]:
dadi.Inference.optimally_scaled_sfs?
In [31]:
func_ex = dadi.Numerics.make_extrap_log_func(split_asym_mig)
popt_asym_mig = [1.5554055 , 3.81288413, 2.34380465, 0.07615248, 0.34730636]
model_asym_mig = func_ex(popt_asym_mig, ns, pts_l)
model_asym_mig = dadi.Inference.optimally_scaled_sfs(model_asym_mig, sfs2d)
In [29]:
#model_asym_mig[10,10] = 300
In [32]:
dadi.Plotting.plot_2d_comp_multinom(data=recent_bottleneck.fold() , model=model_asym_mig.fold(), \
vmin=1, title=['recent bottleneck', 'asym mig'], pop_ids=['ery', 'par'])
In the residual plot, red cells indidcate that the asymmetric migration model predicts more SNP's, blue cells indicate that the recent bottleneck model predicts more SNP's. As can be seen in the plot of the residuals between the recent bottleneck and asymmetric migration model, the asymmetric migration model predicts far more SNP's in frequency categories [0, 1] and [1, 0] than the recent bottleneck model, which is due to the predicted loss of low frequency variants in the latter model. So, the better fit of the recent bottleneck model is mainly due to the better fit to these two SNP frequency categories. I find strange that the recent bottleneck model actually predicts more SNP's in frequency classes ery:0, par: 2-6... . Shouldn't these SNP's also be reduced in frequency after a bottleneck?
In [20]:
resid_rec_bot = sfs2d - recent_bottleneck.fold()
resid_asym_mig = sfs2d - model_asym_mig.fold()
In [24]:
dadi.Plotting.plot_2d_resid?
In [25]:
# plot absolute (non-normalised) residuals for recent bottleneck model
dadi.Plotting.plot_2d_resid(resid_rec_bot)
Blue cells indicate that the recent bottleneck model predicts too many SNP's, red cells indicate that the recent bottleneck predicts too few SNP's as compared to the data. Note that these residuals are not normalised (by the expected Poisson SD of each cell).
In [26]:
# plot absolute (non-normalised) residuals for asymmetric migration model
dadi.Plotting.plot_2d_resid(resid_asym_mig)
The asymmetric migration model predicts far too many SNP's in frequency class [0, 1] and too little in class [0, 2].
I have already fit the recent bottleneck model to the uncorrected spectrum with essentially identical result as above: an extremely recent and severe bottleneck (see 05_2D_models.ipynb). Ludovic's correction should reduce the effect of false homozygote "calling" due to PCR duplicates and X-chromosomal loci on the SFS. In the following I will fit the recent bottleneck model to the corrected spectra with $p=0.3$ and $p=0.39$. I want to look for a correlation of parameter estimates with degree of correction ($p$).
In [9]:
%%px --local
# load spectrum modified with Ludovic's correction
sfs2d_a = dadi.Spectrum.from_file("EryPar_modified_a.2dsfs") # p=0.30
sfs2d_b = dadi.Spectrum.from_file("EryPar_modified_b.2dsfs") # p=0.39
In [38]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d_a # use SFS corrected with p=0.30
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_asym_mig_2epoch_p0.30" # set file name stub for opt. result files
fixed_params = None
In [39]:
p0 = [1.0, 3.0, 2.0, 1.5, 3.8, 0.343805, 0.076152, 0.347306]
In [40]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [41]:
ar_split_asym_mig_2epoch.done()
Out[41]:
In [42]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.30*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [43]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[43]:
In [44]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
popt
Out[44]:
In [45]:
p0 = popt
In [46]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [47]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.30*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [48]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[48]:
In [49]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [50]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [51]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.30*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [52]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[52]:
In [53]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [54]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [55]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.30*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [56]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[56]:
In [58]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [57]:
%%px
fold = 1
In [59]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [60]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.30*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [61]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[61]:
In [62]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [63]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [64]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.30*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [65]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[65]:
In [66]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [67]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [68]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.30*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [69]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[69]:
In [70]:
%%px --local
pts_l = [50, 60, 70]
fold = 1
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
In [71]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [72]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [73]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.30*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [74]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[74]:
This seems to have converged.
In [75]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
popt
Out[75]:
In [77]:
# calculate best-fit model spectrum
model = func_ex(popt, ns, pts_l)
ll_model = dadi.Inference.ll_multinom(model, sfs2d_a)
ll_model
Out[77]:
In [78]:
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d_a)
L = sfs2d_a.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.5f}.".format(theta/L)
In [79]:
mu = 3e-9
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ERY and PAR (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
In [80]:
ery_1, par_1, T1, ery_2, par_2, T2, m1, m2 = popt
In [81]:
print "The ancestral population split apart {0:,} generations ago.".format(int((T1+T2)*2*N_ref)),
print "Immediately after the split the ERY population changed to a size of {0:,} and the PAR population to {1:,}.".format(int(ery_1*N_ref), int(par_1*N_ref)),
print "Since the split of the ancestral population, PAR received 1 individual from ERY every {0:.2f} generations,".format(1.0/(m1*par_1/2)),
print "while ERY received 1 PAR individual every {0:.2f} generations.".format(1.0/(m2*ery_1/2)),
#print "Put another way: The PAR population contained a constant proportion of {0:.2e} of new immigrant alleles each generation".format(popt[3]/2/N_ref/popt[1]),
#print "and the ERY population contained a constant proportion of {0:.2e} of new immigrant alleles each generation.".format(popt[4]/2/N_ref/popt[0])
print "ERY and PAR underwent a second stepwise simultaneous population size change {0:,} generations in the past.".format(int(T2*2*N_ref)),
print "ERY changed to a size of {0:,} individuals and PAR to a size of {1:,} individuals.".format(int(ery_2*N_ref), int(par_2*N_ref))
There is not significant change in the inferred timing and size of the recent bottleneck as compared to the fit with the model corrected with p=0.35.
In [82]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d_b # use SFS corrected with p=0.39
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_asym_mig_2epoch_p0.39" # set file name stub for opt. result files
fixed_params = None
In [83]:
p0 = [1.0, 3.0, 2.0, 1.5, 3.8, 0.343805, 0.076152, 0.347306]
In [84]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [85]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.39*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [86]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[86]:
In [87]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [88]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [89]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.39*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [90]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[90]:
In [91]:
%%px --local
pts_l = [50, 60, 70]
fold = 1
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
In [92]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [93]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [94]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.39*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [95]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[95]:
In [96]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [97]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10))
In [98]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_p0.39*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [99]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[99]:
This looks converged.
In [100]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
popt
Out[100]:
In [101]:
# calculate best-fit model spectrum
model = func_ex(popt, ns, pts_l)
ll_model = dadi.Inference.ll_multinom(model, sfs2d_b)
ll_model
Out[101]:
In [102]:
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d_b)
L = sfs2d_b.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.5f}.".format(theta/L)
In [103]:
mu = 3e-9
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ERY and PAR (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
In [104]:
ery_1, par_1, T1, ery_2, par_2, T2, m1, m2 = popt
In [105]:
print "The ancestral population split apart {0:,} generations ago.".format(int((T1+T2)*2*N_ref)),
print "Immediately after the split the ERY population changed to a size of {0:,} and the PAR population to {1:,}.".format(int(ery_1*N_ref), int(par_1*N_ref)),
print "Since the split of the ancestral population, PAR received 1 individual from ERY every {0:.2f} generations,".format(1.0/(m1*par_1/2)),
print "while ERY received 1 PAR individual every {0:.2f} generations.".format(1.0/(m2*ery_1/2)),
#print "Put another way: The PAR population contained a constant proportion of {0:.2e} of new immigrant alleles each generation".format(popt[3]/2/N_ref/popt[1]),
#print "and the ERY population contained a constant proportion of {0:.2e} of new immigrant alleles each generation.".format(popt[4]/2/N_ref/popt[0])
print "ERY and PAR underwent a second stepwise simultaneous population size change {0:,} generations in the past.".format(int(T2*2*N_ref)),
print "ERY changed to a size of {0:,} individuals and PAR to a size of {1:,} individuals.".format(int(ery_2*N_ref), int(par_2*N_ref))
This corrected spectrum infers the most recent and strongest bottleneck of all versions of the observed SFS's.
The degree of correction has no consistent effect on the parameter estimates of the recent bottleneck model.
| p | T2 | ery_2 | par_2 |
|---|---|---|---|
| 0.0 | 6 | 142 | 74 |
| 0.30 | 10 | 649 | 248 |
| 0.35 | 11 | 892 | 309 |
| 0.39 | 1 | 97 | 29 |
In [107]:
# get the difference between the two corrected spectra
# sfs2d_a: p=0.30
# sfs2d_b: p=0.39
resid = sfs2d_b - sfs2d_a
In [108]:
dadi.Plotting.plot_2d_resid(resid)
Red cells indicate p=0.39 has higher count, blue indicates p=0.30 has higher count. This shows that the two corrected spectra mainly differ in the cells [0, 1] and [0, 2]. I had expected that a greater degree of correction would lead to estimation of a less severe bottleneck as that SFS has a higher count in [0, 1] which I had expected to be reduced by a bottleneck.
In [1]:
from ipyparallel import Client
cl = Client()
cl.ids
Out[1]:
In [14]:
lbview = cl.load_balanced_view()
In [3]:
from glob import glob
import dill
from utility_functions import *
import pandas as pd
# turn on floating point division by default, old behaviour via '//'
from __future__ import division
from itertools import repeat
In [4]:
%matplotlib inline
import pylab
pylab.rcParams['figure.figsize'] = [12, 10]
pylab.rcParams['font.size'] = 14
In [17]:
%%px --local
# run whole cell on all engines a well as in the local IPython session
import numpy as np
import sys
sys.path.insert(0, '/usr/local/lib/python2.7/dist-packages/moments-1.0.0-py2.7.egg')
import moments
In [6]:
%%px --local
# load spectrum modified with Ludovic's correction, p=35
sfs2d = moments.Spectrum.from_file("EryPar_modified.2dsfs")
In [7]:
def split_asym_mig_2epoch(params, ns):
"""
params = (nu1_1,nu2_1,T1,nu1_2,nu2_2,T2,m1,m2)
ns = (n1,n2)
Split into two populations of specified size, with potentially asymmetric migration.
The split coincides with a stepwise size change in the daughter populations. Then,
have a second stepwise size change at some point in time after the split. This is
enforced to happen at the same time for both populations. Migration is assumed to
be the same during both epochs.
nu1: pop size ratio of pop 1 after split (with respect to Na)
nu2: pop size ratio of pop 2 after split (with respect to Na)
T1: Time from split to second size change (in units of 2*Na generations)
nu1_2: pop size ratio of pop 1 after second size change (with respect to Na)
nu2_2: pop size ratio of pop 2 after second size change (with respect to Na)
T2: time in past of second size change (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
n1,n2: Sample sizes of resulting Spectrum
"""
nu1,nu2,T1,nu1_2,nu2_2,T2,m1,m2 = params
# fs for equilibrium ancestral population
sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1])
fs = moments.Spectrum(sts)
# split
fs = moments.Manips.split_1D_to_2D(fs, ns[0], ns[1])
# divergence with potentially asymmetric migration for time T1
fs.integrate([nu1, nu2], T1, m=np.array([[0, m2], [m1, 0]]))
# divergence with potentially asymmetric migration and different pop size for time T2
fs.integrate([nu1_2, nu2_2], T2, m=np.array([[0, m2], [m1, 0]]))
return fs
In [8]:
cl[:].push(dict(split_asym_mig_2epoch=split_asym_mig_2epoch))
Out[8]:
In [9]:
def run_moments(p_init): # for the function to be called with map, it needs to have one input variable
"""
p_init: initial parameter values to run optimisation from
"""
if perturb == True:
p_init = moments.Misc.perturb_params(p_init, fold=fold,
upper_bound=upper_bound, lower_bound=lower_bound)
# note upper_bound and lower_bound variables are expected to be in the namespace of each engine
# run optimisation of paramters
popt = moments_opt_func(p0=p_init, data=sfs, model_func=func, \
lower_bound=lower_bound, upper_bound=upper_bound, \
verbose=verbose, maxiter=maxiter, full_output=full_output, \
fixed_params=fixed_params)
# pickle to file
import dill
name = outname[:] # make copy of file name stub!
for p in p_init:
name += "_%.4f" % (p)
with open(name + ".dill", "w") as fh:
dill.dump((p_init, popt), fh)
return p_init, popt
In [10]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
func = split_asym_mig_2epoch
moments_opt_func = moments.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/moments/split_asym_mig_2epoch" # set file name stub for opt. result files
fixed_params = None
In [11]:
%%px --local
# set lower and upper bounds to nu1, nu2, T, m1, m2
upper_bound = [1e4, 1e4, 6, 1e4, 1e4, 6, 10, 10]
lower_bound = [1e-4, 1e-4, 0, 1e-4, 1e-4, 0, 0, 0]
In [12]:
p0 = [1.0, 3.0, 2.0, 1.5, 3.8, 0.343805, 0.076152, 0.347306]
In [18]:
#ar_split_asym_mig_2epoch = lbview.map(run_moments, repeat(p0, 10), block=False)
In [19]:
ar_split_asym_mig_2epoch.done()
Out[19]:
In [20]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/moments/split_asym_mig_2epoch*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [21]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[21]:
In [22]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
popt
Out[22]:
In [23]:
p0 = popt
In [24]:
#ar_split_asym_mig_2epoch = lbview.map(run_moments, repeat(p0, 10), block=False)
In [25]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/moments/split_asym_mig_2epoch*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [26]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[26]:
In [27]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [30]:
%%px
fold = 1
moments_opt_func = moments.Inference.optimize_log # uses BFGS algorithm
In [29]:
moments.Inference.optimize_log?
In [31]:
#ar_split_asym_mig_2epoch = lbview.map(run_moments, repeat(p0, 10), block=False)
In [32]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/moments/split_asym_mig_2epoch*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [33]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[33]:
In [34]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
p0 = popt
In [35]:
#ar_split_asym_mig_2epoch = lbview.map(run_moments, repeat(p0, 10), block=False)
In [36]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/moments/split_asym_mig_2epoch*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [37]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[37]:
This seems to have converged. The best-fit parameter combination inferred with dadi had a 0logL of 12,034. So the moments model spectrum fits slightly better.
In [38]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 8:16])
popt
Out[38]:
In [39]:
# calculate best-fit model spectrum
model = split_asym_mig_2epoch(popt, ns)
ll_model = moments.Inference.ll_multinom(model, sfs2d)
ll_model
Out[39]:
In [40]:
theta = moments.Inference.optimal_sfs_scaling(model, sfs2d)
L = sfs2d.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.5f}.".format(theta/L)
In [41]:
mu = 3e-9
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ERY and PAR (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
In [42]:
ery_1, par_1, T1, ery_2, par_2, T2, m1, m2 = popt
In [43]:
print "The ancestral population split apart {0:,} generations ago.".format(int((T1+T2)*2*N_ref)),
print "Immediately after the split the ERY population changed to a size of {0:,} and the PAR population to {1:,}.".format(int(ery_1*N_ref), int(par_1*N_ref)),
print "Since the split of the ancestral population, PAR received 1 individual from ERY every {0:.2f} generations,".format(1.0/(m1*par_1/2)),
print "while ERY received 1 PAR individual every {0:.2f} generations.".format(1.0/(m2*ery_1/2)),
#print "Put another way: The PAR population contained a constant proportion of {0:.2e} of new immigrant alleles each generation".format(popt[3]/2/N_ref/popt[1]),
#print "and the ERY population contained a constant proportion of {0:.2e} of new immigrant alleles each generation.".format(popt[4]/2/N_ref/popt[0])
print "ERY and PAR underwent a second stepwise simultaneous population size change {0:,} generations in the past.".format(int(T2*2*N_ref)),
print "ERY changed to a size of {0:,} individuals and PAR to a size of {1:,} individuals.".format(int(ery_2*N_ref), int(par_2*N_ref))
Note the remarkable difference in the inferred timing of the split between ERY and PAR between moments and dadi. dadi had inferred a time of split of 1,084,475 generations. So moments infers a split time twice as old! The bottleneck inferred with moments is even slightly more extreme in time and size than the one inferred with dadi.
In [46]:
moments.Plotting.plot_2d_comp_multinom(model, sfs2d, vmin=1)
In [48]:
def split_asym_mig_2epoch(params, ns, pts):
"""
params = (nu1_1,nu2_1,T1,nu1_2,nu2_2,T2,m1,m2)
ns = (n1,n2)
Split into two populations of specified size, with potentially asymmetric migration.
The split coincides with a stepwise size change in the daughter populations. Then,
have a second stepwise size change at some point in time after the split. This is
enforced to happen at the same time for both populations. Migration is assumed to
be the same during both epochs.
nu1_1: pop size ratio of pop 1 after split (with respect to Na)
nu2_1: pop size ratio of pop 2 after split (with respect to Na)
T1: Time from split to second size change (in units of 2*Na generations)
nu1_2: pop size ratio of pop 1 after second size change (with respect to Na)
nu2_2: pop size ratio of pop 2 after second size change (with respect to Na)
T2: time in past of second size change (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1_1,nu2_1,T1,nu1_2,nu2_2,T2,m1,m2 = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration for time T1
phi = dadi.Integration.two_pops(phi, xx, T1, nu1_1, nu2_1, m12=m2, m21=m1)
# divergence with potentially asymmetric migration and different pop size for time T2
phi = dadi.Integration.two_pops(phi, xx, T2, nu1_2, nu2_2, m12=m2, m21=m1)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [52]:
func_ex = dadi.Numerics.make_extrap_log_func(split_asym_mig_2epoch)
popt_dadi = [1.38930513e+00, 4.23522440e+00, 1.88972686e+00, 3.10953710e-03, 1.07862318e-03, 2.01022845e-05, 7.75847834e-02, 4.16367492e-01]
pts_l = [50, 60, 70]
model_dadi = func_ex(popt_dadi, ns, pts_l)
model_dadi = dadi.Inference.optimally_scaled_sfs(model_dadi, sfs2d)
model_moments = dadi.Inference.optimally_scaled_sfs(model, sfs2d)
In [53]:
dadi.Plotting.plot_2d_comp_multinom(data=model_moments.fold() , model=model_dadi.fold(), \
vmin=1, title=['moments', 'dadi'], pop_ids=['ery', 'par'])
The model spectra from dadi and moments are very similar. I don't understand how the inferred times of split could be such divergent between the two programmes.
Due to extremely low effective coverage, i. e. when ignoring PCR duplicates, it is very likely that rare alleles are underrepresented in the 2D SFS. Although Liu2015 (supplementary notes) caution against the usage and interpretation of SFSs with certain frequency classes excluded, in the following I am going to fit the recent bottleneck model to a 2D spectrum with the singleton classes excluded. This is expected to greatly reduce power to detect recent demographic history. I would like to see whether a recent bottleneck can still be inferred, which would greatly corroborate its existence, or not, which would indicate that the inference of the bottleneck is totally dependent on SNP frequency classes that are probably underrepresented due to failed detection of singletons.
In [11]:
sfs2d.mask[1, :] = True
sfs2d.mask[:, 1] = True
In [12]:
dadi.Plotting.plot_single_2d_sfs(sfs2d, vmin=1, cmap=pylab.cm.jet)
Out[12]:
In [16]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_asym_mig_2epoch_maskSingletons" # set file name stub for opt. result files
fixed_params = None
In [19]:
%%px --local
# set lower and upper bounds to nu1, nu2, T, m1, m2
upper_bound = [1e4, 1e4, 6, 1e4, 1e4, 6, 10, 10]
lower_bound = [1e-4, 1e-4, 0, 1e-4, 1e-4, 0, 0, 0]
In [20]:
p0 = [1.0, 3.0, 2.0, 1.5, 3.8, 0.343805, 0.076152, 0.347306]
In [21]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [22]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_maskSingletons*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [23]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True)
Out[23]:
In [24]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 8:16])
p0 = popt
In [25]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [26]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_maskSingletons*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [28]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[28]:
In [30]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 8:16])
p0 = popt
In [31]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [34]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_maskSingletons*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [35]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[35]:
In [36]:
%%px
fold = 1
In [37]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 8:16])
p0 = popt
In [38]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [39]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_maskSingletons*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [40]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[40]:
In [41]:
%%px --local
pts_l = [50, 60, 70]
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
In [42]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 8:16])
p0 = popt
In [43]:
#ar_split_asym_mig_2epoch = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [44]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_maskSingletons*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [45]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[45]:
It is apparently difficult to achieve convergence. However, the best 10 parameter combinations all do not indicate a recent and strong bottleneck. The inference of a recent and strong bottleneck therefore depends on information in the singleton size class of the 2D spectrum.
In [9]:
def asynchronous_bottleneck(params, ns, pts):
"""
params = (nu1_1,nu2_1,T1,nu1_2,nu2_2,T2,m1,m2)
ns = (n1,n2)
Split into two populations of specified size, with potentially asymmetric migration.
The split coincides with the first stepwise size change in the daughter populations.
Then, ERY has a second stepwise size change at time [T2+T3] after the split followed
by PAR having the 2nd size change at time T3 in the past. Migration is assumed to
be the same during all epochs.
nu1: pop size ratio of pop 1 after split (with respect to Na)
nu2: pop size ratio of pop 2 after split (with respect to Na)
T1: Time from split to 2nd size change of ERY (in units of 2*Na generations)
nu1_2: pop size ratio of ERY after second size change (with respect to Na)
nu2_2: pop size ratio of PAR after second size change (with respect to Na)
T2: time in past of 2nd size change of ERY (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
T3: time in past of 2nd size change of PAR
The split happened T1+T2+T3 2*Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,T1,nu1_2,T2,nu2_2,T3,m1,m2 = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration for time T1
phi = dadi.Integration.two_pops(phi, xx, T1, nu1, nu2, m12=m2, m21=m1)
# divergence with potentially different pop size of ERY for time T2 (+T3)
phi = dadi.Integration.two_pops(phi, xx, T2, nu1_2, nu2, m12=m2, m21=m1)
# divergence with potentially different pop size of PAR for time T3
phi = dadi.Integration.two_pops(phi, xx, T3, nu1_2, nu2_2, m12=m2, m21=m1)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [10]:
cl[:].push(dict(asynchronous_bottleneck=asynchronous_bottleneck))
Out[10]:
In [11]:
%%px --local
func = asynchronous_bottleneck
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [12]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/asynchronous_bottleneck" # set file name stub for opt. result files
fixed_params = None
In [13]:
%%px --local
# set lower and upper bounds to nu1,nu2,T1,nu1_2,T2,nu2_2,T3,m1,m2
upper_bound = [1e4, 1e4, 6, 1e4, 6, 1e4, 6, 10, 10]
lower_bound = [1e-4, 1e-4, 0, 1e-4, 0, 1e-4, 0, 0, 0]
Get optimal parameter values from recent bottleneck model.
In [139]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [140]:
l = 2*8+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
popt = df.sort_values(by='-logL', ascending=True).head(10).iloc[0,8:16]
popt
Out[140]:
In [17]:
p0 = [1.389, 4.235, 1.889, 0.03, 1e-3, 0.01, 1e-3, 0.077, 0.416]
In [142]:
#ar_asynchronous_bottleneck = lbview.map(run_dadi, repeat(p0, 10))
In [143]:
ar_asynchronous_bottleneck = []
for filename in glob("OUT_2D_models/asynchronous_bottleneck*dill"):
ar_asynchronous_bottleneck.append(dill.load(open(filename)))
In [145]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_asynchronous_bottleneck]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','T2_0','par_2_0','T3_0','m1_0','m2_0','ery_1_opt','par_1_opt','T1_opt','ery_2_opt','T2_opt','par_2_opt','T3_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[145]:
In [146]:
popt = df.sort_values(by='-logL', ascending=True).head(10).iloc[0,9:18]
popt
Out[146]:
In [147]:
p0 = np.array(popt)
In [148]:
#ar_asynchronous_bottleneck = lbview.map(run_dadi, repeat(p0, 10))
In [149]:
ar_asynchronous_bottleneck = []
for filename in glob("OUT_2D_models/asynchronous_bottleneck*dill"):
ar_asynchronous_bottleneck.append(dill.load(open(filename)))
In [150]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_asynchronous_bottleneck]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','T2_0','par_2_0','T3_0','m1_0','m2_0','ery_1_opt','par_1_opt','T1_opt','ery_2_opt','T2_opt','par_2_opt','T3_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[150]:
In [151]:
p0 = np.array(df.sort_values(by='-logL', ascending=True).head(10).iloc[0,9:18])
In [152]:
#ar_asynchronous_bottleneck = lbview.map(run_dadi, repeat(p0, 10))
In [153]:
ar_asynchronous_bottleneck = []
for filename in glob("OUT_2D_models/asynchronous_bottleneck*dill"):
ar_asynchronous_bottleneck.append(dill.load(open(filename)))
In [154]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_asynchronous_bottleneck]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','T2_0','par_2_0','T3_0','m1_0','m2_0','ery_1_opt','par_1_opt','T1_opt','ery_2_opt','T2_opt','par_2_opt','T3_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[154]:
In [155]:
p0 = np.array(df.sort_values(by='-logL', ascending=True).head(10).iloc[0,9:18])
In [14]:
%%px --local
pts_l = [50, 60, 70]
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
In [157]:
#ar_asynchronous_bottleneck = lbview.map(run_dadi, repeat(p0, 10))
In [15]:
ar_asynchronous_bottleneck = []
for filename in glob("OUT_2D_models/asynchronous_bottleneck*dill"):
ar_asynchronous_bottleneck.append(dill.load(open(filename)))
In [18]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_asynchronous_bottleneck]
df = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','T2_0','par_2_0','T3_0','m1_0','m2_0','ery_1_opt','par_1_opt','T1_opt','ery_2_opt','T2_opt','par_2_opt','T3_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[18]:
This looks reasonably well converged. The inferred model is qualitatively very similar to the recent bottleneck model. The estimate for the time between the 2nd size change of ERY and the 2nd size change of PAR (T2) is very small and both bittlenecks are inferred to be very recent and very severe.
In [19]:
popt = df.sort_values(by='-logL', ascending=True).iloc[0,9:18]
popt
Out[19]:
| parameter | recent bottleneck | async bottleneck |
|---|---|---|
| ery_1_opt | 1.389305 | 1.384623 |
| par_1_opt | 4.235224 | 4.220191 |
| T1_opt | 1.889727 | 1.854568 |
| ery_2_opt | 0.003110 | 0.005099 |
| par_2_opt | 0.001079 | 0.000129 |
| T2_opt | 0.000020 | 0.000035 |
| T3_opt | 0.000002 | |
| m1_opt | 0.077585 | 0.077798 |
| m2_opt | 0.416367 | 0.418371 |
| -logL | 12,034 | 12,033 |
Allowing for different times of the 2nd population size change for ERY and PAR does not improve the fit to the data significantly. It also does not lead to more sensible parameter values for the inferred time and new population size.
So far, the split coincided with a stepwise population size change.
In [11]:
def IM(params, ns, pts):
"""
ns = (n1,n2)
params = (s,nu1,nu2,T,m1,m2)
Isolation-with-migration model with exponential pop growth or decline.
s: Size of pop 1 after split. (Pop 2 has size 1-s.)
nu1: Final population size ratio of pop 1 (with respect to Na)
nu2: Final population size ratio of pop 2 (with respect to Na)
T: Time in the past of split (in units of 2*Na generations)
m1: Migration from pop 1 to pop 2
m2: Migration from pop 2 to pop 1 (2*Na*m12)
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
s,nu1,nu2,T,m1,m2 = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
nu1_func = lambda t: s * (nu1/s)**(t/T)
nu2_func = lambda t: (1-s) * (nu2/(1-s))**(t/T)
phi = dadi.Integration.two_pops(phi, xx, T, nu1_func, nu2_func, m12=m2, m21=m1)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [12]:
cl[:].push(dict(IM=IM))
Out[12]:
In [13]:
%%px --local
func = IM
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [14]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 300 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/IM" # set file name stub for opt. result files
fixed_params = None # do not fix any parameters, optimise all
In [142]:
%%px --local
# set lower and upper bounds to s, nu1, nu2, T, m1, m2
upper_bound = [1, 1e4, 1e4, 4, 10, 10]
lower_bound = [0, 1e-4, 1e-4, 0, 0, 0]
In [15]:
# using optimal parameter values from `split_asym_mig` model as starting values to perturb
p0 = [0.5, 1.5554055, 3.81288413, 2.34380465, 0.07615248, 0.34730636]
In [ ]:
#ar_IM = lbview.map(run_dadi, repeat(p0, 10))
In [ ]:
# using optimal parameter values from `split_asym_mig` model as starting values to perturb
p0 = [0.5, 15.554055, 38.1288413, 2.34380465, 0.07615248, 0.34730636]
# note I have add 10 to the previous nu1_0 and nu2_0
In [52]:
%%px --local
# set lower and upper bounds to s, nu1, nu2, T, m1, m2
upper_bound = [1, 1e4, 1e4, 6, 10, 10] # note, I have increased the upper bound for T
lower_bound = [0, 1e-4, 1e-4, 0, 0, 0]
In [ ]:
#ar_IM = lbview.map(run_dadi, repeat(p0, 10))
In [37]:
% ll OUT_2D_models/IM_[!l]*dill
In [40]:
ar_IM = []
for filename in glob("OUT_2D_models/IM_[!l]*dill"):
ar_IM.append(dill.load((open(filename))))
In [147]:
get_flag_count(ar_IM, NM=True)
In [41]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_IM]
df = pd.DataFrame(data=returned, columns=['s_0' ,'nu1_0','nu2_0', 'T_0', 'm1_0', 'm2_0', 's_opt' ,'nu1_opt', 'nu2_opt', 'T_opt', 'm1_opt', 'm2_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Out[41]:
None of the parameter combinations has higher likelihood than the best-fit parameter combinations of the previous model (lowest neg. logL 12,359). Also the three best-fit parameter combinations are very close to the upper limit of parameter bound that I set for the time parameter $T$ and the population size ratio for PAR implies a current effective population size of several millions (something I think is unlikely given the dispersal rate and the patchiness of the habitat of Chorthippus parallelus).
In [50]:
%%px
fold = 1
In [51]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 6:12])
popt
Out[51]:
In [53]:
p0 = popt
#ar_IM = lbview.map(run_dadi, repeat(p0, 10))
In [16]:
ar_IM = []
for filename in glob("OUT_2D_models/IM_[!l]*dill"):
ar_IM.append(dill.load((open(filename))))
In [18]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_IM]
df = pd.DataFrame(data=returned, columns=['s_0' ,'nu1_0','nu2_0', 'T_0', 'm1_0', 'm2_0', 's_opt' ,'nu1_opt', 'nu2_opt', 'T_opt', 'm1_opt', 'm2_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[18]:
No better parameter combination could be found. I think it is not worth extending the search.
This model changes the population size growth or decline from exponential to linear.
In [19]:
def IM_l(params, ns, pts):
"""
ns = (n1,n2)
params = (s,nu1,nu2,T,m1,m2)
Isolation-with-migration model with LINEAR pop growth or decline.
s: Size of pop 1 after split. (Pop 2 has size 1-s.)
nu1: Final size of pop 1.
nu2: Final size of pop 2.
T: Time in the past of split (in units of 2*Na generations)
m1: Migration from pop 1 to pop 2
m2: Migration from pop 2 to pop 1 (2*Na*m12)
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
s,nu1,nu2,T,m1,m2 = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# specify linear growth from s to nu
nu1_func = lambda t: s + (nu1-s) * (t/T)
nu2_func = lambda t: (1-s) + (nu2-s) * (t/T)
# divergence for time T with linear growth and asymmetric gene flow
phi = dadi.Integration.two_pops(phi, xx, T, nu1_func, nu2_func,
m12=m2, m21=m1)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [150]:
# test specification of linear growth
T = 2
t = np.linspace(0, 2, 100)
nu = 4
N_t = map(lambda t: 1 + (nu-1) * t/T, t)
pylab.plot(t, N_t)
Out[150]:
In [20]:
cl[:].push(dict(IM_l=IM_l))
Out[20]:
In [21]:
%%px --local
func = IM_l
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [22]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 300 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/IM_l" # set file name stub for opt. result files
fixed_params = None # do not fix any parameters, optimise all
In [23]:
%%px --local
# set lower and upper bounds to s, nu1, nu2, T, m1, m2
upper_bound = [1, 1e4, 1e4, 4, 10, 10]
lower_bound = [0, 1e-4, 1e-4, 0, 0, 0]
In [24]:
# using optimal parameter values from `split_asym_mig` model as starting values to perturb
p0 = [0.5, 1.5554055, 3.81288413, 2.34380465, 0.07615248, 0.34730636]
In [ ]:
#ar_IM_l = lbview.map(run_dadi, repeat(p0, 10))
In [25]:
ar_IM_l = []
for filename in glob("OUT_2D_models/IM_l*dill"):
ar_IM_l.append(dill.load((open(filename))))
In [26]:
get_flag_count(ar_IM_l, NM=True)
In [27]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_IM_l]
df = pd.DataFrame(data=returned, columns=['s_0' ,'nu1_0','nu2_0', 'T_0', 'm1_0', 'm2_0', 's_opt' ,'nu1_opt', 'nu2_opt', 'T_opt', 'm1_opt', 'm2_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Out[27]:
This improves the fit significantly as compared to the same model with exponential growth or decline (best logL 12,888.96). However, the split_asym_mig model from section 3 above still provided a better fit (best logL 12,359.63). This model incorporated a stepwise size change (followed by constant population size) that coincided with the population split. Also, note that the best fit parameter combination here contains an exceedingly small proportion of the ancestral population size ($1-s_{opt}$) allocated to PAR of 0.0004.
In [31]:
%%px --local
# set lower and upper bounds to s, nu1, nu2, T, m1, m2
upper_bound = [1, 1e4, 1e4, 6, 10, 10]
lower_bound = [0, 1e-4, 1e-4, 0, 0, 0]
In [32]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[1, 6:12])
popt
Out[32]:
In [33]:
p0 = popt
#ar_IM_l = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [36]:
ar_IM_l = []
for filename in glob("OUT_2D_models/IM_l*dill"):
ar_IM_l.append(dill.load((open(filename))))
In [37]:
l = 2*len(p0)+1
returned = [flatten(out)[:l] for out in ar_IM_l]
df = pd.DataFrame(data=returned, columns=['s_0' ,'nu1_0','nu2_0', 'T_0', 'm1_0', 'm2_0', 's_opt' ,'nu1_opt', 'nu2_opt', 'T_opt', 'm1_opt', 'm2_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[37]:
There is no good convergence, the inferred $s$ is very high and the inferred T is close to the upper boundary that I set. Furthermore, the best -logL is still higher than the best -logL achieved with the asymmetric migration model (12,359). I don't think it is worth continuing the search.
In [77]:
def split_asym_mig_iso(params, ns, pts):
"""
params = (nu1,nu2,Tc,m1,m2,Ti)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time Tc followed by complete isolation until present.
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
Tc: Time of gene flow after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
Ti: Time of isolation after cessation of gene flow
The split lies Tc+Ti * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,Tc,m1,m2,Ti = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, Tc, nu1, nu2, m12=m2, m21=m1)
# divergence without gene flow
phi = dadi.Integration.two_pops(phi, xx, Ti, nu1, nu2, m12=0, m21=0)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [78]:
cl[:].push(dict(split_asym_mig_iso=split_asym_mig_iso))
Out[78]:
In [15]:
%%px --local
func = split_asym_mig_iso
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [16]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 300 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_asym_mig_iso" # set file name stub for opt. result files
fixed_params = None
In [17]:
%%px --local
# set lower and upper bounds to nu1, nu2, Tc, m1, m2 and Ti
upper_bound = [1e4, 1e4, 6, 10, 10, 6] # note, I have increased the upper bound for T
lower_bound = [1e-4, 1e-4, 0, 0, 0, 0]
In [20]:
# using the optimal parameters from the `split_asym_mig` model
p0 = [1.555405, 3.812884, 2.343805/2, 0.076152, 0.347306, 2.343805/2]
In [168]:
#ar_split_asym_mig_iso = lbview.map(run_dadi, repeat(p0, 10))
In [17]:
ar_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/split_asym_mig_iso*dill"):
ar_split_asym_mig_iso.append(dill.load(open(filename)))
In [18]:
get_flag_count(ar_split_asym_mig_iso, NM=True)
In [19]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Out[19]:
In [20]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,6:12])
popt
Out[20]:
In [22]:
%%px --local
pts_l = [50, 60, 70] # make finer grid
In [23]:
p0 = popt
#ar_split_asym_mig_iso = lbview.map(run_dadi, repeat(p0, 10))
In [24]:
get_flag_count(ar_split_asym_mig_iso, NM=True)
In [25]:
ar_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/split_asym_mig_iso*dill"):
ar_split_asym_mig_iso.append(dill.load(open(filename)))
In [26]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Out[26]:
The best-fit parameter combination has a much higher likelihood than the second best. It infers a very large population size ratio for PAR and a very large $T_c$, i. e. of the time of divergence with gene flow before the beginning of the time of divergence without gene flow. I think I need to further explore this region of the parameter space.
In [27]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,6:12])
popt
Out[27]:
In [28]:
%%px --local
pts_l = [60, 70, 80] # make finer grid
fold = 1
maxiter = 100
In [29]:
p0 = popt
#ar_split_asym_mig_iso = lbview.map(run_dadi, repeat(p0, 10))
In [30]:
get_flag_count(ar_split_asym_mig_iso, NM=True)
In [70]:
ar_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/split_asym_mig_iso*dill"):
ar_split_asym_mig_iso.append(dill.load(open(filename)))
In [73]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[73]:
This looks like decent convergence. Let's try to refine it even more.
In [68]:
%%px --local
# set lower and upper bounds to nu1, nu2, Tc, m1, m2 and Ti
upper_bound = [1e4, 1e4, 8, 10, 10, 6] # note, I have increased the upper bound for T
lower_bound = [1e-4, 1e-4, 0, 0, 0, 0]
In [25]:
%%px --local
pts_l = [60, 70, 80] # make finer grid
fold = 1
maxiter = 300
In [74]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,6:12])
popt
Out[74]:
In [75]:
p0 = popt
In [77]:
#ar_split_asym_mig_iso = lbview.map(run_dadi, repeat(p0, 10))
In [54]:
% ll OUT_2D_models/split_asym_mig_iso_[0-9]*dill
In [55]:
ar_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/split_asym_mig_iso_[0-9]*dill"):
ar_split_asym_mig_iso.append(dill.load(open(filename)))
In [56]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[56]:
This looks like good convergence.
I would like to compare this model with the asymmetric migration model.
In [42]:
ll_s = 12359
ll_c = 12156
D = 2 * (ll_s - ll_c)
D
Out[42]:
With the uncorrected spectrum D was greater than 1000, so the improvement in fit when introducing recent cessation of gene flow was stronger.
In [43]:
# calculate p-value for Chi-square dist.
# the weights specify a weighted sum of chi^2 distributions with 0 and 1 d.o.f
# this is because Ti is 0 in the split_asym_mig model and at the boundary of the parameter space
p = dadi.Godambe.sum_chi2_ppf(D, weights=(0.5, 0.5))
p
Out[43]:
In [65]:
dadi.Godambe.sum_chi2_ppf?
Allowing for a recent period without gene flow significantly improves the fit as compared to the asymmetric migration model. I doubt that there is enough linkage in the data to compromise this result.
In [44]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,6:12])
popt
Out[44]:
In [45]:
model = func_ex(popt, ns, pts_l)
In [46]:
ll_model = dadi.Inference.ll_multinom(model, sfs2d)
ll_model
Out[46]:
In [47]:
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
L = sfs2d.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.5f}.".format(theta/L)
In [48]:
mu = 3e-9
L = sfs2d.data.sum() # this sums over all entries in the spectrum, including masked ones, i. e. also contains invariable sites
print "The total sequence length for the 2D spectrum is {0:,}.".format(int(L))
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ery and par (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
This quite a small effective population size.
In [30]:
nu1, nu2, Tc, m1, m2, Ti = popt
In [49]:
print "The ancestral population split apart {0:,} generations ago.".format(int((Tc+Ti)*2*N_ref)),
print "Immediately after the split the ERY population changed to a size of {0:,} and the PAR population to {1:,}.".format(int(nu1*N_ref), int(nu2*N_ref)),
print "Since the split of the ancestral population, PAR received 1 individual from ERY every {0:.2f} generations,".format(1.0/(m1/2*nu2)),
print "while ERY received 1 PAR individual every {0:.2f} generations.".format(1.0/(m2/2*nu1)),
print "Put another way: The PAR population contained a constant proportion of {0:.2e} of new immigrant alleles each generation".format(m1/2/N_ref/nu2),
print "and the ERY population contained a constant proportion of {0:.2e} of new immigrant alleles each generation.".format(m2/2/N_ref/nu1),
print "ERY and PAR remained in contact for {0:,} generations.".format(int(Tc*2*N_ref)),
print "{0:,} generations ago gene flow between ERY and PAR had ceased.".format(int(Ti*2*N_ref))
The following table compares the inferred parameters for the ancient migration model from the uncorrected and corrected spectrum in their absolute units:
| parameter | uncorrected | corrected |
|---|---|---|
| $N_{a}$ | 125,977 | 113,522 |
| $N_{ery}$ | 363,349 | 389,237 |
| $N_{par}$ | 886,821 | 950,409 |
| $T_c$ | 1,748,505 | 1,525,640 |
| $p_{ery->par}$ | 2.21e-08 | 2.14e-08 |
| $p_{par->ery}$ | 3.27e-07 | 2.60e-07 |
| $T_i$ | 43,628 | 26,970 |
| D | 1044 | 406 |
$N_x$ have unit individuals, $T_x$ has unit generations, $p_x$ are proportions of new immigrant alleles per generation and D is two times the likelihood ratio between the best fit parameters of this model with the hitherto best performing asymmetric migration model.
Note that with the uncorrected spectrum, I could not achieve satisfactory convergence of parameter estimations. The parameters inferred with the corrected spectrum are therefore more reliable. The inferred population sizes for ERY and PAR aren't very different between corrected and uncorrected spectra and so is the inferred time of contact ($T_c$). The migration rates are slightly reduced with the corrected spectrum and the time since cessation of gene flow ($T_i$) is inferred to be more recent than with the uncorrected spectrum.
The addition of $T_i$ improves the model by 203 log likelihood units as compared to the model without a period of complete isolation (split_asym_mig above). If the better fit is not just due to fitting noise and bias in the data, then this means that the final period of complete isolation is significantly greater than 0. I think the fact that both the raw and the corrected spectrum lead to very similar parameter estimates provides some support for the interpretation that $T_i$ is actually greater than zero and that this is not just due to fitting noise or bias in the data.
In [38]:
dadi.Plotting.plot_2d_comp_multinom(model=model, data=sfs2d, vmin=1)
In [45]:
func_ex = dadi.Numerics.make_extrap_log_func(split_asym_mig)
popt_asym_mig = [1.5554055 , 3.81288413, 2.34380465, 0.07615248, 0.34730636]
model_asym_mig = func_ex(popt_asym_mig, ns, pts_l)
model_asym_mig = dadi.Inference.optimally_scaled_sfs(model_asym_mig, sfs2d)
model = dadi.Inference.optimally_scaled_sfs(model, sfs2d)
In [46]:
dadi.Plotting.plot_2d_comp_multinom(data=model.fold() , model=model_asym_mig.fold(), \
vmin=1, title=['ancient mig', 'asym mig'], pop_ids=['ery', 'par'])
I am going to try and allow different rates of gene flow in the two epochs.m
In [79]:
def two_epoch_migration(params, ns, pts):
"""
params = (nu1,nu2,T1,m1,m2,T2,m1_2,m2_2)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time T1 followed by a second epoch T2 of potentially changed rates of gene flow
between the two populations.
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
T1: Time of first epoch after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
T2: Time of second epoch
m1_2: Migration rate from ery into par in 2nd epoch
m2_2: Migration rate from par into ery in 2nd epoch
The split lies T1+T2 * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,T1,m1,m2,T2,m1_2,m2_2 = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, T1, nu1, nu2, m12=m2, m21=m1)
# divergence with potentially different rate of gene flow
phi = dadi.Integration.two_pops(phi, xx, T2, nu1, nu2, m12=m2_2, m21=m1_2)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [82]:
cl[:].push(dict(two_epoch_migration=two_epoch_migration))
Out[82]:
In [83]:
%%px --local
func = two_epoch_migration
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [84]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 300 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/two_epoch_migration" # set file name stub for opt. result files
fixed_params = None
In [85]:
%%px --local
# set lower and upper bounds to nu1,nu2,T1,m1,m2,T2,m1_2,m2_2
upper_bound = [1e4, 1e4, 6, 10, 10, 6, 10, 10]
lower_bound = [1e-4, 1e-4, 0, 0, 0, 0, 0, 0]
Get the best fit parameter values from the ancient migration model.
In [90]:
ar_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/split_asym_mig_iso_[0-9]*dill"):
ar_split_asym_mig_iso.append(dill.load(open(filename)))
In [92]:
l = 2*6+1
success = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df_ami = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
popt_ami = df_ami.sort_values(by='-logL', ascending=True).head(10).iloc[0,6:12]
popt_ami
Out[92]:
For starting values, I am going to set migrations rates equal in both epochs.
In [93]:
popt_ami[-3:-1]
Out[93]:
In [102]:
p0 = np.array(popt_ami)
p0
Out[102]:
In [103]:
# extend array size
p0 = np.resize(p0, 8)
p0
Out[103]:
In [105]:
# set starting values of migration rates in 2nd epoch equal to 1st epoch
p0[-2:] = popt_ami[-3:-1]
p0
Out[105]:
In [106]:
#ar_two_epoch_migration = lbview.map(run_dadi, repeat(p0, 10))
In [113]:
ar_two_epoch_migration = []
for filename in glob("OUT_2D_models/two_epoch_migration_[0-9]*dill"):
ar_two_epoch_migration.append(dill.load(open(filename)))
In [116]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_two_epoch_migration]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'T1_0', 'm1_0', 'm2_0', 'T2_0', 'm1_2_0', 'm2_2_0', 'nu1_opt', 'nu2_opt', 'T1_opt', 'm1_opt', 'm2_opt', 'T2_opt', 'm1_2_opt', 'm2_2_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[116]:
In [118]:
popt = df.sort_values(by='-logL', ascending=True).iloc[0,8:16]
popt
Out[118]:
In [119]:
p0 = np.array(popt)
In [120]:
#ar_two_epoch_migration = lbview.map(run_dadi, repeat(p0, 10))
In [121]:
ar_two_epoch_migration = []
for filename in glob("OUT_2D_models/two_epoch_migration_[0-9]*dill"):
ar_two_epoch_migration.append(dill.load(open(filename)))
In [122]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_two_epoch_migration]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'T1_0', 'm1_0', 'm2_0', 'T2_0', 'm1_2_0', 'm2_2_0', 'nu1_opt', 'nu2_opt', 'T1_opt', 'm1_opt', 'm2_opt', 'T2_opt', 'm1_2_opt', 'm2_2_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[122]:
In [123]:
%%px --local
pts_l = [50, 60, 70] # make finer grid
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
In [124]:
popt = df.sort_values(by='-logL', ascending=True).iloc[0,8:16]
popt
Out[124]:
In [125]:
p0 = np.array(popt)
In [126]:
#ar_two_epoch_migration = lbview.map(run_dadi, repeat(p0, 10))
In [127]:
ar_two_epoch_migration = []
for filename in glob("OUT_2D_models/two_epoch_migration_[0-9]*dill"):
ar_two_epoch_migration.append(dill.load(open(filename)))
In [128]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_two_epoch_migration]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'T1_0', 'm1_0', 'm2_0', 'T2_0', 'm1_2_0', 'm2_2_0', 'nu1_opt', 'nu2_opt', 'T1_opt', 'm1_opt', 'm2_opt', 'T2_opt', 'm1_2_opt', 'm2_2_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[128]:
This looks reasonably converged. There is no good convergence on $m1_2$, the migration rate from ERY into PAR in the 2nd epoch, but all values indicate a drastic reduction. The best-fit parameter combination for the ancient migration model had -logL of 12,156, so only 2 log likelihood units worse.
In [129]:
ll_s = -12156
ll_c = -12154
D = 2*(ll_c - ll_s)
D
Out[129]:
In [130]:
# calculate p-value for Chi-square dist. with 2 d.o.f
# two parameters are at the boundary of the parameter space
# therefore, D should be distributed by a mixture of Chi-square dist.'s
# with different degrees of freedom, but I cannot calculate the mixing probabilities (see Self1987)
# The following test should be conservative.
p = dadi.Godambe.sum_chi2_ppf(D, weights=(0, 0, 1))
p
Out[130]:
I don't think that this model provides a significant improvement over the ancient migration model. All best-fit parameter combinations indicate a drastic reduction in migration from ERY into PAR in the 2nd epoch. The reduction in gene flow from PAR to ERY is inferred to be around 75-80%. The time for the second epoch (T2), is very similar to the time Ti inferred in the ancient migration model. So we cannot rule out that gene flow has continued in the 2nd epoch, especially in the direction from PAR into ERY, but there doesn't seem to be enough power to show this.
I would like to combine the ancient migration and recent bottleneck model. To keep the model simple, I will first run a model that enforces the coincidence of the second size change with the cessation of gene flow. Then I will run a model that allows the cessation of gene flow to occur earlier than the bottleneck (i. e. in the 1st epoch after the split) and a model that allows the cessation of gene flow to occur after the bottleneck (i. e. in the second epoch after the split).
In [36]:
def ancMig_recBotIso(params, ns, pts):
"""
params = (nu1,nu2,Tc,m1,m2,nu1_2,nu2_2,Ti)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time Tc followed by cessation of gene flow and coincident potential
second size change. The split coincides with the first potential pop size change.
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
Tc: Time of gene flow after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
Tbi: Time of isolation after cessation of gene flow and second size change
nu1_2: pop size ratio of pop 1 after second size change
nu2_2: pop size ratio of pop 2 after second size change
The split lies Tc+Tbi * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,Tc,m1,m2,nu1_2,nu2_2,Tbi = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, Tc, nu1, nu2, m12=m2, m21=m1)
# divergence without gene flow and potential second pop size change
phi = dadi.Integration.two_pops(phi, xx, Tbi, nu1_2, nu2_2, m12=0, m21=0)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [38]:
cl[:].push(dict(ancMig_recBotIso=ancMig_recBotIso))
Out[38]:
In [41]:
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(ancMig_recBotIso)
In [40]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/ancMig_recBotIso" # set file name stub for opt. result files
fixed_params = None
In [47]:
%%px --local
# set lower and upper bounds to nu1, nu2, Tc, m1, m2, nu1_2, nu2_2, Tbi
upper_bound = [1e4, 1e4, 6, 10, 10, 1e4, 1e4, 6]
lower_bound = [1e-4, 1e-4, 0, 0, 0, 1e-4, 1e-4, 0]
ancient migration best fit parameter values
nu1,nu2,Tc,m1,m2,Ti [ 3.42873459, 8.372007 , 6.71956526, 0.04062737, 0.20232337, 0.11878969]
recent bottleneck best fit parameter values
nu1,nu2,T1,nu1_2,nu2_2,T2,m1,m2 [ 1.38930513e+00, 4.23522440e+00, 1.88972686e+00, 3.10953710e-03, 1.07862318e-03, 2.01022845e-05, 7.75847834e-02, 4.16367492e-01]
| param | recent bottleneck | ancient mig |
|---|---|---|
| nu1 | 1.389 | 3.428 |
| nu2 | 4.235 | 8.372 |
| Tc | 1.889 | 6.719 |
| m1 | 0.077 | 0.040 |
| m2 | 0.416 | 0.202 |
| nu1_2 | 3.109e-03 | - |
| nu2_2 | 1.078e-03 | - |
| Tb | 2.010e-05 | - |
| Ti | - | 0.118 |
In [44]:
# nu1,nu2,Tc,m1,m2,nu1_2,nu2_2,Tbi
p0 = [1.5, 5, 2, 0.06, 0.3, 3e-1, 1e-1, 0.01]
In [48]:
#ar_ancMig_recBotIso = lbview.map(run_dadi, repeat(p0, 10))
In [49]:
ar_ancMig_recBotIso.done()
Out[49]:
In [54]:
% ll OUT_2D_models/ancMig_recBotIso_[0-9]*dill
In [55]:
ar_ancMig_recBotIso = []
for filename in glob("OUT_2D_models/ancMig_recBotIso_[0-9]*dill"):
ar_ancMig_recBotIso.append(dill.load(open(filename)))
In [57]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_recBotIso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'nu1_2_0', 'nu2_2_0', 'Tbi_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tbi_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[57]:
In [58]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
popt
Out[58]:
In [59]:
p0 = popt
In [60]:
#ar_ancMig_recBotIso = lbview.map(run_dadi, repeat(p0, 10))
In [61]:
ar_ancMig_recBotIso.done()
Out[61]:
In [62]:
ar_ancMig_recBotIso = []
for filename in glob("OUT_2D_models/ancMig_recBotIso_[0-9]*dill"):
ar_ancMig_recBotIso.append(dill.load(open(filename)))
In [63]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_recBotIso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'nu1_2_0', 'nu2_2_0', 'Tbi_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tbi_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[63]:
In [64]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
popt
Out[64]:
In [66]:
p0 = popt
In [67]:
#ar_ancMig_recBotIso = lbview.map(run_dadi, repeat(p0, 10))
In [68]:
ar_ancMig_recBotIso = []
for filename in glob("OUT_2D_models/ancMig_recBotIso_[0-9]*dill"):
ar_ancMig_recBotIso.append(dill.load(open(filename)))
In [69]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_recBotIso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'nu1_2_0', 'nu2_2_0', 'Tbi_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tbi_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[69]:
In [70]:
%%px --local
pts_l = [50, 60, 70] # make finer grid
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
In [71]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
popt
Out[71]:
In [72]:
p0 = popt
In [73]:
#ar_ancMig_recBotIso = lbview.map(run_dadi, repeat(p0, 10))
In [75]:
ar_ancMig_recBotIso = []
for filename in glob("OUT_2D_models/ancMig_recBotIso_[0-9]*dill"):
ar_ancMig_recBotIso.append(dill.load(open(filename)))
In [76]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_recBotIso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'nu1_2_0', 'nu2_2_0', 'Tbi_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tbi_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[76]:
This looks like convergence! The best-fit ancient migration model (6 parameters) had -logL of 12,156 and the best-fit recent bottleneck model (8 parameters) had -logL of 12,034. Note that this model infers a bottleneck in the second epoch only for PAR, not for ERY (compare $\nu_{1_{opt}}$ with $\nu_{1_{2opt}}$ and $\nu_{2_{opt}}$ with $\nu_{2_{2opt}}$).
This ancMig recBotIso model can be turned into the ancient migration model by setting $\nu_{1_2} = \nu_1$ and $\nu_{2_2} = \nu_2$, i. e. no population size change at time $T_{bi}$, just cessation of gene flow. The ancient migration model is therefore nested within this ancMig recBotIso model.
In [77]:
ll_c = -11975
ll_s = -12156
D = 2*(ll_c - ll_s)
D
Out[77]:
In [78]:
# calculate p-value for Chi-square dist. with 2 d.o.f.
# the weights specify a weighted sum of chi^2 distributions with 0, 1 and 2 d.o.f
p = dadi.Godambe.sum_chi2_ppf(D, weights=(0, 0, 1))
p
Out[78]:
This confirms that the ancMig recBotIso model provides a significantly better fit to the data than the ancient migration model despite requiring two more parameters.
The difference between this ancMig recBotIso model and the recent bottleneck model is that the recent bottleneck model enforces continued and constant migration through the second epoch until present while this ancMig recBotIso model enforces complete isolation during the second epoch. Both models have equal number of parameters and one is not nested within the other. A likelihood ratio test can therefore not be used to compare both models. The ancMig recBotIso model provides a better fit to the data by 59 logL units and is therefore clearly the better model. Thus, there seems to be evidence in the data for both a recent "bottleneck" in PAR as well as a recent cessation (or at least reduction) of gene flow.
In [79]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,8:16])
popt
Out[79]:
In [80]:
model = func_ex(popt, ns, pts_l)
In [81]:
ll_model = dadi.Inference.ll_multinom(model, sfs2d)
ll_model
Out[81]:
In [82]:
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
L = sfs2d.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.5f}.".format(theta/L)
In [83]:
mu = 3e-9
L = sfs2d.data.sum() # this sums over all entries in the spectrum, including masked ones, i. e. also contains invariable sites
print "The total sequence length for the 2D spectrum is {0:,}.".format(int(L))
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ery and par (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
In [84]:
nu1,nu2,Tc,m1,m2,nu1_2,nu2_2,Tbi = popt
In [91]:
print "The ancestral population split apart {0:,} generations ago.".format(int((Tc+Tbi)*2*N_ref)),
print "Immediately after the split the ERY population changed to a size of {0:,} and the PAR population to {1:,}.".format(int(nu1*N_ref), int(nu2*N_ref)),
print "Since the split of the ancestral population, PAR received 1 individual from ERY every {0:.2f} generations,".format(1.0/(m1/2*nu2)),
print "while ERY received 1 PAR individual every {0:.2f} generations.".format(1.0/(m2/2*nu1)),
print "Put another way: The PAR population contained a constant proportion of {0:.2e} of new immigrant alleles each generation".format(m1/2/N_ref/nu2),
print "and the ERY population contained a constant proportion of {0:.2e} of new immigrant alleles each generation.".format(m2/2/N_ref/nu1),
print "ERY and PAR remained in contact for {0:,} generations.".format(int(Tc*2*N_ref)),
print "{0:,} generations ago gene flow between ERY and PAR had ceased.".format(int(Tbi*2*N_ref)),
print "At the same time, ERY and PAR underwent a second stepwise simultaneous population size change.",
print "ERY changed to a size of {0:,} individuals and PAR to a size of {1:,} individuals.".format(int(nu1_2*N_ref), int(nu2_2*N_ref))
Note that the bottleneck population size of PAR is still greater than the inferred ancestral population size of ERY and PAR.
In [92]:
dadi.Plotting.plot_2d_comp_multinom(data=sfs2d, model=model, vmin=1)
In [10]:
def ancMig_iso_recBot(params, ns, pts):
"""
params = (nu1,nu2,Tc,m1,m2,Ti,nu1_2,nu2_2,Tb)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time Tc followed by cessation of gene flow at time [Ti+Tb]. The split coincides
with the first potential pop size change. At time Tb there is a potential second size change.
Note, that this model does not allow the cessation of gene flow to occur after the second
pop size change.
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
Tc: Time of gene flow after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
Ti: Time of complete isolation before second size change
nu1_2: pop size ratio of pop 1 after second size change
nu2_2: pop size ratio of pop 2 after second size change
Tb: Time since second pop size change (and continued complete isolation)
The split lies Tc+Ti+Tb * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,Tc,m1,m2,Ti,nu1_2,nu2_2,Tb = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, Tc, nu1, nu2, m12=m2, m21=m1)
# divergence without gene flow
phi = dadi.Integration.two_pops(phi, xx, Ti, nu1, nu2, m12=0, m21=0)
# potential second pop size change
phi = dadi.Integration.two_pops(phi, xx, Tb, nu1_2, nu2_2, m12=0, m21=0)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [11]:
cl[:].push(dict(ancMig_iso_recBot=ancMig_iso_recBot))
Out[11]:
In [12]:
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(ancMig_iso_recBot)
In [13]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/ancMig_iso_recBot" # set file name stub for opt. result files
fixed_params = None
In [14]:
%%px --local
# set lower and upper bounds to nu1,nu2,Tc,m1,m2,Ti,nu1_2,nu2_2,Tb
upper_bound = [1e4, 1e4, 6, 10, 10, 6, 1e4, 1e4, 6]
lower_bound = [1e-4, 1e-4, 0, 0, 0, 0, 1e-4, 1e-4, 0]
In [98]:
popt = df.sort_values(by='-logL', ascending=True).iloc[0,8:16]
popt
Out[98]:
In [18]:
p0 = [1.442, 4.862, 2.376, 0.066266, 0.460058, 1e-2, 1.782, 1.45, 0.039]
In [101]:
#ar_ancMig_iso_recBot = lbview.map(run_dadi, repeat(p0, 10))
In [102]:
ar_ancMig_iso_recBot.done()
Out[102]:
In [103]:
ar_ancMig_iso_recBot = []
for filename in glob("OUT_2D_models/ancMig_iso_recBot_[0-9]*dill"):
ar_ancMig_iso_recBot.append(dill.load(open(filename)))
In [104]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_iso_recBot]
df = pd.DataFrame(data=success, \
columns=['nu1_0', 'nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_2_0', 'nu2_2_0', 'Tb_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tb_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[104]:
In [105]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,9:18])
popt
Out[105]:
In [106]:
p0 = popt
In [107]:
#ar_ancMig_iso_recBot = lbview.map(run_dadi, repeat(p0, 10))
In [108]:
ar_ancMig_iso_recBot.done()
Out[108]:
In [109]:
ar_ancMig_iso_recBot = []
for filename in glob("OUT_2D_models/ancMig_iso_recBot_[0-9]*dill"):
ar_ancMig_iso_recBot.append(dill.load(open(filename)))
In [114]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_iso_recBot]
df = pd.DataFrame(data=success, \
columns=['nu1_0', 'nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_2_0', 'nu2_2_0', 'Tb_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tb_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[114]:
In [15]:
%%px --local
pts_l = [50, 60, 70] # make finer grid
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
In [113]:
#ar_ancMig_iso_recBot = lbview.map(run_dadi, repeat(p0, 10))
In [16]:
ar_ancMig_iso_recBot = []
for filename in glob("OUT_2D_models/ancMig_iso_recBot_[0-9]*dill"):
ar_ancMig_iso_recBot.append(dill.load(open(filename)))
In [19]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_iso_recBot]
df = pd.DataFrame(data=success, \
columns=['nu1_0', 'nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_2_0', 'nu2_2_0', 'Tb_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tb_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[19]:
In [119]:
%%px --local
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
In [120]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,9:18])
popt
Out[120]:
In [121]:
p0 = popt
In [122]:
#ar_ancMig_iso_recBot = lbview.map(run_dadi, repeat(p0, 10))
In [123]:
ar_ancMig_iso_recBot = []
for filename in glob("OUT_2D_models/ancMig_iso_recBot_[0-9]*dill"):
ar_ancMig_iso_recBot.append(dill.load(open(filename)))
In [124]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_iso_recBot]
df = pd.DataFrame(data=success, \
columns=['nu1_0', 'nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_2_0', 'nu2_2_0', 'Tb_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tb_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[124]:
This looks like good convergence, but note that the second population size change has quite divergent parameter value combinations with almost equal likelihood. So its time and magnitude cannot be inferred accurately.
This isolation then bottleneck model can be turned into the the previous ancMig recBotIso model by setting $T_i$ to zero.
In [21]:
ll_s = -11975
ll_c = -11947
D = 2*(ll_c - ll_s)
D
Out[21]:
In [22]:
# calculate p-value for Chi-square dist. with 0 and 1 d.o.f.
# the weights specify a weighted sum of chi^2 distributions with 0 and 1 d.o.f
# this is because in the simple model, Ti is fixed at the boundary of the parameter space
p = dadi.Godambe.sum_chi2_ppf(D, weights=(0.5, 0.5))
p
Out[22]:
A $T_i$ greater than zero seems to be highly significant. I doubt that linkage in the data can compromise this result. This means that allowing for the cessation of gene flow to happen earlier than the 2nd population size change significantly improves the fit to the data.
In [23]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[3,9:18])
popt
Out[23]:
Note, above I have not chosen the most likely parameter combination as popt! Instead, I have chosen the 4th best combination which also has the highest (and therefore most plausible) value for $T_b$ (the time of the recent pop size change). This parameter combination has practically identical likelihood to the most likely parameter combination in the above table.
In [24]:
model = func_ex(popt, ns, pts_l)
ll = dadi.Inference.ll_multinom(model, sfs2d)
ll
Out[24]:
In [25]:
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
L = sfs2d.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.5f}.".format(theta/L)
In [26]:
mu = 3e-9
L = sfs2d.data.sum() # this sums over all entries in the spectrum, including masked ones, i. e. also contains invariable sites
print "The total sequence length for the 2D spectrum is {0:,}.".format(int(L))
N_ref = theta/L/mu/4
print "The effective size of the ancestral population of ery and par (in number of diploid individuals) implied by this theta is:\n {0:,}.".format(int(N_ref))
In [27]:
nu1,nu2,Tc,m1,m2,Ti,nu1_2,nu2_2,Tb = popt
In [34]:
print "The ancestral population split apart {0:,} generations ago.".format(int((Tc+Ti+Tb)*2*N_ref)),
print "Immediately after the split the ERY population changed to a size of {0:,} and the PAR population to {1:,}.".format(int(nu1*N_ref), int(nu2*N_ref)),
print "Since the split of the ancestral population, PAR received 1 individual from ERY every {0:.2f} generations,".format(1.0/(m1/2*nu2)),
print "while ERY received 1 PAR individual every {0:.2f} generations.".format(1.0/(m2/2*nu1)),
print "Put another way: The PAR population contained a constant proportion of {0:.2e} of new immigrant alleles each generation".format(m1/2/N_ref/nu2),
print "and the ERY population contained a constant proportion of {0:.2e} of new immigrant alleles each generation.".format(m2/2/N_ref/nu1),
print "ERY and PAR remained in contact for {0:,} generations.".format(int(Tc*2*N_ref)),
print "{0:,} generations ago gene flow between ERY and PAR had ceased.".format(int(Ti*2*N_ref)),
print "{0:,} generations ago, ERY and PAR underwent a second stepwise simultaneous population size change.".format(int(Tb*2*N_ref)),
print "ERY changed to a size of {0:,} individuals and PAR to a size of {1:,} individuals.".format(int(nu1_2*N_ref), int(nu2_2*N_ref))
It is remarkable that this model gives a time estimate for the cessation of gene flow almost identical to the previous model, ancMig recBotIso that enforced a coincident pop size change, and is also a very similar estimate to the one given by the ancient migration model (26,970). In contrast to the previous model, this model allows the second pop size change to happen at any time (Tb) after the cessation of gene flow (Ti+Tb). Similar to the recent bottleneck model, this model estimates a very recent and very severe population size change for PAR, but in contrast to the recent bottleneck model it infers a substantial population size increase for ERY. I wonder whether this is just an artifact in the data (allele-drop-out?, filtering?) or indicating a real event. For instance, if the environment on the north side of the Pyrenees experienced a drastic change during the last few hundred years. Maybe a regrowth of forest has led to increasingly small and isolated populations of PAR? The question would then be why this is not mirrored on the Spanish side.
In [35]:
dadi.Plotting.plot_2d_comp_multinom(data=sfs2d, model=model, vmin=1)
I would like to compare the best-fit spectrum of the ancMig recBotIso model with the one of this isolation then bottleneck model.
In [37]:
func_ex = dadi.Numerics.make_extrap_log_func(ancMig_recBotIso)
popt_ancMig_recBotIso = [1.44238205, 4.86198073, 2.3766401, 0.06626632, 0.46005763, 1.78231791, 1.4500263, 0.03946682]
model_ancMig_recBotIso = func_ex(popt_ancMig_recBotIso, ns, pts_l)
model_ancMig_recBotIso = dadi.Inference.optimally_scaled_sfs(model_ancMig_recBotIso, sfs2d)
model = dadi.Inference.optimally_scaled_sfs(model, sfs2d)
In [72]:
dadi.Plotting.plot_2d_comp_multinom(data=model.fold() , model=model_ancMig_recBotIso.fold(), \
vmin=1, title=['Iso then Bot', 'Bot+Iso'], pop_ids=['ery', 'par'])
As can be seen in the residual plot, the main difference between these two models is in how their predicted counts in the [0, 1] and [1, 0] SNP classes. Enforcing the pop size change to be as old as the cessation of gene flow leads to higher predictions for those count classes than allowing the pop size changes to be very recent.
In [40]:
def ancMig_bot_iso(params, ns, pts):
"""
params = (nu1,nu2,Tc,m1,m2,nu1_2,nu2_2,Tb,Ti)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time Tc followed by a potential second size change at time [Tb+Ti]. The split coincides
with the first potential pop size change. At time Ti there is a cessation of gene flow.
Note, that this model does not allow the cessation of gene flow to occur earlier than the second
pop size change.
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
Tc: Time of gene flow after split till 2nd pop size change (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
nu1_2: pop size ratio of pop 1 after second size change
nu2_2: pop size ratio of pop 2 after second size change
Tb: Time from second pop size change till cessation of gene flow
Ti: Time of isolation after cessation of gene flow
The split lies Tc+Tb+Ti * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,Tc,m1,m2,nu1_2,nu2_2,Tb,Ti = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, Tc, nu1, nu2, m12=m2, m21=m1)
# potential second pop size change
phi = dadi.Integration.two_pops(phi, xx, Tb, nu1_2, nu2_2, m12=m2, m21=m1)
# divergence without gene flow
phi = dadi.Integration.two_pops(phi, xx, Ti, nu1_2, nu2_2, m12=0, m21=0)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [41]:
cl[:].push(dict(ancMig_bot_iso=ancMig_bot_iso))
Out[41]:
In [42]:
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(ancMig_bot_iso)
In [43]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/ancMig_bot_iso" # set file name stub for opt. result files
fixed_params = None
In [44]:
%%px --local
# set lower and upper bounds to nu1,nu2,Tc,m1,m2,nu1_2,nu2_2,Tb,Ti
upper_bound = [1e4, 1e4, 6, 10, 10, 1e4, 1e4, 6, 6]
lower_bound = [1e-4, 1e-4, 0, 0, 0, 1e-4, 1e-4, 0, 0]
In [45]:
popt = df.sort_values(by='-logL', ascending=True).iloc[3,9:18]
popt
Out[45]:
In [46]:
p0 = [1.492, 4.75, 2.34, 0.069, 0.451, 3.47, 0.032, 0.00056, 0.039]
In [47]:
#ar_ancMig_bot_iso = lbview.map(run_dadi, repeat(p0, 10))
In [48]:
ar_ancMig_bot_iso.done()
Out[48]:
In [49]:
ar_ancMig_bot_iso = []
for filename in glob("OUT_2D_models/ancMig_bot_iso_[0-9]*dill"):
ar_ancMig_bot_iso.append(dill.load(open(filename)))
In [50]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_bot_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0', 'nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'nu1_2_0', 'nu2_2_0', 'Tb_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tb_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[50]:
In [51]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,9:18])
popt
Out[51]:
In [52]:
p0 = popt
In [53]:
#ar_ancMig_bot_iso = lbview.map(run_dadi, repeat(p0, 10))
In [54]:
ar_ancMig_bot_iso = []
for filename in glob("OUT_2D_models/ancMig_bot_iso_[0-9]*dill"):
ar_ancMig_bot_iso.append(dill.load(open(filename)))
In [55]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_bot_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0', 'nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'nu1_2_0', 'nu2_2_0', 'Tb_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tb_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[55]:
In [58]:
%%px --local
pts_l = [50, 60, 70] # make finer grid
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
In [59]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,9:18])
popt
Out[59]:
In [60]:
p0 = popt
In [61]:
#ar_ancMig_bot_iso = lbview.map(run_dadi, repeat(p0, 10))
In [62]:
ar_ancMig_bot_iso.done()
Out[62]:
In [63]:
ar_ancMig_bot_iso = []
for filename in glob("OUT_2D_models/ancMig_bot_iso_[0-9]*dill"):
ar_ancMig_bot_iso.append(dill.load(open(filename)))
In [64]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_ancMig_bot_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0', 'nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'nu1_2_0', 'nu2_2_0', 'Tb_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'nu1_2_opt', 'nu2_2_opt', 'Tb_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[64]:
This seems to have converged on essentially the recent bottleneck best-fit model found previously, which had a best -logL of 12,034. Note that this model can be turned into the recent bottleneck model by setting $T_i$ to zero, i. e. no period without gene flow. This bottleneck then isolation model is therefore not an improvement towards the recent bottleneck model and this indicates that gene flow must have ceased before the recent change in population size.
In [66]:
popt = df.sort_values(by='-logL', ascending=True).iloc[0,9:18]
popt
Out[66]:
Let's compare with the best-fit parameter values of the recent bottleneck model.
In [68]:
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch*dill"):
ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
In [70]:
l = 2*8+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df_rb = pd.DataFrame(data=returned, \
columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'ery_1_opt','par_1_opt','T1_opt','ery_2_opt','par_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
In [71]:
popt_rb = df_rb.sort_values(by='-logL', ascending=True).head(10).iloc[0,8:16]
popt_rb
Out[71]:
The main difference is that the recent bottleneck model estimated a 10 times more severe bottleneck for both populations at a time 10 times more recent.
In [ ]:
In [ ]:
In [ ]:
The two populations are relatively close to the hybrid zone centre at the Col de la Quillane in the Pyrenees. The clines for some characters have been shown to be many kilometers wide. It can therefore be assumed that there is some recent gene flow between the two populations. Can this gene flow be detected from this spectrum?
In [44]:
def split_mig_iso_mig(params, ns, pts):
"""
params = (nu1,nu2,Tc,m1,m2,Ti,Tsc)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time Tc followed by a period of complete isolation Ti which is followed by a restart of
migration until present. Migration rates in the two epochs are assumed to be equal (and constant).
nu1: Size of population 1 after split.
nu2: Size of population 2 after split.
Tc: Time of gene flow after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
Ti: Time of isolation after cessation of gene flow
Tsc: Time of gene flow after period of isolation until present
The split lies Tc+Ti+Tsc * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,Tc,m1,m2,Ti,Tsc = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, Tc, nu1, nu2, m12=m2, m21=m1)
# divergence without gene flow
phi = dadi.Integration.two_pops(phi, xx, Ti, nu1, nu2, m12=0, m21=0)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, Tsc, nu1, nu2, m12=m2, m21=m1)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [46]:
cl[:].push(dict(split_mig_iso_mig=split_mig_iso_mig))
Out[46]:
In [47]:
%%px --local
func = split_mig_iso_mig
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [48]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 300 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/split_mig_iso_mig" # set file name stub for opt. result files
fixed_params = None
In [59]:
%%px --local
# set lower and upper bounds to nu1, nu2, Tc, m1, m2, Ti, Tsc
upper_bound = [1e4, 1e4, 6, 10, 10, 6, 6]
lower_bound = [1e-4, 1e-4, 0, 0, 0, 0, 0]
In [60]:
# using the best fit parameters from the ancient migration model and adding an initial Tsc of 0.05
p0 = [2.76935944, 6.77011784, 5.19634719, 0.04983223, 0.24913079, 0.09420312, 0.05]
In [53]:
#ar_mig_iso_mig = lbview.map(run_dadi, repeat(p0, 10))
In [61]:
% ll OUT_2D_models/split_mig_iso_mig*dill
In [62]:
ar_mig_iso_mig = []
for filename in glob("OUT_2D_models/split_mig_iso_mig*dill"):
ar_mig_iso_mig.append(dill.load(open(filename)))
In [63]:
import pandas as pd
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_mig_iso_mig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'Tsc_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', 'Tsc_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[63]:
There is no good convergence yet.
In [59]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,7:14])
popt
Out[59]:
In [56]:
%%px --local
pts_l = [60, 70, 80]
fold = 1
In [60]:
%%px --local
# set lower and upper bounds to nu1, nu2, Tc, m1, m2, Ti, Tsc
upper_bound = [1e4, 1e4, 8, 10, 10, 6, 6]
lower_bound = [1e-4, 1e-4, 0, 0, 0, 0, 0]
In [61]:
p0 = popt
#ar_mig_iso_mig = lbview.map(run_dadi, repeat(p0, 10))
In [51]:
ar_mig_iso_mig = []
for filename in glob("OUT_2D_models/split_mig_iso_mig_[0-9]*dill"):
ar_mig_iso_mig.append(dill.load(open(filename)))
In [53]:
% ll OUT_2D_models/split_mig_iso_mig_[0-9]*dill
In [54]:
import pandas as pd
l = 2*7+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_mig_iso_mig]
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'Tsc_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', 'Tsc_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[54]:
This looks like good convergence.
The best-fit ancient migration model had a log likelihood of -12156. So only one log likelihood unit worse.
In [64]:
ll_s = 12156.061309
ll_c = 12155.538032
D = 2 * (ll_s - ll_c)
D
Out[64]:
In [65]:
# calculate p-value for Chi-square dist.
# the weights specify a weighted sum of chi^2 distributions with 1 and 2 d.o.f
# this is because Tsc is 0 in the split_asym_mig_iso model and at the boundary of the parameter space
p = dadi.Godambe.sum_chi2_ppf(D, weights=(0.5, 0.5))
p
Out[65]:
A recent restart of gene cannot be detected from this corrected spectrum.
I arrived at the same result with the uncorrected spectrum.
In [55]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,7:14])
popt
Out[55]:
In [56]:
model = func_ex(popt, ns, pts_l)
In [57]:
ll_model = dadi.Inference.ll_multinom(model, sfs2d)
ll_model
Out[57]:
In [63]:
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
L = sfs2d.data.sum()
mu = 3e-9
N_ref = theta/L/mu/4
print "The ancestral population size implied by the theta is is {0:,}.".format(int(N_ref))
print "The optimal time of secondary contact inferred is {0:,} generations.".format(int(popt[-1]*2*N_ref))
I would like to use a module of the programme moments to plot the ancient migration model.
In [1]:
% ll /home/claudius/Downloads/moments/moments
Installed in /usr/local/lib/python2.7/dist-packages/moments-1.0.0-py2.7.egg
In [1]:
import sys
sys.path
Out[1]:
In [2]:
sys.path.insert(0, '/usr/local/lib/python2.7/dist-packages/moments-1.0.0-py2.7.egg')
In [3]:
sys.path
Out[3]:
In [4]:
import moments
In [5]:
from glob import glob
import dill
from utility_functions import *
import pandas as pd
# turn on floating point division by default, old behaviour via '//'
from __future__ import division
import numpy as np
In [6]:
%matplotlib inline
import pylab
pylab.rcParams['figure.figsize'] = [12, 10]
pylab.rcParams['font.size'] = 14
In [7]:
# load spectrum modified with Ludovic's correction, p=35
sfs2d = moments.Spectrum.from_file("EryPar_modified.2dsfs")
In [8]:
ns = sfs2d.sample_sizes
In [17]:
def split_asym_mig_iso(params, ns): # note, I have removed the pts argument from the argument list
"""
params = (nu1,nu2,Tc,m1,m2,Ti)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time Tc followed by complete isolation until present.
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
Tc: Time of gene flow after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
Ti: Time of isolation after cessation of gene flow
The split lies Tc+Ti * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,Tc,m1,m2,Ti = params
# fs for equilibrium ancestral population
sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1])
fs = moments.Spectrum(sts)
# split
fs = moments.Manips.split_1D_to_2D(fs, ns[0], ns[1])
# divergence with potentially asymmetric migration
fs.integrate([nu1, nu2], Tc, m=np.array([[0, m2], [m1, 0]]))
# divergence without gene flow
fs.integrate([nu1, nu2], Ti, m=np.array([[0,0], [0,0]]))
return fs
In [13]:
ar_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/split_asym_mig_iso_[0-9]*dill"):
ar_split_asym_mig_iso.append(dill.load(open(filename)))
In [14]:
l = 2*6+1
success = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[14]:
In [15]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,6:12])
popt
Out[15]:
I am going to use the optimal values from dadi.
In [19]:
# generate ModelPlot object
plot_mod = moments.ModelPlot.generate_model(split_asym_mig_iso, popt, ns)
In [26]:
moments.ModelPlot.plot_model?
In [30]:
# generate the plot for the ancient migration model
# the figure is writen to file, which I am reading in in the following cell using Markdown syntax
moments.ModelPlot.plot_model(plot_mod, save_file='model.png', pop_labels=['ery', 'par'], gen_time=1e-3, \
gen_time_units="KY", fig_title='ancient migration', nref=113522, reverse_timeline=True)
The numbers on the right in the upper plot are a bit misleading. It's a 27 (ky ago) and a zero.
I would like to fit the ancient migration with the programme moments.
In [16]:
from ipyparallel import Client
cl = Client()
cl.ids
Out[16]:
In [31]:
from glob import glob
import dill
from utility_functions import *
import pandas as pd
# turn on floating point division by default, old behaviour via '//'
from __future__ import division
import numpy as np
from itertools import repeat
In [18]:
%matplotlib inline
import pylab
pylab.rcParams['figure.figsize'] = [12, 10]
pylab.rcParams['font.size'] = 14
In [19]:
%%px --local
# run whole cell on all engines a well as in the local IPython session
import sys
sys.path.insert(0, '/usr/local/lib/python2.7/dist-packages/moments-1.0.0-py2.7.egg')
import moments
In [27]:
%%px --local
# load spectrum modified with Ludovic's correction, p=35
sfs2d = moments.Spectrum.from_file("EryPar_modified.2dsfs")
In [20]:
def split_asym_mig_iso(params, ns): # note, I have removed the pts argument from the argument list
"""
params = (nu1,nu2,Tc,m1,m2,Ti)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time Tc followed by complete isolation until present.
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
Tc: Time of gene flow after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
Ti: Time of isolation after cessation of gene flow
The split lies Tc+Ti * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,Tc,m1,m2,Ti = params
# fs for equilibrium ancestral population
sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1])
fs = moments.Spectrum(sts)
# split
fs = moments.Manips.split_1D_to_2D(fs, ns[0], ns[1])
# divergence with potentially asymmetric migration
fs.integrate([nu1, nu2], Tc, m=np.array([[0, m2], [m1, 0]]))
# divergence without gene flow
fs.integrate([nu1, nu2], Ti, m=np.array([[0,0], [0,0]]))
return fs
In [21]:
cl[:].push(dict(split_asym_mig_iso=split_asym_mig_iso))
Out[21]:
In [24]:
def run_moments(p_init): # for the function to be called with map, it needs to have one input variable
"""
p_init: initial parameter values to run optimisation from
"""
if perturb == True:
p_init = moments.Misc.perturb_params(p_init, fold=fold,
upper_bound=upper_bound, lower_bound=lower_bound)
# note upper_bound and lower_bound variables are expected to be in the namespace of each engine
# run optimisation of paramters
popt = moments_opt_func(p0=p_init, data=sfs, model_func=func, \
lower_bound=lower_bound, upper_bound=upper_bound, \
verbose=verbose, maxiter=maxiter, full_output=full_output, \
fixed_params=fixed_params)
# pickle to file
import dill
name = outname[:] # make copy of file name stub!
for p in p_init:
name += "_%.4f" % (p)
with open(name + ".dill", "w") as fh:
dill.dump((p_init, popt), fh)
return p_init, popt
In [12]:
moments.Inference.optimize_log_fmin?
In [28]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
func = split_asym_mig_iso
moments_opt_func = moments.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/moments/split_asym_mig_iso" # set file name stub for opt. result files
fixed_params = None
In [29]:
%%px --local
# set lower and upper bounds to nu1, nu2, Tc, m1, m2 and Ti
upper_bound = [1e3, 1e3, 10, 10, 10, 10] # note, I have increased the upper bound for T
lower_bound = [1e-3, 1e-3, 0, 0, 0, 0]
In [30]:
lbview = cl.load_balanced_view()
In [32]:
p0 = [ 3.42873459, 8.372007, 3.71956526, 0.04062737, 0.20232337, 0.11878969]
#ar_split_asym_mig_iso = lbview.map(run_moments, repeat(p0, 10), block=False)
In [33]:
cl.ids
Out[33]:
In [34]:
ar_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/moments/split_asym_mig_iso_[0-9]*dill"):
ar_split_asym_mig_iso.append(dill.load(open(filename)))
In [35]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[35]:
There is no good convergence yet.
In [36]:
%%px
maxiter = 300
fold = 1
In [37]:
p0 = [2.742542, 6.642431, 4.948508, 0.051360, 0.247458, 0.086862]
#ar_split_asym_mig_iso = lbview.map(run_moments, repeat(p0, 10), block=False)
In [38]:
ar_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/moments/split_asym_mig_iso_[0-9]*dill"):
ar_split_asym_mig_iso.append(dill.load(open(filename)))
In [39]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[39]:
The best three parameter combinations are all hitting the upper bound on the Tc parameter. Also note, that that these three optimisation runs already started with an extremely high Tc. Dadi's best parameter combination had a much lower Tc and slightly better likelihood (12156).
I wonder whether I can improve the fit, if I already let the ancestral population increase in size before the split.
In [41]:
def growth_split_asym_mig_iso(params, ns): # note, I have removed the pts argument from the argument list
"""
params = (nu0,T0,nu1,nu2,Tc,m1,m2,Ti)
ns = (n1,n2)
Let ancestral population grow or shrink for time T0, then split into two populations of specifed size,
with potentially asymmetric migration for a time Tc followed by complete isolation until present.
nu0: population size of ancestral population before split (with respect to Na)
T0: time of increased/reduced pop size before split (in units of 2*Na generations)
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
Tc: Time of gene flow after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
Ti: Time of isolation after cessation of gene flow
The split lies Tc+Ti * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu0,T0,nu1,nu2,Tc,m1,m2,Ti = params
# fs for equilibrium ancestral population
sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1])
fs = moments.Spectrum(sts)
# ancient population size change
fs.integrate([nu0], T0)
# split
fs = moments.Manips.split_1D_to_2D(fs, ns[0], ns[1])
# divergence with potentially asymmetric migration
fs.integrate([nu1, nu2], Tc, m=np.array([[0, m2], [m1, 0]]))
# divergence without gene flow
fs.integrate([nu1, nu2], Ti, m=np.array([[0,0], [0,0]]))
return fs
In [42]:
cl[:].push(dict(growth_split_asym_mig_iso=growth_split_asym_mig_iso))
Out[42]:
In [43]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
func = growth_split_asym_mig_iso
moments_opt_func = moments.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/moments/growth_split_asym_mig_iso" # set file name stub for opt. result files
fixed_params = None
In [44]:
%%px --local
# set lower and upper bounds to nu1, nu2, Tc, m1, m2 and Ti
upper_bound = [1e3, 10, 1e3, 1e3, 10, 10, 10, 10] # note, I have increased the upper bound for T
lower_bound = [1e-3, 0, 1e-3, 1e-3, 0, 0, 0, 0]
In [45]:
p0 = [3, 4, 1.754277, 2.723380, 3.999903, 0.049779, 0.248170, 0.089072]
ar_growth_split_asym_mig_iso = lbview.map(run_moments, repeat(p0, 10), block=False)
In [46]:
ar_growth_split_asym_mig_iso.elapsed/60
Out[46]:
In [52]:
ar_growth_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/moments/growth_split_asym_mig_iso_[0-9]*dill"):
ar_growth_split_asym_mig_iso.append(dill.load(open(filename)))
In [53]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_growth_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu0_0', 'T0_0', 'nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu0_opt', 'T0_opt', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[53]:
Let's refine the search.
In [57]:
popt = np.array(df.sort_values(by='-logL', ascending=True).head(10).iloc[0, 8:16])
popt
Out[57]:
In [58]:
%%px
fold = 1
maxiter = 300
In [59]:
p0 = popt
#ar_growth_split_asym_mig_iso = lbview.map(run_moments, repeat(p0, 10), block=False)
In [60]:
ar_growth_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/moments/growth_split_asym_mig_iso_[0-9]*dill"):
ar_growth_split_asym_mig_iso.append(dill.load(open(filename)))
In [61]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_growth_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu0_0', 'T0_0', 'nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu0_opt', 'T0_opt', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[61]:
Still, no good convergence. Also the optimised parameters values seem to be correlated with the starting values.
Let's try to refine the search one more time.
In [62]:
popt = np.array(df.sort_values(by='-logL', ascending=True).head(10).iloc[0, 8:16])
popt
Out[62]:
In [63]:
p0 = popt
#ar_growth_split_asym_mig_iso = lbview.map(run_moments, repeat(p0, 10), block=False)
In [ ]:
ar_growth_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/moments/growth_split_asym_mig_iso_[0-9]*dill"):
ar_growth_split_asym_mig_iso.append(dill.load(open(filename)))
In [64]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_growth_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu0_0', 'T0_0', 'nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu0_opt', 'T0_opt', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[64]:
There is still no good convergence and the best parameter combination is still worse than the ancient migration model above.
In [1]:
from ipyparallel import Client
cl = Client()
cl.ids
Out[1]:
In [18]:
lbview = cl.load_balanced_view()
In [2]:
%%px --local
# run whole cell on all engines a well as in the local IPython session
import numpy as np
import sys
sys.path.insert(0, '/home/claudius/Downloads/dadi')
import dadi
In [74]:
def growth_split_asym_mig_iso(params, ns, pts):
"""
params = (nu1,nu2,Tc,m1,m2,Ti)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time Tc followed by complete isolation until present.
nu0: relative population size of ancient population (with respect to Na)
T0: length of time of ancient population size change (in 2*Na generations)
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
Tc: Time of gene flow after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
Ti: Time of isolation after cessation of gene flow
The split lies Tc+Ti * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu0,T0,nu1,nu2,Tc,m1,m2,Ti = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# ancient population growth
phi = dadi.Integration.one_pop(phi, xx, T0, nu0)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, Tc, nu1, nu2, m12=m2, m21=m1)
# divergence without gene flow
phi = dadi.Integration.two_pops(phi, xx, Ti, nu1, nu2, m12=0, m21=0)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [65]:
cl[:].push(dict(growth_split_asym_mig_iso=growth_split_asym_mig_iso))
Out[65]:
In [75]:
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(growth_split_asym_mig_iso)
In [67]:
from glob import glob
import dill
from utility_functions import *
import pandas as pd
# turn on floating point division by default, old behaviour via '//'
from __future__ import division
import numpy as np
from itertools import repeat
In [13]:
def run_dadi(p_init): # for the function to be called with map, it needs to have one input variable
"""
p_init: initial parameter values to run optimisation from
"""
if perturb == True:
p_init = dadi.Misc.perturb_params(p_init, fold=fold,
upper_bound=upper_bound, lower_bound=lower_bound)
# note upper_bound and lower_bound variables are expected to be in the namespace of each engine
# run optimisation of paramters
popt = dadi_opt_func(p0=p_init, data=sfs, model_func=func_ex, pts=pts_l, \
lower_bound=lower_bound, upper_bound=upper_bound, \
verbose=verbose, maxiter=maxiter, full_output=full_output, \
fixed_params=fixed_params)
# pickle to file
import dill
name = outname[:] # make copy of file name stub!
for p in p_init:
name += "_%.4f" % (p)
with open(name + ".dill", "w") as fh:
dill.dump((p_init, popt), fh)
return p_init, popt
In [14]:
%%px --local
# load spectrum modified with Ludovic's correction, p=35
sfs2d = dadi.Spectrum.from_file("EryPar_modified.2dsfs")
In [21]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
pts_l = [50, 60, 70]
func_ex = func_ex
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = sfs2d
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/dadi/growth_split_asym_mig_iso" # set file name stub for opt. result files
fixed_params = None
In [16]:
%matplotlib inline
import pylab
pylab.rcParams['figure.figsize'] = [12, 10]
pylab.rcParams['font.size'] = 14
In [17]:
%%px --local
# set lower and upper bounds to nu0, T0, nu1, nu2, Tc, m1, m2 and Ti
upper_bound = [1e3, 10, 1e3, 1e3, 10, 10, 10, 10] # note, I have increased the upper bound for T
lower_bound = [1e-3, 0, 1e-3, 1e-3, 0, 0, 0, 0]
In [22]:
# using best parameter combination from optimisations with moments
p0 = [0.322405, 0.721422, 0.939739, 2.281812, 3.222201, 0.147156, 0.714855, 0.029154]
#ar_growth_split_asym_mig_iso = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [23]:
ar_growth_split_asym_mig_iso.done()
Out[23]:
In [24]:
ar_growth_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/dadi/growth_split_asym_mig_iso_[0-9]*dill"):
ar_growth_split_asym_mig_iso.append(dill.load(open(filename)))
In [25]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_growth_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu0_0', 'T0_0', 'nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu0_opt', 'T0_opt', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[25]:
No good convergence yet. Let's refine the grid and increase maxiter and reduce fold.
In [68]:
%%px
fold = 1
maxiter = 300
pts_l = [50, 60, 70]
In [27]:
popt = np.array(df.sort_values(by='-logL', ascending=True).head(10).iloc[0, 8:16])
popt
Out[27]:
In [28]:
p0 = popt
#ar_growth_split_asym_mig_iso = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [29]:
ar_growth_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/dadi/growth_split_asym_mig_iso_[0-9]*dill"):
ar_growth_split_asym_mig_iso.append(dill.load(open(filename)))
In [31]:
get_flag_count(ar_growth_split_asym_mig_iso, NM=True)
None of the optimisation runs achieved convergence within the limit of 300 iterations.
In [30]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_growth_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu0_0', 'T0_0', 'nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu0_opt', 'T0_opt', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[30]:
It is difficult to achieve convergence. It may already be too complex. It does not achieve a higher likelihood than the ancient migration model (-logL: 12156).
I conclude that the addition of a population size change of the ancestral population of ERY and PAR does not improve the fit. However, it does achieve a similarly good fit with a time since split between ERY and PAR much shorter than inferred with the ancient migration model.
Maybe I can achieve convergence with another optimsation algorithm.
In [32]:
dadi.Inference.optimize?
In [33]:
%%px --local
# set up global variables on engines required for run_dadi function call
ns = sfs2d.sample_sizes # both populations have the same sample size
pts_l = [50, 60, 70]
func_ex = func_ex
dadi_opt_func = dadi.Inference.optimize # uses BFGS algorithm
sfs = sfs2d
perturb = True
fold = 1 # perturb randomly up to `fold` times 2-fold
maxiter = 300 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_2D_models/dadi/growth_split_asym_mig_iso" # set file name stub for opt. result files
fixed_params = None
In [34]:
popt = np.array(df.sort_values(by='-logL', ascending=True).head(10).iloc[0, 8:16])
popt
Out[34]:
In [35]:
p0 = popt
#ar_growth_split_asym_mig_iso = lbview.map(run_dadi, repeat(p0, 10), block=False)
In [37]:
ar_growth_split_asym_mig_iso.elapsed/60
Out[37]:
In [69]:
ar_growth_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/dadi/growth_split_asym_mig_iso_[0-9]*dill"):
ar_growth_split_asym_mig_iso.append(dill.load(open(filename)))
In [70]:
l = 2*len(p0)+1
success = [flatten(out)[:l] for out in ar_growth_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu0_0', 'T0_0', 'nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu0_opt', 'T0_opt', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[70]:
There is still no good convergence. Quite different demographic scenarios have almost identical likelihood. One of which is the ancient migration model from above, which had -logL of 12,156 with it's best fit parameter values. A population size change in the ancestral population to ERY and PAR cannot be inferred from this spectrum.
This illustrates the uncertainty in the inference of the demographic scenario. The most parsimonious model found is the ancient migration model.
In [80]:
def translate_time(i, j):
"""
i: row index in parameter table with best-fit combinations at the top
j: index for parameter in parameter list
"""
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[i,8:16])
model = func_ex(popt, ns, pts_l)
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
L = sfs2d.data.sum()
mu = 3e-9
N_ref = theta/L/mu/4
return int(popt[j]*2*N_ref)
In [85]:
# get translation of time of split parameter Tc for the top 5 parameter combinations
map(translate_time, range(5), repeat(4, 5))
Out[85]:
Note, that despite the much smaller values for Tc (the time of split between ERY and PAR) in genetic units (2Nref) as compared to the ancient migration model, this ancestral size change model infers very similar values for Tc in absolute units (generations) (compare eith ancient).
In [86]:
# get translation of time of isolation for the top 5 parameter combinations
map(translate_time, range(5), repeat(7, 5))
Out[86]:
The inferred time of complete isolation is also very similar to the ancient migration model best-fit parameter values.
In [87]:
# get translation of time of isolation for the top 5 parameter combinations
map(translate_param, range(5), repeat(5, 5))
Out[87]:
In [91]:
def translate_migration(i, j, k):
"""
i: row index in parameter table with best-fit combinations at the top
j: index for migration parameter in parameter list
k: index of pop size parameter in param list
"""
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[i,8:16])
model = func_ex(popt, ns, pts_l)
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
L = sfs2d.data.sum()
mu = 3e-9
N_ref = theta/L/mu/4
return popt[j]/2/N_ref/popt[k]
In [92]:
# get translation of migration parameter ery->par for top 5 param comb
map(translate_migration, range(5), repeat(5, 5), repeat(3, 5))
Out[92]:
In [93]:
# get translation of migration parameter par->ery for top 5 param comb
map(translate_migration, range(5), repeat(6, 5), repeat(2, 5))
Out[93]:
The migration rate in proportion of new immigrant alleles per generation is in the same order of magnitude as for the ancient migration model. Migration from ery into par is again about 10 times weaker than in the other direction. The migration rates with this ancestral population size change model are generally slightly higher than those from the ancient migration model.
In conclusion it can be said that this ancestral population size change model predicts a very similar demographic history as the ancient migration model.
In [1]:
import numpy as np
import sys
sys.path.insert(0, '/home/claudius/Downloads/dadi')
import dadi
In [2]:
from glob import glob
import dill
from utility_functions import *
import pandas as pd
# turn on floating point division by default, old behaviour via '//'
from __future__ import division
In [3]:
# load spectrum modified with Ludovic's correction, p=35
sfs2d = dadi.Spectrum.from_file("EryPar_modified.2dsfs")
In [4]:
def split_asym_mig_iso(params, ns, pts):
"""
params = (nu1,nu2,Tc,m1,m2,Ti)
ns = (n1,n2)
Split into two populations of specifed size, with potentially asymmetric migration
for a time Tc followed by complete isolation until present.
nu1: population size ratio of population 1 after split (with respect to Na)
nu2: population size ratio of population 2 after split (with respect to Na)
Tc: Time of gene flow after split (in units of 2*Na generations)
m1: Migration rate from ery into par (in units of 2*Na ind per generation)
m2: Migration rate from par into ery (in units of 2*Na ind per generation)
Ti: Time of isolation after cessation of gene flow
The split lies Tc+Ti * 2Na generations in the past.
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,Tc,m1,m2,Ti = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
# split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# divergence with potentially asymmetric migration
phi = dadi.Integration.two_pops(phi, xx, Tc, nu1, nu2, m12=m2, m21=m1)
# divergence without gene flow
phi = dadi.Integration.two_pops(phi, xx, Ti, nu1, nu2, m12=0, m21=0)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [5]:
ns = sfs2d.sample_sizes
pts_l = [50, 60, 70]
func_ex = dadi.Numerics.make_extrap_log_func(split_asym_mig_iso)
In [6]:
ar_split_asym_mig_iso = []
for filename in glob("OUT_2D_models/split_asym_mig_iso_[0-9]*dill"):
ar_split_asym_mig_iso.append(dill.load(open(filename)))
In [7]:
l = 2*6+1
success = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'Tc_0', 'm1_0', 'm2_0', 'Ti_0', 'nu1_opt', 'nu2_opt', 'Tc_opt', 'm1_opt', 'm2_opt', 'Ti_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[7]:
In [8]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,6:12])
popt
Out[8]:
This is the definition of FIM_uncert in Godambe.py after I changed it:
def FIM_uncert(func_ex, grid_pts, p0, data, log=False, multinom=True, eps=0.01, return_FIM=False): """ Parameter uncertainties from Fisher Information Matrix
Returns standard deviations of parameter values.
func_ex: Model function
all_boot: List of bootstrap frequency spectra
p0: Best-fit parameters for func_ex
data: Original data frequency spectrum
eps: Fractional stepsize to use when taking finite-difference derivatives.
Note that if eps*param is < 1e-6, then the step size for that parameter
will simply be eps, to avoid numerical issues with small parameter
perturbations.
log: If True, assume log-normal distribution of parameters. Returned values
are then the standard deviations of the *logs* of the parameter values,
which can be interpreted as relative parameter uncertainties.
multinom: If True, assume model is defined without an explicit parameter for
theta. Because uncertainty in theta must be accounted for to get
correct uncertainties for other parameters, this function will
automatically consider theta if multinom=True. In that case, the
final entry of the returned uncertainties will correspond to
theta.
return_FIM: if true, return the Hessian matrix, else return standard devia-
tions of parameters
"""
if multinom:
func_multi = func_ex
model = func_multi(p0, data.sample_sizes, grid_pts)
theta_opt = Inference.optimal_sfs_scaling(model, data)
p0 = list(p0) + [theta_opt]
func_ex = lambda p, ns, pts: p[-1]*func_multi(p[:-1], ns, pts)
H = get_godambe(func_ex, grid_pts, [], p0, data, eps, log, just_hess=True)
#return numpy.sqrt(numpy.diag(numpy.linalg.inv(H)))
uncerts = numpy.sqrt(numpy.diag(numpy.linalg.inv(H)))
if return_FIM:
return H
else:
return uncertsI have changed the FIM_uncert function to return the inverse of the full Hessian, not just it's diagonal.
In [9]:
H = dadi.Godambe.FIM_uncert(func_ex=func_ex, grid_pts=pts_l, p0=popt, data=sfs2d, return_FIM=True)
The returned Hessian matrix is the Fisher Information Matrix CrossValidated.
Thus, the Fisher information may be seen as the curvature of the support curve (the graph of the log-likelihood). Near the maximum likelihood estimate, low Fisher information therefore indicates that the maximum appears "blunt", that is, the maximum is shallow and there are many nearby values with a similar log-likelihood. Conversely, high Fisher information indicates that the maximum is sharp.
from the Wikipedia article on Fisher Information
If your data are all unlinked, you can get estimates of your parameter uncertainties using the Hessian (a.k.a. Fisher Information Matrix), which is the second derivate matrix of the log-likelihood with respect to your parameters. (Essentially you're making a quadratic approximation of the likelihood surface.)
Ryan Gutenkunst on the dadi forum
The inverse (reciprocal) of the FIM is the approximated variance-covariance matrix.
In [10]:
import numpy as np
import pandas as pd
In [11]:
var_covar = np.linalg.inv(H)
In [12]:
df_var_covar = pd.DataFrame(data=var_covar, \
columns=['nu1', 'nu2', 'Tc', 'm1', 'm2', 'Ti', 'theta'], \
index=['nu1', 'nu2', 'Tc', 'm1', 'm2', 'Ti', 'theta'])
df_var_covar
Out[12]:
In [13]:
# standard deviations for the parameters (genetic units)
popt_sd = np.sqrt(np.diag(np.linalg.inv(H)))
popt_sd
Out[13]:
In [14]:
model = func_ex(popt, ns, pts_l)
In [15]:
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
theta
Out[15]:
In [32]:
np.resize?
In [16]:
popt = np.resize(popt, len(popt)+1)
popt
Out[16]:
In [17]:
# add theta to popt
popt[-1] = theta
popt
Out[17]:
In [61]:
round?
In [18]:
print "p_opt SD pSE"
for p, sd, psd in zip(popt, popt_sd, popt_sd/popt*100):
print "{0:.3f} (+-{1:.3f}) (+-{2:2d}%)".format(p, sd, int(round(psd)))
In [19]:
nu1, nu2, Tc, m1, m2, Ti, theta = popt
In [20]:
T_total = Tc + Ti
T_total
Out[20]:
In [35]:
from math import sqrt
In [22]:
sigma_T_total = sqrt( df_var_covar["Tc"]["Tc"]**2 + df_var_covar["Ti"]["Ti"]**2 + 2*df_var_covar["Tc"]["Ti"]**2 )
sigma_T_total
Out[22]:
That is the estimated standard deviation of the split time $T_{total}$ in genetic units (2Nref).
In [23]:
# 95% CI
print T_total - 1.96*sigma_T_total
print T_total + 1.96*sigma_T_total
Now, it would be nice to get the uncertainty for the split time translated to absolute units (generations).
To make things more simple, let's get the standard deviation for $T_c$ first. This should be very close to $T_{total}$.
equation B.12 in Gillespie: Populations genetics - A concise guide. $V()$ is the variance.
equation (4) in:
Gray, G. Covariances in Multiplicative Estimates. Transactions of the American Fisheries Society, 1999, 128, 475-482
Note the sign change for the last two terms as compared to equation (6) in:
Bohrnstedt, G. W. & Goldberger, A. S.: On the Exact Covariance of Products of Random Variables. Journal of the American Statistical Association, 1969, 64, 1439-1442
The sign change makes $\hat{V}(xy)$ an unbiased estimator of $V(xy)$.
$x$ and $y$ are assumed to be bivariate normally distributed. $E(x)$ and $E(y)$ are the maximum likelihood estimates returned by dadi. $C()$ is the covariance.
In [24]:
mu = 3e-9
L = sfs2d.data.sum()
In [25]:
var_theta = df_var_covar['theta']['theta']
var_theta
Out[25]:
In [26]:
theta
Out[26]:
In [27]:
var_Tc = df_var_covar['Tc']['Tc']
var_Tc
Out[27]:
In [28]:
cov_Tc_theta = df_var_covar['Tc']['theta']
cov_Tc_theta
Out[28]:
In [64]:
var_Tc_theta = Tc**2 * var_theta + theta**2 * var_Tc + 2*theta*Tc*cov_Tc_theta - var_Tc*var_theta - cov_Tc_theta**2
var_Tc_theta
Out[64]:
The highly negative $C(T_c, \theta)$ makes this variance negative. Maybe the formula of Gray1999 fails because $T_c$ and $\theta$ are not bivariate normal or the approximation of the covariance via the Hessian is not very good? The assumption of normality is necessary to drop the third moments (e. g. $C(x^2, y^2)$), which are unknown, from equation (5) in Bohrnstedt1969. If I want to calculate confidence intervals, I also have to assume that the parameters are normally distributed.
In [31]:
1/(2*L*mu)**2
Out[31]:
In [65]:
var_Tc_abs = 1/(2*L*mu)**2 * var_Tc_theta
var_Tc_abs
Out[65]:
Instead of the unbiased estimator of the variance of the product of two variables, I am going to use "a conventional asymptotic approximation" (Bohrnstedt1969, equation (14)) and simply ignore the last two terms in $\hat{V}(yx)$ of Gray1999: $V(T_c)V(\theta)$ and $C^2(T_c, \theta)$. This leaves:
$$ \hat{V}(xy) \approx E^2(x) V(y) + E^2(y) V(x) + 2 E(x)E(y) C(x, y) $$Note, that the approximation is conservative, i. e. the approximation will always overestimate the true variance.
In [81]:
var_Tc_theta = Tc**2 * var_theta + theta**2 * var_Tc + 2*theta*Tc*cov_Tc_theta
var_Tc_theta
Out[81]:
In [82]:
var_Tc_abs = 1/(2*L*mu)**2 * var_Tc_theta
var_Tc_abs
Out[82]:
In [83]:
# standard deviation of Tc in generations
int( sqrt(var_Tc_abs) )
Out[83]:
If this is actually close to the true standard deviation of $T_c$, then this suggests that the estimate of $T_c$ in generations is remarkably precise! The optimal $T_c$ in generations was:
In [45]:
Nref = theta/(4*L*mu)
int(Tc * 2 * Nref)
Out[45]:
The estimated correlation coefficient between $T_c$ and $\theta$ is:
In [46]:
corr_Tc_theta = cov_Tc_theta / sqrt(var_theta * var_Tc)
corr_Tc_theta
Out[46]:
$T_c$ and $\theta$ are highly negatively correlated. This explains why the realtively large estimated standard deviation of $T_c$ in genetic units (0.971) translates to such a small standard deviation in absolute units (generations).
Now, what is the estimate for the standard deviation for $T_i$, the time since complete isolation between ERY and PAR?
In [47]:
var_Ti = df_var_covar['Ti']['Ti']
var_Ti
Out[47]:
In [70]:
cov_Ti_theta = df_var_covar['Ti']['theta']
cov_Ti_theta
Out[70]:
In [84]:
var_Ti_abs = 1/(2*L*mu)**2 * (Ti**2 * var_theta + theta**2 * var_Ti + 2*theta*Ti*cov_Ti_theta)
var_Ti_abs
Out[84]:
In [85]:
# standard deviation for Ti in generations
int( sqrt(var_Ti_abs) )
Out[85]:
This again suggests a remarkably precise estimate!
In [51]:
# correlation coefficient between Ti and theta
cov_Ti_theta / sqrt(var_theta * var_Ti)
Out[51]:
$T_i$ and $\theta$ are also highly negatively correlated.
Finally, it would be nice to have an approximate standard deviation for the total time since the split between ERY and PAR.
equation B.12 in Gillespie: Populations genetics - A concise guide
equation B.14 in Gillespie: Populations genetics - A concise guide
equation B.16 in Gillespie: Populations genetics - A concise guide
In [86]:
cov_Tc_Ti = df_var_covar['Tc']['Ti']
cov_Tc_Ti
Out[86]:
In [87]:
var_Tc_plus_Ti = var_Tc + var_Ti + 2 * cov_Tc_Ti
var_T_total_abs = 1/(2*L*mu)**2 * ((Tc+Ti)**2 * var_theta + theta**2 * var_Tc_plus_Ti + 2*theta*(Tc+Ti)*(cov_Tc_theta+cov_Ti_theta))
var_T_total_abs
Out[87]:
In [89]:
# standard deviation of T_total in generations
int(sqrt(var_T_total_abs))
Out[89]:
In [56]:
# correlation coefficient between Tc and Ti
cov_Tc_Ti / sqrt(var_Tc * var_Ti)
Out[56]:
$T_c$ and $T_i$ are highly correlated. That's why the $SD\left(T_{total}^{abs}\right) < SD\left(T_c^{abs}\right) + SD\left(T_i^{abs}\right)$: 34345 < 34012 + 2045.
In [57]:
var_Nref = (1/(4*L*mu))**2 * var_theta
int(sqrt(var_Nref))
Out[57]:
In [58]:
# optimal estimate of Nref
int(theta/(4*L*mu))
Out[58]:
In [93]:
var_nu1 = df_var_covar['nu1']['nu1']
var_nu1
Out[93]:
In [94]:
cov_nu1_theta = df_var_covar['nu1']['theta']
cov_nu1_theta
Out[94]:
In [95]:
var_Nery = (1/(4*L*mu))**2 * (nu1**2 * var_theta + theta**2 * var_nu1 + 2*nu1*theta*cov_nu1_theta)
int(sqrt(var_Nery))
Out[95]:
In [96]:
# correlation coefficient between nu1 and theta
cov_nu1_theta / sqrt(var_nu1 * var_theta)
Out[96]:
The lower variance of $N_{ery}$ as compared to $N_{ref}$ is due to the highly negative correlation between $\nu_1$ and $\theta$.
In [97]:
var_nu2 = df_var_covar['nu2']['nu2']
var_nu2
Out[97]:
In [98]:
cov_nu2_theta = df_var_covar['nu2']['theta']
cov_nu2_theta
Out[98]:
In [99]:
var_Npar = (1/(4*L*mu))**2 * (nu2**2 * var_theta + theta**2 * var_nu2 + 2*nu2*theta*cov_nu2_theta)
int(sqrt(var_Npar))
Out[99]:
In [100]:
# correlation coefficient between nu2 and theta
cov_nu2_theta / sqrt(var_nu2 * var_theta)
Out[100]:
taken from Wikipedia article Propagation of uncertainty and citations therein and this CrossValidated post. I couldn't find a better source so far.
equation (12) in:
Bohrnstedt, G. W. & Goldberger, A. S.: On the Exact Covariance of Products of Random Variables. Journal of the American Statistical Association, 1969, 64, 1439-1442
With the assumption of multivariate normality, the third moment from the original equation could be dropped.
see Bohrnstedt1969 and Gray1999
In [101]:
var_nu2 = df_var_covar['nu2']['nu2']
var_nu2
Out[101]:
In [102]:
var_m1 = df_var_covar['m1']['m1']
var_m1
Out[102]:
In [103]:
cov_nu2_theta = df_var_covar['nu2']['theta']
cov_nu2_theta
Out[103]:
In [104]:
cov_nu2_m1 = df_var_covar['nu2']['m1']
cov_nu2_m1
Out[104]:
In [105]:
cov_theta_m1 = df_var_covar['theta']['m1']
cov_theta_m1
Out[105]:
In [107]:
var_nu2_theta = nu2**2 * var_theta + theta**2 * var_nu2 + 2*nu2*theta*cov_nu2_theta
var_nu2_theta
Out[107]:
In [108]:
cov_theta_nu2_m1 = theta*cov_nu2_m1 + nu2*cov_theta_m1
cov_theta_nu2_m1
Out[108]:
In [122]:
var_m1_over_nu2theta = (m1/(nu2*theta))**2 * (var_m1/m1**2 + var_nu2_theta/(theta*nu2)**2 - 2/(m1*nu2*theta) * cov_theta_nu2_m1)
var_m1_over_nu2theta
Out[122]:
In [123]:
var_p_ery_par = (2*L*mu)**2 * var_m1_over_nu2theta
In [124]:
# standard deviation for the proportion of new gametes from ERY
# in the PAR population each generation
sqrt(var_p_ery_par)
Out[124]:
In [115]:
# MLE of p_ery->par was:
m1/(2*Nref*nu2)
Out[115]:
In [116]:
var_m2 = df_var_covar['m2']['m2']
var_m2
Out[116]:
In [118]:
var_nu1_theta = nu1**2 * var_theta + theta**2 * var_nu1 + 2*nu1*theta*cov_nu1_theta
var_nu1_theta
Out[118]:
In [119]:
cov_nu1_m2 = df_var_covar['nu1']['m2']
cov_nu1_m2
Out[119]:
In [120]:
cov_theta_m2 = df_var_covar['theta']['m2']
cov_theta_m2
Out[120]:
In [121]:
cov_m2_nu1_theta = theta*cov_nu1_m2 + nu1*cov_theta_m2
cov_m2_nu1_theta
Out[121]:
In [125]:
var_m2_over_thetanu1 = (m2/(theta*nu1))**2 * ( var_m2/m2**2 + var_nu1_theta/(nu1*theta)**2 - 2/(m2*theta*nu1) * cov_m2_nu1_theta )
var_m2_over_thetanu1
Out[125]:
In [126]:
var_p_par_ery = (2*L*mu)**2 * var_m2_over_thetanu1
var_p_par_ery
Out[126]:
In [128]:
sqrt(var_p_par_ery)
Out[128]:
In [129]:
m2/(2*Nref*nu1)
Out[129]:
The following table shows the inferred parameters for the ancient migration model from the corrected spectrum in their absolute units together with the estimated standard errors (SD) and relative standard errors SD/MLE in percent:
| parameter | MLE | SD | %SD/MLE |
|---|---|---|---|
| $N_{a}$ | 113,522 | 14,058 | 12% |
| $N_{ery}$ | 389,237 | 4,138 | 1% |
| $N_{par}$ | 950,409 | 7,393 | 0.78% |
| $T_c$ | 1,525,640 | 34,012 | 2.2% |
| $p_{ery->par}$ | 2.14e-08 | 2.68e-09 | 12% |
| $p_{par->ery}$ | 2.60e-07 | 3.14e-08 | 12% |
| $T_i$ | 26,970 | 2,045 | 7.6% |
| $T_c+T_i$ | 1,552,610 | 34,345 | 2.2% |
$N_x$ have unit individuals, $T_x$ has unit generations, $p_x$ are proportions of new immigrant individuals per generation.
The effective population sizes of ERY and PAR seem to be estimated with extremely high precision. The precision for the estimated time parameters also seems to be extremely good. Note that these standard errors are only valid if all SNP's in the spectrum are unlinked. If there is linkage between SNP's then these standard errors will underestimate the true standard errors. Better estimates of variance can only be achieved with bootstraps over contigs and the Godambe Information Matrix. However, acquiring bootstraps over contigs is currently not possible with ANGSD/realSFS.