Simulation

Don’t know how to program in Eidos? Don’t worry! You can use our script generator to create SLiM scripts for use with slime! Hooray for you.

Contents of this section:

Quickstart

SLiM scripts can be created from slime.RecentHistory objects. The following code shows a simple, minimal demographic history of admixture:

>>> ref_pops = [slime.PopulationConfiguration(initial_size=100),
...      slime.PopulationConfiguration(initial_size=100)]
>>> adm_pop = slime.PopulationConfiguration(initial_size=50)
>>> script = slime.RecentHistory(final_gen=100, chrom_length=1e7,
...    recombination=1e-8, reference_configs=ref_pops, adm_configs=adm_pop,
...    prop=[0.3, 0.7])

For example, in the code below, I specify two ancestral populations, each with population size 100 at the time of admixture. Our admixed population has initial population size 50, and is initially created with a single-pulse of admixture between the two reference populations that occurred 100 generations ago. Each reference population contributes migrants to the admixed population with proportions of 0.3 and 0.7 respectively. The simulated chromosome contains 1e7 bases, and recombination occurs at a uniform rate of 1e-8 per base per generation.

At minimum, your history of admixture should involve:

  • some number of discrete ancestral populations
  • a present-day admixed population
  • the age of the admixed population, in generations
  • an admixture proportion
  • the number of bases in the genomic segment being simulated
  • a recombination rate

The sections below explain how more complicated admixtures can be represented with slime.

Generating a SLiM script

Use the dump_script() method:

>>> script.dump_script('recent-history.slim')

You should then have a file called recent-history.slim in your root directory.

You can also view the SLiM script corresponding to your RecentHistory object at any stage using the print_script() method:

>>> script.print_script()

Recombination

Recombination rates are specified per base, per generation. Can be specified for a slime.RecentHistory object via the recombination parameter.

Uniform recombination

Pass the rate to your slime.RecentHistory object via the recombination parameter.

For example, the following code indicates a recombination rate of 1e-8 per base per gen:

>>> script = slime.RecentHistory(final_gen=100, chrom_length=1e7,
...    recombination=1e-8, reference_configs=ref_pops, adm_configs=adm_pop,
...    prop=[0.3, 0.7])

Non-uniform recombination

A path to a space-delimited text file containing positions and rates can also be passed to the recombination parameter. The first column of the file should contain (ordered) integers representing base positions. The second column should contain recombination rates. The rate in the ‘n’th row is the per-base recombination rate for all bases between the base positions given in the `n`th (inclusive) and `n+1`th (not inclusive) row of the file.

Special rows: The position in the final row must be the length of the simulated sequence + n, but note that this row isn’t actually used (it’s a weird quirk of SLiM)… The first row is a header and will be ignored.

Here is an example of a valid recombination map for a sequence of length 10, with a recombination hotspot between bases 6 and 10:

position rate(cM/Mb)
1 0.1
6 0.5
11 0.1

Mutations

To initialise a slime simulation with mutations, you make a MutationTypes object and feed it into your RecentHistory object via the mutations parameter.

Note

If you are appending your SLiM simulation to a msprime simulation, there’s no need to code neutral mutations into your RecentHistory object, as they are more easily generated later using recapitate.

Fixed selection coefficients

Say we wish to specify 2 different types of mutations to occur:

  1. Deleterious mutations with selection coefficient 0.2. These will occur with relative proportion 0.2, and have a dominance coefficient of 0.7.
  2. Weakly beneficial mutations with selection coefficient 0.6. These will occur with a relative proportion 0.8, and have a dominance coefficient of 0.4.

The following code initialises a MutationTypes object that specifies these types of mutations to occur in our forward-in-time simulation.

>>> muts = slime.MutationTypes(mutation_rate=.005,
...     selection_coeffs=[0.2,0.6], proportions=[0.2, 0.8],
...     dominance_coeffs=[0.7, 0.4])

We’ll then use this MutationTypes object to initialize a RecentHistory object via the mutations parameter.

>>> config = slime.PopulationConfiguration(initial_size=10)
>>> script = slime.RecentHistory(final_gen=20, chrom_length=100,
...     reference_configs=[config, config], adm_configs=config,
...     prop=[0.3,0.7], mutations=muts)

Random selection coefficients

Coming soon.

Population size changes

Constant growth

To specify a constant rate of population growth for any of the populations, pass an exponential growth parameter to the corresponding slime.PopulationConfiguration object.

For instance, suppose our admixed population grows at a rate of exp(0.1) in each new generation:

>>> adm_pop = slime.PopulationConfiguration(initial_size=50, growth_rate=0.1)

