Coalescense simulation download

We encode this model as follows:. After we set up our model, we use the DemographyDebugger to check our calculations. Then, we have a demographic event that changes the growth rate to 0. At generation , another event occurs, setting the growth rate for the population to 0. Then, the population size is constant at 20, from generation until the indefinite past.

A more complex example involving a three-population out-of-Africa human model is available in the online documentation. Up to this point we have assumed that all samples are taken at the present time.

However, msprime allows us to specify arbitrary sampling times and locations, allowing us to simulate for example ancient samples. One of the key innovations of msprime is that it makes simulation of the full coalescent with recombination possible at whole-chromosome scale. Adding recombination to a simulation is simple, requiring very minor changes to the methods given above. It is often useful to think of both sequence lengths and recombination rates as defined in units of base-pairs.

Note, however, that these are continuous values, so this correspondence should not be taken too literally. Note also that because msprime assumes an infinite sites mutation model the length parameter is not connected to the number of mutational sites. Thus any number of mutations can occur on a given sequence length, depending on the mutation rate specified. The result of this particular simulation is a tree sequence that contains 82 distinct trees.

Other replicate simulations with different random seeds will usually result in different numbers of trees. Up to this point we have focused on simulations that returned a single tree representing the genealogy of a sample. The inclusion of recombination, however, means that there may be more than one tree relating our samples. The TreeSequence object returned by msprime is a very concise and efficient representation of these highly correlated trees. To process the trees, we simply consider them one at a time, using the trees iterator.

The HapMap genetic map for chromosome 22 blue matches the density of breakpoints for a simulated chromosome green well. Although coordinates are specified in floating point values, msprime uses a discrete loci model when performing simulations.

However, the number of loci is configurable and it is possible to simulate a specific number of discrete loci. Here we simulate the history of two samples in a system with ten loci, each of length 1 with recombination rate of 1 between adjacent loci per generation. In the output, we see that the breakpoints between trees now occur exactly at the integer boundaries between these loci.

This shows that we can also simulate models of recombination with discrete loci in msprime , as well as the more standard continuous genome. In the previous section we showed how to run simulations in msprime , and how to construct population models and demographic histories. In this section we show how to process the results of simulations.

This is not a comprehensive review of the capabilities of the msprime Python API, but concentrates on some useful examples. To illustrate this, we consider a simulation of , samples of ten megabases from a simple two-population model with human-like parameters:.

We are often interested in finding the most recent common ancestor MRCA of a pair or many pairs of samples. For example, identity-by-descent IBD tracts are defined as contiguous stretches of genome in which the MRCA for a pair of samples is the same. Computing IBD segments for a pair of samples is very straightforward:. The distribution of the length of IBD segments for a pair of samples taken from the same or different populations. Of course, we would need to sample many such pairs of samples or longer sequences to get a reasonable approximation of the real distribution of block lengths.

The msprime API provides an extremely efficient way to count the number of samples that are beneath a particular node in a tree. This can be used, for example, to compute allele frequencies efficiently and is the basis for many of the fast algorithms in the API. In this example we iterate over all the trees in the tree sequence, and then iterate over all the sites in each tree. The underlying implementation ensures that this operation requires constant time, and so it is very efficient.

We see that such rare alleles are common. We reiterate that msprime currently generates mutations under the infinitely many sites model so that each mutation occurs at a unique site. Future versions of msprime or other software packages may produce tree sequences with back or recurrent mutations, where this simple approach will not work. To emphasize this point and to ensure that the above code chunk is not accidentally applied in such situations we have included an assert statement.

We use asserts in a similar way in later code chunks. A powerful feature of this sample-counting approach is that we can perform the same operation over an arbitrary subset of the samples. For example, suppose we wished to count the number of sites that are private to a specific population:. This example is very similar, except we provide an extra argument to ts.

Here we indicate that we are interested in tracking the set of samples within the population in question. Again, we iterate over all trees and over all sites within trees. If the total count is equal to the within-population count, we know that this mutation is private to the population. In some situations it is useful to analyze data for different subsets of the samples separately.

This is possible using the simplify method:. Here we extract the tree sequence representing the history of a tiny subset of the original samples, with IDs 1, 3, 5, and 7. The subset tree sequence contains all the genealogical information relevant to the subsamples, but no more.

Concretely, both coalescences that are not ancestral to the subsample and coalescences that predate the MRCA of the subsample are excluded.

Thus, the number of distinct trees is greatly reduced. By default, we also remove any sites that have no mutations within these subtrees i. Node IDs in the simplified tree sequence are not the same as in the original. Tree of a subset of the samples in a large simulation. Node IDs in the subset and full tree sequences are shown. While it is nearly always more efficient to work with mutations in terms of their context within the trees, it is sometimes more convenient to work with the allelic states of the samples.

