In [1]:
import numpy
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
import numpy as np
# turn on floating point division by default, old behaviour via '//'
from __future__ import division
In [3]:
# import 2D unfolded spectrum
sfs2d_unfolded = dadi.Spectrum.from_file('dadiExercises/EryPar.unfolded.2dsfs.dadi_format')
# add population labels
sfs2d_unfolded.pop_ids = ["ery", "par"]
In [4]:
ns = sfs2d_unfolded.sample_sizes
pts_l = [40, 50, 60]
In [5]:
# print number of segregating sites in the SFS
sfs2d_unfolded.S()
Out[5]:
The 2D spectrum contains counts from 60k sites that are variable in par or ery or both.
For the estimation of the 2D SFS, realSFS has only taken sites that had data from at least 9 individuals in each population (see assembly.sh, lines 1430 onwards).
In [6]:
import pylab
%matplotlib inline
pylab.rcParams['font.size'] = 14.0
pylab.rcParams['figure.figsize'] = [14.0, 12.0]
dadi.Plotting.plot_single_2d_sfs(sfs2d_unfolded, vmin=1, cmap=pylab.cm.jet)
Out[6]:
Cells with counts below 1 are masked and appear white in the upper plot.
In [7]:
sfs2d = sfs2d_unfolded.fold()
In [8]:
# plot the folded GLOBAL minor allele frequency spectrum
dadi.Plotting.plot_single_2d_sfs(sfs2d, vmin=1, cmap=pylab.cm.jet)
Out[8]:
In [9]:
ns = sfs2d.sample_sizes
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
In [10]:
def split_nomig(params, ns, pts):
"""
params = (nu1,nu2,T)
ns = (n1,n2)
Split into two populations of specifed size, no migration.
nu1: Size of population 1 after split.
nu2: Size of population 2 after split.
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
This model specifies a split of the ancestral population $T \times 2N_{ref}$ generations ago. The two daughter populations have a constant population size until the present of $\nu_1 \times N_{ref}$ and $\nu_2 \times N_{ref}$, respectively.
In [11]:
# load optimisation results from file
ar_split_nomig = []
for filename in glob("OUT_2D_models/split_nomig*dill"):
ar_split_nomig.append(dill.load(open(filename)))
In [12]:
success = [flatten(out)[:7] for out in ar_split_nomig if out[1][4] == 0]
df = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'T_0', 'nu1_opt', 'nu2_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Out[12]:
The optimal parameter values look very robust: very different starting values converge to almost identical optimal parameter values.
In [13]:
# create link
func = split_nomig
# create extrapolating version of the model function
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [14]:
# calculate best-fit model spectrum
model_spectrum = func_ex((0.636945, 1.229606, 0.290465), ns, pts_l)
In [15]:
theta = dadi.Inference.optimal_sfs_scaling(model_spectrum, sfs2d)
In [16]:
L = sfs2d.data.sum() # this sums over all entries in the spectrum, including masked ones, i. e. also contains invariable sites
print "The optimal value of theta per site for the ancestral population is {0:.4f}.".format(theta/L)
In [17]:
mu = 3e-9
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: {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 [18]:
print "The ancestral population of ery and par split apart {0:,} generations ago.".format(int(0.290465*2*N_ref)),
print "Since then until present the ery population had a size of {0:,} and the par population of {1:,}.".format(int(0.636945*N_ref)
, int(1.229606*N_ref))
Note, that the fitting of 1D models to the spectra of each population had indicated a population size reduction for both populations and more so for parallelus, which seems to contradict what this 2D model says about the two populations.
As mentioned above, the population size history inferred from the 2D spectrum contradicts what was inferred from the individual 1D spectra of each population. The 2D spectrum is not computed from the same number of sites as the two 1D spectra. While inclusion of a site in each 1D spectrum required read data from at least 9 individuals, the 2D spectrum required read data from at least 9 individuals in both populations. That is because each population's SAF files were created with minInd 9 (see lines 1426-1435 in assembly.sh). When estimating the 2D SFS with realSFS only sites that occurred in each population's SAF files were included (see lines 1661-1664 in assembly.sh).
In [19]:
# get marginal 1D spectra from the 2D spectrum
fsm_ery = sfs2d.marginalize([1])
fsm_par = sfs2d.marginalize([0])
In [20]:
# get original 1D spectra
# import 1D unfolded spectrum of ery
fs_ery_unfolded = dadi.Spectrum.from_file('ERY.unfolded.sfs')
fs_ery = fs_ery_unfolded.fold()
# import 1D unfolded spectrum of par
fs_par_unfolded = dadi.Spectrum.from_file('PAR.unfolded.sfs')
fs_par = fs_par_unfolded.fold()
In [21]:
pylab.plot(fs_par, 'go-', label='individual', linewidth=2.5, markersize=12)
pylab.plot(fsm_par, 'g>--', label='marginal', linewidth=2.5, markersize=12)
pylab.grid()
pylab.title('PAR')
pylab.legend()
Out[21]:
The individual 1D spectrum for PAR and the 1D spectrum from marginalising over ERY in the 2D spectrum do not look drastically different.
In [22]:
pylab.plot(fs_ery, 'ro-', label='individual', linewidth=2.5, markersize=12)
pylab.plot(fsm_ery, 'r>--', label='marginal', linewidth=2.5, markersize=12)
pylab.grid()
pylab.title('ERY')
pylab.legend()
Out[22]:
The marginal 1D spectrum for ERY looks quite different. It mostly looks shifted down. This makes sense, when a smaller number of sites (those that are present in the SAF files of both ERY and PAR) go into the 2D spectrum than into the individual 1D spectrum.
Since the individual and marginal 1D spectra have similar shape, they should lead to similar demographic inference. The difference in the inference of population sizes for the two populations from 1D and 2D models can therefore not be explained by sytematic differences between 1D and 2D spectra.
In [73]:
dadi.Plotting.plot_2d_comp_multinom(model_spectrum, sfs2d, vmin=1)
The residuals show that the best-fit model has a deficit of fixed divergent SNP's as well as of low frequency shared SNP's. Maybe the fit can be improved by introducing gene flow.
In [23]:
# use dadi's built-in model function
func = dadi.Demographics2D.split_mig
func_ex = dadi.Numerics.make_extrap_log_func(func)
This model adds a single migration rate parameter $m$ . It is assumed that migration is of equal size $m$ in each direction.
In [32]:
# load optimisation results from file
ar_split_mig = []
for filename in glob("OUT_2D_models/split_mig_[0-9]*dill"):
ar_split_mig.append(dill.load(open(filename)))
In [33]:
success = [flatten(out)[:9] for out in ar_split_mig if out[1][4] == 0]
df = pd.DataFrame(data=success, \
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)
Out[33]:
From many different initial parameter value combinations the optimisations converge on essentially the same set of parameters that includes a migration rate that is not 0. This optimal parameter combination is much more likely than the other combinations. From many initial parameter values the optimisations converge to a model that is essentially the same as the split-no-migration model from above. However, these parameter combinations have much lower likelihood.
In [27]:
popt = numpy.array(df.sort_values(by='-logL', ascending=True).iloc[0, 4:8])
# calculate best-fit model spectrum
model_spectrum = func_ex(popt, ns, pts_l)
theta = dadi.Inference.optimal_sfs_scaling(model_spectrum, sfs2d)
print "The optimal value of theta per site for the ancestral population is {0:.4f}.".format(theta/L)
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 individuals) implied by this theta is: {0:,}.".format(int(N_ref))
In [78]:
dadi.Plotting.plot_2d_comp_multinom(model_spectrum, sfs2d, vmin=1)
This has reduced the residuals for low frequency shared polymorphisms along the diagonal, but the model still cannot fit divergently fixed SNP's well.
The dadi manual says about the migration rate parameter (on page 11):
The migration parameter m12 specifies the rate of migration from pop 2 into pop 1. It is equal to the fraction of individuals each generation in pop 1 that are new migrants from pop 2, times the 2Nref.
So the fraction of new immigrant alleles each generation from $pop_2$ in $pop_1$ is:
$$ p_{m12} = \frac{m_{12}}{2N_{ref}} $$Since dadi's migration rate is scaled by 2Nref, the migration rate parameter as reported by dadi has units of immigrant gametes per generation. However, whenever pop 1 does not have the same size as the ancestral population ($N_{ref}$), dadi's $m_{12}$ needs to be rescaled by multiplication with $\nu_1$ to get the number of immigrant gametes per generation ($M_{12}$):
$$ \begin{align} p_{m12} &= \frac{m_{12}}{2N_{ref}\nu_1} \\[5pt] M_{12} &= p_{m12} \times 2N_{ref}\nu_1 \\[5pt] M_{12} &= \nu_1 m_{12} \end{align} $$See this confirm units of dadi thread on the dadi forum.
What this model fit says in plain English:
In [79]:
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 {0:.2e} of the population size of PAR were made up each generation of new immigrant individuals from 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])
Is the model with migration significantly better than the model without migration?
From the dadi manual:
The full likelihood ratio test statistic is equal to $D = 2(ll_c - ll_s)$ where $ll_c$ and $ll_s$ are the likelihoods of the complex and simple model, respectively.
$D$ should be $\chi^2$ distributed
with degrees of freedom equal to the difference in number of parameters between the simple and complex model.
Not all SNP's in the spectrum will be unlinked with each other. Therefore, the likelihoods calculated are composite likelihoods. Dadi can calculate an adjustment factor for $D$ that corrects for linkage. However, this requires bootstrapped data sets, which I haven't created yet. In the following calculation I am assuming that all SNP's can be assumed independent data points.
In [80]:
# calculate logL for split-migration model
ll_c = dadi.Inference.ll_multinom(model_spectrum, sfs2d)
func = split_nomig
# create extrapolating version of the model function
func_ex = dadi.Numerics.make_extrap_log_func(func)
# calculate best-fit model spectrum with optimal parameter values
model_spectrum = func_ex((0.636945, 1.229606, 0.290465), ns, pts_l)
# calculate logL for split-no-migration model
ll_s = dadi.Inference.ll_multinom(model_spectrum, sfs2d)
In [81]:
D = 2*(ll_c - ll_s)
D
Out[81]:
From the dadi example script YRI_CEY.py:
Because this is a test of a parameter on the boundary of the parameter space (m cannot be less than zero), our null distribution is an even proportion of chi^2 distributions with 0 and 1 d.o.f. To evaluate the p-value, we use the point percent function for a weighted sum of chi^2 dists.
See Self2007 for more complex cases.
In [82]:
# calculate p-value for Chi-square dist. with 1 degree of freedom
p = dadi.Godambe.sum_chi2_ppf(D, weights=(0.5, 0.5))
p
Out[82]:
Gene flow significantly improves the fit to the observed spectrum. I doubt that there could be enough linkage in the data to compromise this result.
Phylogeographic studies as well as paleoclimatic data strongly suggest that the two subspecies, erythropus and parallelus, diverged in isolation during at least the last ice age and only recently (~10,000 YBP) came into secondary contact. Clines of several morphological as well as molecular markers have been shown to be very wide, which indicates substantial gene flow across the hybrid zone into adjacent pure populations. The two population samples come from such "pure" populations in close proximity to the hybrid zone. They may therefore show evidence of recent gene flow between the two subspecies.
I am now defining a model function that is a mixture of the split_nomig and split_mig models. It specifies a split, divergence in isolation for a time $T_i$, then secondary contact of the two diverged populations since $T_c$ times $2N_{ref}$ generations in the past.
In [28]:
def secondary_contact(params, ns, pts):
"""
params = (nu1,nu2,Ti,Tc,m)
ns = (n1,n2)
Split into two populations of specified size,
then diverge in isolation,
finally come into secondary contact.
nu1: Size of population 1 after split.
nu2: Size of population 2 after split.
Ti: Time of divergence in isolation (in units of 2*Na generations)
Tc: Time in the past when secondary contact began (and continues till present)
m: symmetrical migration rate (in units of 2*Na per generation)
n1,n2: Sample sizes of resulting Spectrum
pts: Number of grid points to use in integration.
"""
nu1,nu2,Ti,Tc,m = params
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi) # split
# divergence in isolation
phi = dadi.Integration.two_pops(phi, xx, Ti, nu1, nu2, m12=0, m21=0)
# secondary contact
phi = dadi.Integration.two_pops(phi, xx, Tc, nu1, nu2, m12=m, m21=m)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
In [37]:
ar_sec_contact = []
for filename in glob("OUT_2D_models/sec_contact*dill"):
ar_sec_contact.append(dill.load(open(filename)))
In [38]:
# get parameters from all runs into a table
returned = ([flatten(out)[:11] for out in ar_sec_contact])
df = pd.DataFrame(data=returned, \
columns=['nu1_0','nu2_0', 'Ti_0', 'Tc_0', 'm_0', 'nu1_opt', 'nu2_opt', 'Ti_opt', 'Tc_opt', 'm_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Out[38]:
The most likely parameter combination has the same likelihood as the best-fit for the split_mig model (i. e. the model with $T_i = 0$). The best-fit parameters for this secondary contact model suggest only a short time of complete isolation followed by a long time of secondary contact, which is qualitatively equivalent to a model without complete isolation. A $T_i$ greater 0 does not improve the likelihood, so I think it is not significant.
In [34]:
func = secondary_contact
# create extrapolating version of the model function
func_ex = dadi.Numerics.make_extrap_log_func(func)
In [35]:
import numpy as np
In [39]:
# determine theta
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 5:10])
# get the unscaled optimal secondary contact model spectrum
model_spectrum = func_ex(popt, ns, pts_l)
# get theta
theta = dadi.Inference.optimal_sfs_scaling(model_spectrum.fold(), sfs2d)
theta
Out[39]:
In [40]:
dadi.Plotting.plot_2d_comp_multinom(model_spectrum, sfs2d, vmin=1)
In [41]:
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 Nref of ery and par (in number of diploid individuals) implied by this theta is: {0:,}.".format(int(N_ref))
In [42]:
print "After the split of the ancestral population {0:,} generations ago, the daughter populations increased in size (instantaneously).".format(int((popt[2]+popt[3])*2*N_ref)),
print "ERY increased to about the same size as the ancestral population (Nref), while PAR increased to {0:.2f} times of Nref.".format(popt[1]),
print "After the split ERY and PAR diverged in isolation for {0:,} generations.".format(int(popt[2]*2*N_ref)),
print "The two populations came into secondary contact {0:,} generations ago.".format(int(popt[3]*2*N_ref)),
print "ERY and PAR stayed in continuous secondary contact until the present and",
print "during that time exchanged one individual from par to ery every {0:.1f} generations".format(1.0/(popt[4]/2*popt[0])),
print "and one individual from ery into par every {0:.1f} generations.".format(1.0/(popt[4]/2*popt[1]))
I would like to compare expected spectra from a recent secondary contact model with expected spectra from an isolation with continuous migration model (previous section).
In [43]:
success = [flatten(out)[:9] for out in ar_split_mig if out[1][4] == 0]
df_split_mig = pd.DataFrame(data=success, \
columns=['nu1_0','nu2_0', 'T_0', 'm_0', 'nu1_opt', 'nu2_opt', 'T_opt', 'm_opt', '-logL'])
popt_split_mig = np.array(df_split_mig.sort_values(by='-logL', ascending=True).iloc[0, 4:8])
popt_split_mig
Out[43]:
In [44]:
# get split_mig model function (built-in)
func_ex_split_mig = dadi.Numerics.make_extrap_log_func(dadi.Demographics2D.split_mig)
# get optimally scaled split-migration model spectrum
model_spectrum_split_mig = dadi.Inference.optimally_scaled_sfs(func_ex_split_mig(popt_split_mig, ns, pts_l), sfs2d)
In [45]:
# parameter combination that specifies a recent secondary contact model
p_sec_contact = [0.67858642, 1.29659005, 0.35015044, 0.00483189, 1.92454244]
These parameters says the following about the history of PAR and ERY:
In [46]:
print "After the split of the ancestral population {0:,} generations ago, the daughter populations increased in size (instantaneously).".format(int((p_sec_contact[2]+p_sec_contact[3])*2*N_ref)),
print "ERY decreased to about {0:.2f} the size of the ancestral population (Nref), while PAR increased to {1:.2f} times of Nref.".format(p_sec_contact[0],p_sec_contact[1]),
print "After the split ERY and PAR diverged in isolation for {0:,} generations.".format(int(p_sec_contact[2]*2*N_ref)),
print "The two populations came into secondary contact {0:,} generations ago.".format(int(p_sec_contact[3]*2*N_ref)),
print "ERY and PAR stayed in continuous secondary contact until present and",
print "during that time contained a proportion {0:.2e} new immigrant alleles each generations.".format(p_sec_contact[4]/2/N_ref)
In [47]:
# get the optimally scaled secondary contact model spectrum
func_ex = dadi.Numerics.make_extrap_log_func(secondary_contact)
SC_model_spectrum = dadi.Inference.optimally_scaled_sfs(func_ex(p_sec_contact, ns, pts_l), sfs2d)
In [48]:
# get the log likelihood of this RSC parameterisation
ll_RSC = dadi.Inference.ll_multinom(SC_model_spectrum, sfs2d)
ll_RSC
Out[48]:
In [49]:
18574-19711
Out[49]:
This RSC parameterisation reduces the fit to the data by 1137 log likelihood units cmpared to the best-fit parameters which suggest practically no time of complete isolation.
In [50]:
# I have edited the functions `dadi.Plotting.plot_2d_comp_multinom` and `dadi.Plotting.plot_2d_comp_Poisson`
# in `Plotting.py` in the dadi repo to add the argument `title`
dadi.Plotting.plot_2d_comp_multinom(data=SC_model_spectrum.fold(), model=model_spectrum_split_mig.fold(), \
vmin=1, title=['RSC', 'continuous migration'], pop_ids = ['ery', 'par'])
The recent secondary contact model (with the parameters given above) predicts far fewer low frequency shared polymorphisms along the diagonal. Despite the long time of isolation it also predicts fewer highly diverged SNP's (red in the residual plot).
On the other hand, recent secondary contact predicts more high frequency shared polymorphisms than isolation with continuous migration (blue in the residual plot).
Define an IM model that allows for unequal migration rates between the two populations.
In [44]:
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: Size of population 1 after split.
nu2: Size of population 2 after split.
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 [99]:
# load optimal parameters from files
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 [100]:
# create table
l = 2*5+1
success = [flatten(out)[:l] for out in ar_split_asym_mig]
df = pd.DataFrame(data=success, \
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).head(10)
Out[100]:
These parameter combinations consistently indicate a much smaller gene flow from erythropus into parallelus than in the other direction. This seems plausible when considering the observed expansion of the range of parallelus towards the south of the Pyrenees. Shapes of neutral clines across the hybrid zone should also be less steep on the erythropus side than on the parallelus side. A more trivial explanation is that the sampling location for ERY (Greixer) is closer to the cline centre than the sampling location for PAR (Aunat). A quick look at the map does not seem to support this, but migration barriers, also those that do not exist anymore, like forests, will have had an effect on the past gene flow.
In [101]:
# get optimal parameter combination
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 5:10])
popt
Out[101]:
In [39]:
func = split_asym_mig
func_ex = dadi.Numerics.make_extrap_log_func(func)
# calculate best-fit model spectrum
model_spectrum = func_ex(popt, ns, pts_l)
# optimally scale model spectrum
model_spectrum_asym_mig_scaled = dadi.Inference.optimally_scaled_sfs(model_spectrum, sfs2d)
theta = dadi.Inference.optimal_sfs_scaling(model_spectrum, sfs2d)
theta
Out[39]:
In [40]:
dadi.Plotting.plot_2d_comp_multinom(model_spectrum_asym_mig_scaled, sfs2d, vmin=1)
In [50]:
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 optimal value of theta per site for the ancestral population is {0:.5f}.".format(theta/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: {0:,}.".format(int(N_ref))
In [49]:
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])
This model therefore indicates a much stronger introgression from parallelus into erythropus than in the opposite direction. Is the asymmetrical gene flow significant?
In [85]:
ll_c = dadi.Inference.ll_multinom(model_spectrum_asym_mig_scaled, sfs2d)
ll_c
Out[85]:
In [86]:
# get maximum log likelihood of best fit parameter combination for the model
# with symmetrical migration rates
func = dadi.Demographics2D.split_mig
func_ex = dadi.Numerics.make_extrap_log_func(func)
# best-fit paramter combination for the split_mig model (see section 3 above)
popt_split_mig = [0.994185, 1.766127, 0.922632, 0.250688]
model_spectrum_sym_mig = func_ex(popt_split_mig, ns, pts_l)
model_spectrum_sym_mig_scaled = dadi.Inference.optimally_scaled_sfs(model_spectrum_sym_mig, sfs2d)
ll_s = dadi.Inference.ll_multinom(model_spectrum_sym_mig, sfs2d)
ll_s
Out[86]:
In [87]:
D = 2*(ll_c - ll_s)
D
Out[87]:
In [88]:
# calculate p-value for Chi-square dist. with 1 degree of freedom
p = dadi.Godambe.sum_chi2_ppf(D)
p
Out[88]:
It looks as though allowing for asymmetrical gene flow significantly improves the fit to the data. Linkage is unlikely to push the p-value above a significance value of 1e-2, for instance.
In [89]:
ar_split_asym_mig_nom1 = []
for filename in glob("OUT_2D_models/split_asym_mig_nom1*dill"):
ar_split_asym_mig_nom1.append(dill.load(open(filename)))
In [90]:
l = 2*5+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_nom1]
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).head(10)
Out[90]:
In [91]:
ll_s = -18340.9
D = 2*(ll_c - ll_s)
D
Out[91]:
In [92]:
# calculate p-value for Chi-square dist.
# the weights specify a weighted sum of a chi^2 distributions with 1 and 2 d.o.f
#
p = dadi.Godambe.sum_chi2_ppf(D, weights=(0.5, 0.5))
p
Out[92]:
The gene flow from ery into par is highly significant. It is unlikely that this result can be compromised by linkage between the SNP's.
In [93]:
# get maximum log likelihood of best fit parameter combination for the model
# with symmetrical migration rates
func = dadi.Demographics2D.split_mig
func_ex = dadi.Numerics.make_extrap_log_func(func)
# best-fit paramter combination for the split_mig model (see section 2 above)
popt_split_mig = [0.994185, 1.766127, 0.922632, 0.250688]
model_spectrum_sym_mig = func_ex(popt_split_mig, ns, pts_l)
model_spectrum_sym_mig_scaled = dadi.Inference.optimally_scaled_sfs(model_spectrum_sym_mig, sfs2d)
In [94]:
dadi.Plotting.plot_2d_comp_multinom(model=model_spectrum_asym_mig_scaled.fold(), \
data=model_spectrum_sym_mig_scaled.fold(), vmin=1,\
title=['sym mig', 'asym mig'], pop_ids=['ery', 'par'])
The model spectrum for asymmetric gene flow has an excess (brown) of low to medium MAF SNP's in PAR that occur at low frequency in ERY. It has fewer medium to high frequency SNP's that are private to PAR and more medium to high frequency SNP's that are private to ERY. It has far fewer low frequency SNP's that are private to ERY.
In [95]:
popt_split_mig = popt_split_mig + [popt_split_mig[-1]]
popt_split_mig
Out[95]:
In [102]:
popt
Out[102]:
In [103]:
pd.DataFrame(data=np.array([popt_split_mig, popt]), index=['sym_mig', 'asym_mig'], columns=['nu1', 'nu2', 'T', 'm1', 'm2'])
Out[103]:
The excess of SNP's in MAF categroy 1 that did not introgress from PAR into ERY could be due to the greater population size increase inferred for PAR in the asymmetric gene flow model, although this model also inferred a greater time since the population size increase which should lead to greater equilibration of the initial rise in low frequency variants.
I would now like to add a gradual population size change after the split of the ancestral population. So far, the split coincided with a stepwise population size change.
Dadi's built-in IM model function contains what I would like to model now. To avoid confusion of migration rate parameters, I am anyway going to specify a slightly modified version of this model that uses the same naming convention and order of m as used above.
In [6]:
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 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)
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 [105]:
% ll OUT_2D_models/IM*dill
In [106]:
ar_IM = []
for filename in glob("OUT_2D_models/IM*dill"):
ar_IM.append(dill.load((open(filename))))
In [107]:
l = 2*6+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[107]:
This model is hard to fit to the data. The best-fit parameter combinations contain an extremely long divergence time T of around 3, i. e. $6 \times N_{ref}$ generations. No parameter combination with higher likelihood could be found than in the previous runs and they are all much lower than from the split_asym_mig model (-18104), which imposed a stepwise size change coincident with the population split and otherwise constant population sizes.
In [108]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 6:12])
popt
Out[108]:
In [109]:
func_ex = dadi.Numerics.make_extrap_log_func(IM)
model_spectrum = func_ex(popt, ns, pts_l)
In [110]:
dadi.Plotting.plot_2d_comp_multinom(model_spectrum, sfs2d, vmin=1)
In [111]:
# optimally scale the best-fit IM model spectrum,
# by multiplication with optimal theta
model_spectrum_scaled = dadi.Inference.optimally_scaled_sfs(model_spectrum, sfs2d)
In [112]:
# make extraplating function from `split_asym_mig` model function
func_ex = dadi.Numerics.make_extrap_log_func(split_asym_mig)
# these were the optimal parameters for the `split_asym_mig` model
popt_asym_mig = [0.92091098, 2.2386172 , 1.24468148, 0.08414741, 0.49491818]
# get (unscaled) expected SFS for optimal `split_asym_mig` model
model_spectrum_asym_mig = func_ex(popt_asym_mig, ns, pts_l)
# optimally scale spectrum
model_spectrum_asym_mig = dadi.Inference.optimally_scaled_sfs(model_spectrum_asym_mig, sfs2d)
In [113]:
dadi.Plotting.plot_2d_comp_multinom(data=model_spectrum_asym_mig.fold(), model=model_spectrum_scaled.fold(), \
vmin=1, title=['asym_mig', 'IM'], pop_ids=['ery', 'par'])
The IM model with gradual population size changes predicts far fewer low frequency SNP's that are private to PAR and far more high frequency SNP's that are private to ERY.
So far, a model with an ancient split and low but constant gene flow (more from par into ery than in the other direction) has fit the data best. I would like to know whether I can achieve an improvement if I let gene flow cease completely at some point in time.
In [47]:
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: 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
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
This is a piecewise model. I could make gene flow fade out gradually, but that would require four more parameters that specify the initial and final levels of gene flow. Also, it would make the model less comparable to the split_asym_mig model, which is also piecewise.
In [48]:
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 [49]:
l = 2*6+1
# get all optimisation results
returned = [flatten(out)[:l] for out in ar_split_asym_mig_iso]
df = pd.DataFrame(data=returned, \
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[49]:
It was difficult to achieve convergence. Parameter combinations with quite different time of primary contact ($T_C$) have very similar likelihood.
The log likelihood of the split_asym_mig model with best-fit parameters was: -18104.96.
In [50]:
D = 2*(-17582.88 - -18104.96)
D
Out[50]:
In [51]:
# 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 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[51]:
This ancient migration model fits the data significantly better than the split asym mig model. That also means that the recent period of complete isolation ($T_i$) is significant. Assuming that I had found the best-fit parameters for that model, this means that a model that allows for a very ancient split (>1 million generations ago) and that specifies a recent cessation of gene flow fits the data better than a model with a less ancient split and constant migration until the present.
In [52]:
popt_mig_iso = np.array(df.sort_values(by='-logL', ascending=True).iloc[0,6:12])
popt_mig_iso
Out[52]:
In [53]:
# create unscaled model spectrum
func_ex = dadi.Numerics.make_extrap_log_func(split_asym_mig_iso)
model_mig_iso = func_ex(popt_mig_iso, ns, pts_l)
In [54]:
dadi.Plotting.plot_2d_comp_multinom(model=model_mig_iso, data=sfs2d, vmin=1)
In [55]:
theta = dadi.Inference.optimal_sfs_scaling(model_mig_iso, sfs2d)
theta
Out[55]:
In [56]:
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: {0:,}.".format(int(N_ref))
In [57]:
nu1, nu2, Tc, m1, m2, Ti = popt_mig_iso
In [58]:
popt = popt_mig_iso
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))
This is another indication that a recent secondary contact model is not compatible with the data. If the better fit of this ancient migration model (compared to constant migration) is not just due to fitting noise or biases in the data, then it would also mean that there is no detectable gene flow since the last ice age, despite the fact that the two populations are just a few kilometers apart and the clines of several markers are very wide across this transect of the hybrid zone.
In [42]:
model_mig_iso_scaled = dadi.Inference.optimally_scaled_sfs(model_mig_iso, sfs2d)
In [45]:
func_ex = dadi.Numerics.make_extrap_log_func(split_asym_mig)
popt_asym_mig = [0.92091098, 2.2386172 , 1.24468148, 0.08414741, 0.49491818]
model_asym_mig = func_ex(popt_asym_mig, ns, pts_l)
model_asym_mig_scaled = dadi.Inference.optimally_scaled_sfs(model_asym_mig, sfs2d)
In [46]:
dadi.Plotting.plot_2d_comp_multinom(data=model_mig_iso_scaled.fold() , model=model_asym_mig_scaled.fold(), vmin=1, title=['ancient mig', 'continuous mig'], pop_ids=['ery', 'par'])
Apparently the continuous migration model predicts many more variants at frequency 1 in erythropus that are shared with parallelus than the ancient migration model.
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 [59]:
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 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
This model just adds another epoch to the ancient migration model in the previous section. The period without gene flow is now followed by another period where gene flow recommences and lasts until the present. Note that migration rates are assumed to be the same in both epochs of gene flow.
In [60]:
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 [61]:
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[61]:
The best parameter combination from the ancient migration model had a likelihood of 17582.88. So the addition of $T_{SC}$ to the model had not improved the fit significantly. A time of secondary contact cannot be inferred from this spectrum.
In [62]:
popt = np.array(df.sort_values(by='-logL', ascending=True).iloc[1,7:14])
popt
Out[62]:
In [65]:
func_ex = dadi.Numerics.make_extrap_log_func(split_mig_iso_mig)
model = func_ex(popt, ns, pts_l)
theta = dadi.Inference.optimal_sfs_scaling(model, sfs2d)
mu = 3e-9
L = sfs2d.data.sum() # this sums over all entries in the spectrum, including masked ones, i. e. also contains invariable sites
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: {0:,}.".format(int(N_ref))
print "A Tsc of {0:.6f} corresponds to {1:,} generations.".format(popt[-1], int(popt[-1]*2*N_ref))
Such a low time since secondary contact is obviously nonsensical.
The split between ery and par is ancient, maybe >1 Mya. A much higher effective population size is inferred for parallelus than for erythropus. There is a strong signal of gene flow in the data, which is inferred to have been much higher in the direction par-->ery than the other way round. The gene flow has not begun after a long period of complete isolation (i. e. no recent secondary contact). Rather, there seems to be a strong signal of ancient gene flow that recently has ceased completely. A model that allowed for exponential (instead of instantaneous) population size change did not fit the data better.