Skip to content

integration

Integration of phi for triallelic diffusion These methods include ones for the injection of density for new triallelic sites and integration forward in time

advance(phi, x, T, y1, y2, nu=1.0, sig1=0.0, sig2=0.0, theta1=1.0, theta2=1.0, dt=0.001)

Integrate phi, y1, and y2 forward in time

lam - proportion of mutations that occur from simulateous mutation model (Hodgkinson/Eyre-Walker 2010) Args: phi (array-like): density function for triallelic sites x (array-like): grid T (float): amount of time to integrate, scaled by 2N generations y1 (array-like): density of biallelic background sites, integrated forward alongside phi y2 (array-like): density of biallelic background sites, integrated forward alongside phi nu (float): relative size of population to ancestral size sig1 (float): population scaled selection coefficient for the first derived allele. sig2 (float): population scaled selection coefficients for the second derived allele. theta1 (float): population scaled mutation rates for the first derived allele. theta2 (float): population scaled mutation rates for the second derived allele. dt (float): time step for integration

Source code in dadi/Triallele/integration.py 
 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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
def advance(phi, x, T, y1, y2, nu=1., sig1=0., sig2=0., theta1=1., theta2=1., dt=0.001):
    """
    Integrate phi, y1, and y2 forward in time

    lam - proportion of mutations that occur from simulateous mutation model (Hodgkinson/Eyre-Walker 2010)
    Args:
        phi (array-like): density function for triallelic sites
        x (array-like): grid
        T (float): amount of time to integrate, scaled by 2N generations
        y1 (array-like): density of biallelic background sites, integrated forward alongside phi
        y2 (array-like): density of biallelic background sites, integrated forward alongside phi
        nu (float): relative size of population to ancestral size
        sig1 (float): population scaled selection coefficient for the first derived allele.
        sig2 (float): population scaled selection coefficients for the second derived allele.
        theta1 (float): population scaled mutation rates for the first derived allele.
        theta2 (float): population scaled mutation rates for the second derived allele.
        dt (float): time step for integration
    """
    dx = numerics.grid_dx(x)
    U01 = numerics.domain(x)

    C_base = numerics.transition12(x,dx,U01)
    V1_base,M1_base = numerics.transition1(x,dx,U01,sig1,sig2)
    V2_base,M2_base = numerics.transition2(x,dx,U01,sig1,sig2)
    V1D1_base,M1D1_base = numerics.transition1D(x,dx,sig1)
    V1D2_base,M1D2_base = numerics.transition1D(x,dx,sig2)
    sig_line = sig1-sig2
    Vline_base,Mline_base = numerics.transition1D(x,dx,sig_line)

    if np.isscalar(nu):
        C = identity(len(x)**2) + dt/nu*C_base
        P1 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*(V1_base/nu+M1_base)
        P2 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*(V2_base/nu+M2_base)
        P1D1 = np.eye(len(x)) + dt*(V1D1_base/nu+M1D1_base)
        P1D2 = np.eye(len(x)) + dt*(V1D2_base/nu+M1D2_base)
        Pline = np.eye(len(x)) + dt*(Vline_base/nu+Mline_base)
        P = numerics.remove_diag_density_weights_nonneutral(x,dt,nu,sig1,sig2)

        for ii in range(int(T/dt)):
            y1[1] += dt/dx[1]/x[1]/2 * theta1
            y1 = numerics.advance1D(y1,P1D1)
            y2[1] += dt/dx[1]/x[1]/2 * theta2
            y2 = numerics.advance1D(y2,P1D2)
            phi = inject_mutations_1(phi, dt, x, dx, y2, theta1)
            phi = inject_mutations_2(phi, dt, x, dx, y1, theta2)
            phi = numerics.advance_adi(phi,U01,P1,P2,x,ii)
            phi = numerics.advance_cov(phi,C,x,dx)
            #phi *= 1-P
            # move density to diagonal boundary and integrate it
            phi = numerics.move_density_to_bdry(x,phi,P)
            phi = numerics.advance_line(x,phi,Pline)

        T_elapsed = int(T/dt)*dt
        if T - T_elapsed > 1e-8:
            # adjust dt and integrate last time step
            dt = T-T_elapsed
            C = identity(len(x)**2) + dt/nu*C_base
            P1 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*(V1_base/nu+M1_base)
            P2 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*(V2_base/nu+M2_base)
            P1D1 = np.eye(len(x)) + dt*(V1D1_base/nu+M1D1_base)
            P1D2 = np.eye(len(x)) + dt*(V1D2_base/nu+M1D2_base)
            Pline = np.eye(len(x)) + dt*(Vline_base/nu+Mline_base)
            P = numerics.remove_diag_density_weights_nonneutral(x,dt,nu,sig1,sig2)

            y1[1] += dt/dx[1]/x[1]/2 * theta1
            y1 = numerics.advance1D(y1,P1D1)
            y2[1] += dt/dx[1]/x[1]/2 * theta2
            y2 = numerics.advance1D(y2,P1D2)
            phi = inject_mutations_1(phi, dt, x, dx, y2, theta1)
            phi = inject_mutations_2(phi, dt, x, dx, y1, theta2)
            phi = numerics.advance_adi(phi,U01,P1,P2,x,0)
            phi = numerics.advance_cov(phi,C,x,dx)
            #phi *= 1-P
            # move density to diagonal boundary and integrate it
            phi = numerics.move_density_to_bdry(x,phi,P)
            phi = numerics.advance_line(x,phi,Pline)
    else:
        Ts = np.concatenate(( np.linspace(0,np.floor(T/dt)*dt,np.floor(T/dt)+1), np.array([T]) ))

        for ii in range(int(T/dt)):
            nu_current = nu(Ts[ii])
            C = identity(len(x)**2) + dt/nu_current*C_base
            P1 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*(V1_base/nu_current+M1_base)
            P2 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*(V2_base/nu_current+M2_base)
            P1D1 = np.eye(len(x)) + dt*(V1D1_base/nu_current+M1D1_base)
            P1D2 = np.eye(len(x)) + dt*(V1D2_base/nu_current+M1D2_base)
            Pline = np.eye(len(x)) + dt*(Vline_base/nu_current+Mline_base)
            P = numerics.remove_diag_density_weights_nonneutral(x,dt,nu_current,sig1,sig2)                

            y1[1] += dt/dx[1]/x[1]/2 * theta1
            y1 = numerics.advance1D(y1,P1D1)
            y2[1] += dt/dx[1]/x[1]/2 * theta2
            y2 = numerics.advance1D(y2,P1D2)
            phi = inject_mutations_1(phi, dt, x, dx, y2, theta1)
            phi = inject_mutations_2(phi, dt, x, dx, y1, theta2)
            phi = numerics.advance_adi(phi,U01,P1,P2,x,ii)
            phi = numerics.advance_cov(phi,C,x,dx)
            #phi *= 1-P
            # move density to diagonal boundary and integrate it
            phi = numerics.move_density_to_bdry(x,phi,P)
            phi = numerics.advance_line(x,phi,Pline)

        T_elapsed = int(T/dt)*dt
        if T - T_elapsed > 1e-8:
            # adjust dt and integrate last time step
            dt = T-T_elapsed
            nu_current = nu(Ts[-1])
            C = identity(len(x)**2) + dt/nu_current*C_base
            P1 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*(V1_base/nu_current+M1_base)
            P2 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*(V2_base/nu_current+M2_base)
            P1D1 = np.eye(len(x)) + dt*(V1D1_base/nu_current+M1D1_base)
            P1D2 = np.eye(len(x)) + dt*(V1D2_base/nu_current+M1D2_base)
            Pline = np.eye(len(x)) + dt*(Vline_base/nu_current+Mline_base)
            P = numerics.remove_diag_density_weights_nonneutral(x,dt,nu_current,sig1,sig2)                

            y1[1] += dt/dx[1]/x[1]/2 * theta1
            y1 = numerics.advance1D(y1,P1D1)
            y2[1] += dt/dx[1]/x[1]/2 * theta2
            y2 = numerics.advance1D(y2,P1D2)
            phi = inject_mutations_1(phi, dt, x, dx, y2, theta1)
            phi = inject_mutations_2(phi, dt, x, dx, y1, theta2)
            phi = numerics.advance_adi(phi,U01,P1,P2,x,0)
            phi = numerics.advance_cov(phi,C,x,dx)
            #phi *= 1-P
            # move density to diagonal boundary and integrate it
            phi = numerics.move_density_to_bdry(x,phi,P)
            phi = numerics.advance_line(x,phi,Pline)

    return phi,y1,y2

advance_old(phi, x, T, y1, y2, nu=1.0, sig1=0.0, sig2=0.0, theta1=1.0, theta2=1.0, dt=0.001)

Integrate phi, y1, and y2 forward in time

lam - proportion of mutations that occur from simulateous mutation model (Hodgkinson/Eyre-Walker 2010) Args: phi (array-like): density function for triallelic sites x (array-like): grid T (float): amount of time to integrate, scaled by 2N generations y1 (array-like): density of biallelic background sites, integrated forward alongside phi y2 (array-like): density of biallelic background sites, integrated forward alongside phi nu (float): relative size of population to ancestral size sig1 (float): population scaled selection coefficient for the first derived allele. sig2 (float): population scaled selection coefficients for the second derived allele. theta1 (float): population scaled mutation rates for the first derived allele. theta2 (float): population scaled mutation rates for the second derived allele. dt (float): time step for integration

Source code in dadi/Triallele/integration.py 
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
def advance_old(phi, x, T, y1, y2, nu=1., sig1=0., sig2=0., theta1=1., theta2=1., dt=0.001):
    """
    Integrate phi, y1, and y2 forward in time

    lam - proportion of mutations that occur from simulateous mutation model (Hodgkinson/Eyre-Walker 2010)
    Args:
        phi (array-like): density function for triallelic sites
        x (array-like): grid
        T (float): amount of time to integrate, scaled by 2N generations
        y1 (array-like): density of biallelic background sites, integrated forward alongside phi
        y2 (array-like): density of biallelic background sites, integrated forward alongside phi
        nu (float): relative size of population to ancestral size
        sig1 (float): population scaled selection coefficient for the first derived allele.
        sig2 (float): population scaled selection coefficients for the second derived allele.
        theta1 (float): population scaled mutation rates for the first derived allele.
        theta2 (float): population scaled mutation rates for the second derived allele.
        dt (float): time step for integration
    """
    dx = numerics.grid_dx(x)
    U01 = numerics.domain(x)
    C = identity(len(x)**2) + dt/nu*numerics.transition12(x,dx,U01)
    P1 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*numerics.transition1(x,dx,U01,sig1,sig2,nu) 
    P2 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*numerics.transition2(x,dx,U01,sig1,sig2,nu)
    P1D1 = numerics.transition1D(x,dx,dt,sig1,nu)
    P1D2 = numerics.transition1D(x,dx,dt,sig2,nu)

    sig_line = sig1-sig2
    Pline = numerics.transition1D(x,dx,dt,sig_line,nu)
    P = numerics.remove_diag_density_weights_nonneutral(x,dt,nu,sig1,sig2)

    for ii in range(int(T/dt)):
        y1[1] += dt/dx[1]/x[1]/2 * theta1
        y1 = numerics.advance1D(y1,P1D1)
        y2[1] += dt/dx[1]/x[1]/2 * theta2
        y2 = numerics.advance1D(y2,P1D2)
        phi = inject_mutations_1(phi, dt, x, dx, y2, theta1)
        phi = inject_mutations_2(phi, dt, x, dx, y1, theta2)
        phi = numerics.advance_adi(phi,U01,P1,P2,x,ii)
        phi = numerics.advance_cov(phi,C,x,dx)
        #phi *= 1-P
        # move density to diagonal boundary and integrate it
        phi = numerics.move_density_to_bdry(x,phi,P)
        phi = numerics.advance_line(x,phi,Pline)

    T_elapsed = int(T/dt)*dt
    if T - T_elapsed > 1e-8:
        # adjust dt and integrate last time step
        dt = T-T_elapsed
        C = identity(len(x)**2) + dt/nu*numerics.transition12(x,dx,U01) # covariance term
        P1 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*numerics.transition1(x,dx,U01,sig1,sig2,nu) 
        P2 = np.outer(np.array([0,1,0]),np.ones(len(x))) + dt*numerics.transition2(x,dx,U01,sig1,sig2,nu)
        P1D1 = numerics.transition1D(x,dx,dt,sig1,nu)
        P1D2 = numerics.transition1D(x,dx,dt,sig2,nu)

        sig_line = sig1-sig2
        Pline = numerics.transition1D(x,dx,dt,sig_line,nu)
        P = numerics.remove_diag_density_weights_nonneutral(x,dt,nu,sig1,sig2)

        y1[1] += dt/dx[1]/x[1]/2 * theta1
        y1 = numerics.advance1D(y1,P1D1)
        y2[1] += dt/dx[1]/x[1]/2 * theta2
        y2 = numerics.advance1D(y2,P1D2)
        phi = inject_mutations_1(phi, dt, x, dx, y2, theta1)
        phi = inject_mutations_2(phi, dt, x, dx, y1, theta2)
        phi = numerics.advance_adi(phi,U01,P1,P2,x,0)
        phi = numerics.advance_cov(phi,C,x,dx)
        #phi *= 1-P
        # move density to diagonal boundary and integrate it
        phi = numerics.move_density_to_bdry(x,phi,P)
        phi = numerics.advance_line(x,phi,Pline)

    return phi,y1,y2

