Skip to main content
Access keys NCBI Homepage MyNCBI Homepage Main Content Main Navigation
eLife. 2020; 9: e54967.
Published online 2020 Jun 23. doi: 10.7554/eLife.54967
PMCID: PMC7438115
PMID: 32573438

A community-maintained standard library of population genetic models

Jeffrey R Adrion,#1,†§ Christopher B Cole,#2,†§ Noah Dukler,#3,†§ Jared G Galloway,#1,†§ Ariella L Gladstein,#4,†§ Graham Gower,#5,†§ Christopher C Kyriazis,#6,†§ Aaron P Ragsdale,#7,†§ Georgia Tsambos,#8,†§ Franz Baumdicker,9 Jedidiah Carlson,10 Reed A Cartwright,11 Arun Durvasula,12 Ilan Gronau,13 Bernard Y Kim,14 Patrick McKenzie,15 Philipp W Messer,16 Ekaterina Noskova,17 Diego Ortega-Del Vecchyo,18 Fernando Racimo,5 Travis J Struck,19 Simon Gravel,#7,‡# Ryan N Gutenkunst,#19,‡# Kirk E Lohmueller,#6,12,‡# Peter L Ralph,#1,20,‡# Daniel R Schrider,#4,‡# Adam Siepel,#3,‡# Jerome Kelleher,#21,‡ and Andrew D Kern#1,‡
Graham Coop, Reviewing Editor and Patricia J Wittkopp, Senior Editor
Graham Coop, University of California, Davis, United States;
Author information Article notes Copyright and License information PMC Disclaimer

Associated Data

Supplementary Materials
Data Availability Statement

Abstract

The explosion in population genomic data demands ever more complex modes of analysis, and increasingly, these analyses depend on sophisticated simulations. Recent advances in population genetic simulation have made it possible to simulate large and complex models, but specifying such models for a particular simulation engine remains a difficult and error-prone task. Computational genetics researchers currently re-implement simulation models independently, leading to inconsistency and duplication of effort. This situation presents a major barrier to empirical researchers seeking to use simulations for power analyses of upcoming studies or sanity checks on existing genomic data. Population genetics, as a field, also lacks standard benchmarks by which new tools for inference might be measured. Here, we describe a new resource, stdpopsim, that attempts to rectify this situation. Stdpopsim is a community-driven open source project, which provides easy access to a growing catalog of published simulation models from a range of organisms and supports multiple simulation engine backends. This resource is available as a well-documented python library with a simple command-line interface. We share some examples demonstrating how stdpopsim can be used to systematically compare demographic inference methods, and we encourage a broader community of developers to contribute to this growing resource.

Research organism: Human

Introduction

While population genetics has always used statistical methods to make inferences from data, the degree of sophistication of the questions, models, data, and computational approaches used have all increased over the past two decades. Currently, there exist a myriad of computational methods that can infer the histories of populations (Gutenkunst et al., 2009; Li and Durbin, 2011; Excoffier et al., 2013; Schiffels and Durbin, 2014; Terhorst et al., 2017; Ragsdale and Gravel, 2019), the distribution of fitness effects (Boyko et al., 2008; Kim et al., 2017; Tataru et al., 2017; Fortier et al., 2019; Huang and Siepel, 2019; Vecchyo et al., 2019), recombination rates (McVean et al., 2004; Chan et al., 2012; Lin et al., 2013; Adrion et al., 2020; V Barroso et al., 2019), and the extent of positive selection in genome sequence data (Kim and Stephan, 2002; Eyre-Walker and Keightley, 2009; Alachiotis et al., 2012; Garud et al., 2015; DeGiorgio et al., 2016; Kern and Schrider, 2018; Sugden et al., 2018). While these methods have undoubtedly increased our understanding of genetic and evolutionary processes, very little has been done to systematically benchmark the quality of these inferences or their robustness to deviations from their underlying assumptions. As large databases of population genetic variation begin to be used to inform public health procedures, the accuracy and quality of these inferences is becoming ever more important.

Assessing the accuracy of inference methods for population genetics is challenging in large part because the ‘ground-truth’ in question generally comes not from direct empirical observations, as the relevant historical processes can rarely be observed, but instead from simulations. Population genetic simulations are therefore critically important to the field, yet there has been no systematic attempt to establish community standards or best practices for executing them. Instead, the general modus operandi to date has been for individual groups to validate their own methods using simulations coded from scratch. Often these simulations are more useful to showcase a novel method than to rigorously compare it with competing methods. Moreover, this situation results in a great deal of duplicated effort, and contributes to decreased reproducibility and transparency across the entire field. It is also a barrier to entry to the field, because new researchers can struggle with the many steps involved in implementing a state-of-the-art population genetics simulation, including identifying appropriate demographic models from the literature, translating them into input for a simulator, and choosing appropriate values for key population genetic parameters, such as the mutation and recombination rates.

A related issue is that it has been challenging to assess the degree to which modeling assumptions and choices of data summaries can affect population genetic inferences. Standardized simulations would enable these questions to be systematically examined. Importantly, there are clear examples of different methods yielding fundamentally different conclusions. For example, Markovian coalescent methods applied to human genomes have suggested large ancient (> 100,000 years ago) ancestral population sizes and bottlenecks that have not been detected by other methods based on allele frequency spectra (see Beichman et al., 2017). These distinct methods differ in how they model, summarize, and optimize fit to genetic variation data, suggesting that such design choices can greatly affect the performance of the inference. Furthermore, some methods are likely to perform better than others under certain scenarios, but researchers lack principled guidelines for selecting the best method for addressing their particular questions. The need for guidance from simulated data will only increase as researchers seek to apply population genetic methods to a growing collection of non-model taxa.

For these reasons, we have generated a standardized, community-driven resource for simulating published demographic models from a number of popular study systems. This resource, which we call stdpopsim, makes running realistic simulations for population genetic analysis a simple matter of choosing pre-implemented models from a community-maintained catalog. The stdpopsim catalog currently contains six species: humans, Pongo abelii, Canis familiaris, Drosophila melanogaster, Arabidopsis thaliana, and Escherichia coli. For each species, the catalog contains curated information on our current understanding of the physical organization of its genome, inferred genetic maps, population-level parameters (e.g. mutation rate and generation time estimates), and published demographic models. These models and parameters are meant to represent the field’s current understanding, and we intend for this resource to evolve as new results become available, and other existing models are added to stdpopsim by the community. We have implemented both a command line interface and a simple Python API that can be used to simulate genomic data from a choice of organism, genetic map, chromosome, and demographic history. In this way, stdpopsim will lower the barrier to high-quality simulation for exploratory analyses, enable rigorous evaluation of population genetic software, and contribute to increased reliability of population genetic inferences.

The stdpopsim library has been developed by the PopSim Consortium using a distributed open source model, with strong procedures in place to continue its growth and maintain quality. Importantly, we developed rigorous quality control methods to ensure that we have correctly implemented the models as described in their original publication and provided documented methods for others to contribute new models. We invite new collaborators to join our community: those interested should visit our developer documentation at https://stdpopsim.readthedocs.io/en/latest/development.html. Below we describe the resource and give examples of how it can be used to benchmark demographic inference methods.

Results

