Running a WF simulation in shadie¶

We will start with a simple example using the Wright-Fisher model.

shadie provides a standard Wright-Fisher Model, as well as some alternatives (e.g. the Moran model). These models are implemented as SLiM nonWF models.

Running the Wright-Fisher¶

The classic Wright-Fisher Model makes the following assumptions:

1. Discrete, non-overlapping generations
2. Random mating between hermaphroditic individuals
3. Constant population size (of diploid individuals; 2N haplotypes)
4. No fitness effects; simulation is neutral (no selection)
5. No recombination
   • Note that the default recombination rate in shadie is NOT 0.0, so to meet this requirement you must be sure to set it in the initialize function
6. Infinite sites

shadie uses context manager to check model parameters before running a simulation. The syntax requires the user to define a model using a with statement and then defining all parameters. The model needs a .initialize() function, which takes general parameters such as the chromosome, mutation rate, recombination rate, and file names. You can also specify what kind of reproduction the model will use with .reproduction. followed by the appropriate function (a list of options can be seen in Jupyter Notebooks by pressing tab after typing .reproduction..

See the example below:

In [2]:

import shadie
default_chrom = shadie.chromosome.default()

with shadie.Model() as WF_model:
    WF_model.initialize(
        chromosome=default_chrom, 
        recomb_rate=0.0,
        sim_time=1000,
        file_out="WF.trees",
    )
    WF_model.reproduction.wright_fisher(pop_size=1000)

import shadie default_chrom = shadie.chromosome.default() with shadie.Model() as WF_model: WF_model.initialize( chromosome=default_chrom, recomb_rate=0.0, sim_time=1000, file_out="WF.trees", ) WF_model.reproduction.wright_fisher(pop_size=1000)

Now, your model exists as the model object you defined, in this case WF_model. You can print the model script by running:

In [3]:

print(WF_model.script)

print(WF_model.script)

initialize() {

  // model type
  initializeSLiMModelType("nonWF");

  // config
  initializeRecombinationRate(0.0, 10000);
  initializeMutationRate(1e-08);
  initializeTreeSeq(simplificationInterval=NULL);

  // MutationType init
  initializeMutationType('m2', 0.1, 'g', -3.0, 1.5);
  initializeMutationType('m3', 0.8, 'e', 0.04);
  c(m3, m2).haploidDominanceCoeff = 1.0;
  c(m3, m2).convertToSubstitution = T;

  

  // ElementType init
  initializeGenomicElementType('g1', c(m2,m3), c(8,0.1));
  initializeGenomicElementType('g2', m2, 1);

  // Chromosome (GenomicElement init)
  types = rep(g1, asInteger(floor(1999/3)));
  starts = 2001 + seqLen(integerDiv(1999, 3)) * 3;
  ends = starts + 1;
  initializeGenomicElement(types, starts, ends);

  types = rep(g2, asInteger(floor(1999/3)));
  starts = 4001 + seqLen(integerDiv(1999, 3)) * 3;
  ends = starts + 1;
  initializeGenomicElement(types, starts, ends);

  types = rep(g1, asInteger(floor(1999/3)));
  starts = 6001 + seqLen(integerDiv(1999, 3)) * 3;
  ends = starts + 1;
  initializeGenomicElement(types, starts, ends);


  // constants (Population sizes and others)
  defineConstant('K', 1000);

  // globals (metadata dictionary)
  defineGlobal('METADATA', Dictionary('file_in', 'None', 'file_out', 'WF.trees', 'mutation_rate', '1e-08', 'recomb_rate', '0.0', 'model', 'shadie WF', 'length', '1000', 'spo_pop_size', '1000', 'gam_pop_size', 'NA', 'spo_mutation_rate', '1e-08', 'recombination_rate', '0.0', 'selection', 'none', 'gens_per_lifecycle', '1'));

  // extra scripts (Optional)
  
}


// executes before offspring are generated
// define starting diploid population.

1 first() {
    sim.addSubpop('p1', K);
}


// generates offspring
// WF model with no selection; random mating
reproduction(p1 ) {
    subpop.addCrossed(individual, subpop.sampleIndividuals(1));
}


// implements survival adjustments
// non-overlapping generations
survival() {
    return (individual.age == 0);
}


// executes after selection occurs
// end of sim; save .trees file

1001 late() {
    sim.treeSeqRememberIndividuals(sim.subpopulations.individuals);
    sim.treeSeqOutput('WF.trees', metadata = METADATA);
}

If desired, this script can be copied and pasted into the SLiMgui, or saved as a file and run using the command line.

To run the model in shadie, simply call the run() function:

In [4]:

WF_model.run()

WF_model.run()

The output of the SLiM simulations, a .trees file with the indicated file_out name, should now be located in your working directory.

Access the Tree Sequence¶

You can inspect the output of the simulation using the postsim module.

In [5]:

postsim = shadie.postsim.OneSim("WF.trees", default_chrom)
postsim.tree_sequence

postsim = shadie.postsim.OneSim("WF.trees", default_chrom) postsim.tree_sequence

🌿 INFO | one_sim.py | shadie assumes sims were run without a burn-in and without neutral mutations and will recapitate and mutate (add neutral mutations) any loaded sims by default. Current settings:
Recapitate: True
Add neutral mutations: True

🌿 INFO | one_sim.py | Simulation currently has 1 existing mutations of type(s) {2}.
🌿 INFO | one_sim.py | Mutating using mutation rate: 1e-08. Keeping 6 existing mutations of type(s) {0, 2}.

Out[5]:

Tree Sequence
Trees	1
Sequence Length	10001.0
Time Units	ticks
Sample Nodes	7513
Total Size	783.9 KiB
Metadata	dict SLiM: dict cycle: 1001 file_version: 0.8 model_type: nonWF name: sim nucleotide_based: False separate_sexes: False spatial_dimensionality: spatial_periodicity: stage: late tick: 1001 user_metadata: dict file_in: list None file_out: list WF.trees gam_pop_size: list NA gens_per_lifecycle: list 1 length: list 1000 model: list shadie WF mutation_rate: list 1e-08 recomb_rate: list 0.0 recombination_rate: list 0.0 selection: list none spo_mutation_rate: list 1e-08 spo_pop_size: list 1000

Table	Rows	Size	Has Metadata
Edges	7520	235.0 KiB
Individuals	2000	197.1 KiB	✅
Migrations	0	8 Bytes
Mutations	6	1.5 KiB	✅
Nodes	7521	279.7 KiB	✅
Populations	2	2.4 KiB	✅
Provenances	4	7.0 KiB
Sites	6	160 Bytes

Plot the tree sequence using the draw_tree_sequence() function. This function samples 10 haplotypes by default:

In [6]:

postsim.draw_tree_sequence();

postsim.draw_tree_sequence();

NOTE: The default Wright Fisher model in shadie is the classic model, which has no selection. Even if non-neutral mutation types are defined and occur, they do not affect the simulation outcome.

Breaking the Assumptions of Wright-Fisher¶

There is a parameter in the shadie Wright-Fisher models, selection, that controls what kind of selection can occur in the model. Adding selection will break the neutrality assumption of the Wright-Fisher.

There are two options for adding selection:

• soft selection: fitness affects mating success only
• hard selection: fitness affects survival only

Soft and hard selection are implemented in the examples below. These examples also break the no recombination assumption, simply by passing a non-zero value to the recomb_rate parameter.

In [4]:

with shadie.Model() as WF_soft_model:
    WF_soft_model.initialize(
        chromosome=default_chrom, 
        recomb_rate=1e-9,
        file_out="WF_soft.trees",
    )
    WF_soft_model.reproduction.wright_fisher(
        pop_size=1000,
        selection= "soft",
    )

with shadie.Model() as WF_soft_model: WF_soft_model.initialize( chromosome=default_chrom, recomb_rate=1e-9, file_out="WF_soft.trees", ) WF_soft_model.reproduction.wright_fisher( pop_size=1000, selection= "soft", )

In [5]:

with shadie.Model() as WF_hard_model:
    WF_hard_model.initialize(
        chromosome=default_chrom, 
        recomb_rate=1e-9,
        file_out="WF_hard.trees",
    )
    WF_hard_model.reproduction.wright_fisher(
        pop_size=1000,
        selection= "hard",
    )

with shadie.Model() as WF_hard_model: WF_hard_model.initialize( chromosome=default_chrom, recomb_rate=1e-9, file_out="WF_hard.trees", ) WF_hard_model.reproduction.wright_fisher( pop_size=1000, selection= "hard", )

WF models incorporating selection may be more appropriate to use as control simulations to compare against simulations with complex life cycle models.

Alternative Models in shadie based on Wright-Fisher¶

Because the life cycle models in shadie incorporate alternations of generations and both haploid and diploid individuals, there are further modifications to the Wright-Fisher model that may be useful. The same three options for selection ("none", "soft" and "hard) are available for each of them.

Moran Model¶

The Moran model is a simple alternative to Wright-Fisher that allows overlapping generations.

In [6]:

with shadie.Model() as moran_model:
    moran_model.initialize(
        chromosome=default_chrom, 
        recomb_rate= 0.0,
        file_out="moran.trees",
    )
    moran_model.reproduction.moran(pop_size=1000)

with shadie.Model() as moran_model: moran_model.initialize( chromosome=default_chrom, recomb_rate= 0.0, file_out="moran.trees", ) moran_model.reproduction.moran(pop_size=1000)

Haploid Models¶

The default Wright-Fisher model usually assumes diploid individuals. There are two haploid Wright-Fisher models in shadie:

• clonal: A single parent is chosen N times with replacement; each time a parents is chosen it produces one offspring
• sexual: The parent genomes undergo reocmbination and one offspring is produced per pair. This requires recombination rate != 0.0 Both models have all three selection options: "none", "soft", and "hard".

In [8]:

with shadie.Model() as WF_1N_clonal_model:
    WF_1N_clonal_model.initialize(
        chromosome=default_chrom, 
        file_out="1N_clonal.trees",
    )
    WF_1N_clonal_model.reproduction.wright_fisher_haploid_clonal(pop_size=1000)

with shadie.Model() as WF_1N_clonal_model: WF_1N_clonal_model.initialize( chromosome=default_chrom, file_out="1N_clonal.trees", ) WF_1N_clonal_model.reproduction.wright_fisher_haploid_clonal(pop_size=1000)

In [3]:

with shadie.Model() as WF_1N_sexual_model:
    WF_1N_sexual_model.initialize(
        chromosome=default_chrom,
        recomb_rate=1e-9,
        file_out="1N_sexual.trees",
    )
    WF_1N_sexual_model.reproduction.wright_fisher_haploid_sexual(
        pop_size=1000,
        selection="soft"
    )

with shadie.Model() as WF_1N_sexual_model: WF_1N_sexual_model.initialize( chromosome=default_chrom, recomb_rate=1e-9, file_out="1N_sexual.trees", ) WF_1N_sexual_model.reproduction.wright_fisher_haploid_sexual( pop_size=1000, selection="soft" )

Alternation-of-Generations "Wright-Fisher"¶

This is a highly altered Wright-Fisher model that incorporates alternation of generations, which may be useful as a control or comparison to more complex models with alternation-of-generations.

In [3]:

with shadie.Model() as WF_altgen_model:
    WF_altgen_model.initialize(
        chromosome=default_chrom, 
        recomb_rate= 0.0,
        file_out="moran.trees",
    )
    WF_altgen_model.reproduction.wright_fisher_altgen(
        spo_pop_size=1000,
        gam_pop_size=1000,
        selection="hard"
    )

with shadie.Model() as WF_altgen_model: WF_altgen_model.initialize( chromosome=default_chrom, recomb_rate= 0.0, file_out="moran.trees", ) WF_altgen_model.reproduction.wright_fisher_altgen( spo_pop_size=1000, gam_pop_size=1000, selection="hard" )

In [4]:

print(WF_altgen_model.script)

print(WF_altgen_model.script)

initialize() {

  // model type
  initializeSLiMModelType("nonWF");

  // config
  initializeRecombinationRate(0.0, 10000);
  initializeMutationRate(1e-08);
  initializeTreeSeq(simplificationInterval=NULL);

  // MutationType init
  initializeMutationType('m2', 0.1, 'g', -3.0, 1.5);
  initializeMutationType('m3', 0.8, 'e', 0.04);
  c(m3, m2).haploidDominanceCoeff = 1.0;
  c(m3, m2).convertToSubstitution = T;

  

  // ElementType init
  initializeGenomicElementType('g1', c(m2,m3), c(8,0.1));
  initializeGenomicElementType('g2', m2, 1);

  // Chromosome (GenomicElement init)
  types = rep(g1, asInteger(floor(1999/3)));
  starts = 2001 + seqLen(integerDiv(1999, 3)) * 3;
  ends = starts + 1;
  initializeGenomicElement(types, starts, ends);

  types = rep(g2, asInteger(floor(1999/3)));
  starts = 4001 + seqLen(integerDiv(1999, 3)) * 3;
  ends = starts + 1;
  initializeGenomicElement(types, starts, ends);

  types = rep(g1, asInteger(floor(1999/3)));
  starts = 6001 + seqLen(integerDiv(1999, 3)) * 3;
  ends = starts + 1;
  initializeGenomicElement(types, starts, ends);


  // constants (Population sizes and others)
  defineConstant('SPO_POP_SIZE', 1000);
  defineConstant('GAM_POP_SIZE', 1000);

  // globals (metadata dictionary)
  defineGlobal('METADATA', Dictionary('model', 'shadie haploid clonal WF', 'length', '1000', 'spo_pop_size', '1000', 'gam_pop_size', '1000', 'spo_mutation_rate', '1e-08', 'selection', 'hard', 'recombination_rate', '0.0'));

  // extra scripts (Optional)
  
}


// executes before offspring are generated
// define starting haploid population.

1 first() {
    sim.addSubpop('p1', SPO_POP_SIZE);
sim.addSubpop('p0', 0);
}


// generates offspring
// clonal random mating; mating success weighted by fitness.
reproduction(p1 ) {
    // parents are chosen randomly
    fitness = p1.cachedFitness(NULL);
    parents = sample(p1.individuals, GAM_POP_SIZE, replace=T);
    for (i in seqLen(GAM_POP_SIZE)){
        breaks = sim.chromosome.drawBreakpoints(parents[i]);
        p0.addRecombinant(parents.genome1[i], parents.genome2[i], breaks, NULL, NULL, NULL);
    }
    self.active = 0;
}


// generates offspring
// clonal random mating; mating success weighted by fitness.
reproduction(p0 ) {
    // parents are chosen randomly
    fitness = p0.cachedFitness(NULL);
    parents1 = sample(p0.individuals, SPO_POP_SIZE, replace=T);
    parents2 = sample(p0.individuals, SPO_POP_SIZE, replace=T);
    for (i in seqLen(SPO_POP_SIZE))
        p1.addRecombinant(parents1.genome1[i], NULL, NULL, parents2.genome1[i], NULL, NULL);
    self.active = 0;
}


// executes after offspring are generated
// Fitness scaling for hard selection

early() {
    // diploids (p1) just generated haploid into p0
    if (community.tick % 2 == 0) {
        // fitness affects haploid survival
        p0.fitnessScaling = GAM_POP_SIZE / p0.individualCount;
        }
    // haploids (p0) just generated diploids into p1
    else {
        //fitness affects diploid survival
        p1.fitnessScaling = SPO_POP_SIZE / p1.individualCount;
    }
}


// implements survival adjustments
// non-overlapping generations
survival() {
    return (individual.age == 0);
}


// executes after selection occurs
// end of sim; save .trees file

2001 late() {
    sim.treeSeqRememberIndividuals(sim.subpopulations.individuals);
    sim.treeSeqOutput('moran.trees', metadata = METADATA);
}