alt_mut_mech_sample_spectrum(ns)

alternate mutation mechanism, mutations inserted at [1,1] turns out that changing population size does not effect the distribution of mutations entering the population this way we implement Jenkins et al (2014) exact solution this is for neutral spectrum only, for selected spectrum, integrate as above with lam = 1 ns - number of sampled individuals from the population

Source code in dadi/Triallele/integration.py 
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
def alt_mut_mech_sample_spectrum(ns):
    """
    alternate mutation mechanism, mutations inserted at [1,1]
    turns out that changing population size does not effect the distribution of mutations entering the population this way
    we implement Jenkins et al (2014) exact solution
    this is for neutral spectrum only, for selected spectrum, integrate as above with lam = 1
    ns - number of sampled individuals from the population
    """
    fs = np.zeros((ns+1,ns+1))
    for ii in range(ns)[1:]:
        for jj in range(ns)[1:]:
            if ii + jj < ns:
                na = ns - ii - jj
                fs[ii,jj] = 2*ns/(ns-2) * 1./((ns-na-1)*(ns-na)*(ns-na+1))
    fs = dadi.Spectrum(fs)
    fs[:,0].mask = True
    fs[0,:].mask = True
    for ii in range(len(fs)):
        fs.mask[ii,ns-ii:] = True
    return fs