This information is obtained in msprime using the variants iterator, which returns a Variant object for each site in the tree sequence. A Variant consists of: a a reference to the Site in question; b the alleles at this site the strings representing the actual states ; and c the genotypes representing the observed state for each sample. The genotypes are encoded in a NumPy array, such that variant.

The values in the genotypes array are therefore indexes into the alleles list. The ancestral state at a given site is guaranteed to be the first element in the alleles list, but no other assumptions about ordering of the alleles list should be made.

For biallelic sites, working with genotypes is straightforward as the genotypes array can only contain 0 and 1 values, which correspond to the ancestral and derived states, respectively. The genotypes values are returned as a NumPy array, and so the full NumPy library is available for efficient processing. Using the genotypes in this way is convenient, as complex patterns of back and recurrent mutations can be handled without difficulty.

This code is straightforward, as we simply iterate over all variants and count the number of one values in the genotypes array. Using the np. Generating all the genotypes for , samples at , sites, however, is an expensive operation and the overall calculation takes about 1. In the case of infinite sites mutations, we can recast this operation to use the efficient sample counting methods described in Section 3.

This approach is far more efficient, requiring less than 2 s to compute the same value. A powerful property of the tree sequence representation is that we can efficiently find the differences between adjacent trees. This is very useful when we have some value that we wish to compute that changes in a simple way between trees. Here we use it to keep a running track of the total branch length of our trees, without needing to perform a full traversal each time.

This function returns the total branch length value for each tree in the sequence as a NumPy array. It works by keeping track of the total branch length as we proceed from left to right, and storing this value in the output array for each tree.

Computing the current value for the total branch length is then simply a case of subtracting the branch lengths for all outgoing edges and adding the branch lengths for all incoming edges. This is extremely efficient because, after the first tree has been constructed there is at most four incoming and outgoing edges [ 23 ]. Thus, each tree transition costs constant time. In contrast, if we compute the total branch length by performing a full traversal for each tree, each tree transition is very costly when we have a large sample size.

In this example, computing the array of branch lengths using the incremental approach given here took 8 s. Computing the same array using the tree. This is because msprime currently implements this operation by a full traversal in Python; in future, this may change to using the algorithm given here.

The simplify method is useful here if you wish to export data from a subset of the simulated samples. However, it is worth noting that for large sample sizes, exporting genotype data may require a great deal of memory and take some time. One of the advantages of the msprime API is that we do not need to explicitly generate genotypes in order to compute many statistics of interest. In this section we show some examples of validating simple analytic predictions from coalescent theory using simulations.

The number of segregating sites is the total number of mutations that occurred in the history of the sample assuming the infinite sites mutation model. Since mutations happen as a Poisson process along the branches of the tree, what we are really interested in is the distribution of the total branch length of the tree. The results in this section are well-known classical results from coalescent theory; this section is intended as a demonstration of how to proceed when comparing analytic results to simulations.

We show some idiomatic examples for integrating with the state-of-the-art data analysis packages such as Pandas [ 30 ] and Seaborn [ 44 ]. All analytic predictions are taken from [ 43 ]. The first properties we are interested in are the mean and the variance of the total branch length of coalescent trees. We first create an array of the six different n values that we wish to simulate, and then create arrays to hold the results of the simulations.

Because we are running 10, replicates for each sample size, we allocate arrays to hold 60, values. This approach of storing the data in arrays is convenient because it allows us to use Pandas dataframes in an idiomatic fashion. We then iterate over all of our sample sizes and run 10, replicates of each. For each simulation, we simply store the sample size value and the total branch length in a Pandas dataframe.

This gives us access to many powerful data analysis tools including the Seaborn library, which we use for visualization here. After we have created our simulation data, we define our analytic predictions and plot the data. Comparisons of the distribution of simulated total branch lengths with analytic results. Ideally, we wish to capture the full distribution analytically.

In the following code chunk we define the analytic prediction for the total branch length distribution, and compare it with the simulated distribution for a sample of size The results are shown in Fig.

We can see an excellent agreement between the smoothed kernel density estimate produced by Seaborn and the theoretical prediction. Since we cannot directly observe branch lengths, we are usually more interested in mutations when working with data.

The mutation process is intimately related to the distribution of branch lengths, since mutations occur randomly along tree branches. One simple summary of the mutational process is the total number of segregating sites, that is, the number of sites at which we observe variation. We can obtain this very easily from simulations simply by specifying a mutation rate parameter.

Simulations of the number of segregating sites, and comparisons with analytic predictions. In the previous section we saw how to run simulations to generate trees under the assumptions of the single-locus coalescent and compare these with analytic predictions. This assumes that our data is not affected by recombination, which is often unrealistic. Here we show how to compute empirical distributions of equivalent quantities, and compare these with classical results from the literature.