The stdpopsim library is a community-maintained collection of empirical genome data and population genetics simulation models, illustrated in Figure 1. The package (https://github.com/popsim-consortium/stdpopsim) centers on a catalog of genomic information and demographic models for a growing list of species (Figure 1A), and software resources to facilitate efficient simulations (Figure 1B–C). Given the genome data and simulation model descriptions defined within the library, it is straightforward to run standardized simulations across a range of organisms. Stdpopsim has a Python API and a user-friendly command line interface, allowing users with minimal experience direct access to state-of-the-art simulations. Simulations are output in the ‘succinct tree sequence’ format (Kelleher et al., 2016; Kelleher et al., 2018; Kelleher et al., 2019), which contains complete genealogical information about the simulated samples, is extremely compact, and can be processed efficiently using the tskit library (Kelleher et al., 2016; Kelleher et al., 2018). The tree sequence format could also be converted to other formats (e.g., VCF) by the user if desired.

An external file that holds a picture, illustration, etc. Object name is elife-54967-fig1.jpg
Structure of stdpopsim.

(A) The hierarchical organization of the stdpopsim catalog contains all model simulation information within individual species (expanded information shown here for H. sapiens only). Each species is associated with a representation of the physical genome, and one or more genetic maps and demographic models. Dotted lines indicate that only a subset of these categories is shown. At right we show example code to specify and simulate models using (B) the python API or (C) the command line interface.

The species catalog

The central feature of stdpopsim is the species catalog, a systematic organization of the key quantitative data needed to simulate a given species. Data are currently available for humans, P. abelii, C. familiaris, D. melanogaster, A. thaliana, and E. coli. A species definition consists of two key elements. Firstly, the library defines some basic information about our current understanding of each species’ genome, including information about chromosome lengths, average mutation rate estimates, and generation times. We also provide access to detailed empirical information such as inferred genetic maps, which model observed heterogeneity in recombination rate along chromosomes. Such maps are often large, so we do not distribute them directly with the software, but make them available for download in a standard format. When a simulation using such a map is requested by the user, stdpopsim will transparently download the map data into a local cache, where it can be quickly retrieved for subsequent simulations. In the initial version of stdpopsim, we support the HapMapII (Frazer et al., 2007) and deCODE (Kong et al., 2010) genetic maps for humans; the Nater et al., 2017 maps for P. abelii; the Campbell et al., 2016 map for C. familiaris; the Salomé et al., 2012 map for A. thaliana; and the Comeron et al., 2012 map for D. melanogaster. Adding further maps to the library is straightforward. The second key element of a species description within stdpopsim is a set of carefully curated population genetic model descriptions from the literature, which allow simulation under specific historical scenarios that have been fit to present-day patterns of genetic variation (see the Materials and methods for a description of the community development and quality-control process for these models.)

The current demographic models in the stdpopsim catalog are shown in Table 1. Homo sapiens currently has the richest selection of population models. These include: a simplified version of the Tennessen et al., 2012 model with only the African population specified (expansion from the ancestral population and recent growth; Africa_1T12); the three-population model of Gutenkunst et al., 2009, which specifies the out-of-Africa bottleneck as well as the subsequent divergence of the European and Asian populations (OutOfAfrica_3G09); the Tennessen et al., 2012 two-population variant of the Gutenkunst et al. model, which does not include Asian populations but more explicitly models recent rapid human population growth in Europe (OutOfAfrica_2T12); the Browning et al., 2018 admixture model for American populations, which specifies ancestral African, European, and Asian population components (AmericanAdmixture_4B11); a three-population out-of-Africa model from Ragsdale and Gravel, 2019, which includes archaic admixture (OutOfAfricaArchaicAdmixture_5R19); a complex model of ancient Eurasian admixture from Kamm et al., 2019 (AncientEurasia_9K19); and a synthetic model of oscillating population size from Schiffels and Durbin, 2014 (Zigzag_1S14).

Table 1.

Initial set of demographic models in the catalog and summary of computing resources needed for simulation.

For each model, we report the CPU time, maximum memory usage and the size of the output tskit file, as simulated using the msprime simulation engine (version 0.7.4). In each case, we simulate 100 samples drawn from the first population, for the shortest chromosome of that species and a constant chromosome-specific recombination rate. The times reported are for a single run on an Intel i5-7600K CPU. Computing resources required will vary widely depending on sample sizes, chromosome length, recombination rates and other factors.

Model IDCitationCPU(s)Ram(MB)File(MB)
HomSap (Homo sapiens)
 Africa_1T12Tennessen et al., 201210.0194.223.3
 Zigzag_1S14Schiffels and Durbin, 20143.3106.17.9
 AshkSub_7G19Gladstein and Hammer, 201913.8216.326.4
 OutOfAfrica_3G09Gutenkunst et al., 200910.2182.021.1
 OutOfAfrica_2T12Tennessen et al., 201210.7198.424.1
 AncientEurasia_9K19Kamm et al., 201963.1304.441.2
 AmericanAdmixture_4B11Browning et al., 201810.6188.122.3
 PapuansOutOfAfrica_10J19Jacobs et al., 2019204.5524.777.8
 OutOfAfricaArchaicAdmixture_5R19Ragsdale and Gravel, 20198.8185.421.7
DroMel (Drosophila melanogaster)
 OutOfAfrica_2L06Li and Stephan, 2006252.8678.0106.7
 African3Epoch_1S16Sheehan and Song, 20163.0123.911.5
AraTha (Arabidopsis thaliana)
 African2Epoch_1H18Huber et al., 20184.3220.516.5
 African3Epoch_1H18Huber et al., 20182.6241.318.4
PonAbe (Pongo abelii) 
 TwoSpecies_2L11Locke et al., 20117.2171.914.7

For D. melanogaster, we have implemented the three-epoch model estimated by Sheehan and Song, 2016 from an African sample (African3Epoch_1S16), as well as the out-of-Africa divergence and associated bottleneck model of Li and Stephan, 2006, which jointly models African and European populations (OutOfAfrica_2L06). For A. thaliana, we implemented the model in Durvasula et al., 2017 inferred using MSMC. This model includes a continuous change in population size over time, rather than pre-specified epochs of different population sizes (SouthMiddleAtlas_1D17). We have also implemented a two-epoch and a three-epoch model estimated from African samples of A. thaliana in Huber et al., 2018 (African2Epoch_1H18 and African3Epoch_1H18).

In addition to organism-specific models, stdpopsim also includes a generic piecewise constant size model and isolation with migration (IM) model which can be used with any genome and genetic map. Together, these models contain many features believed to affect observed patterns of polymorphism (e.g. bottlenecks, population growth, admixture) and therefore provide useful benchmarks for method development.

To guarantee reproducibility, we have standardized naming conventions for species, genetic maps, and demographic models that will enable long-term stability of unique identifiers used throughout stdpopsim, as described in our documentation (https://stdpopsim.readthedocs.io/en/latest/development.html#naming-conventions).

Simulation engines

Currently, stdpopsim uses the msprime coalescent simulator (Kelleher et al., 2016) as the default simulation engine. Coalescent simulations, while highly efficient, are limited in their ability to model continuous geography or complex selection scenarios, such as recurrent sweeps and background selection. For these reasons, we have also implemented the forward-time simulator, SLiM (Haller and Messer, 2019; Haller and Messer, 2019), as an alternative backend engine to stdpopsim, allowing for the simulation of processes that cannot be modeled under the coalescent. However, as forward-time simulators explicitly model all individuals in a population, simulating large population sizes can be highly demanding of computational resources. One common practice used to address this challenge is to simulate a smaller population, but to rescale resulting times, mutation rates, recombination rates, and selection coefficients so that the intensity of mutation, recombination, and allele frequency change due to selection per unit time remains the same (see the SLiM manual and Uricchio and Hernandez, 2014). Our implementation of the SLiM backend allows easy use of this rescaling through a single ‘scaling factor’ argument. Such down-scaled simulations are not completely equivalent to simulating all individuals in the population, and may lead to subtle differences, especially in the presence of selection. However, since many sequence-based measures of population diversity remain nearly unchanged when rescaling in this fashion, this practice is effective for many purposes and widely employed.

We validated our implementation of the SLiM engine by comparing estimates of several population genetic summary statistics for neutral simulations generated by both SLiM and msprime. Examples of this validation for the AncientEurasia_9K19 model (Kamm et al., 2019) are shown in Appendix 1—figure 1 and Appendix 1—figure 2. For this model, down-scaling factors of up to 10 produce patterns of both diversity and linkage disequilibrium that are indistinguishable from those observed under the coalescent (i.e. msprime). Scaling down by a factor of 50 does appear to modify the distribution of these sequence statistics. Interestingly, the apparent difference between distributions is somewhat larger when simulating using a uniform recombination rate (Appendix 1—figure 2), likely due to the lower variation in the values of these statistics. Importantly, both comparisons validate the equivalence of SLiM and msprime when no down-scaling is applied. The results are also optimistic about the rescaling strategy to reduce computational burden, but the possible effects are not well-understood, so results relying on rescaled simulations should be carefully validated.

Documentation and reproducibility

The stdpopsim command-line interface, by default, outputs citation information for the models, genetic maps, and simulation engines used in any particular run. We hope that this feature will encourage users to appropriately acknowledge the resources used in published work, and encourage authors publishing demographic models to contribute to our ongoing community-driven development process. Together with the stdpopsim version number and the long-term stable identifiers for population models and genetic maps, this citation information will result in well-documented and reproducible simulation workflows. The individual tree sequence files produced by stdpopsim also contain complete provenance information including the command line arguments, operating system environment and versions of key libraries used.

Use case: comparing methods of demographic inference

As an example of the utility of stdpopsim, we demonstrate how it can be easily used to perform a fair comparison of popular demographic inference methods. Although we present comparison of results from several methods, our aim at this stage is not to provide an exhaustive evaluation or ranking of these methods. Our hope is instead to demonstrate how stdpopsim will facilitate more detailed future explorations of the strengths and weaknesses of the numerous inference methods that are available to the population genetics community (see Discussion).

We start by comparing popular methods for estimating population size histories of single populations and subsequently show simple examples of multi-population inference. To reproducibly evaluate and compare the performance of inference methods, we developed workflows using snakemake (Köster and Rahmann, 2012), available from https://github.com/popsim-consortium/analysis, that allow efficient computing in multicore or cluster environments. Our workflow generates R replicates of C chromosomes, producing n population samples in each of a total of R×C simulations for each demographic model. After simulation, the workflow prepares input files for each inference method by grouping all n×R×C simulated chromosomes into a single file. Each file is then converted into an input file appropriate for each inference method (such that all inference methods run on the same simulation replicates). Each of the inference programs are then run in parallel, and finally, estimates of population size history from each program are plotted.

Single-population demographic models

For single-population demographic models, we compared MSMC (Schiffels and Durbin, 2014), SMC++ (Terhorst et al., 2017), and stairway plot (Liu and Fu, 2015) on simulated genomes sampled from a single population, under several of the demographic models described above. However, these experiments raise the question of what to use as the ‘true’ population sizes in the case of multi-population models with migration. In particular, a simple single-population model that is fit to data simulated under a multi-population model, is not expected to recover the actual simulated population sizes because of model misspecification. Instead, we argue that the best one may expect in such a scenario is to infer a model that accurately reflects the coalescence time distribution of the simulated model. Under a multi-population model, the coalescence time distribution is influenced by migration between the target population and populations not analyzed in inference, as well as by the ancestral effective population sizes. The inverse coalescence rate is commonly interpreted as the effective population size, since these are equal in a single-population model with random mating. We thus analytically computed inverse coalescence rates in msprime for each simulated model, and used them as benchmarks for the ‘true’ effective population sizes. See the Appendix for a precise definition and description of the inverse coalescence rate computation.

Figure 2 presents the results from simulations under OutOfAfricaArchaicAdmixture_5R19, a model of human migration out of Africa that includes archaic admixture (Ragsdale and Gravel, 2019), along with an empirical genetic map.In each column of this figure we show the inferred population size history (denoted N⁢(t)) from samples taken from each of the three extant populations in the model. In each row we show comparisons among the methods (including two sample sizes for MSMC). Blue lines show estimates from each of three replicate whole genome simulations, and black lines indicate the ‘true’ values depicted by the inverse coalescence rates (although in this specific model the inverse coalescence rates are very close to the simulated population sizes; Appendix 1—figure 3). While there is variation in accuracy among methods, populations, and individual replicates, the methods generally produce a good estimate of the true effective population sizes of the simulations, with inferred values mostly within a factor of two of the truth, and most methods inferring a bottleneck at approximately the correct time.

An external file that holds a picture, illustration, etc. Object name is elife-54967-fig2.jpg
Comparing estimates of N⁢(t) in humans.

Here we show estimates of population size over time (N⁢(t)) inferred using four different methods: smc++, stairway plot, and MSMC with n=2 and n=8 samples. Data were generated by simulating replicate human genomes under the OutOfAfricaArchaicAdmixture_5R19 model (Ragsdale and Gravel, 2019) and using the HapMapII_GRCh37 genetic map (Frazer et al., 2007). From top to bottom, we show estimates for each of the three populations in the model (YRI, CEU, and CHB). In shades of blue we show the estimated N⁢(t) trajectories for each of three replicates. As a proxy for the ‘truth’, in black we show inverse coalescence rates as calculated from the demographic model used for simulation (see text).

Using stdpopsim, we can readily compare performance on this benchmark to that based on a different model of human history. In Appendix 1—figure 4, we show estimates of N⁢(t) from simulations using the same physical and genetic maps, but from the OutOfAfrica_3G09 demographic model that does not include archaic admixture. Again we see that each of the methods is capturing relevant parts of the population history, although the accuracy varies across time. In comparing inferences between the models it is interesting to note that N⁢(t) estimates for the CHB and CEU simulated populations are generally better across methods than estimates from the YRI simulated population.

We can also see how well methods might do at recovering the population history of a constant-sized population, with human genome architecture and genetic map. We show results of such an experiment in Appendix 1—figure 5. All methods recover population size within a factor of two of the simulated values; however, SMC-based methods tend to infer sinusoidal patterns of population size even though no such change is present.

As most method development for population genetics has been focused on human data, it is important to ask how such methods might perform in non-human genomes. Figure 3 shows parameter estimates from the African3Epoch_1S16 model, originally estimated from an African sample of D. melanogaster (Sheehan and Song, 2016), and Appendix 1—figure 6 shows estimates from simulations of A. thaliana under the African2Epoch_1H18 model originally inferred by Huber et al., 2018. In both cases, as with humans, we use stdpopsim to simulate replicate genomes using an empirically-derived genetic map, and try to infer back parameters of the simulation model. Accuracy is mixed among methods when doing inference on simulated data from these D. melanogaster and A. thaliana models, and generally worse than what we observe for simulations of the human genome.

An external file that holds a picture, illustration, etc. Object name is elife-54967-fig3.jpg
Comparing estimates of N⁢(t) in Drosophila.

Population size over time (N⁢(t)) estimated from an African population sample. Data were generated by simulating replicate D. melanogaster genomes under the African3Epoch_1S16 model (Sheehan and Song, 2016) with the genetic map of Comeron et al., 2012. In shades of blue we show the estimated N⁢(t) trajectories for each replicate. As a proxy for the ‘truth’, in black we show inverse coalescence rates as calculated from the demographic model used for simulation (see text).

Multi-population demographic models

As stdpopsim implements multi-population demographic models, we also explored parameter estimation of population divergence parameters. In particular, we simulated data under multi-population models for humans and D. melanogaster and then inferred parameters using ∂⁡a⁢∂⁡i, fastsimcoal2, and smc++. For simplicity, we conducted inference in ∂⁡a⁢∂⁡i and fastsimcoal2 by fitting an isolation with migration (IM) model with constant population sizes and bi-directional migration (Hey and Nielsen, 2004). Our motivation for fitting this simple IM model was to mimic the typical approach of two population inference on empirical data, where the user is not aware of the ‘true’ underlying demography and the inference model is often misspecified. For human models with more than two populations (e.g. Gutenkunst et al., 2009) this limitation means that users are inferring parameters for a model that does not match the model from which the data were generated (Figure 4A and B). However, since the model used for inference also allows gene flow between populations, we directly compare estimated effective population sizes to the values used in simulations (black line in Figure 4C) and not the inverse coalescence rates.

An external file that holds a picture, illustration, etc. Object name is elife-54967-fig4.jpg
Parameters estimated using a multi-population human model.

Here we show estimates of N⁢(t) inferred using ∂⁡a⁢∂⁡i, fastsimcoal2, and smc++. (A) Data were generated by simulating replicate human genomes under the OutOfAfrica_3G09 model and using the HapMapII_GRCh37 genetic map inferred in Frazer et al., 2007. (B) For ∂⁡a⁢∂⁡i and fastsimcoal2 we show parameters inferred by fitting the depicted IM model, which includes population sizes, migration rates, and a split time between CEU and YRI samples. (C) Population size estimates for each population (rows) from ∂⁡a⁢∂⁡i, fastsimcoal2, and smc++ (columns). In shades of blue we show N⁢(t) trajectories estimated from each simulation, and in black simulated population sizes for the respective population. The population split time, TD⁢I⁢V, is shown at the bottom (simulated value in black and inferred values in blue), with a common x-axis to the population size panels.

In Figure 4C, we show estimates of population sizes and divergence time, for each of the inference methods, using samples drawn from African and European populations simulated under the OutOfAfrica_3G09 model. Our results highlight many of the strengths and weaknesses of the different methods. For instance, the SFS-based approaches with simple IM models do not capture recent exponential growth in the CEU population, but do consistently recover the simulated YRI population size history. Moreover, these approaches allow migration rates to be estimated (Appendix 1—figure 7), and lead to more accurate inferences of divergence times. However, these migration rate estimates are somewhat biased. In contrast, smc++ is much better at capturing the recent exponential growth in the CEU population, though it consistently underestimates divergence times because it assumes no migration between populations (Figure 4C).

Again, we can extend this analysis to other taxa and examine the performance of these methods for a two-population model of D. melanogaster. Appendix 1—figure 8 shows inference results using data simulated under the OutOfAfrica_2L06 model. This model includes an ancestral population in Africa from which a European population splits off following a bottleneck, with no post-divergence gene flow between the African and European population (Appendix 1—figure 8A). Here again, we find that ∂⁡a⁢∂⁡i and fastsimcoal2 infer more consistent histories, but they do not detect the brief bottleneck in Europe, due to the inference model not allowing for population size changes after the population split. In addition, ∂⁡a⁢∂⁡i and fastsimcoal2 both do reasonably well at correctly inferring the absence of migration (Appendix 1—figure 9). In contrast, the inferred demographic parameters from smc++ are more noisy, though in some cases better capture the short bottleneck in the European population.

Although these results do not represent an exhaustive benchmarking, we have begun to highlight some of the strengths and weaknesses of these methods. Future work should build on these results and undertake more in-depth comparisons under a wider range of simulated demographic models.

Discussion

Here, we have described the first major product from the PopSim Consortium: the stdpopsim library. We have founded the Consortium with a number of specific goals in mind: standardization of simulation within the population genetics community, increased reproducibility and ease of use of complex simulations, community-based development and decision making guiding best practices in population genetics, and benchmarking of inference methods.

The stdpopsim library allows for rigorous standardization of complex population genetic simulations. Population genetics, as a field, has yet to coalesce around a set of standards for the crucial task of method evaluation, which in our discipline hinges on simulation. In contrast, other fields such as structural biology (Moult et al., 1995) and machine learning Russakovsky et al., 2015 have a long track record of standardized method testing. We hope that our efforts represent the beginning of what will prove to be an equally longstanding and valuable tradition in population genetics.

Besides being a resource for developers of computational methods, we aim for stdpopsim to be a resource for empirical researchers using genomic data. For instance, stdpopsim could be used in power analyses to determine adequate sample sizes, or in sanity checks to see if observed data (e.g. levels of divergence or the allele frequency spectrum) are roughly consistent with the hypothesized scenario. Currently, many studies would benefit from such simulation-based checks. However, there are major barriers to implementation, since individual research groups must reimplement complex, previously published demographic models, a task made especially daunting by additional layers of realism (e.g. recombination maps).

Benchmarking population size inference

We have illustrated in this paper how stdpopsim can be used for direct comparisons of inferential methods on a common set of simulations. Our benchmarking comparisons have been limited, but nevertheless reveal some informative features. For example, at the task of estimating population size histories for simulated human populations, we find that the sequence-based methods (MSMC and smc++) perform somewhat better overall—at least for moderate times in the past—than the site frequency spectrum-based method (stairway plot), which tends to over-estimate the sizes of oscillations (Figure 2 and Appendix 1—figure 4). In contrast, stairway plot outperforms the sequence-based methods on simulations of D. melanogaster or A. thaliana populations, in which linkage disequilibrium is reduced (Figure 3 and Appendix 1—figure 6). In simulations of two human populations (Figure 4), ∂⁡a⁢∂⁡i and fastsimcoal2 do reasonably well at reconstructing the simulated YRI history and estimating divergence times, but struggle with the more complex simulated CEU history, in large part because the methods assume constant population sizes. On the other hand, smc++ does not have the same restrictions on its inferred history, and as a result does much better with the CEU history but tends to underestimate divergence times due to the assumption of no migration. The results for the two-population D. melanogaster model (Appendix 1—figure 8) are generally similar. In these comparisons, fastsimcoal2 and ∂⁡a⁢∂⁡i perform almost identically, which is expected because they fit the same models to the same summaries of the data, differing only in how they calculate model expectations and optimize parameters.

All methods for inferring demographic history have strengths and weaknesses (as recently reviewed by Beichman et al., 2018). We compared inferences from simulated whole genome data, but many factors affect choice of methodology. Markovian coalescent methods (MSMC and smc++) require long contiguous stretches of sequence data. In contrast, frequency spectrum methods (stairway plot, ∂⁡a⁢∂⁡i, and fastsimcoal2) can use reduced-representation sequencing data, such as RADseq (Andrews et al., 2016). ∂⁡a⁢∂⁡i and fastsimcoal2 require a pre-specified parametric model, unlike MSMC, smc++, and stairway plot. Using a parametric approach yields less noisy results, but a model that is too simple may not capture important demographic events (Figure 4 and Appendix 1—figure 8), and other forms of model misspecification may also produce undesirable behavior. From a software engineering perspective, methods also differ in their ease of installation and use. We hope our workflows will assist in the application of all the methods we have considered.

Altogether, these preliminary experiments highlight the utility of stdpopsim for comparing a variety of inference methods on the same footing, under a variety of different demographic models. In addition, the ability of stdpopsim to generate data with and without significant features, such as a genetic map or population-size changes (e.g., Appendix 1—figure 5), allows investigation of the failure modes of popular methods. Moreover the comparison of methods across the various genome organizations, genetic maps, and demographic histories of different organisms, provides valuable information about how methods might perform on non-human systems. Finally, comparison of results across methods or simulation runs provides an estimate of inference uncertainty, analogous to parametric bootstrapping, especially when different methods are vulnerable to model misspecification in different ways.

Next steps

Stdpopsim is intended to be a fully open, community-developed project. Our implementations of genome representations and genetic maps for the some of the most common study systems in computational genetics—humans, Drosophila, and Arabidopsis (among others)—are only intended to be a starting point for future development. Researchers are invited to contribute to the resource by adding their organisms and models of choice. The stdpopsim resource is accompanied by clearly documented standard operating procedures that are intended to minimize barriers to entry for new developers. In this way, we expect the resource to expand and adapt to meet the evolving needs of the population genomics community.

One of our goals is to engage research communities studying other taxa, so as to expand the resource to many more species. Although we have included demographic models and recombination maps, there are many biological processes that we do not model. Some of the additions that we are enthusiastic to add are: selection (including distributions of fitness effects, maps of functional elements, both single and recurrent hitchhiking events, and selection on polygenic traits), gene conversion, mutation models (rate heterogeneity), more realistic demography (overlapping generations, separate sexes, mortality/fecundity schedules), geographic population structure, and downstream aspects of data quality (genotyping and mapping error). Moreover, an in-depth investigation into the effects of population-size rescaling under many of the above scenarios is warranted, given our preliminary findings using neutral simulations (Appendix 1—figures 1 and ​and2).2). Some other important processes are more challenging to model with current simulation software, such as structural variation, changing recombination maps over time, transposable elements, and context-dependent mutation.

We wish to emphasize that although the included demographic histories are some of the most widely used models for our current set of species, we anticipate the set of available models to expand as new methods and new modeling frameworks are developed. For instance, the current models all describe a small set of discrete, randomly mating populations, which are likely good approximations for deep-time population history, but may be less useful for methods describing dynamics of contemporary populations. Stdpopsim’s framework is sufficiently general that more realistic population models will be easily incorporated, as they are published. Additional aspects of the framework, such as genome builds, will also continue to change as improvements are made to our understanding of genome structure.

Materials and methods

Model quality control

As a consortium we have agreed to a standardized procedure for model inclusion into stdpopsim that allows for rigorous quality control. Imagine Developer A wants to introduce a new model into stdpopsim. Developer A implements the demographic model for the relevant organism along with clear documentation of the model parameters and populations. This model is submitted as a ‘pull request’, where it is evaluated by a reviewer and then included as ‘preliminary’, but is not linked to the online documentation nor the command line interface. Developer A submits a quality control (QC) issue, after which a second developer, Developer B (perhaps found by requesting review from the broader Consortium), then independently reimplements the model from the relevant primary sources and adds an automatic unit test for equality between the QC implementation and the preliminary production model. If the two implementations are equivalent, the original model is included in stdpopsim. If not, we move to an arbitration process whereby A and B first try to work out the details of what went wrong. If that fails, the original authors of the published model must be contacted to resolve ambiguities. Further details of our QC process can be found in our developer documentation (https://stdpopsim.readthedocs.io/en/latest/development.html).

The possibility for error and the importance of careful qualty control was illustrated very clearly during our own development process: while carrying out the final revisions of this paper, we noticed that the OutOfAfrica_3G09 model (Gutenkunst et al., 2009) had not gone through our QC process. The subsequent QC revealed that our implementation was in fact slightly wrong—migration rates had not been set to zero to the European population in the most ancient time period when there should have only been a single population. This error was propagated from the msprime documentation, where the model was presented as an illustrative example. A number of studies have been published using copies of this erroneous example code.

Workflow for analysis of simulated data

To demonstrate the utility of stdpopsim we created Snakemake workflows (Köster and Rahmann, 2012) that perform demographic inference on tree sequence output from our package using a few common software packages (see Appendix 1—figure 10 for an example workflow). Our choice of Snakemake allows complete reproducibility of the analyses shown, and all code is available from https://github.com/popsim-consortium/analysis.

We performed two types of demographic inference. Our first task was to infer effective population size over time (denoted N⁢(t)). This was done using three software packages: stairway plot, which uses site frequency spectrum information only (Liu and Fu, 2015); MSMC (Schiffels and Durbin, 2014), which is based on the sequentially Markovian coalescent (SMC), run with two different sample sizes (n=2,8); and smc++ (Terhorst et al., 2017), which combines information from the site frequency spectrum with recombination information as in SMC-based methods. No attempt was made at trying to optimize the analysis from any particular software package, as our goal was not to benchmark performance of methods but instead show how such benchmarking could be easily done using the stdpopsim resource. In this spirit, we ran each software package as near to default parameters as possible. For stairway plot, we set the parameters numRuns = 1 and dimFactor = 5000. For smc++ we used the ‘estimate’ run mode to infer N⁢(t) with all other parameters set to their default values. For MSMC, we used the --fixedRecombination option and used the default number of iterations.

For the single-population task, we ran human (HomSap) simulations using a variety of models (see Table 1): OutOfAfricaArchaicAdmixture_5R19, OutOfAfrica_3G09, and a constant-sized generic model. Each simulation used the HapmapII_GRCh37 genetic map. For D. melanogaster we estimated N⁢(t) from an African sample simulated under the DroMel, African3Epoch_1S16 model using the Comeron2012_dm6 map. Finally, we ran simulations of A. thaliana genomes using the AraTha African2Epoch_1H18 model under the Salome2012_TAIR7 map. For each model, three replicate whole genomes were simulated and the population size estimated from those data. In all cases, we set the sample size of the focal population to N=50 chromosomes.

Following simulation, low-recombination portions of chromosomes were masked from the analysis in a manner that reflects the ‘accessible’ subset of sites used in empirical population genomic studies (e.g. Danecek et al., 2011; Langley et al., 2012). Specifically we masked all regions of 1 cM or greater in the lowest 5th percentile of the empirical distribution of recombination, regions which are nearly uniformly absent for empirical analysis. This approach to masking was chosen to prevent marginal trees with low or no recombination from biasing the comparisons of demographic inference methods. It should be noted that masking is not implemented within stdpopsim proper; tree sequences generated by stdpopsim are always raw and unmasked. This allows users the flexibility to implement masking approaches that are specific to their needs for downstream analysis.

Our second task was to explore inference with two-population models using some of the multi-population demographic models implemented in stdpopsim. For HomSap, we used the OutOfAfrica_3G09 model with the HapmapII_GRCh37 genetic map, and for DroMel we used the OutOfAfrica_2L06 model with the Comeron2012_dm6 map. The HomSap model is a three population model (Africa, Europe, and Asia) including post-divergence migration and exponential growth (Figure 4C), whereas the DroMel model is a two population model (Africa and Europe) with no post-divergence migration and constant population sizes (Appendix 1—figure 8).

To conduct inference on these models, we applied three commonly used methods: ∂⁡a⁢∂⁡i(Gutenkunst et al., 2009), fastsimcoal2 (Excoffier et al., 2013), and smc++ (Terhorst et al., 2017). As above, these methods were used generally with default settings and we did not attempt to optimize their performance or fit parameter-rich demographic models.

For both ∂⁡a⁢∂⁡i and fastsimcoal2, we fit a two population isolation-with-migration (IM) model with constant population sizes. This IM model contains six parameters: the ancestral population size, the sizes of each population after the split, the divergence time, and two migration rate parameters. Importantly, this meant that for both species, the fitted model did not match the simulated model (Figure 4 and Appendix 1—figure 8). In the HomSap case, we therefore performed inference solely on the Africa and Europe populations, meaning that the Asia population functioned as a ‘ghost’ population that was ignored by our inference. To validate our inference approach, we also conducted inference on a generic IM model that was identical to the model used for inference (Appendix 1—figure 11).

From HomSap simulations, we took 20 whole genome samples each from the Europe and Africa populations from each replicate. Runtimes of DroMel simulations were prohibitively slow when simulating whole genomes with the Comeron2012_dm6 map due to large effective population sizes leading to high effective recombination rates. For this reason, we present only data from 50 samples of a 3 MB region of chromosome 2R from simulations under OutOfAfrica_2L06. For the generic IM simulations, we used the HomSap genome along with the HapmapII_GRCh37 genetic map and sampled 20 individuals from each population.

Following simulation, we output tree sequences and masked low-recombination regions using the same approach described for the single population workflow above. We converted tree sequences into a two-dimensional site frequency spectrum for all chromosomes in the appropriate format for ∂⁡a⁢∂⁡i and fastsimcoal2. For each simulation replicate, we performed 10 runs of ∂⁡a⁢∂⁡i and fastsimcoal2, checking to ensure that each method reached convergence.

Detailed settings for ∂⁡a⁢∂⁡i and fastsimcoal2 can be found in the Snakefile on our git repository (https://github.com/popsim-consortium/analysis). Estimates from the highest log-likelihood (out of 10 runs) for each simulation replicate are shown in Figure 4C and Appendix 1—figure 8C.

For smc++, we converted the tree sequences into VCF format and performed inference with default settings. Importantly, smc++ assumes no migration post-divergence, deviating from the simulated model. However, because smc++ allows for continuous population size changes, it is better equipped to capture many of the more complex aspects of the simulated demographic models (e.g. exponential growth).

To visualize our results, we plotted the inferred population size trajectories for each simulation replicate alongside the simulated population sizes (Figure 4C and Appendix 1—figure 8C). Here, unlike the single-population workflow, we compare our inferred population sizes only to the simulated population sizes and not the inverse coalescence rates.

Resource availability

The stdpopsim package is available for download on the Python Package Index: https://pypi.org/project/stdpopsim/. Documentation for the project can be found here: https://stdpopsim.readthedocs.io/en/latest/.

Acknowledgements

We thank the Probabilistic Modeling in Genomics conference organizers for making this collaboration possible, and the Simons Center for Quantitative Biology at Cold Spring Harbor Laboratory for sponsoring the first workshop. Early on in the project we were encouraged by many people including Patrick Phillips, Richard Durbin, Dmitri Petrov, and Sohini Ramachandran. In addition, we thank NESCENT and Matt Hahn, Victoria Sork, and Michael Whitlock for organizing a 2014 catalysis meeting in which many of the goals of this effort were first laid out. CCK and KEL were funded under NIH Award R35GM119856. JRA and ADK were funded under NIH Award R01GM117241. TJS and RNG were funded under NIH Award R01GM127348. ALG and DRS were funded under NIH award R00HG008696. ND and AS were supported in part by NIH Awards R01HG010346 and R35GM127070. FR and GG were supported by a Villum Young Investigator award (project no. 00025300). DODV is funded by a UC MEXUS-CONACYT Collaborative Grant and a DGAPA-PAPIIT grant (PAPIIT-IA200620). JK is supported by the Robertson Foundation.

Appendix 1

Calculating coalescence rates

In population genetics, the ‘effective population size’ of a population model with constant (census) size is often defined to be the number of diploids in a Wright-Fisher population that would have the same coalescence rate (or, equivalently, genetic drift) as the population in question (reviewed in Crow and Denniston, 1988). One reason the concept is useful is because theory predicts that genetic data from distinct populations with the same effective population size will look similar in many ways: for instance, their mean coalescence times will be the same. Conversely, this implies that effective population size should be easier to infer from genomic data than aspects of population demography that do not affect effective population size. An analogous observation holds for populations of changing size, if we define the ‘coalescence rate’ of a given demographic model at a particular point back in time to be the rate of coalescence of remaining lineages and define the ‘coalescence effective size’ at that time, denoted Ne⁢(t), so that the coalescence rate at time t in the past is 1/(2⁢Ne⁢(t)). With these definitions, any two models with the same effective population size trajectory (Ne⁢(t)) will have the same distribution of coalescence times. For this reason, we might guess that if we apply an inference method that assumes a Wright-Fisher population with changing size through time to a different population model, the inferred demographic history will match the ‘effective population size history’ defined in this way. These observations and the following calculations are standard in coalescent theory (see e.g. Wakeley, 2005), but they are provided here for completeness.

We compute the coalescence rate of a collection of samples in a given demographic model at a particular point back in time as the expected number of coalescences happening at that time per unit of time and per pair of as-yet-uncoalesced lineages. More concretely, let p⁢(t) denote the probability that the lineages of a randomly chosen pair of samples have not yet coalesced t units of time ago, let p⁢(z,t) denote the probability that those lineages have not yet coalesced and are furthermore both in location z, and let Ne⁢(z,t) be the (effective) diploid population size in location z at the time, so that 1/(2⁢Ne⁢(z,t)) is the rate of coalescence there. Then, we compute the mean coalescence rate as

r⁢(t)=1p⁢(t)⁢∑zp⁢(z,t)2⁢Ne⁢(z,t).

This follows because if we have m diploid samples, and hence (2m2) lineages, the expected number of coalescences in location z between times t and t+d⁢t ago is

(2m2)p(z,t)dt2Ne(z,t),

and the expected number of pairs of uncoalesced lineages at that time is

(2m2)p(t).

The expression for r⁢(t) is a ratio of these two quantities; to obtain it we need to compute p⁢(t) and p⁢(z,t). This is relatively straightforward using the general theory of Markov chains (e.g,. Kemeny et al., 2012), and is implemented in msprime.

Note that since these quantities are per pair of lineages, this definition depends on the locations of the samples. The coalescence rate also has the intuitive interpretation that it is the average between-lineage coalescence rate, averaged over where uncoalesced lineages might be. Since the local coalescence rate is the inverse of the population size, 1/r⁢(t) (as shown for instance in Figure 2) is a weighted harmonic mean of the census sizes of the different populations present at that time. This is as expected: suppose that we have two populations, one big and one small, connected by migration. If all our samples are from the big population, the number of recent coalescences should be small, reflecting the large population size, while in the long run, the coalescence rate approaches an intermediate rate. On the other hand, more recent coalescences are expected if all samples are from the small population, A method that fits a single, time-varying population size to the data might be expected to find a population size trajectory to match these time-varying rates of coalescence.

We use the same computations to analytically compute mean coalescence times: since for any nonnegative random variable T, the mean value is 𝔼[T]=∫0∞ℙ{T>t}dt, we can obtain the mean coalescence time as

∫0∞p⁢(t)⁢𝑑t,

where p⁢(t) is defined above.

The coalescence rate trajectories can be computed from a model in msprime using the coalescence_rate_trajectory method of the Demography Debugger class, which can be obtained from a stdpopsim model using the model.get_demography_debugger() method.

Appendix 1—figure 1.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig1.jpg
Validating the SLiM engine backend under a genetic map.

Here, we validate our integration of the SLiM (Haller et al., 2019; Haller and Messer, 2019) engine backend. We show quantile-quantile plots between SLiM and msprime engines for three population genetic summary statistics: r2, Tajima’s π, and Tajima’s D. Additionally, we show runtimes for generating each simulation replicate. Data were generated by simulating 100 replicates of human chromosome 22 under the AncientEurasia_9K19 model (Kamm et al., 2019) using the HapMapII_GRCh37 genetic map (Frazer et al., 2007). 12 samples were drawn from each population (excluding basal Eurasians). From top to bottom, we show results using three scaling factors for the population sizes: Q = 1, Q = 10, and Q = 50. Kolmogorov-Smirnov two-sample test statistics (D) and p-values are shown, testing the null hypothesis that the quantiles were drawn from the same continuous distribution.

Appendix 1—figure 2.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig2.jpg
Validating the SLiM engine backend under uniform recombination.

Here, we validate our integration of the SLiM (Haller et al., 2019; Haller and Messer, 2019) engine backend. We show quantile-quantile plots between SLiM and msprime engines for three population genetic summary statistics: r2, Tajima’s π, and Tajima’s D. Additionally, we show runtimes for generating each simulation replicate. Data were generated by simulating 100 replicates of human chromosome 22 under the AncientEurasia_9K19 model (Kamm et al., 2019) using a uniform rate of recombination across the chromosome. 12 samples were drawn from each population (excluding basal Eurasians). From top to bottom, we show results using three scaling factors for the population sizes: Q = 1, Q = 10, and Q = 50. Kolmogorov-Smirnov two-sample test statistics (D) and p-values are shown, testing the null hypothesis that the quantiles were drawn from the same continuous distribution.

Appendix 1—figure 3.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig3.jpg
Comparing simulated population sizes and inverse coalescence rates in humans.

Data are shown from human genomes under the OutOfAfricaArchaicAdmixture_5R19 model (Ragsdale and Gravel, 2019) and using the HapMapII_GRCh37 genetic map (Frazer et al., 2007). From left to right, we show sizes for each of the three populations in the model: YRI, CEU, and CHB. We plot the simulated sizes for each population in black, and in red we plot inverse coalescence rates as calculated from the demographic model used for simulation (see text). In this specific model, these two measures are near identical, but in other models with higher migration rates we expect to see a larger departure between the two.

Appendix 1—figure 4.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig4.jpg
Comparing estimates of N⁢(t) in humans.

Estimates of population size over time (N⁢(t)) inferred using four different methods, smc++, stairway plot, and MSMC with n=2 and n=8. Data were generated by simulating replicate human genomes under the OutOfAfrica_3G09 model (Gutenkunst et al., 2009) and using the HapMapII_GRCh37 genetic map (Frazer et al., 2007). From top to bottom, we show estimates for each of the three populations in the model: YRI, CEU, and CHB. In shades of blue, we show the estimated N⁢(t) trajectories for each replicate. As a proxy for the ‘truth’, in black we show inverse coalescence rates as calculated from the demographic model used for simulation (see text).

Appendix 1—figure 5.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig5.jpg
Comparing estimates of N⁢(t) in humans.

Here, we show estimates of population size over time (N⁢(t)) inferred using fourdifferent methods, smc+, and stairway plot, and MSMC with n=2 and n=8. Data were generated by simulating replicate human genomes under a constant sized population model with N=104 and using the HapMapII_GRCh37 genetic map (Frazer et al., 2007). As a proxy for the ‘truth’, in black we show inverse coalescence rates as calculated from the demographic model used for simulation (see text).

Appendix 1—figure 6.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig6.jpg
Comparing estimates of N⁢(t) in A. thaliana.

Here, we show estimates of population size over time (N⁢(t)) inferred using four different methods, smc++, and stairway plot, and MSMC with n=2 and n=8. Data were generated by simulating replicate A. thaliana genomes under the African2Epoch_1H18 model (Durvasula et al., 2017) and using the SalomeAveraged_TAIR7 genetic map (Salomé et al., 2012). As a proxy for the ‘truth’, in black we show inverse coalescence rates as calculated from the demographic model used for simulation (see text).

Appendix 1—figure 7.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig7.jpg
Migration rate estimates for the human Gutenkunst model.

Here, we show inferred migration rates from ∂⁡a⁢∂⁡i and fastsimcoal2. Data were generated by simulating replicate human genomes under the Gutenkunst et al., 2009 model and using the genetic map inferred in Frazer et al., 2007. Directional migration from Europe to Africa is represented as M⁢I⁢G⁢_⁢A⁢F⁢_⁢E⁢U and migration from Africa to Europe is represented as M⁢I⁢G⁢_⁢E⁢U⁢_⁢A⁢F. Note that the x-axis coordinates are arbitrary.

Appendix 1—figure 8.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig8.jpg
Parameters estimated using a two-population Drosophila model.

Here, we show estimates of N⁢(t) inferred using ∂⁡a⁢∂⁡i, fastsimcoal2, and smc++. Data were generated by simulating replicate Drosophila genomes under the Li and Stephan, 2006 model and using the genetic map inferred in Comeron et al., 2012. See legend of Figure 4 for details. In shades of blue, we show the estimated N⁢(t) trajectories for each replicate. In black we show the simulated population sizes.

Appendix 1—figure 9.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig9.jpg
Migration rate parameters estimated under a two-population Drosophila model.

Here, we show inferred migration rates from ∂⁡a⁢∂⁡i and fastsimcoal2. Data were generated by simulating replicate Drosophila genomes under the Li and Stephan, 2006 model and using the genetic map inferred in Comeron et al., 2012. Directional migration from Europe to Africa is represented as M⁢I⁢G⁢_⁢A⁢F⁢_⁢E⁢U and migration from Africa to Europe is represented as M⁢I⁢G⁢_⁢E⁢U⁢_⁢A⁢F. Note that the x-axis coordinates are arbitrary.

Appendix 1—figure 10.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig10.jpg
Workflow for our N(t) inference methods comparison.

Here, we show single replicate for two chromosomes, chr22 and chrX, simulated under the HomSap OutOfAfrica_3G09 demographic model, with a HapmapII_GRCh37 genetic map. Note that the data used as input by all inference methods smc++, MSMC, and stairway plot, come from the same set of simulations.

Appendix 1—figure 11.

An external file that holds a picture, illustration, etc. Object name is elife-54967-app1-fig11.jpg
Parameters estimated from a generic IM model Here we show estimates of N⁢(t) inferred using ∂⁡a⁢∂⁡i, fastsimcoal2, and smc++.

Data were generated by simulating under a generic IM model with a human genome and Frazer et al., 2007 genetic map. In shades of blue we show the estimated N⁢(t) trajectories for each replicate. In black we show the simulated population sizes.

Funding Statement

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Contributor Information

Graham Coop, University of California, Davis, United States.

Patricia J Wittkopp, University of Michigan, United States.

Funding Information

This paper was supported by the following grants:

  • National Institute of General Medical Sciences R35GM119856 to Christopher C Kyriazis, Kirk E Lohmueller.
  • National Institute of General Medical Sciences R01GM117241 to Jeffrey R Adrion, Andrew D Kern.
  • National Institute of General Medical Sciences R01GM127348 to Travis J Struck, Ryan N Gutenkunst.
  • National Institute of General Medical Sciences R00HG008696 to Ariella L Gladstein, Daniel R Schrider.
  • National Institute of General Medical Sciences R35GM127070 to Noah Dukler, Adam Siepel.
  • National Human Genome Research Institute R01HG010346 to Noah Dukler, Adam Siepel.
  • Villum Fonden 00025300 to Graham Gower, Fernando Racimo.
  • University of California Institute for Mexico and the United States UC MEXUS-CONACYT Collaborative Grant to Diego Ortega Del Vecchyo.
  • Consejo Nacional de Ciencia y Tecnología UC MEXUS-CONACYT Collaborative Grant to Diego Ortega Del Vecchyo.
  • Dirección General de Asuntos del Personal Académico, Universidad Nacional Autónoma de México PAPIIT-IA200620 to Diego Ortega Del Vecchyo.
  • Robertson Foundation to Jerome Kelleher.

Additional information

Competing interests

Reviewing editor, eLife.

No competing interests declared.

Author contributions

Major contribution to stdpopsim, documentation or analysis.

Major contribution to stdpopsim, documentation or analysis.

Major contribution to stdpopsim, documentation or analysis.

Major contribution to stdpopsim, documentation or analysis.

Major contribution to stdpopsim, documentation or analysis.

Major contribution to stdpopsim, documentation or analysis.

Major contribution to stdpopsim, documentation or analysis.

Major contribution to stdpopsim, documentation or analysis.

Major contribution to stdpopsim, documentation or analysis.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Contribution to software/significant community contribution.

Conceptualization, Methodology, Software, Validation, Formal Analysis, Resources, Data Curation, Writing - Original Draft Preparation, Writing - Review & Editing, Supervision, Project Administration.

Conceptualization, Methodology, Software, Validation, Formal Analysis, Resources, Data Curation, Writing - Original Draft Preparation, Writing - Review & Editing, Supervision, Project Administration.

Conceptualization, Methodology, Software, Validation, Formal Analysis, Resources, Data Curation, Writing - Original Draft Preparation, Writing - Review & Editing, Supervision, Project Administration.

Conceptualization, Methodology, Software, Validation, Formal Analysis, Resources, Data Curation, Writing - Original Draft Preparation, Writing - Review & Editing, Supervision, Project Administration.

Conceptualization, Methodology, Software, Validation, Formal Analysis, Resources, Data Curation, Writing - Original Draft Preparation, Writing - Review & Editing, Supervision, Project Administration.

Conceptualization, Methodology, Software, Validation, Formal Analysis, Resources, Data Curation, Writing - Original Draft Preparation, Writing - Review & Editing, Supervision, Project Administration.

Conceptualization, Methodology, Software, Validation, Formal Analysis, Resources, Data Curation, Writing - Original Draft Preparation, Writing - Review & Editing, Supervision, Project Administration.

Conceptualization, Methodology, Software, Validation, Formal Analysis, Resources, Data Curation, Writing - Original Draft Preparation, Writing - Review & Editing, Supervision, Project Administration.

Additional files

Transparent reporting form

Data availability

All resources are available from https://github.com/popsim-consortium/stdpopsim (copy archived at https://github.com/elifesciences-publications/stdpopsim).

References

  • Adrion JR, Galloway JG, Kern AD. Predicting the landscape of recombination using deep learning. Molecular Biology and Evolution. 2020;37:1790–1808. doi: 10.1093/molbev/msaa038. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Alachiotis N, Stamatakis A, Pavlidis P. OmegaPlus: a scalable tool for rapid detection of selective sweeps in whole-genome datasets. Bioinformatics. 2012;28:2274–2275. doi: 10.1093/bioinformatics/bts419. [PubMed] [CrossRef] [Google Scholar]
  • Andrews KR, Good JM, Miller MR, Luikart G, Hohenlohe PA. Harnessing the power of RADseq for ecological and evolutionary genomics. Nature Reviews Genetics. 2016;17:81–92. doi: 10.1038/nrg.2015.28. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Beichman AC, Phung TN, Lohmueller KE. Comparison of single genome and allele frequency data reveals discordant demographic histories. G3: Genes, Genomes, Genetics. 2017;7:3605–3620. doi: 10.1534/g3.117.300259. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Beichman AC, Huerta-Sanchez E, Lohmueller KE. Using genomic data to infer historic population dynamics of nonmodel organisms. Annual Review of Ecology, Evolution, and Systematics. 2018;49:433–456. doi: 10.1146/annurev-ecolsys-110617-062431. [CrossRef] [Google Scholar]
  • Boyko AR, Williamson SH, Indap AR, Degenhardt JD, Hernandez RD, Lohmueller KE, Adams MD, Schmidt S, Sninsky JJ, Sunyaev SR, White TJ, Nielsen R, Clark AG, Bustamante CD. Assessing the evolutionary impact of amino acid mutations in the human genome. PLOS Genetics. 2008;4:e1000083. doi: 10.1371/journal.pgen.1000083. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Browning SR, Browning BL, Daviglus ML, Durazo-Arvizu RA, Schneiderman N, Kaplan RC, Laurie CC. Ancestry-specific recent effective population size in the americas. PLOS Genetics. 2018;14:e1007385. doi: 10.1371/journal.pgen.1007385. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Campbell CL, Bhérer C, Morrow BE, Boyko AR, Auton A. A Pedigree-Based map of recombination in the domestic dog genome. G3: Genes, Genomes, Genetics. 2016;6:3517–3524. doi: 10.1534/g3.116.034678. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Chan AH, Jenkins PA, Song YS. Genome-wide fine-scale recombination rate variation in Drosophila melanogaster. PLOS Genetics. 2012;8:e1003090. doi: 10.1371/journal.pgen.1003090. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Comeron JM, Ratnappan R, Bailin S. The many landscapes of recombination in Drosophila melanogaster. PLOS Genetics. 2012;8:e1002905. doi: 10.1371/journal.pgen.1002905. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Crow JF, Denniston C. Inbreeding and variance effective population numbers. Evolution. 1988;42:482–495. doi: 10.1111/j.1558-5646.1988.tb04154.x. [PubMed] [CrossRef] [Google Scholar]
  • Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, McVean G, Durbin R, 1000 Genomes Project Analysis Group The variant call format and VCFtools. Bioinformatics. 2011;27:2156–2158. doi: 10.1093/bioinformatics/btr330. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • DeGiorgio M, Huber CD, Hubisz MJ, Hellmann I, Nielsen R. SweepFinder2: increased sensitivity, robustness and flexibility. Bioinformatics. 2016;32:1895–1897. doi: 10.1093/bioinformatics/btw051. [PubMed] [CrossRef] [Google Scholar]
  • Durvasula A, Fulgione A, Gutaker RM, Alacakaptan SI, Flood PJ, Neto C, Tsuchimatsu T, Burbano HA, Picó FX, Alonso-Blanco C, Hancock AM. African genomes illuminate the early history and transition to selfing in Arabidopsis thaliana. PNAS. 2017;114:5213–5218. doi: 10.1073/pnas.1616736114. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Excoffier L, Dupanloup I, Huerta-Sánchez E, Sousa VC, Foll M. Robust demographic inference from genomic and SNP data. PLOS Genetics. 2013;9:e1003905. doi: 10.1371/journal.pgen.1003905. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Eyre-Walker A, Keightley PD. Estimating the rate of adaptive molecular evolution in the presence of slightly deleterious mutations and population size change. Molecular Biology and Evolution. 2009;26:2097–2108. doi: 10.1093/molbev/msp119. [PubMed] [CrossRef] [Google Scholar]
  • Fortier AL, Coffman AJ, Struck TJ, Irby MN, Burguete JEL, Ragsdale AP, Gutenkunst RN. DFEnitely different: genome-wide characterization of differences in mutation fitness effects between populations. bioRxiv. 2019 doi: 10.1101/703918. [CrossRef]
  • Frazer KA, Ballinger DG, Cox DR, Hinds DA, Stuve LL, Gibbs RA, Belmont JW, Boudreau A, Hardenbol P, Leal SM, Pasternak S, Wheeler DA, Willis TD, Yu F, Yang H, Zeng C, Gao Y, Hu H, Hu W, Li C, Lin W, Liu S, Pan H, Tang X, Wang J, Wang W, Yu J, Zhang B, Zhang Q, Zhao H, Zhao H, Zhou J, Gabriel SB, Barry R, Blumenstiel B, Camargo A, Defelice M, Faggart M, Goyette M, Gupta S, Moore J, Nguyen H, Onofrio RC, Parkin M, Roy J, Stahl E, Winchester E, Ziaugra L, Altshuler D, Shen Y, Yao Z, Huang W, Chu X, He Y, Jin L, Liu Y, Shen Y, Sun W, Wang H, Wang Y, Wang Y, Xiong X, Xu L, Waye MM, Tsui SK, Xue H, Wong JT, Galver LM, Fan JB, Gunderson K, Murray SS, Oliphant AR, Chee MS, Montpetit A, Chagnon F, Ferretti V, Leboeuf M, Olivier JF, Phillips MS, Roumy S, Sallée C, Verner A, Hudson TJ, Kwok PY, Cai D, Koboldt DC, Miller RD, Pawlikowska L, Taillon-Miller P, Xiao M, Tsui LC, Mak W, Song YQ, Tam PK, Nakamura Y, Kawaguchi T, Kitamoto T, Morizono T, Nagashima A, Ohnishi Y, Sekine A, Tanaka T, Tsunoda T, Deloukas P, Bird CP, Delgado M, Dermitzakis ET, Gwilliam R, Hunt S, Morrison J, Powell D, Stranger BE, Whittaker P, Bentley DR, Daly MJ, de Bakker PI, Barrett J, Chretien YR, Maller J, McCarroll S, Patterson N, Pe'er I, Price A, Purcell S, Richter DJ, Sabeti P, Saxena R, Schaffner SF, Sham PC, Varilly P, Altshuler D, Stein LD, Krishnan L, Smith AV, Tello-Ruiz MK, Thorisson GA, Chakravarti A, Chen PE, Cutler DJ, Kashuk CS, Lin S, Abecasis GR, Guan W, Li Y, Munro HM, Qin ZS, Thomas DJ, McVean G, Auton A, Bottolo L, Cardin N, Eyheramendy S, Freeman C, Marchini J, Myers S, Spencer C, Stephens M, Donnelly P, Cardon LR, Clarke G, Evans DM, Morris AP, Weir BS, Tsunoda T, Mullikin JC, Sherry ST, Feolo M, Skol A, Zhang H, Zeng C, Zhao H, Matsuda I, Fukushima Y, Macer DR, Suda E, Rotimi CN, Adebamowo CA, Ajayi I, Aniagwu T, Marshall PA, Nkwodimmah C, Royal CD, Leppert MF, Dixon M, Peiffer A, Qiu R, Kent A, Kato K, Niikawa N, Adewole IF, Knoppers BM, Foster MW, Clayton EW, Watkin J, Gibbs RA, Belmont JW, Muzny D, Nazareth L, Sodergren E, Weinstock GM, Wheeler DA, Yakub I, Gabriel SB, Onofrio RC, Richter DJ, Ziaugra L, Birren BW, Daly MJ, Altshuler D, Wilson RK, Fulton LL, Rogers J, Burton J, Carter NP, Clee CM, Griffiths M, Jones MC, McLay K, Plumb RW, Ross MT, Sims SK, Willey DL, Chen Z, Han H, Kang L, Godbout M, Wallenburg JC, L'Archevêque P, Bellemare G, Saeki K, Wang H, An D, Fu H, Li Q, Wang Z, Wang R, Holden AL, Brooks LD, McEwen JE, Guyer MS, Wang VO, Peterson JL, Shi M, Spiegel J, Sung LM, Zacharia LF, Collins FS, Kennedy K, Jamieson R, Stewart J, International HapMap Consortium A second generation human haplotype map of over 3.1 million SNPs. Nature. 2007;449:851–861. doi: 10.1038/nature06258. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Garud NR, Messer PW, Buzbas EO, Petrov DA. Recent selective sweeps in north american Drosophila melanogaster show signatures of soft sweeps. PLOS Genetics. 2015;11:e1005004. doi: 10.1371/journal.pgen.1005004. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Gladstein AL, Hammer MF. Substructured population growth in the ashkenazi jews inferred with approximate bayesian computation. Molecular Biology and Evolution. 2019;36:1162–1171. doi: 10.1093/molbev/msz047. [PubMed] [CrossRef] [Google Scholar]
  • Gutenkunst RN, Hernandez RD, Williamson SH, Bustamante CD. Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLOS Genetics. 2009;5:e1000695. doi: 10.1371/journal.pgen.1000695. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Haller BC, Galloway J, Kelleher J, Messer PW, Ralph PL. Tree-sequence recording in SLiM opens new horizons for forward-time simulation of whole genomes. Molecular Ecology Resources. 2019;19:552–566. doi: 10.1111/1755-0998.12968. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Haller BC, Messer PW. SLiM 3: forward genetic simulations beyond the Wright-Fisher model. Molecular Biology and Evolution. 2019;36:632–637. doi: 10.1093/molbev/msy228. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Hey J, Nielsen R. Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics. 2004;167:747–760. doi: 10.1534/genetics.103.024182. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Huang YF, Siepel A. Estimation of allele-specific fitness effects across human protein-coding sequences and implications for disease. Genome Research. 2019;29:1310–1321. doi: 10.1101/gr.245522.118. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Huber CD, Durvasula A, Hancock AM, Lohmueller KE. Gene expression drives the evolution of dominance. Nature Communications. 2018;9:7. doi: 10.1038/s41467-018-05281-7. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Jacobs GS, Hudjashov G, Saag L, Kusuma P, Darusallam CC, Lawson DJ, Mondal M, Pagani L, Ricaut FX, Stoneking M, Metspalu M, Sudoyo H, Lansing JS, Cox MP. Multiple deeply divergent denisovan ancestries in papuans. Cell. 2019;177:1010–1021. doi: 10.1016/j.cell.2019.02.035. [PubMed] [CrossRef] [Google Scholar]
  • Kamm J, Terhorst J, Durbin R, Song YS. Efficiently inferring the demographic history of many populations with allele count data. Journal of the American Statistical Association. 2019;155:1–16. doi: 10.1080/01621459.2019.1635482. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Kelleher J, Etheridge AM, McVean G. Efficient coalescent simulation and genealogical analysis for large sample sizes. PLOS Computational Biology. 2016;12:e1004842. doi: 10.1371/journal.pcbi.1004842. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Kelleher J, Thornton KR, Ashander J, Ralph PL. Efficient pedigree recording for fast population genetics simulation. PLOS Computational Biology. 2018;14:e1006581. doi: 10.1371/journal.pcbi.1006581. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Kelleher J, Wong Y, Wohns AW, Fadil C, Albers PK, McVean G. Inferring whole-genome histories in large population datasets. Nature Genetics. 2019;51:1330–1338. doi: 10.1038/s41588-019-0483-y. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Kemeny JG, Snell JL, Knapp AW. Denumerable Markov Chains. Springer Science & Business Media; 2012. [CrossRef] [Google Scholar]
  • Kern AD, Schrider DR. diploS/HIC: an updated approach to classifying selective sweeps. G3: Genes, Genomes, Genetics. 2018;8:1959–1970. doi: 10.1534/g3.118.200262. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Kim BY, Huber CD, Lohmueller KE. Inference of the distribution of selection coefficients for new nonsynonymous mutations using large samples. Genetics. 2017;206:345–361. doi: 10.1534/genetics.116.197145. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Kim Y, Stephan W. Detecting a local signature of genetic hitchhiking along a recombining chromosome. Genetics. 2002;160:765–777. [PMC free article] [PubMed] [Google Scholar]
  • Kong A, Thorleifsson G, Gudbjartsson DF, Masson G, Sigurdsson A, Jonasdottir A, Walters GB, Jonasdottir A, Gylfason A, Kristinsson KT, Gudjonsson SA, Frigge ML, Helgason A, Thorsteinsdottir U, Stefansson K. Fine-scale recombination rate differences between sexes, populations and individuals. Nature. 2010;467:1099–1103. doi: 10.1038/nature09525. [PubMed] [CrossRef] [Google Scholar]
  • Köster J, Rahmann S. Snakemake--a scalable bioinformatics workflow engine. Bioinformatics. 2012;28:2520–2522. doi: 10.1093/bioinformatics/bts480. [PubMed] [CrossRef] [Google Scholar]
  • Langley CH, Stevens K, Cardeno C, Lee YC, Schrider DR, Pool JE, Langley SA, Suarez C, Corbett-Detig RB, Kolaczkowski B, Fang S, Nista PM, Holloway AK, Kern AD, Dewey CN, Song YS, Hahn MW, Begun DJ. Genomic variation in natural populations of Drosophila melanogaster. Genetics. 2012;192:533–598. doi: 10.1534/genetics.112.142018. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Li H, Durbin R. Inference of human population history from individual whole-genome sequences. Nature. 2011;475:493–496. doi: 10.1038/nature10231. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Li H, Stephan W. Inferring the demographic history and rate of adaptive substitution in Drosophila. PLOS Genetics. 2006;2:e166. doi: 10.1371/journal.pgen.0020166. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Lin K, Futschik A, Li H. A fast estimate for the population recombination rate based on regression. Genetics. 2013;194:473–484. doi: 10.1534/genetics.113.150201. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Liu X, Fu YX. Exploring population size changes using SNP frequency spectra. Nature Genetics. 2015;47:555–559. doi: 10.1038/ng.3254. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Locke DP, Hillier LW, Warren WC, Worley KC, Nazareth LV, Muzny DM, Yang SP, Wang Z, Chinwalla AT, Minx P, Mitreva M, Cook L, Delehaunty KD, Fronick C, Schmidt H, Fulton LA, Fulton RS, Nelson JO, Magrini V, Pohl C, Graves TA, Markovic C, Cree A, Dinh HH, Hume J, Kovar CL, Fowler GR, Lunter G, Meader S, Heger A, Ponting CP, Marques-Bonet T, Alkan C, Chen L, Cheng Z, Kidd JM, Eichler EE, White S, Searle S, Vilella AJ, Chen Y, Flicek P, Ma J, Raney B, Suh B, Burhans R, Herrero J, Haussler D, Faria R, Fernando O, Darré F, Farré D, Gazave E, Oliva M, Navarro A, Roberto R, Capozzi O, Archidiacono N, Della Valle G, Purgato S, Rocchi M, Konkel MK, Walker JA, Ullmer B, Batzer MA, Smit AF, Hubley R, Casola C, Schrider DR, Hahn MW, Quesada V, Puente XS, Ordoñez GR, López-Otín C, Vinar T, Brejova B, Ratan A, Harris RS, Miller W, Kosiol C, Lawson HA, Taliwal V, Martins AL, Siepel A, Roychoudhury A, Ma X, Degenhardt J, Bustamante CD, Gutenkunst RN, Mailund T, Dutheil JY, Hobolth A, Schierup MH, Ryder OA, Yoshinaga Y, de Jong PJ, Weinstock GM, Rogers J, Mardis ER, Gibbs RA, Wilson RK. Comparative and demographic analysis of orang-utan genomes. Nature. 2011;469:529–533. doi: 10.1038/nature09687. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • McVean GA, Myers SR, Hunt S, Deloukas P, Bentley DR, Donnelly P. The fine-scale structure of recombination rate variation in the human genome. Science. 2004;304:581–584. doi: 10.1126/science.1092500. [PubMed] [CrossRef] [Google Scholar]
  • Moult J, Pedersen JT, Judson R, Fidelis K. A large-scale experiment to assess protein structure prediction methods. Proteins: Structure, Function, and Genetics. 1995;23:303. doi: 10.1002/prot.340230303. [PubMed] [CrossRef] [Google Scholar]
  • Nater A, Mattle-Greminger MP, Nurcahyo A, Nowak MG, de Manuel M, Desai T, Groves C, Pybus M, Sonay TB, Roos C, Lameira AR, Wich SA, Askew J, Davila-Ross M, Fredriksson G, de Valles G, Casals F, Prado-Martinez J, Goossens B, Verschoor EJ, Warren KS, Singleton I, Marques DA, Pamungkas J, Perwitasari-Farajallah D, Rianti P, Tuuga A, Gut IG, Gut M, Orozco-terWengel P, van Schaik CP, Bertranpetit J, Anisimova M, Scally A, Marques-Bonet T, Meijaard E, Krützen M. Morphometric, behavioral, and genomic evidence for a new orangutan species. Current Biology. 2017;27:3487–3498. doi: 10.1016/j.cub.2017.09.047. [PubMed] [CrossRef] [Google Scholar]
  • Ragsdale AP, Gravel S. Models of archaic admixture and recent history from two-locus statistics. PLOS Genetics. 2019;15:e1008204. doi: 10.1371/journal.pgen.1008204. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Russakovsky O, Deng J, Su H, Krause J, Satheesh S, Ma S, Huang Z, Karpathy A, Khosla A, Bernstein M, Berg AC, Fei-Fei L. ImageNet large scale visual recognition challenge. International Journal of Computer Vision. 2015;115:211–252. doi: 10.1007/s11263-015-0816-y. [CrossRef] [Google Scholar]
  • Salomé PA, Bomblies K, Fitz J, Laitinen RA, Warthmann N, Yant L, Weigel D. The recombination landscape in Arabidopsis thaliana F2 populations. Heredity. 2012;108:447–455. doi: 10.1038/hdy.2011.95. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Schiffels S, Durbin R. Inferring human population size and separation history from multiple genome sequences. Nature Genetics. 2014;46:919–925. doi: 10.1038/ng.3015. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Sheehan S, Song YS. Deep learning for population genetic inference. PLOS Computational Biology. 2016;12:e1004845. doi: 10.1371/journal.pcbi.1004845. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Sugden LA, Atkinson EG, Fischer AP, Rong S, Henn BM, Ramachandran S. Localization of adaptive variants in human genomes using averaged one-dependence estimation. Nature Communications. 2018;9:703. doi: 10.1038/s41467-018-03100-7. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Tataru P, Mollion M, Glémin S, Bataillon T. Inference of distribution of fitness effects and proportion of adaptive substitutions from polymorphism data. Genetics. 2017;207:1103–1119. doi: 10.1534/genetics.117.300323. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Tennessen JA, Bigham AW, O'Connor TD, Fu W, Kenny EE, Gravel S, McGee S, Do R, Liu X, Jun G, Kang HM, Jordan D, Leal SM, Gabriel S, Rieder MJ, Abecasis G, Altshuler D, Nickerson DA, Boerwinkle E, Sunyaev S, Bustamante CD, Bamshad MJ, Akey JM, Broad GO, Seattle GO, NHLBI Exome Sequencing Project Evolution and functional impact of rare coding variation from deep sequencing of human exomes. Science. 2012;337:64–69. doi: 10.1126/science.1219240. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Terhorst J, Kamm JA, Song YS. Robust and scalable inference of population history from hundreds of unphased whole genomes. Nature Genetics. 2017;49:303–309. doi: 10.1038/ng.3748. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Uricchio LH, Hernandez RD. Robust forward simulations of recurrent hitchhiking. Genetics. 2014;197:221–236. doi: 10.1534/genetics.113.156935. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • V Barroso G, Puzović N, Dutheil JY. Inference of recombination maps from a single pair of genomes and its application to ancient samples. PLOS Genetics. 2019;15:e1008449. doi: 10.1371/journal.pgen.1008449. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
  • Vecchyo DOD, Lohmueller KE, Novembre J. Haplotype-based inference of the distribution of fitness effects. bioRxiv. 2019 doi: 10.1101/770966. [PMC free article] [PubMed] [CrossRef]
  • Wakeley J. Coalescent Theory, an Introduction. Roberts and Company; 2005. http://www.coalescentheory.com/ [Google Scholar]
2020; 9: e54967.
Published online 2020 Jun 23. doi: 10.7554/eLife.54967.sa1

Decision letter

Graham Coop, Reviewing Editor
Graham Coop, University of California, Davis, United States;
Reviewed by John Novembre, Reviewer, Arun Sethuraman, Reviewer, and Sara Mathieson, Reviewer
Author information Copyright and License information PMC Disclaimer

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

Simulations have long been central to population genetics. Population genomics has, in turn, become central to numerous areas of evolutionary biology and human genetics, and a vast array of different statistical methods have been developed. However, these methods are rarely ground-truthed in any truly reproducible manner. This paper takes an important set of steps in developing a flexible framework to standardize testing of population genetics software.

Decision letter after peer review:

Thank you for submitting your article "A community-maintained standard library of population genetic models" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Patricia Wittkopp as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: John Novembre (Reviewer #1); Arun Sethuraman (Reviewer #2); Sara Mathieson (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Simulations have long been central to population genetics. Population genomics has, in turn, become central to numerous areas of evolutionary biology and human genetics, and a vast array of different statistical methods have been developed. However, these methods are rarely ground-truthed in any truly reproducible manner. This paper takes an important set of steps in developing a flexible framework to standardize testing of population genetics software.

The reviewers' comments are all reasonably straight forward, often hit on similar points, and so most should be easily addressable. Please respond to each point in your response to reviews. In our online discussions with the reviewers three points came the fore:

1) Ensuring that the language about the models and parameters conveys the correct sense of uncertainty. Obviously the reviewers and the authors know that these just represent best guesses at the moment, but as this platform catches on we don't want these numbers and models to be seen as gospel.

2) We would like to see the SLIM integration demonstrated with an application that MsPrime doesn't cover. This could be, for example, a figure of the average value some summary statistics surrounding a selective under the human demographic model. The reviewers didn't want this to be a lot of work, e.g. it would be outside of the scope of the current paper to demonstrate the power of a range of selection methods. However, we felt that a simple demonstration of the functionality would substantially increase the scenarios easily available and imaginable by readers of the paper.

3) We felt that the Discussion would benefit from a discussion of the details not yet incorporated into the platform, e.g. sex-specific recombination, gene conversion, assembly error, structural variants etc. Such an addition might spur future development, and also give the authors an opportunity to discuss general ways forward on these issues.

4) We would be interested to see more of a discussion of the strengths and weakness of the various demographic inference methods.

