In [48]:
import numpy as np
from cycler import cycler
import matplotlib.pyplot as plt
from pylab import rcParams
numIndividuals = 100
def generation (genolist):
newGenos = list()
length = len(genolist)
selection = np.random.choice(length, length)
for i in range(length):
newGenos.append(genolist[selection[i]])
return(newGenos)
#############################
def simulate (numIndividuals, numGenerations, numPopulations):
# Two alleles: true and false
init = [True]*int(numIndividuals/2) + [False] * int(numIndividuals/2)
allFreqs = list()
for pop in range(numPopulations):
freqs = list()
freqs.append(sum(init)/ numIndividuals)
current = init
for gen in range(numGenerations):
new = generation(current)
freqs.append(sum(new) / numIndividuals)
current = new
allFreqs.append(freqs)
return allFreqs
def plotFreqs (allFreqs):
%matplotlib inline
rcParams['figure.figsize'] = 18, 8
rcParams['xtick.labelsize'] = 20
rcParams['ytick.labelsize'] = 20
plt.ylim([0,1.0])
plt.suptitle('Allele frequency over generations', fontsize=20, fontweight='bold')
plt.rc('lines', linewidth=2)
plt.rc('axes', prop_cycle=(cycler('color', ['r', 'g', 'b', 'y'])))
# + cycler('linestyle', ['-', '--', ':', '-.'])))
for i in range(len(allFreqs)):
plt.plot (allFreqs[i]) #(position, coverage, 'r-')
plt.xlabel ('Generation', fontsize=18)
plt.ylabel ('Allele frequency', fontsize=18)
#plt.tight_layout()
In [49]:
plotFreqs(simulate (20, 50, 10))
In [51]:
plotFreqs(simulate (200, 50, 10))
In [53]:
plotFreqs(simulate (2000, 50, 10))
In [63]:
def generationSelection (genolist):
newGenos = list()
length = len(genolist)
weights = list()
for i in range(length):
if genolist[i]:
weights.append(1.05)
else:
weights.append(1.0)
probas = [weights[i]/sum(weights) for i in range(len(weights))]
newGenos = np.random.choice(genolist, length, p=probas)
return(newGenos)
#############################
def simulateSelection (numIndividuals, numGenerations, numPopulations):
# Two alleles: true and false
init = [True]*int(numIndividuals/2) + [False] * int(numIndividuals/2)
allFreqs = list()
for pop in range(numPopulations):
freqs = list()
freqs.append(sum(init)/ numIndividuals)
current = init
for gen in range(numGenerations):
new = generationSelection(current)
freqs.append(sum(new) / numIndividuals)
current = new
allFreqs.append(freqs)
return allFreqs
In [64]:
plotFreqs(simulateSelection (20, 50, 10))
In [65]:
plotFreqs(simulateSelection (200, 50, 10))
In [68]:
plotFreqs(simulateSelection (2000, 50, 10))
In [ ]: