from __future__ import division
from collections import defaultdict, OrderedDict
from copy import deepcopy
import seaborn as sns
sns.set_style('white')
import matplotlib.pyplot as plt
import simuPOP as sp
from simuPOP import demography
%matplotlib inline
simuPOP Version 1.1.4 : Copyright (c) 2004-2011 Bo Peng Revision 4950 (Nov 11 2014) for Python 2.7.8 (64bit, 1thread) Random Number Generator is set to mt19937 with random seed 0x7e4566e8103c3a64. This is the standard short allele version with 256 maximum allelic states. For more information, please visit http://simupop.sourceforge.net, or email simupop-list@lists.sourceforge.net (subscription required).
num_loci = 10
pop_size = 50
num_gens = 101
num_pops = 10
migs = [0, 0.005, 0.01, 0.02, 0.05, 0.1]
init_ops = OrderedDict()
pre_ops = OrderedDict()
post_ops = OrderedDict()
def init_accumulators(pop, param):
accumulators = param
for accumulator in accumulators:
if accumulator.endswith('_sp'):
pop.vars()[accumulator] = defaultdict(list)
else:
pop.vars()[accumulator] = []
return True
def update_accumulator(pop, param):
accumulator, var = param
if var.endswith('_sp'):
for sp in range(pop.numSubPop()):
pop.vars()[accumulator][sp].append(deepcopy(pop.vars(sp)[var[:-3]]))
else:
pop.vars()[accumulator].append(deepcopy(pop.vars()[var]))
return True
pops = sp.Population([pop_size] * num_pops, loci=[1] * num_loci, infoFields=['migrate_to'])
init_ops['accumulators'] = sp.PyOperator(init_accumulators, param=['fst'])
init_ops['Sex'] = sp.InitSex()
init_ops['Freq'] = sp.InitGenotype(freq=[0.5, 0.5])
for i, mig in enumerate(migs):
post_ops['mig-%d' % i] = sp.Migrator(demography.migrIslandRates(mig, num_pops), reps=[i])
post_ops['Stat-fst'] = sp.Stat(structure=sp.ALL_AVAIL)
post_ops['fst_accumulation'] = sp.PyOperator(update_accumulator, param=('fst', 'F_st'))
mating_scheme = sp.RandomMating()
sim = sp.Simulator(pops, rep=len(migs))
sim.evolve(initOps=init_ops.values(), preOps=pre_ops.values(), postOps=post_ops.values(),
matingScheme=mating_scheme, gen=num_gens)
(101L, 101L, 101L, 101L, 101L, 101L)
fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(111)
for i, pop in enumerate(sim.populations()):
ax.plot(pop.dvars().fst, label='mig rate %.4f' % migs[i])
ax.legend(loc=2)
ax.set_ylabel('FST')
ax.set_xlabel('Generation')
<matplotlib.text.Text at 0x7f46ddc14350>
num_gens = 400
num_loci = 5
init_ops = OrderedDict()
pre_ops = OrderedDict()
post_ops = OrderedDict()
init_ops['Sex'] = sp.InitSex()
init_ops['Freq'] = sp.InitGenotype(freq=[0.5, 0.5])
post_ops['Stat-freq'] = sp.Stat(alleleFreq=sp.ALL_AVAIL, vars=['alleleFreq', 'alleleFreq_sp'])
init_ops['accumulators'] = sp.PyOperator(init_accumulators, param=['allele_freq', 'allele_freq_sp'])
post_ops['freq_accumulation'] = sp.PyOperator(update_accumulator, param=('allele_freq', 'alleleFreq'))
post_ops['freq_sp_accumulation'] = sp.PyOperator(update_accumulator,
param=('allele_freq_sp', 'alleleFreq_sp'))
for i, mig in enumerate(migs):
post_ops['mig-%d' % i] = sp.Migrator(demography.migrSteppingStoneRates(mig, num_pops), reps=[i])
#two 2 - talk
pops = sp.Population([pop_size] * num_pops, loci=[1] * num_loci, infoFields=['migrate_to'])
sim = sp.Simulator(pops, rep=len(migs))
sim.evolve(initOps=init_ops.values(), preOps=pre_ops.values(), postOps=post_ops.values(),
matingScheme=mating_scheme, gen=num_gens)
(400L, 400L, 400L, 400L, 400L, 400L)
def get_maf(var):
locus_data = [gen[locus] for gen in var]
maf = [min(freq.values()) for freq in locus_data]
maf = [v if v != 1 else 0 for v in maf]
return maf
fig, axs = plt.subplots(3, num_pops // 2 + 1, figsize=(16, 9), sharex=True, sharey=True, squeeze=False)
fig.suptitle('Minimum allele frequency at the meta-population and 5 demes', fontsize='xx-large')
for line, pop in enumerate([sim.population(0), sim.population(1), sim.population(len(migs) - 1)]):
for locus in range(num_loci):
maf = get_maf(pop.dvars().allele_freq)
axs[line, 0].plot(maf)
axs[line, 0].set_axis_bgcolor('black')
for nsp in range(num_pops // 2):
for locus in range(num_loci):
maf = get_maf(pop.dvars().allele_freq_sp[nsp * 2])
axs[line, nsp + 1].plot(maf)
fig.subplots_adjust(hspace=0, wspace=0)