Beyond testing methods, another major use axis that I can personally see for this platform is for empirical researchers to see whether patterns they see in their data are consistent with known demographic models. For example, demonstrating whether the frequency spectrum in larger samples lines up with previously inferred models. Having all of these models implemented in a central place will significantly lower the barrier of entry of empirical researchers into rigorous simulation frameworks.

Overall we are very excited to see this big move forward and welcome the authors' careful work in this important area.

Reviewer #1:

The manuscript describes a community effort to standardize population genetic simulations and it presents an example of the resource's utility for method's testing. This development answers a long-standing need in the population genetics community for greater standardization and ease-of-implementation of simulation protocols. As such, I expect it will be a well-used resource and it represents an important advance for the field's practices. From the test implementation included in the paper, it was remarkable to see how highly variable (and often poor) the performance of the demographic inference methods was, and those results add to the scientific contribution of this manuscript.

Some of the specific features that are nice about the work are the incorporation of inferred genetic maps where possible, the automated output of citation information, the use of long-term stable identifiers for models and genetic maps, and the outputting of provenance information in the tree sequence files.

1) Emphasize parameters are current best-guesses:

a) I am most worried about the fragility of the models over time and the misperception one might have that these models are "accurate" for a species. Genome assemblies change, recombination rates improve, mutation rates change, etc. I would like to see that uncertainty reflected more in the language used in the paper so that it's clear the catalog is a collection of inferred parameters that are subject to change over time. This is a subtle and cosmetic, but important I think. For example, "the library defines basic information about each species genome, including information about chromosome lengths, average mutation rates, and generation times." Yet, this is information we don't have – each of these parameters is actively under revision/discussion for even the best studied species. It would be great to hammer home these are all inferred values. (Also see sentence on "details on the physical organization of its genome"; the same comment applies, and an improvement might be something like "details on the physical organization of the latest reference genome assembled for the species").