equilibrium_neutral_exact(x)

With thetas = 1 nu = 1 sig1 = sig2 = 0

Source code in dadi/Triallele/integration.py 
49
50
51
52
53
54
55
56
57
58
def equilibrium_neutral_exact(x):
    """
    With thetas = 1
    nu = 1
    sig1 = sig2 = 0
    """
    phi = np.zeros((len(x),len(x)))
    for ii in range(len(phi))[1:]:
        phi[ii,1:-ii-1] = 1./x[ii]/x[1:-ii-1]
    return phi

inject_mutations_1(phi, dt, x, dx, y2, theta1)

new mutations injected along phi[1,:] against a background given by y2 Args: phi (array-like): numerical density function dt (float): given time step x (array-like): one dimensional grid dx (array-like): grid spacing y2 (array-like): the biallelic density function theta1 (float): population scaled mutation rate for mutation 1

Source code in dadi/Triallele/integration.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
def inject_mutations_1(phi, dt, x, dx, y2, theta1):
    """
    new mutations injected along phi[1,:] against a background given by y2
    Args:
        phi (array-like): numerical density function
        dt (float): given time step 
        x (array-like): one dimensional grid
        dx (array-like): grid spacing
        y2 (array-like): the biallelic density function
        theta1 (float): population scaled mutation rate for mutation 1
    """
    phi[1,1:-1] += y2[1:-1] / dx[1] * 1./x[1] * dt * theta1/2
    return phi

inject_mutations_2(phi, dt, x, dx, y1, theta2)

new mutations injected along phi[:,1] against a background given by y1 Args: phi (array-like): numerical density function dt (float): given time step x (array-like): one dimensional grid dx (array-like): grid spacing y1 (array-like): the biallelic density function theta2 (float): population scaled mutation rate for mutation 2

Source code in dadi/Triallele/integration.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
def inject_mutations_2(phi, dt, x, dx, y1, theta2):
    """
    new mutations injected along phi[:,1] against a background given by y1
    Args:
        phi (array-like): numerical density function
        dt (float): given time step
        x (array-like): one dimensional grid
        dx (array-like): grid spacing
        y1 (array-like): the biallelic density function
        theta2 (float): population scaled mutation rate for mutation 2
    """
    phi[1:-1,1] += y1[1:-1] / dx[1] * 1./x[1] * dt * theta2/2
    return phi

inject_simultaneous_muts(phi, dt, x, dx, theta)

simultaneous mutation model - see Hodgkinson and Eyre-Walker 2010, injected at (Delta,Delta)

Source code in dadi/Triallele/integration.py 
42
43
44
45
46
47
def inject_simultaneous_muts(phi, dt, x, dx, theta):
    """
    simultaneous mutation model - see Hodgkinson and Eyre-Walker 2010, injected at (Delta,Delta)
    """
    phi[1,1] += 1. / x[1] / x[1] / dx[1] / dx[1] * dt * theta
    return phi