Instant population size changes

Create a msprime.PopulationParametersChange object and add it to your RecentHistory object using the add_size_change method.

For instance, suppose our admixed population has size 25 from generation 11 onwards:

>>> size_change = msprime.PopulationParametersChange(time=11, initial_size=25,
...             population_id=2)
>>> script = slime.RecentHistory(final_gen=gens, chrom_length=1e7,
...    recombination=1e-8, reference_configs=ref_pops, adm_configs=adm_pop,
...    prop = [0.5, 0.7])
>>> script.add_size_change(PopulationParametersChange=size_change)

Read more about msprime.PopulationParametersChange objects here.

Changes in growth rate

Note

Fill this out.

Migration

Constant migration

By default, slime assumes that the admixed populations is created by a single pulse of migration from the reference populations. This can be changed by specifying migration rates from each reference population using the add_migration_rate method:

>>> config = slime.PopulationConfiguration(initial_size=100)
>>> script = slime.RecentHistory(final_gen=20, chrom_length=10,
...         recombination=.01, reference_configs=[config, config],
...         adm_configs=config, prop=[0.3,0.7])
>>> script.add_migration_rate(rates=[.01,.02])

To change migration rates at a particular time in the simulation, use the optional argument time. For example, adding this line to the above code will stop all migration into the admixed population at the end of the 5th generation:

>>> script.add_migration_rate(rates=[0, 0], time=(5, 'late'))

Mass migration

You can also add one-off migrations into the admixed population using the add_mass_migration method.

>>> script.add_mass_migration(prop=[.1, .1], time=(9, 'late'))

Note

At the moment, the add_mass_migration method will set migration rates in the following generation to 0 (not to what they were in earlier generations of the simulation). You can set migration rates in generation n+1 to what they were before using add_migration_rate.

API

class scriptgenerator.EventList(start_time, end_time=None)

Holds a start and and end time, and a list of events that happens between those times.

class scriptgenerator.MutationTypes(mutation_rate, selection_coeffs, proportions, dominance_coeffs=None)

Holds information about mutations to put into the simulation.

class scriptgenerator.PopulationConfiguration(initial_size, sample_size=None, growth_rate=0)

Holds information about the populations in the simulation. The analogue of msprime.PopulationConfiguration objects.

Variables:
  • initial_size (int) – The number of haploids in the population at the start of the simulation.
  • sample_size – Think – IS THIS NEEDED???
  • growth_rate (float) – The initial growth rate. (GIVEN AS EXPONENT??)
class scriptgenerator.RecentHistory(final_gen, chrom_length, reference_configs, adm_configs, prop, recombination=0, ref_labels=None, adm_label=None, reference_populations=None, admixed_population=None, mutations=None, outfile='recent-history.trees', model_type='WF', scriptfile='recent-history.slim')

An object holding demographic history that can be converted into a SLiM script.

Variables:
  • final_gen (int) – The number of generations to simulate.
  • chrom_length (int) – The number of bases in the simulated segment.
  • reference_configs (slime.PopulationConfiguration) – Configurations of the ancestral populations.
  • adm_configs (slime.PopulationConfiguration) – Configurations of the admixed population.
  • prop (list, dtype=float) – The initial proportion of admixture from each reference population.
  • recombination (float or recombination map (FORMAT?)) – Recombination rate or recombination map.
  • ref_labels
  • adm_label
  • reference_populations
  • admixed_population
  • mutations
  • outfile – The filename to which to write the outputted tree sequence.
  • model_type – To be removed.
  • scriptfile – The filename to which to write the SLiM script.
add_continuous_process(new_range, event, early=True)

Adds an event that occurs over an interval of generations. Breaks up existing events if necessary.

add_mass_migration(prop, time)

Specifies a mass migration event from reference populations into admixed population at a specified time. This is like changing the migration rate, but the migration rate is changed to 0 afterwards. # Later: would be useful if it is changed back to original migration rate?

add_migration_rate(rates, time=(2, 'late'))

Specifies rates of migration from references to admixed population.

add_new_single_time(time)

Creates a new code block for the supplied single time.

add_population_parameters_change(paramChanges)

Takes a list of msprime.PopulationParametersChange objects and ‘translates’ them into SLiM commands. https://msprime.readthedocs.io/en/stable/api.html#demographic-events There are two types of changes: exponential growth rate changes and instantaneous population size changes.

find_event_index(time=None, start=False, initialization=False)

Finds the index of the script at which a new event should be inserted.

find_time_index(time)

Finds the index of the script at which a new generation and time should be inserted.