1 SLiM recipes

SLiM recipes can be run with wrapper pgip-slim to generate multiple simulations.

Code

pgip-slim --help
## Usage: pgip-slim [OPTIONS] SLIM
## 
##   pgip slim simulator CLI.
## 
##   Wrapper that runs slim on SLIM controlfile.
## 
## Options:
##   -o, --outdir PATH               Output directory
##   -N, --population-size INTEGER   Population size
##   -r, --recombination_rate FLOAT  Recombination rate
##   -m, --mutation_rate FLOAT       Mutation rate
##   -l, --sequence_length INTEGER   Sequence length
##   -n, --repetitions INTEGER RANGE
##                                   Number of repetitions  [x>=1]
##   --seed TEXT                     Random seed
##   --threads INTEGER RANGE         Number of parallel threads to run  [x>=1]
##   --prefix TEXT                   File output prefix
##   --no-recapitate                 Don't do recapitation
##   --debug                         Print debugging info
##   --help                          Show this message and exit.

As an example, the following command can be used to generate 10 simulations using the selective sweep recipe.

Code

pgip-slim --seed 42 -n 10 -r 1e-6 -m 1e-7 --threads 12 \
    selective_sweep.slim -l 1000000 --outdir results/slim

1.1 Selective sweep

slim code to simulate a selective sweep at a locus on a 1Mbp chromosome, with mutation rate \(\mu=10^{-8}\), recombination rate \(\rho=10^{-8}\), and \(N_e=5,000\).

initialize() {
  defineConstant("Nref", 5e3);
  if (!exists("seed"))
    defineConstant("seed", getSeed());
  defineConstant("SimID", seed);
  if (!exists("N"))
    defineConstant("N", Nref);
  if (!exists("outdir"))
    defineConstant("outdir", tempdir());
  if (!exists("mu"))
    defineConstant("mu", 1e-8);
  if (!exists("rho"))
    defineConstant("rho", 1e-8);
  if (!exists("seqlength"))
    defineConstant("seqlength", 1e6);
  if (!exists("position"))
    defineConstant("position", asInteger(seqlength / 2));
  if (!exists("outfile"))
    defineConstant("outfile", outdir + "slim_" + SimID + ".rho_" + rho + ".mu_" + mu + ".N_" + N + ".trees");

  setSeed(seed);

  initializeTreeSeq(simplificationRatio=INF);
  initializeMutationRate(0);
  initializeMutationType("m1", 0.5, "f", 0.0);
  initializeMutationType("m2", 1.0, "f", 0.5);  // introduced mutation
  initializeGenomicElementType("g1", m1, 1.0);
  initializeGenomicElement(g1, 0, seqlength - 1);
  initializeRecombinationRate(rho);
}
1 early() {
  // Chapter 14.9: SLiM models diploid individuals that contain two
  // haploid genomes; this is, at present, a design constraint in SLiM
  // that cannot be modified. IOW: N individuals, 2*N genomes
  sim.addSubpop("p1", N);
}
1000 late() {
  // introduce the sweep mutation
  target = sample(p1.genomes, 1);
  target.addNewDrawnMutation(m2, position);
  // save the state of the simulation
  sim.treeSeqOutput(tempdir() + "slim_" + SimID + ".trees");
}
1000:100000 late() {
  if (sim.countOfMutationsOfType(m2) == 0)
    {
      fixed = (sum(sim.substitutions.mutationType == m2) == 1);

      if (fixed)
    {
      cat(SimID + ": FIXED\n");
      cat("Writing output file : " + outfile + "\n");
      sim.treeSeqOutput(outfile);
      sim.simulationFinished();
    }
      else
    {
      cat(SimID + ": LOST - RESTARTING\n");

      // go back to tick 1000
      sim.readFromPopulationFile(tempdir() + "slim_" + SimID + ".trees");
      setSeed(rdunif(1, 0, asInteger(2^62) - 1));
    }
    }
}

Recipe to simulate a selective sweep.

1.2 Neutral simulation

slim code to simulate a neutrally evolving locus on a 100kbp chromosome, with mutation rate \(\mu=10^{-8}\), recombination rate \(\rho=10^{-8}\), and \(N_e=5,000\).

initialize() {
  defineConstant("Nref", 5e3);
  if (!exists("seed"))
    defineConstant("seed", getSeed());
  defineConstant("SimID", seed);
  if (!exists("N"))
    defineConstant("N", Nref);
  if (!exists("outdir"))
    defineConstant("outdir", tempdir());
  if (!exists("mu"))
    defineConstant("mu", 1e-7);
  if (!exists("rho"))
    defineConstant("rho", 1e-8);
  if (!exists("seqlength"))
    defineConstant("seqlength", 1e5);
  if (!exists("outfile"))
    defineConstant("outfile", outdir + "slim_" + SimID + ".rho_" + rho + ".mu_" + mu + ".N_" + N + ".trees");
  setSeed(seed);

  initializeTreeSeq(simplificationRatio=INF);
  initializeMutationRate(0);
  initializeMutationType("m1", 0.5, "f", 0.0);
  initializeGenomicElementType("g1", m1, 1.0);
  initializeGenomicElement(g1, 0, seqlength - 1);
  initializeRecombinationRate(rho);
}
1 early() {
  sim.addSubpop("p1", N);
}
10000 early() {
      cat("Writing output file : " + outfile + "\n");
      sim.treeSeqOutput(outfile);
      sim.simulationFinished();
}

Recipe to simulate a neutrally evolving locus.