BasicsFFPopSim can be imported into Python as a module, as explained in the FFPopSim Documentation. To set up a population and evolve it for a few generation, nothing more than the following is needed:
This script is included in the scripts folder under the name ffpopsim_intro.py. To execute the script in an interactive mode, open IPython by typing
import FFPopSim #load the FFPopSim library L = 256 # simulate 256 loci pop = h.haploid_highd(L) # instance of haploid_highd with L loci pop.carrying_capacity = 5000 # set the average population size to 5000 pop.outcrossing_rate = 1 # make the species obligate outcrossing pop.crossover_rate = 0.02 / pop.L # set the segment map to 0.02 pop.mutation_rate = 0.001 # per locus mutation rate equal to 0.1/N pop.evolve(100) # evolve for 100 generations
on the command line. To run the script type
where "path_to_script" has to be replaced with the location of the script file. The population can now be inspected interactively in the Python shell. The flexibility of a scripting language allows the user to put together more complicated scenarios with little effort. In particular, the simulation can be augmented by plotting and analysis. This requires two additional libraries numpy and matplotlib. The following example sets up a sexual population with a few loci that initially are at frequency 0.5. The script tracks the allele frequencies and plots how they change through genetic drift.
The result from the script should look more or less like this: All allele frequencies start at 0.5, which is how the population was initialized. They then change over the time through genetic drift. Open the script in your favorite text editor, change parameters, and rerun.
import numpy as np #numerical library from matplotlib import pyplot as plt #plotting import matplotlib.cm as cm #colormaps import FFPopSim # specify parameters L = 256 # simulate 256 loci # set up population pop = FFPopSim.haploid_highd(L) #instance of haploid_highd with L loci pop.carrying_capacity = 2000 # set the average population size to 2000 pop.outcrossing_rate = 1 # make the species obligate outcrossing pop.crossover_rate = 0.02 / pop.L # set the map length of the segment to 0.02 pop.mutation_rate = 0.1 / pop.carrying_capacity # per locus mutation rate equal to 0.1/N # initialize the population in linkage equilibrium with the specified allele frequencies initial_allele_frequencies = 0.5*np.ones(pop.L) # initial allele frequencies at 1/2 pop.set_allele_frequencies(initial_allele_frequencies, pop.carrying_capacity) # save initial allele frequencies and starting time point allele_frequencies = [pop.get_allele_frequencies()] tp = [pop.generation] # evolve and track the allele frequencies while pop.generation < 500: pop.evolve(10) # save allele frequencies and time allele_frequencies.append(pop.get_allele_frequencies()) tp.append(pop.generation) # convert to an array to enable slicing allele_frequencies = np.array(allele_frequencies) # plot the result plt.figure() for locus in xrange(5,pop.L,50): # plot a few neutral trajectories plt.plot(tp, allele_frequencies[:,locus], c=cm.cool(locus), lw=2)