Simulation-Based Model Validation#

We used SLiM 3 to perform forward simulations as a means to validate the diffusion approximation for collapsed polyploid spectra in dadi. Below we give details on how to run each of the validation scripts. The scripts can be found in the validation/ folder in the GitHub repo. To run multiple replicates, each call to SLiM can be wrapped with a Bash for loop with the loop index variable being passed to the rep parameter, followed by any other parameters needed by the model and then the name of the SLiM script:

Example run for SLiM simulation.#
for i in {1..100}
do
    slim -d "rep=${i}" [other parameters] <script-name>
done

When run, each script generates a single realization of the demographic model and outputs a site frequency spectrum that is the format required by dadi. The corresponding Python code for comparing results from the diffusion approximation in dadi can be found in the supp/SupplementalSfsData/plot_slim_dadi_comp.py script.

SLiM scripts#

autotetraploid_snm.slim#

The standard neutral model for an autotetraploid.

# Parameters: rep
slim -d "rep=1" autotetraploid_snm.slim

autotetraploid_bottleneck.slim#

Simulates a bottleneck to population size nuBot times the original population size.

# Parameters: rep, nuBot
slim -d "rep=1" -d "nuBot=0.5" autotetraploid_bottleneck.slim

allotetraploid_iso.slim#

Simulates an equilibrium population that splits at time T1 in the past before being sampled.

# Parameters: rep, T1
slim -d "rep=1" -d "T1=0.5" allotetraploid_iso.slim

allotetraploid_bottleneck.slim#

Simulates an equilibrium population that splits at time T1 + T2 in the past before being sampled, with a bottleneck of size nuBot occurring at time T2 in the past.

# Parameters: rep, T1, T2, nuBot
slim -d "rep=1" -d "T1=0.5" -d "T2=0.25" -d "nuBot=0.5" allotetraploid_iso.slim

segtetraploid_iso.slim#

Simulates an equilibrium population that splits at time T1 + T2 in the past before being sampled, with homoeologous exchange occurring at a rate of dij starting at time T2 in the past. Here dij is equivalent to the exchange parameter \(e_{i \leftrightarrow j}\) described in the manuscript. The use of iso to describe this model is a bit of a misnomer since the subgenomes are not isolated, but we kept it for naming consistency with the allotetraploid models.

# Parameters: rep, T1, T2, dij
slim -d "rep=1" -d "T1=0.5" -d "T2=0.25" -d "dij=0.001" segtetraploid_iso.slim

segtetraploid_bottleneck.slim#

Simulates an equilibrium population that splits at time T1 + T2 in the past before being sampled, with a bottleneck of size nuBot and homoeologous exchange occurring at a rate of dij starting at time T2 in the past. Here dij is equivalent to the exchange parameter \(e_{i \leftrightarrow j}\) described in the manuscript.

# Parameters: rep, T1, T2, dij, nuBot
slim -d "rep=1" -d "T1=0.5" -d "T2=0.25" -d "dij=0.001" -d "nuBot=0.5" segtetraploid_iso.slim

Inference and plotting scripts#

In the inference/optimization/ folder, there are Python and R scripts to simulate data under the allotetraploid and segmental allotetraploid models, run inference on the simulated data with dadi, and then collate and plot the maximum likelihood parameter estimates across all replicates for each combination of parameters. The subfolders for each model contain the simulated frequency spectra and optimization results (in CSV format) for all combinations of parameters tested.

References

Haller, B. C. and P. W. Messer. 2019. SLiM 3: Forward genetic simulations beyond the Wright–Fisher model. Molecular Biology and Evolution 36:632–-637.