forked from timbeissinger/Maize-Teo-Scripts
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrealisticBottleneck.py~
111 lines (75 loc) · 3.31 KB
/
realisticBottleneck.py~
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
"""
Realistic bottleneck model, no migration.
"""
"""
In words:
At time TF in the past, an equilibrium teosinte population will split into
a maize population and a teosinte population. The maize population will immediately go
through a bottleneck of size nuMB. Then, it will grow exponentially at
time TF. The Teosinte population will not go through a bottleneck. It will grow linearly to
end at size nuTeoF.
"""
import numpy
import dadi
def maizeBottleneck (params, ns, pts):
nuMB , nuMF , nuTeoF , TF , m12, m21 = params
# define the grid
xx = dadi.Numerics.default_grid(pts)
# phi for the equilibrium ancestral population
phi = dadi.PhiManip.phi_1D(xx)
# now do the population split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# Send maize population through a bottleneck
#phi = dadi.Integration.two_pops(phi, xx, T = TMB , nu1 = 1, nu2 = nuMB)
# define the maize exponential growth function
nu_funcM = lambda t: nuMB*(nuMF/nuMB) ** (t/TF)
# Option 1: define the teosinte exponential growth function
nu_funcTeo = lambda t: 1*(nuTeoF/1) ** (t/TF)
# Recover maize population exponentially, allow migration
phi = dadi.Integration.two_pops(phi, xx, T = TF, nu1 = nu_funcTeo, nu2 = nu_funcM, m12 = m12, m21 = m21)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs
def oneD_expansion ((nu , T ), ns, pts):
# define the grid
xx = dadi.Numerics.default_grid(pts)
# phi for the equilibrium ancestral population
phi = dadi.PhiManip.phi_1D(xx)
# define the exponential growth function
nu_func = lambda t: numpy.exp(numpy.log(nu) * t/T)
# Grow population
phi = dadi.Integration.one_pop(phi, xx, T, nu_func)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,))
return fs
def oneD_bottleneck ((nuB , nuF , TB, TF ), ns, pts):
# define the grid
xx = dadi.Numerics.default_grid(pts)
# phi for the equilibrium ancestral population
phi = dadi.PhiManip.phi_1D(xx)
# send maize through bottleneck
# phi = dadi.Integration.one_pop(phi, xx, TB, nuB)
# define the exponential growth function
nu_func = lambda t: nuB*(nuF/nuB) ** (t/TF)
# Grow population
phi = dadi.Integration.one_pop(phi, xx, TF, nu_func)
# phi = dadi.Integration.one_pop(phi, xx, TF, nuF)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,))
return fs
def maizeBottleneck_selection (params, ns, pts):
nuMB , nuMF , nuTeoF , TF , m12, m21, sM = params
# define the grid
xx = dadi.Numerics.default_grid(pts)
# phi for the equilibrium ancestral population
phi = dadi.PhiManip.phi_1D(xx)
# now do the population split
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# Send maize population through a bottleneck
#phi = dadi.Integration.two_pops(phi, xx, T = TMB , nu1 = 1, nu2 = nuMB)
# define the maize exponential growth function
nu_funcM = lambda t: nuMB*(nuMF/nuMB) ** (t/TF)
# define the teosinte exponential growth function
# nu_funcTeo = lambda t: numpy.exp(numpy.log(nuTeoF) * t/T)
nu_funcTeo = lambda t: 1*(nuTeoF/1) ** (t/TF)
# Recover maize population exponentially, allow migration and selection
phi = dadi.Integration.two_pops(phi, xx, T = TF, nu1 = nu_funcTeo, nu2 = nu_funcM, m12 = m12, m21 = m21, gamma1 = 0, gamma2 =sM)
fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
return fs