Since analytic results for many quantities are generally unknown for the case of recombination along a linear sequence, we limit ourselves to the pairwise samples. In this code chunk we again run 10 4 replicate simulations for a range of input parameters, and store the results in a Pandas data frame. After defining our analytic predictions for the mean and variance of the number of segregating sites, we then plot the observed and predicted values in Fig.

Comparing the simulated results to analytic predictions we see excellent agreement. The mean number of segregating sites is not affected by recombination, but recombination does substantially reduce the variance. The analytical challenges of deriving likelihood functions even under highly idealized models of population structure and history have led to the development of likelihood-free inference methods, in particular Approximate Bayesian Computation ABC [ 2 ].

ABC approximates the posterior distribution of model parameters by drawing from simulations. Because of its flexibility ABC has become a standard inference tool in statistical population genetics see ref.

We will demonstrate how msprime can be used to set up an ABC inference by means of a simple toy example. We stress that this is meant as an illustration rather than an inference tool for practical use. However, given the flexibility of msprime , it should be relatively straightforward to implement more a realistic framework focused on specific inference applications. We assume that data for loci or sequence blocks these could be RAD loci in practice for a single diploid individual have been generated from each of two populations.

We would like to infer the amount of gene flow between the two populations. For the sake of simplicity, we will assume the simplest possible model of population structure; that is, two populations, of the same effective size exchanging migrants at a constant rate of m migrants per generation. Note that higher level population genetic summaries, e. Since msprime simulates rooted trees, the columns and rows of the unfolded jSFS correspond to the frequency of derived mutations in each population and the entries of the jSFS are simply mutation counts.

For example, for the first locus we have:. One could base inference on the bSFS [ 4 , 28 ], but we will for the sake of simplicity use a simpler and lossy summary of the data: the average jSFS across loci. To illustrate a simple ABC inference, we will focus on a single parameter of interest, the migration rate M. ABC measures the fit of data simulated under the prior to the observed data via a vector of summary statistics.

We assume an exponential distribution, a common choice of prior [ 13 ]. Here we run simulation replicates for each of the 10, m values drawn from the prior, giving a total of one million individual simulations. We use the multiprocessing module to distribute these computations over the available CPU cores. ABC results. This bias is in fact expected given that our prior is also strongly biased towards low m. We can check the effect the acceptance threshold on the inference and get a sense of the expected information about m using a cross-validation procedure: we repeat the inference on pseudo-observed data sets PODS simulated under a known truth.

Since we can re-use the same set of replicates simulated under the prior for inference, such cross-validation is computationally efficient. Figure 13 b shows the mean and the root mean square error RMSE of m estimates across PODS against the acceptance threshold and confirms that both the downward bias in m estimates and the associated RMSE increase with larger acceptance thresholds.

While this toy example illustrates the principle of ABC inference, sampling only a small fraction of simulations generated under the prior is clearly computationally inefficient and more efficient sampling strategies for ABC inference have been developed [ 2 ]. In practice, we are generally interested in fitting parameter-rich models and it would be straightforward to implement ABC inference for complex model of population structure and demography in msprime. In this chapter we have focused on the usage of msprime as a coalescent simulator, and illustrated its flexibility through concrete examples.

While many examples discuss how to create and run the simulations themselves, others are concerned with how we analyze the output of these simulations. We have shown particularly in Section 3 that these methods can be very efficient, allowing us to easily analyze chromosome scale data for hundreds of thousands of samples.

The data structures and APIs used in msprime are currently being developed to increase their generality and applicability. Recent work [ 11 , 24 ] has shown that forward-time simulations can also benefit from these methods. By recording all genealogical information for the simulated population in the form of a succinct tree sequence, we avoid the need to generate and carry forward neutral mutations; by definition, they do not affect the genealogies, and can therefore be placed on them afterwards.

Through the use of a well-documented interchange API and thoroughly specified data formats, forward-time simulators can output data that is compatible with the msprime API, and precisely the same techniques described here can be used to analyze the results. Thus, code written to analyze coalescent simulations can equally be applied to analyze forwards simulations. There is currently a great deal of activity from a growing community around msprime.

We plan to separate the tree sequence processing code from the simulator and create a library, provisionally known as tskit. This standalone library C and Python interfaces are planned will greatly facilitate integration with forwards-time simulators, allowing them to easily offload tree sequence processing to tskit. Algorithms for efficiently calculating statistics using the incremental techniques outlined in Section 3. Also in development are methods to estimate the tree sequence data structure from real data, which would allow us to use these efficient algorithms on observed as well as simulated data.