b) On a related note, there is wording regarding ensuring "implemented models are accurate" – I think what is meant is that the "implemented models are faithful to the source publications from which they derive". As the second half of this paper makes clear, because of errors in inference, many of these models will not be accurate, in the sense of representing the true history of the species.

2) It would be wonderful to have a comments section for the models that would be either curated by a set of editors or crowd-sourced. I say this because overtime, models will proliferate, and some will come to be regarded as out-of-date. One can imagine those approaching the field afresh will be overwhelmed by the possible selections and possible implement models that become outdated. If the goal is standardization, how do we communicate that standards change? A comment system (or even star-rating system?) may be wise to implement now. Assuming it can be layered on top of the existing structure, it may be enough for this publication to note this as a future challenge that needs development/addressing.

3) In terms of the maturity of the examples developed for the initial release, I would have liked to see at least one simulation model with a selective sweep, one with background selection, and one with spatial stepping-stone structure. Each of these would be helpful test cases to implement to be sure that the existing catalog framework has the breadth/flexibility necessary to accommodate future use cases. I do not think this is a requirement for publication, but it would add great value to this initial release of the resource.

4) The approach of masking "low-recombination" portions of the chromosomes seems like an incomplete/indirect attempt to model the inherent limitations of sequencing to an "accessible" genome.

