<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Table of Contents
<p><div class="lev1 toc-item"><a href="#Plot-the-data" data-toc-modified-id="Plot-the-data-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Plot the data</a></div><div class="lev2 toc-item"><a href="#To-fold-or-not-to-fold-by-ANGSD" data-toc-modified-id="To-fold-or-not-to-fold-by-ANGSD-11"><span class="toc-item-num">1.1&nbsp;&nbsp;</span>To fold or not to fold by ANGSD</a></div><div class="lev1 toc-item"><a href="#$\theta$-and-implied-$N_{ref}$" data-toc-modified-id="$\theta$-and-implied-$N_{ref}$-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>$\theta$ and implied $N_{ref}$</a></div><div class="lev2 toc-item"><a href="#standard-neutral-model" data-toc-modified-id="standard-neutral-model-21"><span class="toc-item-num">2.1&nbsp;&nbsp;</span>standard neutral model</a></div><div class="lev2 toc-item"><a href="#reducing-noise-and-bias" data-toc-modified-id="reducing-noise-and-bias-22"><span class="toc-item-num">2.2&nbsp;&nbsp;</span>reducing noise and bias</a></div><div class="lev1 toc-item"><a href="#exponential-growth-model" data-toc-modified-id="exponential-growth-model-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>exponential growth model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-31"><span class="toc-item-num">3.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-32"><span class="toc-item-num">3.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#two-epoch-model" data-toc-modified-id="two-epoch-model-4"><span class="toc-item-num">4&nbsp;&nbsp;</span>two epoch model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-41"><span class="toc-item-num">4.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-42"><span class="toc-item-num">4.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#bottleneck-then-exponential-growth-model" data-toc-modified-id="bottleneck-then-exponential-growth-model-5"><span class="toc-item-num">5&nbsp;&nbsp;</span>bottleneck then exponential growth model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-51"><span class="toc-item-num">5.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-52"><span class="toc-item-num">5.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#three-epoch-model" data-toc-modified-id="three-epoch-model-6"><span class="toc-item-num">6&nbsp;&nbsp;</span>three epoch model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-61"><span class="toc-item-num">6.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-62"><span class="toc-item-num">6.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#exponential-growth-then-bottleneck-model" data-toc-modified-id="exponential-growth-then-bottleneck-model-7"><span class="toc-item-num">7&nbsp;&nbsp;</span>exponential growth then bottleneck model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-71"><span class="toc-item-num">7.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-72"><span class="toc-item-num">7.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#Conclusion" data-toc-modified-id="Conclusion-8"><span class="toc-item-num">8&nbsp;&nbsp;</span>Conclusion</a></div><div class="lev1 toc-item"><a href="#Can-gene-flow-explain-the-difficulty-of-fitting-1D-models?" data-toc-modified-id="Can-gene-flow-explain-the-difficulty-of-fitting-1D-models?-9"><span class="toc-item-num">9&nbsp;&nbsp;</span>Can gene flow explain the difficulty of fitting 1D models?</a></div>
Step1:  Plot the data
Step2:  To fold or not to fold by ANGSD
Does estimating an unfolded spectrum with ANGSD and then folding yield a sensible folded SFS when the sites are not polarised with respect to an ancestral allele but with respect to the reference allele? Matteo Fumagalli thinks that this is sensible.
Step3:  For parallelus, the spectrum folded by dadi looks better than the spectrum folded by ANGSD. I am therefore going to use spectra folded in dadi.
Step4:  $\theta$ and implied $N_{ref}$
standard neutral model
Step5:  What effective ancestral population size would that imply?
According to section 5.2 in the dadi manual
Step6:  These effective population sizes are consistent with those reported for other insect species (Lynch2016, fig. 3b).
In dadi, times are given in units of $2N_{ref}$ generations.
reducing noise and bias
Step8:  The lower plot (green line) is for the scaled Poisson residuals. 
$$
residuals = (model - data)/\sqrt{model}
$$
The model is the expected counts in each frequency class. If these counts are Poisson distributed, then their variance is equal to their expectation. The differences between model and data are therefore scaled by the expected standard deviation of the model counts.
Although a standard neutral model is unrealistic for either population, the observed counts deviate by up to 90 standard deviations from this model!
What could be done about this? 
The greatest deviations are seen for the first two frequency classes, the ones that should provide the greatest amount of information for theta (Fu1994) and therefore probably also other parameters. I think masking out just one of the first two frequency classes will lead to highly biased inferences. Masking both frequency classes will reduce a lot of the power to infer the demographic history. I therefore think that masking will most likely not lead to better estimates.
Toni has suggested that the doubleton class is inflated due to "miscalling" heterozygotes as homozygotes. When they contain a singleton they will be "called" as homozygote and therefore contribute to the doubleton count. This is aggravated by the fact that the sequenced individuals are all male which only possess one X chromosome. The X chromosome is the fourth largest of the 9 chromosomes of these grasshoppers (8 autosomes + X) (see Gosalvez1988, fig. 2). That is, about 1/9th of the sequenced RAD loci are haploid but ANGSD assumes all loci to be diploid. The genotype likelihoods it calculates are all referring to diploid genotypes.
I think one potential reason for the extreme deviations is that the genotype likelihoods are generally biased toward homozygote genotypes (i. e. also for autosomal loci) due to PCR duplicates (see eq. 1 in Nielsen2012). So, one potential improvement would be to remove PCR duplicates. 
Another potential improvement could be found by subsampling 8/9th to 8/10th of the contigs in the SAF files and estimating an SFS from these. Given enough subsamples, one should eventually be found that maximally excludes loci from the X chromosome. This subsample is expected to produce the least squared deviations from an expected SFS under the standard neutral model. However, one could argue that this attempt to exclude problematic loci could also inadvertently remove loci that strongly deviate from neutral expectations due to non-neutral evolution, again reducing power to detect deviations from the standard neutral model. I think one could also just apply the selection criterion of the second MAF class to be lower than the first and just save all contig subsamples and SFS's that fulfill that criterioin, since that should be true for all demographic scenarios. This approach, however, is compute intensive and difficult to implement.
Ludovic has suggested optimising for a fraction $p$ of all SNP's that lie on the X chromosome. SNP's on the X chromosome are counted twice in the SFS. This fraction $p$ should therefore be subtracted from even frequency classes and added to the respective frequency class that contains SNP's that are 1/2 as frequent, e. g. from class 2 --> 1 or from 8 --> 4. One could optimise for the minimum deviation from a neutral spectrum.
Finally, it may be possible to acquire the the Locusta genome contigs scaffolded or annotated by linkage group (see Wang2014). One could then blast-annotate the RAD contigs with the homologous Locusta linkage group and create subsets of RAD contigs for each linkage group that exclude that linkage group. There may be a subset (presumably the one that excluded RAD contigs from the X chromosome) that has markedly reduced residuals when compared to a neutral spectrum.
exponential growth model
This model assumes that exponential growth (or decline) started some time $T$ in the past and the current effective population size is a multiple $\nu$ of the ancient populations size, i. e. before exponential growth began. So this model just estimates two parameters. If $\nu$ is estimated to be 1, this indicates that the population size hasn't changed (although see Myers2008). If it is below one, this indicates exponential population decline (how realistic is that?). If it is above 1, this indicates exponential population growth.
erythropus
Step9:  This optimal parameter set was difficult to find. It can only be found from parameter values that are already quite close to the optimal values. The model enforces exponential population decline. This may not be very realistic. 
parallelus
A sweep through the parameter space did not lead to convergence of parameters for the parallelus spectrum. Optimisations starting from up to 6-fold randomly perturbed neutral parameter values (1, 1) also did not lead to successful optimisations. The exponential growth/decline model could therefore not be fit to the par spectrum.
two epoch model
erythropus
The built-in two epoch model defines a piecewise history, where at some time $T$ in the past the ancestral population instantaneously changed in size to the contemporary population size, which has ratio of $\nu$ to the ancestral population size.
The parameters values with the highest likelihood ($\nu=0.0001$ and $T=0.000002$) are hitting the lower bound of parameter space that I've set. The time parameter $T$ that dadi infers with the highest likelihood corresponds to 3 generations and therefore does not make sense. The two epoch model could therefore not be fit to the erythropus spectrum.
parallelus
Step10:  Very different population size histories have the same likelihood. The most likely parameter combinations are hitting the upper parameter bound for the time parameter (6).
Step11:  From the dadi manual
Step12:  These optimal parameter values were found by a broad sweep across the parameter space (data not shown).
The optimal parameters suggest the following
Step13:  Can a shortterm effective population size of 21 million reasonably be assumed?
Step14:  parallelus
There was no convergence on a set of optimal parameter values.
three epoch model
The built-in three_epoch model specifies a piecewise history (with only instantaneous population size changes instead of gradual changes). At time $TF+TB$ the ancient population underwent a size change, stayed at this size ($\nu_B \times N_{ref}$) for $TB \times 2N_{ref}$ generations and then underwent a second size size change $TF \times 2N_{ref}$ generations in the past to the contemporary population size ($\nu_F \times N_{ref}$). The model has therefore two population size parameters, $\nu_B$ and $\nu_F$ as well as two time parameters, $TB$ and $TF$.
With 4 parameters, this model is already so complex that a sweep through the parameter space (i. e. starting optimisations from many different positions) is prohibitively time-consuming. I have therefore started optimisations with initial parameter values that were generated by randomly perturbing a neutral parameter combination, $\nu_B = \nu_F = TB = TF = 1$, by up to 4-fold. I can therefore not be certain whether the optimal parameter combinations thus found correspond to the global or a just a local optimum of the likelihood function.
erythropus
Step15:  Note, all but the last parameter combination have the same likelihood. Fairly different parameter combinations have practically identical likelihood. A reduction of the contemporary population size to 1/4 of the ancient population size is quite conistently inferred ($\nu_F$). The ancient population size change ($\nu_B$) is not inferred consistently. It ranges from a 186-fold increase to a reduction to 1/3 of the ancient population size.
Step16:  This plot visualises the range of stepwise population size histories that the above parameter combinations imply (all with likelihood equal to 2168.11186). Most parameter combinations infer an ancient population size increase followed by a drastic population size collapse to less than 1/3 of the ancient population size that happened more than 400,000 generations ago. 
Are the strength of population size expansion, $(\nu_B)^{TB}$, and the strength population size reduction, $(\frac{1}{\nu_F})^{TF}$, correlated with each other?
Step17:  If a long period of increased population size would be correlated with a long period of decreased population size, this plot should show a positive correlation. This does not seem to be the case.
parallelus
Step18:  As can be seen, extremely different population size histories have the same likelihood.
Step19:  All parameter combinations consistently show a reduction in the contemporary population size with respect to the ancient population size ($\nu_F$). The ancient population size change is much less clear and ranges from a 300-fold increase to 100-fold reduction. Also this event is inferred to have occurred quite distant in the past, generally much more than 1 million generations ago.
exponential growth then bottleneck model
This model specifies exponential growth/decline toward $\nu_B \times N_{ref}$ during a time period of $TB \times 2N_{ref}$ generations, after which (at time $TF \times 2N_{ref}$) the population size undergoes an instantaneous size change to the contemporary size ($\nu_F \times N_{ref}$). 
erythropus
Step20:  Except for the last two, all optimal parameter combinations have the same likelihood and they do not deviate much from their respective initial values. In all successful optimisations above, $\nu_F$, the ratio of contemporary population size to ancient population size, converges to a value below 1/3. I think an ancient size increase or decrease cannot be inferred.
Step21:  There is a clear linear correlation between the size of the inferred contemporary population size ($\nu_F \times N_{ref}$) and the time ($TF$) for which the population has been at this size.
parallelus
Step22:  Often, but not always are optimal parameter values close to their respective initial values. All optimal parameter combinations have the same likelihood. Only a recent (stepwise) population decrease is consistently inferred ($\nu_F$).
Step23:  There could be a linear correlation between the duration of population size reduction and the size of the reduction, but it's weaker than for erythropus.
Conclusion
Not much can be inferred from the 1D spectra of both erythropus and parallelus. However, the dadi optimisations have consistently suggested a population size reduction more than 400,000 generations ago to less than 1/3 of the ancient population size. The reduction seems to have been much stronger in parallelus than in erythropus.
Can gene flow explain the difficulty of fitting 1D models?
The fitting of 2D models was substantially easier than fitting 1D models and led to conistent and robust parameter estimates.
Could it be that the inability to model gene flow from an external source is responsible for the issues with fitting 1D models?
I have saved to files the results of each optimisation run (i. e. for each combination of starting parameter values).
Step24:  These are the successful parameter optimisations I have for the model that specifies a split at time $T$ in the past into two daughter populations with population size ratios (relative to $N_{ref}$) of $\nu_1$ for ery and $\nu_2$ for par. $m$ is the rate of migration, which is assumed to be symmetrical. The migration rate is the fraction of individuals from the other population times $2N_{ref}$, i. e. number of immigrant alleles.
The best paramter combination, that I got from these optimisation runs is
<ASSISTANT_TASK:>
Python Code:
# load dadi module
import sys
sys.path.insert(0, '/home/claudius/Downloads/dadi')
import dadi
%ll
%ll dadiExercises
# import 1D spectrum of ery
fs_ery = dadi.Spectrum.from_file('dadiExercises/ERY.FOLDED.sfs.dadi_format')
# import 1D spectrum of ery
fs_par = dadi.Spectrum.from_file('dadiExercises/PAR.FOLDED.sfs.dadi_format')
Explanation: Table of Contents
<p><div class="lev1 toc-item"><a href="#Plot-the-data" data-toc-modified-id="Plot-the-data-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Plot the data</a></div><div class="lev2 toc-item"><a href="#To-fold-or-not-to-fold-by-ANGSD" data-toc-modified-id="To-fold-or-not-to-fold-by-ANGSD-11"><span class="toc-item-num">1.1&nbsp;&nbsp;</span>To fold or not to fold by ANGSD</a></div><div class="lev1 toc-item"><a href="#$\theta$-and-implied-$N_{ref}$" data-toc-modified-id="$\theta$-and-implied-$N_{ref}$-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>$\theta$ and implied $N_{ref}$</a></div><div class="lev2 toc-item"><a href="#standard-neutral-model" data-toc-modified-id="standard-neutral-model-21"><span class="toc-item-num">2.1&nbsp;&nbsp;</span>standard neutral model</a></div><div class="lev2 toc-item"><a href="#reducing-noise-and-bias" data-toc-modified-id="reducing-noise-and-bias-22"><span class="toc-item-num">2.2&nbsp;&nbsp;</span>reducing noise and bias</a></div><div class="lev1 toc-item"><a href="#exponential-growth-model" data-toc-modified-id="exponential-growth-model-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>exponential growth model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-31"><span class="toc-item-num">3.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-32"><span class="toc-item-num">3.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#two-epoch-model" data-toc-modified-id="two-epoch-model-4"><span class="toc-item-num">4&nbsp;&nbsp;</span>two epoch model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-41"><span class="toc-item-num">4.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-42"><span class="toc-item-num">4.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#bottleneck-then-exponential-growth-model" data-toc-modified-id="bottleneck-then-exponential-growth-model-5"><span class="toc-item-num">5&nbsp;&nbsp;</span>bottleneck then exponential growth model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-51"><span class="toc-item-num">5.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-52"><span class="toc-item-num">5.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#three-epoch-model" data-toc-modified-id="three-epoch-model-6"><span class="toc-item-num">6&nbsp;&nbsp;</span>three epoch model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-61"><span class="toc-item-num">6.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-62"><span class="toc-item-num">6.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#exponential-growth-then-bottleneck-model" data-toc-modified-id="exponential-growth-then-bottleneck-model-7"><span class="toc-item-num">7&nbsp;&nbsp;</span>exponential growth then bottleneck model</a></div><div class="lev2 toc-item"><a href="#erythropus" data-toc-modified-id="erythropus-71"><span class="toc-item-num">7.1&nbsp;&nbsp;</span><em>erythropus</em></a></div><div class="lev2 toc-item"><a href="#parallelus" data-toc-modified-id="parallelus-72"><span class="toc-item-num">7.2&nbsp;&nbsp;</span><em>parallelus</em></a></div><div class="lev1 toc-item"><a href="#Conclusion" data-toc-modified-id="Conclusion-8"><span class="toc-item-num">8&nbsp;&nbsp;</span>Conclusion</a></div><div class="lev1 toc-item"><a href="#Can-gene-flow-explain-the-difficulty-of-fitting-1D-models?" data-toc-modified-id="Can-gene-flow-explain-the-difficulty-of-fitting-1D-models?-9"><span class="toc-item-num">9&nbsp;&nbsp;</span>Can gene flow explain the difficulty of fitting 1D models?</a></div>
End of explanation
import pylab
%matplotlib inline
pylab.rcParams['figure.figsize'] = [12.0, 10.0]
pylab.plot(fs_par, 'go-', label='par')
pylab.plot(fs_ery, 'rs-', label='ery')
pylab.xlabel('minor allele frequency')
pylab.ylabel('SNP count')
pylab.legend()
Explanation: Plot the data
End of explanation
! cat ERY.unfolded.sfs
! cat PAR.unfolded.sfs
# import 1D unfolded spectrum of ery
fs_ery_unfolded = dadi.Spectrum.from_file('ERY.unfolded.sfs')
# import 1D unfolded spectrum of par
fs_par_unfolded = dadi.Spectrum.from_file('PAR.unfolded.sfs')
pylab.plot(fs_ery, 'ro-', label='folded by ANGSD')
pylab.plot(fs_ery_unfolded.fold(), 'bs-', label='folded by ' + r'$\delta$a$\delta$i')
pylab.xlabel('minor allele frequency')
pylab.ylabel('SNP count')
pylab.legend()
pylab.plot(fs_par, 'go-', label='folded by ANGSD')
pylab.plot(fs_par_unfolded.fold(), 'bs-', label='folded by ' + r'$\delta$a$\delta$i')
pylab.xlabel('minor allele frequency')
pylab.ylabel('SNP count')
pylab.legend()
Explanation: To fold or not to fold by ANGSD
Does estimating an unfolded spectrum with ANGSD and then folding yield a sensible folded SFS when the sites are not polarised with respect to an ancestral allele but with respect to the reference allele? Matteo Fumagalli thinks that this is sensible.
End of explanation
fs_ery = fs_ery_unfolded.fold()
fs_par = fs_par_unfolded.fold()
Explanation: For parallelus, the spectrum folded by dadi looks better than the spectrum folded by ANGSD. I am therefore going to use spectra folded in dadi.
End of explanation
# create link to built-in model function
func = dadi.Demographics1D.snm
# make the extrapolating version of the demographic model function
func_ex = dadi.Numerics.make_extrap_log_func(func)
# setting the smallest grid size slightly larger than the largest population sample size
pts_l = [40, 50, 60]
ns = fs_ery.sample_sizes
ns
# calculate unfolded AFS under standard neutral model (up to a scaling factor theta)
neutral_model = func_ex(0, ns, pts_l)
neutral_model
theta_ery = dadi.Inference.optimal_sfs_scaling(neutral_model, fs_ery)
theta_ery
theta_par = dadi.Inference.optimal_sfs_scaling(neutral_model, fs_par)
theta_par
Explanation: $\theta$ and implied $N_{ref}$
standard neutral model
End of explanation
mu = 3e-9
L = fs_ery.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 ery spectrum is {0:,}.".format(int(L))
N_ref_ery = theta_ery/L/mu/4
print "The effective ancestral population size of ery (in number of diploid individuals) implied by this theta is: {0:,}.".format(int(N_ref_ery))
mu = 3e-9
L = fs_par.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 par spectrum is {0:,}.".format(int(L))
N_ref_par = theta_par/L/mu/4
print "The effective ancestral population size of par (in number of diploid individuals) implied by this theta is: {0:,}.".format(int(N_ref_par))
Explanation: What effective ancestral population size would that imply?
According to section 5.2 in the dadi manual:
$$
\theta = 4 N_{ref} \mu_{L} \qquad \text{L: sequence length}
$$
Then
$$
\mu_{L} = \mu_{site} \times L
$$
So
$$
\theta = 4 N_{ref} \mu_{site} \times L
$$
and
$$
N_{ref} = \frac{\theta}{4 \mu_{site} L}
$$
Let's assume the mutation rate per nucleotide site per generation is $3\times 10^{-9}$ (see e. g. Liu2017).
End of explanation
# compare neutral model prediction with ery spectrum
dadi.Plotting.plot_1d_comp_multinom(neutral_model.fold()[:19], fs_ery[:19], residual='linear')
# compare neutral model prediction with par spectrum
dadi.Plotting.plot_1d_comp_multinom(neutral_model.fold()[:19], fs_par[:19], residual='linear')
Explanation: These effective population sizes are consistent with those reported for other insect species (Lynch2016, fig. 3b).
In dadi, times are given in units of $2N_{ref}$ generations.
reducing noise and bias
End of explanation
import dill
# loading optimisation results from previous analysis
growth_ery = dill.load(open('OUT_exp_growth_model/ERY_perturb_ar_ery.dill'))
def flatten(array):
    
        Returns a list of flattened elements of every inner lists (or tuples)
        ****RECURSIVE****
    
    import numpy
    res = []
    for el in array:
        if isinstance(el, (list, tuple, numpy.ndarray)):
            res.extend(flatten(el))
            continue
        res.append(el)
    return list( res )