New features are being added to the msprime simulator also, with support for a discrete time Wright-Fisher model and a family of multiple-merger coalescent models in development. We hope that in the coming years a diverse ecosystem of tools and applications using these APIs and data structures will emerge.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. Skip to main content Skip to sections. This service is more advanced with JavaScript available.

Advertisement Hide. Coalescent Simulation with msprime. Open Access. First Online: 24 January Key words Population genetics Coalescent theory Simulation Python. Download protocol PDF. For example, here we simulate a history for a sample of three chromosomes: Open image in new window This code chunk illustrates the basic approach required to draw a tree in a Jupyter notebook.

The output of one random realization of this process is shown in Fig. The resulting tree has five nodes: nodes 0, 1, and 2 are leaves , and represent our samples. Node 3 is an internal node, and is the parent of 0 and 2. Node 4 is also an internal node, and is the root of the tree. In msprime , we always refer to nodes by their integer IDs and obtain information about these nodes by calling methods on the tree object.

For example, the code tree. Similarly, tree. Here we have two mutations, shown by the red squares. Mutations occur above a given node in the tree, and all samples beneath this node will inherit the mutation. The infinite sites mutations used here are simple binary mutations, that is, the ancestral state is 0 and the derived state is 1. Thus, if G is the genotype matrix, G [ j , k ] is the state of the k th sample at the j th site.

In our example above, the site 0 has a mutation over node 3, and site 1 has a mutation over node 1, and so we get the following matrix: Open image in new window.

We create our model by first making a list of two PopulationConfiguration objects. This results in samples being allocated sequentially to the populations when simulate is called: 0 and 1 are placed in population 0, and samples 2 and 3 are placed in population 1.

We then declare our migration matrix, which is asymmetric in this example. This is illustrated in Fig. Each node has been colored by its population red is population 0 and blue population 1.

Thus, the leaf nodes 0 and 1 are both from population 0, and 2 and 3 are both from population 2 as explained above. As we go up the tree, the first event that occurs is 2 and 3 coalescing in population 1, creating node 4. After this, 4 coalesces with node 0, which has at some point before this migrated into deme 1, creating node 5. Node 1 also migrates into deme 1, where it coalesces with 5. Because migration is asymmetric here, the MRCA of the four samples must occur within deme 1.

This code produces the plot in Fig. We can see that node 0 experienced very few migration events before it ended up in deme 2, where it coalesced with 4 which never migrated.

Node 2, on the other hand, migrated 30 times before it finally coalesced with 7 in deme 0. Note that there are many more migration events than nodes here, implying that most migration events are not identifiable from a genealogy in real data [ 38 ]. The samples 0 and 1, and 2 and 3 coalesce quickly within their own populations. However, because the migration rate between the populations is zero these lineages are isolated and would never coalesce without some change in demography.

The migration rate change event happens at time 20, resulting in node 5 migrating to deme 1 soon afterwards. The lineages then coalesce at time Epoch: 0 -- Events generation Epoch: All of the trees that we previously considered had leaf nodes at time zero. In this case, the samples 0, 1, and 2 are taken at time 0 in population 0, but node 3 is sampled at time 0.

Note that in this case we used the samples parameter to simulate to specify our samples. This is the most general approach to assigning samples, and allows samples to be assigned to arbitrary populations and at arbitrary times.

This code generates the plot in Fig. Mar 19, Feb 14, Feb 11, Feb 9, Feb 5, Feb 4, Feb 3, Jan 10, Jan 8, Jan 7, Jan 5, Jan 3, Dec 10, Dec 6, Dec 5, Nov 13, Nov 8, Oct 3, Oct 2, Sep 29, Sep 12, Sep 11, Sep 5, Sep 3, Aug 24, Aug 23, Aug 22, Aug 19, Aug 13, Jul 23, Jul 22, Jul 15, Jun 29, Jun 25, Jun 24, Jun 22, Jun 21, Jun 20, Download the file for your platform.

If you're not sure which to choose, learn more about installing packages. Warning Some features may not work without JavaScript. Please try enabling it if you encounter problems. Search PyPI Search. Latest version Released: Dec 9, Navigation Project description Release history Download files. Project links Homepage Docs Bug Tracker. Maintainers thompsonsed. Project description Project details Release history Download files Project description A package for coalescence-based spatially explicit neutral ecology simulations Introduction pycoalescence is a Python package for spatially explicit coalescence neutral simulations.

Installation Usage of conda is recommended to aid handling installation of dependencies. The SQLite library available here comes included with Python. Numerical Python numpy package pip install numpy. The proj library for converting between coordinate systems. Recommended Scipy package for generating fragmented landscapes pip install scipy. Matplotlib package for plotting fragmented landscapes pip install matplotlib. Optional For work involving large csv files, the fast-cpp-csv-parser by Ben Strasser, available here can be used.