a) Shouldn't the approach instead be to drop "low complexity" regions (e.g., as defined by an excessive number of "N"'s in the reference, low mapability scores, or via tools like RepeatMasker?). This part of the pipeline seems open to refinement.

b) Are the "masks" a separate configuration file for the simulations? It seems that it would be preferable for them to be separate from the recombination rate files – right now it reads as if the mask applied is a function of the genetic map file, but this seems too inflexible for users who prefer an alternative approach to masking.

Reviewer #2:

The authors in this manuscript describe the implementation of a publicly curated, open source simulation package called stdpopsim – equipped with commonly utilized population genetic demographic models in humans, Drosophila melanogaster, and Arabidopsis thaliana. I am in awe of what these authors have achieved, in terms of benchmarking these standard models in an effort to avoid duplication of effort, and possibility of erroneous inference. The package is currently equipped with several "in-built" models that allow the simulation of trees with msprime and SLiM. The authors explain the application of these models by simulating and benchmarking estimates of Ne under a couple of scenarios. The manuscript is also well written, and use popular tools like dadi and smc++ to estimate and benchmark the simulations under a variety of models. Across all simulated scenarios (except under more complex models), the simulations seem relatively accurate.

Despite a little hiccup with python version mismatch prior to me successfully installing stdpopsim, I was able to successfully get the tutorials running within minutes after. I have one recommendation however for the tutorials – it would definitely help if the CLI versus python tutorials were kept separate. I found it a little confusing since they are all listed on the same page (https://stdpopsim.readthedocs.io/en/latest/tutorial.html). The simulations, testing models, calculating divergences, plotting ran without a hitch, and I am impressed and excited to play around with more models in coming days. Having also developed similar libraries/pipelines, I have also found it extremely useful for developers to provide some more detailed documentation/tutorials via Jupyter notebooks, or some similar platform. I did however notice that the authors have provided their analyses as Snakemake files in the interest of replication. I did not replicate their analyses, but I trust that the documentation for these analyses are detailed enough to aid readers/users in establishing similar analysis pipelines for stdpopsim simulated data.

I thoroughly enjoyed reading this manuscript, and learning of all the new features that have been written into stdpopsim, and believe that this will be an invaluable contribution to the population genomics community.

Reviewer #3:

The authors present a standardized framework for creating reproducible population genetic simulations. This resource will allow researchers to create models for new species/scenarios, and easily compare methods on the same dataset. The authors are correct that the current state of benchmarking in population genetics causes inconsistency, duplicated effort, and confusion about the performance of different methods. The authors highlight that creating a realistic and meaningful simulation study is a barrier to entering this field, and I absolutely agree. I will be passing this resource on to my students and look forward to using it myself. Stdpopsim is a crucial step forward. The manuscript itself is well-written and describes an extensive comparison of demography methods in a variety of species/scenarios. I have a few comments.

1) The authors mention that SLiM can be used as an alternative backend, which would presumably allow for simulations with selection. Although I don't think an extensive comparison of selection methods is necessary for this paper, it would be ideal if the authors can give some idea of how this would work (example command line, etc). There are also a myriad of methods for detecting/quantifying selection, and these simulations are not consistent either.