import pandas as pd
i = 4 # where to find flag, 6 for BFGS, 4 for Nelder-Mead
successfull_popt_ery = [flatten(out)[:5] for out in growth_ery if out[1][i] == 0]
df = pd.DataFrame(data=successfull_popt_ery, \
                  columns=['nu_0', 'T_0', 'nu_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True) # smaller is better
print "The most likely parameter combination for this model suggests that the ery population size started to exponentially shrink at {0:,} generations ago to a current population size of {1:,}.".format(int(0.00769*2*N_ref_ery), int(0.14929*N_ref_ery))
Explanation: The lower plot (green line) is for the scaled Poisson residuals. 
$$
residuals = (model - data)/\sqrt{model}
$$
The model is the expected counts in each frequency class. If these counts are Poisson distributed, then their variance is equal to their expectation. The differences between model and data are therefore scaled by the expected standard deviation of the model counts.
Although a standard neutral model is unrealistic for either population, the observed counts deviate by up to 90 standard deviations from this model!
What could be done about this? 
The greatest deviations are seen for the first two frequency classes, the ones that should provide the greatest amount of information for theta (Fu1994) and therefore probably also other parameters. I think masking out just one of the first two frequency classes will lead to highly biased inferences. Masking both frequency classes will reduce a lot of the power to infer the demographic history. I therefore think that masking will most likely not lead to better estimates.
Toni has suggested that the doubleton class is inflated due to "miscalling" heterozygotes as homozygotes. When they contain a singleton they will be "called" as homozygote and therefore contribute to the doubleton count. This is aggravated by the fact that the sequenced individuals are all male which only possess one X chromosome. The X chromosome is the fourth largest of the 9 chromosomes of these grasshoppers (8 autosomes + X) (see Gosalvez1988, fig. 2). That is, about 1/9th of the sequenced RAD loci are haploid but ANGSD assumes all loci to be diploid. The genotype likelihoods it calculates are all referring to diploid genotypes.
I think one potential reason for the extreme deviations is that the genotype likelihoods are generally biased toward homozygote genotypes (i. e. also for autosomal loci) due to PCR duplicates (see eq. 1 in Nielsen2012). So, one potential improvement would be to remove PCR duplicates. 
Another potential improvement could be found by subsampling 8/9th to 8/10th of the contigs in the SAF files and estimating an SFS from these. Given enough subsamples, one should eventually be found that maximally excludes loci from the X chromosome. This subsample is expected to produce the least squared deviations from an expected SFS under the standard neutral model. However, one could argue that this attempt to exclude problematic loci could also inadvertently remove loci that strongly deviate from neutral expectations due to non-neutral evolution, again reducing power to detect deviations from the standard neutral model. I think one could also just apply the selection criterion of the second MAF class to be lower than the first and just save all contig subsamples and SFS's that fulfill that criterioin, since that should be true for all demographic scenarios. This approach, however, is compute intensive and difficult to implement.
Ludovic has suggested optimising for a fraction $p$ of all SNP's that lie on the X chromosome. SNP's on the X chromosome are counted twice in the SFS. This fraction $p$ should therefore be subtracted from even frequency classes and added to the respective frequency class that contains SNP's that are 1/2 as frequent, e. g. from class 2 --> 1 or from 8 --> 4. One could optimise for the minimum deviation from a neutral spectrum.
Finally, it may be possible to acquire the the Locusta genome contigs scaffolded or annotated by linkage group (see Wang2014). One could then blast-annotate the RAD contigs with the homologous Locusta linkage group and create subsets of RAD contigs for each linkage group that exclude that linkage group. There may be a subset (presumably the one that excluded RAD contigs from the X chromosome) that has markedly reduced residuals when compared to a neutral spectrum.
exponential growth model
This model assumes that exponential growth (or decline) started some time $T$ in the past and the current effective population size is a multiple $\nu$ of the ancient populations size, i. e. before exponential growth began. So this model just estimates two parameters. If $\nu$ is estimated to be 1, this indicates that the population size hasn't changed (although see Myers2008). If it is below one, this indicates exponential population decline (how realistic is that?). If it is above 1, this indicates exponential population growth.
erythropus
End of explanation
ar_par_te = dill.load(open("OUT_two_epoch/PAR_perturb_ar_par.dill"))
i = 4 # index of flag with NM algorithm
successfull_popt = [flatten(out)[:5] for out in ar_par_te if out[1][i] == 0]
df = pd.DataFrame(data=successfull_popt, \
                  columns=['nu_0','T_0', 'nu_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Explanation: This optimal parameter set was difficult to find. It can only be found from parameter values that are already quite close to the optimal values. The model enforces exponential population decline. This may not be very realistic. 
parallelus
A sweep through the parameter space did not lead to convergence of parameters for the parallelus spectrum. Optimisations starting from up to 6-fold randomly perturbed neutral parameter values (1, 1) also did not lead to successful optimisations. The exponential growth/decline model could therefore not be fit to the par spectrum.
two epoch model
erythropus
The built-in two epoch model defines a piecewise history, where at some time $T$ in the past the ancestral population instantaneously changed in size to the contemporary population size, which has ratio of $\nu$ to the ancestral population size.
The parameters values with the highest likelihood ($\nu=0.0001$ and $T=0.000002$) are hitting the lower bound of parameter space that I've set. The time parameter $T$ that dadi infers with the highest likelihood corresponds to 3 generations and therefore does not make sense. The two epoch model could therefore not be fit to the erythropus spectrum.
parallelus
End of explanation
print "A T of 6 corresponds to {0:,} generations.".format(6*2*N_ref_par)
Explanation: Very different population size histories have the same likelihood. The most likely parameter combinations are hitting the upper parameter bound for the time parameter (6).
End of explanation
ar_ery = dill.load(open("OUT_bottlegrowth/ERY_perturb_ar_ery.dill"))
successfull_popt = [flatten(out)[:7] for out in ar_ery if out[1][4] == 0]
df = pd.DataFrame(data=successfull_popt, columns=['nuB_0', 'nuF_0', 'T_0', 'nuB_opt', 'nuF_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=False)
Explanation: From the dadi manual:
If your fits often push the bounds of your parameter space (i.e., results are often at the bounds of one or more parameters), this indicates a problem. It may be that your bounds are too conservative, so try widening them. It may also be that your model is misspecified or that there are unaccounted biases in your data.
Ryan Gutenkunst in a dadi forum thread on parameters hitting the boundary:
This indicates that dadi is having trouble fitting your data. One possibility is that the history of the population includes important events that aren’t in your models. Another possibility is that your data is biased in ways that aren’t in your models. For example, maybe your missing calls for rare alleles.
I am not sure whether an even higher $T$ can reasonably be assumed. The expected time to the most recent common ancestor (MRCA) in a neutral genealogy is:
$$
E\big[T_{MRCA}\big] = 2 \Big(1-\frac{1}{n}\Big)
$$
measured in $N_e$ generations (for coalescent time unit see p. 66, 92 in Wakeley2009). Note, that $T_{MRCA}$ is close to its large sample size limit of 2 already for moderate sample sizes. In order to put dadi's time parameter on the coalescent time scale it needs to be multiplied by 2! A T of 6 from dadi therefore corresponds to a T of 12 on the coalescent time scale, which should be completely beyond the possible $T_{MRCA}$ of any genealogy.
See figure 3.4, p. 79, in Wakeley2009 for the distribution of the $T_{MRCA}$. For $n=36$ it has a mode close to 1.2 and an expected value of 1.94. Values for $T_{MRCA}$ greater than 4 are very unlikely given a standard coalescent model, but may be more likely under models including population expansion or gene flow from another population.
I think a two epoch model cannot be fit to the parallelus spectrum.
bottleneck then exponential growth model
The built-in bottlegrowth model specifies an instantaneous size change $T \times 2N_{ref}$ generations in the past immediately followed by the start of exponential growth (or decline) toward the contemporary population size. The model has three parameters: 
- ratio of population size ($\nu_B$) after instantaneous size change with respect to the ancient population size ($N_{ref}$)
- time of instantaneous size change ($T$) in $2N_{ref}$ generations in the past
- ratio of contemporary to ancient population size ($\nu_F$)
erythropus
End of explanation
print "At time {0:,} generations ago, the ery population size instantaneously increased by almost 40-fold (to {1:,}).".format(int(0.36794*2*N_ref_ery), int(39.2*N_ref_ery))
Explanation: These optimal parameter values were found by a broad sweep across the parameter space (data not shown).
The optimal parameters suggest the following:
End of explanation
print "This was followed by exponential decline towards a current population size of {0:,}.".format(int(0.4576*N_ref_ery))
Explanation: Can a shortterm effective population size of 21 million reasonably be assumed?
End of explanation
# load optimisation results from file
ar_ery = []
import glob
for filename in glob.glob("OUT_three_epoch/ERY_perturb_ar_ery*.dill"):
    ar_ery.extend(dill.load(open(filename)))
# extract successful (i. e. converged) optimisation
successfull_popt = [flatten(out)[:9] for out in ar_ery if out[1][4] == 0]
# create data frame
df = pd.DataFrame(data=successfull_popt, \
                  columns=['nuB_0', 'nuF_0', 'TB_0', 'TF_0', 'nuB_opt', 'nuF_opt', 'TB_opt', 'TF_opt', '-logL'])
# sort data frame by negative log likelihood
df.sort_values(by='-logL', ascending=True) # smaller is better
Explanation: parallelus
There was no convergence on a set of optimal parameter values.
three epoch model
The built-in three_epoch model specifies a piecewise history (with only instantaneous population size changes instead of gradual changes). At time $TF+TB$ the ancient population underwent a size change, stayed at this size ($\nu_B \times N_{ref}$) for $TB \times 2N_{ref}$ generations and then underwent a second size size change $TF \times 2N_{ref}$ generations in the past to the contemporary population size ($\nu_F \times N_{ref}$). The model has therefore two population size parameters, $\nu_B$ and $\nu_F$ as well as two time parameters, $TB$ and $TF$.
With 4 parameters, this model is already so complex that a sweep through the parameter space (i. e. starting optimisations from many different positions) is prohibitively time-consuming. I have therefore started optimisations with initial parameter values that were generated by randomly perturbing a neutral parameter combination, $\nu_B = \nu_F = TB = TF = 1$, by up to 4-fold. I can therefore not be certain whether the optimal parameter combinations thus found correspond to the global or a just a local optimum of the likelihood function.
erythropus
End of explanation
# remove one parameter combination with slightly lower logL than the others
df = df.sort_values(by='-logL').head(-1)
# the time of the ancient pop size change is TB+TF
df['TB+TF'] = pd.Series(df['TB_opt']+df['TF_opt'])
# extract columns from table
nuB = df.loc[:,'nuB_opt']
nuF = df.loc[:, 'nuF_opt']
Tb_Tf = df.loc[:, 'TB+TF']
TF = df.loc[:, 'TF_opt']
# turn nu (a ratio) into absolute Ne and T into generations
nuB_n = nuB*N_ref_ery
nuF_n = nuF*N_ref_ery
Tb_Tf_g = Tb_Tf*2*N_ref_ery
TF_g = TF*2*N_ref_ery
anc = [N_ref_ery] * len(nuB) # ancestral pop size
pres = [1] * len(nuB) # 1 generation in the past
past = [max(Tb_Tf_g)+1000] * len(nuB) # furthest time point in the past
pylab.rcParams['figure.figsize'] = [12.0, 10.0]
pylab.rcParams['font.size'] = 12.0
# plot the best 21 parameter combinations
for x, y in zip(zip(past, Tb_Tf_g, Tb_Tf_g, TF_g, TF_g, pres), zip(anc, anc, nuB_n, nuB_n, nuF_n, nuF_n)):
    pylab.semilogy(x, y)
pylab.xlabel('generations in the past')
pylab.ylabel('effective population size')
Explanation: Note, all but the last parameter combination have the same likelihood. Fairly different parameter combinations have practically identical likelihood. A reduction of the contemporary population size to 1/4 of the ancient population size is quite conistently inferred ($\nu_F$). The ancient population size change ($\nu_B$) is not inferred consistently. It ranges from a 186-fold increase to a reduction to 1/3 of the ancient population size.
End of explanation
pylab.plot(df['nuB_opt']**df['TB_opt'], (1.0/df['nuF_opt'])**df['TF_opt'], 'bo')
pylab.xlabel('strength of first size change')
pylab.ylabel('strength of second size change')
Explanation: This plot visualises the range of stepwise population size histories that the above parameter combinations imply (all with likelihood equal to 2168.11186). Most parameter combinations infer an ancient population size increase followed by a drastic population size collapse to less than 1/3 of the ancient population size that happened more than 400,000 generations ago. 
Are the strength of population size expansion, $(\nu_B)^{TB}$, and the strength population size reduction, $(\frac{1}{\nu_F})^{TF}$, correlated with each other?
End of explanation
import dill
ar_par = dill.load(open("OUT_three_epoch/PAR_perturb_ar_par.dill"))
ar_par_extreme = dill.load(open("OUT_three_epoch/PAR_perturb_extreme_ar_par.dill"))
# add new output to previous output
successfull_popt = [flatten(out)[:9] for out in ar_par if out[1][4] == 0]
successfull_popt.extend([flatten(out)[:9] for out in ar_par_extreme if out[1][4] == 0])
# create data frame
df = pd.DataFrame(data=successfull_popt, \
                  columns=['nuB_0', 'nuF_0', 'TB_0', 'TF_0', 'nuB_opt', 'nuF_opt', 'TB_opt', 'TF_opt', '-logL'])
# sort data frame by negative log likelihood
df.sort_values(by='nuB_opt', ascending=True)
Explanation: If a long period of increased population size would be correlated with a long period of decreased population size, this plot should show a positive correlation. This does not seem to be the case.
parallelus
End of explanation
# add time for ancient size change
df['TB+TF'] = pd.Series(df['TB_opt']+df['TF_opt'])
# extract columns from table
nuB = df.loc[:,'nuB_opt']
nuF = df.loc[:, 'nuF_opt']
Tb_Tf = df.loc[:, 'TB+TF']
TF = df.loc[:, 'TF_opt']
# turn nu into absolute Ne and T into generations
nuB_n = nuB*N_ref_par
nuF_n = nuF*N_ref_par
Tb_Tf_g = Tb_Tf*2*N_ref_par
TF_g = TF*2*N_ref_par
# auxilliary
anc = [N_ref_par] * len(nuB)
pres = [1] * len(nuB)
past = [max(Tb_Tf_g)+1000] * len(nuB)
pylab.rcParams['figure.figsize'] = [12.0, 10.0]
pylab.rcParams['font.size'] = 12.0
for x, y in zip(zip(past, Tb_Tf_g, Tb_Tf_g, TF_g, TF_g, pres), zip(anc, anc, nuB_n, nuB_n, nuF_n, nuF_n)):
    pylab.semilogy(x, y)
pylab.xlabel('generations in the past')
pylab.ylabel('effective population size')
Explanation: As can be seen, extremely different population size histories have the same likelihood.
End of explanation
ar_ery = dill.load(open("OUT_expGrowth_bottleneck/ERY_perturb_ar_ery.dill"))
success = [flatten(out)[:9] for out in ar_ery if out[1][4] == 0]
df = pd.DataFrame(data=success, \
                  columns=['nuB_0','TB_0', 'nuF_0', 'TF_0', 'nuB_opt', 'TB_opt', 'nuF_opt', 'TF_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Explanation: All parameter combinations consistently show a reduction in the contemporary population size with respect to the ancient population size ($\nu_F$). The ancient population size change is much less clear and ranges from a 300-fold increase to 100-fold reduction. Also this event is inferred to have occurred quite distant in the past, generally much more than 1 million generations ago.
exponential growth then bottleneck model
This model specifies exponential growth/decline toward $\nu_B \times N_{ref}$ during a time period of $TB \times 2N_{ref}$ generations, after which (at time $TF \times 2N_{ref}$) the population size undergoes an instantaneous size change to the contemporary size ($\nu_F \times N_{ref}$). 
erythropus
End of explanation
F = df[df['-logL'] < 2172].loc[:,['nuF_opt', 'TF_opt']]
pylab.plot(F['TF_opt']*2*N_ref_ery, F['nuF_opt']*N_ref_ery, 'bo')
pylab.xlabel('generations')
pylab.ylabel('contemporary pop. size')
Explanation: Except for the last two, all optimal parameter combinations have the same likelihood and they do not deviate much from their respective initial values. In all successful optimisations above, $\nu_F$, the ratio of contemporary population size to ancient population size, converges to a value below 1/3. I think an ancient size increase or decrease cannot be inferred.
End of explanation
ar_par = dill.load(open("OUT_expGrowth_bottleneck/PAR_perturb_ar_ery.dill"))
success = [flatten(out)[:9] for out in ar_par if out[1][4] == 0]
df = pd.DataFrame(data=success, \
                  columns=['nuB_0','TB_0', 'nuF_0', 'TF_0', 'nuB_opt', 'TB_opt', 'nuF_opt', 'TF_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Explanation: There is a clear linear correlation between the size of the inferred contemporary population size ($\nu_F \times N_{ref}$) and the time ($TF$) for which the population has been at this size.
parallelus
End of explanation
F = df[df['TF_opt'] < 3].loc[:,['nuF_opt', 'TF_opt']]
pylab.plot(F['TF_opt']*2*N_ref_par, F['nuF_opt']*N_ref_par, 'bo')
pylab.xlabel('TF generations')
pylab.ylabel('contemporary pop. size')
Explanation: Often, but not always are optimal parameter values close to their respective initial values. All optimal parameter combinations have the same likelihood. Only a recent (stepwise) population decrease is consistently inferred ($\nu_F$).
End of explanation
# load results of fitting isolation with migration model to 2D sepctrum
import dill, glob
split_mig_res = []
for filename in glob.glob("OUT_2D_models/split_mig*dill"):
    split_mig_res.append(dill.load(open(filename)))
from utility_functions import *
import pandas as pd
success = [flatten(out)[:9] for out in split_mig_res 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)
Explanation: There could be a linear correlation between the duration of population size reduction and the size of the reduction, but it's weaker than for erythropus.
Conclusion
Not much can be inferred from the 1D spectra of both erythropus and parallelus. However, the dadi optimisations have consistently suggested a population size reduction more than 400,000 generations ago to less than 1/3 of the ancient population size. The reduction seems to have been much stronger in parallelus than in erythropus.
Can gene flow explain the difficulty of fitting 1D models?
The fitting of 2D models was substantially easier than fitting 1D models and led to conistent and robust parameter estimates.
Could it be that the inability to model gene flow from an external source is responsible for the issues with fitting 1D models?
I have saved to files the results of each optimisation run (i. e. for each combination of starting parameter values).
End of explanation
popt = df.sort_values(by='-logL', ascending=True).iloc[0, 4:8]
popt
import numpy as np
popt = np.array(popt)
func = dadi.Demographics2D.split_mig # divergence-with-gene-flow model function (built-in)
func_ex = dadi.Numerics.make_extrap_log_func(func)
ns = np.array([36, 36])
pts_l = [40, 50, 60]
# get optimal model spectrum from 2D model 'divergence-with-migration' model
split_mig_2D_best_fit_model_spectrum = func_ex(popt, ns, pts_l)
split_mig_2D_best_fit_model_spectrum
# 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"]
# fold the spectrum
sfs2d = sfs2d_unfolded.fold()
# need to scale model spectrum by optimal theta, which depends on the number of sites in the data
model_spectrum = dadi.Inference.optimally_scaled_sfs(split_mig_2D_best_fit_model_spectrum, sfs2d)
model_spectrum
import pylab
%matplotlib inline
pylab.rcParams['figure.figsize'] = [12, 10]
pylab.rcParams['font.size'] = 14
model_spectrum.pop_ids = sfs2d.pop_ids
dadi.Plotting.plot_single_2d_sfs(model_spectrum.fold(), vmin=1, cmap=pylab.cm.jet)
# get the marginal model spectra for ery and par from the optimal 2D model spectrum
fs_ery_model = model_spectrum.marginalize([1])
fs_par_model = model_spectrum.marginalize([0])
pylab.plot(fs_par_model.fold(), 'go-', label='par')
pylab.plot(fs_ery_model.fold(), 'rs-', label='ery')
pylab.xlabel('minor allele frequency')
pylab.ylabel('SNP count')
pylab.legend()
# import 1D unfolded spectrum of ery
fs_ery_unfolded = dadi.Spectrum.from_file('ERY.unfolded.sfs')
# import 1D unfolded spectrum of par
fs_par_unfolded = dadi.Spectrum.from_file('PAR.unfolded.sfs')
fs_ery = fs_ery_unfolded.fold()
fs_par = fs_par_unfolded.fold()
dadi.Plotting.plot_1d_comp_multinom(fs_ery_model.fold()[:19], fs_ery[:19])
dadi.Plotting.plot_1d_comp_multinom(fs_par_model.fold()[:19], fs_par[:19])
# log likelihoods
ll_model_ery = dadi.Inference.ll_multinom(fs_ery_model.fold(), fs_ery)
ll_model_ery
ll_model_par = dadi.Inference.ll_multinom(fs_par_model.fold(), fs_par)
ll_model_par
Explanation: These are the successful parameter optimisations I have for the model that specifies a split at time $T$ in the past into two daughter populations with population size ratios (relative to $N_{ref}$) of $\nu_1$ for ery and $\nu_2$ for par. $m$ is the rate of migration, which is assumed to be symmetrical. The migration rate is the fraction of individuals from the other population times $2N_{ref}$, i. e. number of immigrant alleles.
The best paramter combination, that I got from these optimisation runs is:
End of explanation
<END_TASK>