2) I like the inclusion of the "zigzag" history, as well as generic piecewise constant models and IM models (subsection “The Species Catalog”). I wonder if these could be included in a separate section (not organism specific) in the documentation and software (and then in Table 1). Right now the zigzag model is under humans in the catalog.

3) In the subsection “Use case: comparing methods of demographic inference”, the authors set up notation for the number of replicates (R), number of chromosomes (C), and sample size (n), but don't seem to use it afterward (or use it inconsistently). It would be helpful if all the figure legends and main text included this notation (I am guessing the number of replicates is 3 based on the images, but this should be clarified). The authors use N in the Materials and methods (i.e. subsection “Workflow for analysis of simulated data”) to refer to population size (which makes sense), but then also say "In all cases we set the sample size of the focal population to N = 50 chromosomes." For MSMC, the sample size was set to n=2,8 which suggests haploid samples, but the "Calculating coalescence rates" section says that n is the diploid sample size.

4) "Calculating coalescence rates" section needs a read through. Reword first sentence and add some citations (especially regarding computing p(t) and p(z,t)). It was unclear to me how the "mean coalescence times" were used (the rate was used to compute the ground truth over time). This section is also referred to as the Appendix in the main text.

Contributor Information

John Novembre, University of Chicago, United States.

Arun Sethuraman,

Sara Mathieson,

2020; 9: e54967.
Published online 2020 Jun 23. doi: 10.7554/eLife.54967.sa2

Author response

Copyright and License information PMC Disclaimer

[…] In our online discussions with the reviewers three points came the fore:

1) Ensuring that the language about the models and parameters conveys the correct sense of uncertainty. Obviously the reviewers and the authors know that these just represent best guesses at the moment, but as this platform catches on we don't want these numbers and models to be seen as gospel.

We fully agree that the language in the original manuscript was not careful enough in conveying uncertainty about the models and parameters we discuss. We have made edited the text throughout the manuscript to address this. See our response to reviewer 1 below for examples. We have also added a bit more to the Discussion (below “next steps”) on this point.

2) We would like to see the SLIM integration demonstrated with an application that MsPrime doesn't cover. This could be, for example, a figure of the average value some summary statistics surrounding a selective under the human demographic model. The reviewers didn't want this to be a lot of work, e.g. it would be outside of the scope of the current paper to demonstrate the power of a range of selection methods. However, we felt that a simple demonstration of the functionality would substantially increase the scenarios easily available and imaginable by readers of the paper.

We would like to see this also and agree that it would be of great interest. However, we felt that performing adequate simulations with selection would require a substantial amount of additional work and would significantly expand the scope of the manuscript. Furthermore, discussing selection would distract from the current relatively simple focus on inference of demographic history. We definitely intend to include selection (and expand the available generic models) in future work, as we now discuss in greater detail in the Discussion (“next steps”).

Instead, we have decided to add a section validating the use of SLiM in stdpopsim neutral simulations (see “simulation engines” section in the Results and Appendix—figures 1 and 2). Our focus there was to show that under neutral simulations, SLiM produced data that were consistent with the coalescent simulations of msprime. To expand on this fundamental point, we also examined the influence of population size down-scaling, which is a common practice used in forward simulations, and indeed often crucial for tractable compute times. While this is by no means a comprehensive investigation into this issue, we do believe that it demonstrates the use of SLiM in future applications that will also simulate selection.

3) We felt that the Discussion would benefit from a discussion of the details not yet incorporated into the platform, e.g. sex-specific recombination, gene conversion, assembly error, structural variants etc. Such an addition might spur future development, and also give the authors an opportunity to discuss general ways forward on these issues.

Thank you for this good idea – we have expanded the section of the Discussion

(“next steps”) with more detail on our future plans. While that is so we are wary of promising too much at this point.

4) We would be interested to see more of a discussion of the strengths and weakness of the various demographic inference methods.

We have added an additional paragraph to the Discussion exploring various factors that may affect choice of inference method. These include the data required by the method, the type of model, and the implementation. We also provide a qualitative comparison between the performance of methods, as shown in our limited analysis. We believe that a more comprehensive comparison between methods is beyond the scope of the current manuscript, which focuses on the resource itself, rather than its application.

Beyond testing methods, another major use axis that I can personally see for this platform is for empirical researchers to see whether patterns they see in their data are consistent with known demographic models. For example, demonstrating whether the frequency spectrum in larger samples lines up with previously inferred models. Having all of these models implemented in a central place will significantly lower the barrier of entry of empirical researchers into rigorous simulation frameworks.

This is an important point – we have (naturally) focused on methods benchmarking, but our impact could well be greater if stdpopsim becomes widely used among empirical researchers (for whom on average running realistic simulations presents a greater barrier than to developers of computational methods). We’ve added a paragraph about this to the Discussion and additional words to the Abstract.

Reviewer #1:

[…] 1) Emphasize parameters are current best-guesses:

a) I am most worried about the fragility of the models over time and the misperception one might have that these models are "accurate" for a species. Genome assemblies change, recombination rates improve, mutation rates change, etc. I would like to see that uncertainty reflected more in the language used in the paper so that it's clear the catalog is a collection of inferred parameters that are subject to change over time. This is a subtle and cosmetic, but important I think. For example, "the library defines basic information about each species genome, including information about chromosome lengths, average mutation rates, and generation times." Yet, this is information we don't have – each of these parameters is actively under revision/discussion for even the best studied species. It would be great to hammer home these are all inferred values. (Also see sentence on "details on the physical organization of its genome"; the same comment applies, and an improvement might be something like "details on the physical organization of the latest reference genome assembled for the species").

We thank the reviewer for bringing up this important point and fully agree that the language in the manuscript needs to better emphasize the uncertainty in our current parameter estimates. To address this, we have made several changes throughout the manuscript, including the passages cited above:

The first quoted passage now reads: “For each species, the catalog contains curated information on our current understanding of the physical organization of its genome, inferred genetic maps, population-level parameters (e.g., mutation rate and generation time estimates), and published demographic models. These models and parameters are meant to represent the field’s current understanding, and we intend for this resource to evolve as new results become available, and other existing models are added to stdpopsim by the community.”

The second passage now reads: “Firstly, the library defines some basic information about our current understanding of each species’ genome, including information about chromosome lengths, average mutation rate estimates, and generation times. We also provide access to detailed empirical information such as inferred genetic maps, which model observed heterogeneity in recombination rate along chromosomes.”

b) On a related note, there is wording regarding ensuring "implemented models are accurate" – I think what is meant is that the "implemented models are faithful to the source publications from which they derive". As the second half of this paper makes clear, because of errors in inference, many of these models will not be accurate, in the sense of representing the true history of the species.

We made the suggested change: “Importantly, we developed rigorous quality control methods to ensure that we have correctly implemented the models as described in their original publication and provided documented methods for others to contribute new models.”

2) It would be wonderful to have a comments section for the models that would be either curated by a set of editors or crowd-sourced. I say this because overtime, models will proliferate, and some will come to be regarded as out-of-date. One can imagine those approaching the field afresh will be overwhelmed by the possible selections and possible implement models that become outdated. If the goal is standardization, how do we communicate that standards change? A comment system (or even star-rating system?) may be wise to implement now. Assuming it can be layered on top of the existing structure, it may be enough for this publication to note this as a future challenge that needs development/addressing.

Thank you for this forward-thinking comment. We considered various ways to do this (e.g., by enabling the wiki associated with the github repository: https://github.com/popsim-consortium/stdpopsim/issues/415)

3) In terms of the maturity of the examples developed for the initial release, I would have liked to see at least one simulation model with a selective sweep, one with background selection, and one with spatial stepping-stone structure. Each of these would be helpful test cases to implement to be sure that the existing catalog framework has the breadth/flexibility necessary to accommodate future use cases. I do not think this is a requirement for publication, but it would add great value to this initial release of the resource.

As mentioned in our response to the editor’s comments, we agree that this is an important avenue to pursue in the near future. However, performing adequate simulations with selection would require a substantial amount of additional work and would likely distract from the current relatively simple focus on inference of demographic history. We definitely intend to include selection (and expand the available generic models) in future work, as we now discuss in greater detail in the Discussion (“next steps”).

4) The approach of masking "low-recombination" portions of the chromosomes seems like an incomplete/indirect attempt to model the inherent limitations of sequencing to an "accessible" genome.

a) Shouldn't the approach instead be to drop "low complexity" regions (e.g., as defined by an excessive number of "N"'s in the reference, low mapability scores, or via tools like RepeatMasker?). This part of the pipeline seems open to refinement.

Our initial motivation for masking was to reduce the overrepresentation of marginal trees with little to no recombination from biasing patterns of diversity in such a way that demographic inference methods would be misled. While we agree with the reviewer that different masking approaches might better reflect the masking done on real genomic data, at this time the best way to mask remains an open question, and we feel that a nuanced analysis of masking is outside the scope of this paper. For these reasons, we feel that masking using a simple recombination rate threshold is a reasonable approach. We note that there is however a substantial correlation between recombination rate and mappability (by any of the measures mentioned above) wherein the notoriously difficult portions of genomes to assemble are generally in regions of low recombination. We agree that there may indeed be better ways to mask, which is why we allow users to mask their simulations as they see fit.

b) Are the "masks" a separate configuration file for the simulations? It seems that it would be preferable for them to be separate from the recombination rate files – right now it reads as if the mask applied is a function of the genetic map file, but this seems too inflexible for users who prefer an alternative approach to masking.

Mask files are not currently a component of stdpopsim proper, rather they were implemented separately from running stdpopsim, for the sole purpose of comparing demographic inference methods. All of the masks that we have used are available on the analysis repository that is associated with this manuscript. Users who download stdpopsim will always be simulating raw and unmasked tree sequence files, to which they can apply any variety of masks post hoc, if they so choose. We are currently making plans, however, to incorporate some version of this process into a future stdpopsim release.

Reviewer #2:

[…] Despite a little hiccup with python version mismatch prior to me successfully installing stdpopsim, I was able to successfully get the tutorials running within minutes after. I have one recommendation however for the tutorials – it would definitely help if the CLI versus python tutorials were kept separate. I found it a little confusing since they are all listed on the same page (https://stdpopsim.readthedocs.io/en/latest/tutorial.html). The simulations, testing models, calculating divergences, plotting ran without a hitch, and I am impressed and excited to play around with more models in coming days. Having also developed similar libraries/pipelines, I have also found it extremely useful for developers to provide some more detailed documentation/tutorials via Jupyter notebooks, or some similar platform. I did however notice that the authors have provided their analyses as Snakemake files in the interest of replication. I did not replicate their analyses, but I trust that the documentation for these analyses are detailed enough to aid readers/users in establishing similar analysis pipelines for stdpopsim simulated data.

Good idea – we have reorganized the Tutorials and added to them, and have also added in a basic stdpopsim API and CLI example in a Jupyter Notebook that can be accessed and used interactively via Binder and encourage users to try out the tutorials there. We have linked to the Binder in the README on the GitHub. See https://mybinder.org/v2/gh/popsimconsortium/stdpopsim/master?filepath=stdpopsim_example.ipynb

Reviewer #3:

[…] I have a few comments.

1) The authors mention that SLiM can be used as an alternative backend, which would presumably allow for simulations with selection. Although I don't think an extensive comparison of selection methods is necessary for this paper, it would be ideal if the authors can give some idea of how this would work (example command line, etc). There are also a myriad of methods for detecting/quantifying selection, and these simulations are not consistent either.

As mentioned in our response to the editor’s comment, we added a demonstration of the use of SLiM to the Results (see “Simulation engines” section and Appendix—figures 2 and 2). See also our response to related comments made by reviewers #1 and #2.

2) I like the inclusion of the "zigzag" history, as well as generic piecewise constant models and IM models (subsection “The Species Catalog”). I wonder if these could be included in a separate section (not organism specific) in the documentation and software (and then in Table 1). Right now the zigzag model is under humans in the catalog.

We considered both of these suggestions, and here’s why we’ve left it the way it is. Most of the columns of the table don’t apply to generic models, so it seems strange to include them there. The “zigzag” model could definitely be a generic model, and indeed was initially implemented as such. However, the reason we settled on the zigzag model is being defined as a human model (see discussion at https://github.com/popsim-consortium/stdpopsim/issues/106) is that its effective population size values are taken from (or at least inspired by) values inferred from human genomes: as implemented, it is not as “generic” as it might seem.

3) In the subsection “Use case: comparing methods of demographic inference”, the authors set up notation for the number of replicates (R), number of chromosomes (C), and sample size (n), but don't seem to use it afterward (or use it inconsistently). It would be helpful if all the figure legends and main text included this notation (I am guessing the number of replicates is 3 based on the images, but this should be clarified). The authors use N in the Materials and methods (i.e. subsection “Workflow for analysis of simulated data”) to refer to population size (which makes sense), but then also say "In all cases we set the sample size of the focal population to N = 50 chromosomes." For MSMC, the sample size was set to n=2,8 which suggests haploid samples, but the "Calculating coalescence rates" section says that n is the diploid sample size.

This is a good catch! We see why this is confusing. We introduced this notation

(R, C, and n) to make it completely clear in this paragraph what exactly was being simulated, but we feel that continuing to use this notation elsewhere would actually obscure things, since we don’t do any calculations with these quantities. We have changed the “n” in the Appendix to an “m”.

4) "Calculating coalescence rates" section needs a read through. Reword first sentence and add some citations (especially regarding computing p(t) and p(z,t)). It was unclear to me how the "mean coalescence times" were used (the rate was used to compute the ground truth over time). This section is also referred to as the Appendix in the main text.

We’ve expanded the first sentence substantially, and added a few more citations. We apologize that we don’t know of a paper to cite that does precisely the same calculations (but don’t doubt that such a paper exists), and instead refer only vaguely to “general theory of Markov chains” (but now with a citation). The section is now explicitly labeled “Appendix”. Hopefully this is more clear now.

Articles from eLife are provided here courtesy of eLife Sciences Publications, Ltd
-