Inference for one-step beneficial mutations using next generation sequencing

Appendix

A Simulation of population dynamics for a single first-step beneficial mutation

Below we present a procedure of simulation of a population dynamics in which one beneficial mutation competes against the wild type. The relative frequency of a mutation is recorded in consecutive generations.

1. Let N_wt(t) denote the count of the wild type and let N_mut(t) denote the count of the mutant in generation t. Let N_wt=N and N_mut=0 at t=0.

2. The wild type individuals produce Gwt(t)=Nwt(t)2wwt progeny and mutants produce Gmut(t)=Nmut(t)2wwt(1+si) progeny.

3. The probability that an offspring of the wild type is a mutant equals μ. Let M(t) denote the number of mutants among the progeny of the wild type, then M(t)∼Bin(G_wt(t), μ).

4. The population size is constant and equal to N. Therefore, at the end of each generation the number of progeny needs to be reduced to N by binomial sampling. Let

   Nmut(t+1)~ Bin (N,Gmut(t)+M(t)Gwt(t)+Gmut(t))

   and

   Nwt(t+1)=N−Nmut(t+1).

5. Record current relative frequency of a mutant. Go back to step 2.

B Solution to the difference equation for the mean proportion of a mutation

Here we provide a solution to the following difference equation

π(t+1)=π(t)(qs−μ)+μπ(t)(qs−1)+1.

Let a=q^s–μ, b=μ, c=q^s–1, and d=1. Since

(d−a)2+4bc=((1−μ)−qs)2

is always positive, we can apply the following method. Let

ρ=(d−a)2+4bc=(1−qs+μ)2+4μ(qs−1)=((1−μ)−qs)2=qs+μ−1.

Note that in our case ρ is always positive since q^s>1–μ, because s≥0, w_wt≥0, and 0≤μ≤1. Let

η=d−a+ρ2c=1−qs+μ+qs−1+μ2(qs−1)=μqs−1.

Then

(4)π(t)=1x(t)−η. (4)

To calculate x(t) the following first-order difference equation has to be solved

x(t+1)=d−ηcηc+ax(t)+cηc+a.

Let

α=d−ηcηc+aβ=cηc+a

then

x(t)=αtx(0)+β(1−αt)1−α,

where

x(0)=1η+π(0).

Therefore,

x(t)=αtη+β(1−αt)1−α=(ηc+a−d−cπ(0))(d−ηc)t+c(ηc+a)t+c(ηc+a)t(η+π(0))(ηc+a)t(η+π(0))(2ηc+a−d).

Next, using equation (4) we obtained the final formula for π(t)

π(t)=(ηc+a)t(η+π(0))(2ηc+a−d)(ηc+a−d−cπ(0))(d−ηc)t+c(ηc+a)t+c(ηc+a)t(η+π(0))−η=1qs−1(qst(qs+μ−1)(μ+π(0)(qs−1))(qs−1)(1−π(0))(1−μ)t+qst(μ+π(0)(qs−1))−μ).

C Likelihood function ofrandδ

Result: The log likelihood function of r and δ is given by:

lnL(r,δ)=ln(rY)+∑a=1Y∑b=1nln(kaxa,b)+(r−Y)ln∫0∞(1−g(s))nkf(s)ds+∑i=1Yln∫0∞g(si)∑j=1nxi,j(1−g(si))nki−∑j=1nxi,jf(si)dsi.

Proof: Since all entries in the matrix X are assumed to be independent random variables distributed according to a binomial distribution, the probability of observing a particular matrix X is equal to the product of binomial probabilities of cells in the matrix multiplied by the combinatorial term that accounts for different possible order of observations. This probability is given by the following formula

P(X1,1=x1,1, X1,2=x1,2,...,XY,n=xY,n, XY+1,1=0,...,Xr,n=0|s1,...,sr)=(rY)P(X1,1=x1,1)P(X1,2=x1,2)...P(XY,n=xY,n)P(XY+1,1=0)...P(Xr,n=0)=(rY)∏a=1Y∏b=1n(kaxa,b)∏i=1Yg(si)∑j=1nxi,j(1−g(si))nki−∑j=1nxi,j∏l=Y+1r(1−g(sl))nk,

where g(s_i) is given by equation (1). Here we assume that selection coefficients s₁, s₂,…, s_r are unobserved random variables with probability density function f(s_i). Therefore, we define the joint probability of a matrix X and a vector of selection coefficients as

P(X1,1=x1,1, X1,2=x1,2,...,XY,n=xY,n, XY+1,1=0,...,Xr,n=0,s1,...,sr)=(rY)∏a=1Y∏b=1n(kaxa,b)∏i=1Yg(si)∑j=1nxi,j(1−g(si))nki−∑j=1nxi,jf(si)∏l=Y+1r(1−g(sl))nkf(sl).

To obtain the likelihood function of r and δ we integrated out selection coefficients. Here we present the log likelihood function

lnL(r,δ)=ln(rY)+∑a=1Y∑b=1nln(kaxa,b)+(r−Y)ln∫0∞(1−g(s))nkf(s)ds+∑i=1Yln∫0∞g(si)∑j=1nxi,j(1−g(si))nki−∑j=1nxi,jf(si)dsi.

D Expected number of singletons

Result: The expected number of singletons is given by the following formula

E(Z)=rnk∫0∞g(s)(1−g(s))nk−1f(s)ds.

Proof: Denote the number singletons by Z, then

Z=∑i=1rIi,

where I_i is an indicator function such that

Ii=(1if the i−th mutant was observed only once in one replicate0otherwise.

Then

E(Ii|si)=P(Ii=1)=P(∑j=1nXi,j=1).

Since

∑j=1nXi,j~ Bin (nk, g(si)),

where g(s_i) is given by formula (1),

E(Ii|si)=nkg(si)(1−g(si))nk−1.

To obtain the expected value of I_i, the selection coefficient s_i has to be integrated out from the above formula. Therefore,

E(Ii)=nk∫0∞g(si)(1−g(si))nk−1f(si)dsi,

where f(s_i) is the probability density function of s_i. Based on this result we obtained the expected value of the number of singletons

E(Z)=∑i=1rE(Ii)=rnk∫0∞g(s)(1−g(s))nk−1f(s)ds.

E Expected number of neutral mutations observed more than once

Result: The expected number of neutral mutations observed more than once is given by the following formula

E(V)=rneu(1−(1−tμ−π(0))nk−nk(tμ+π(0))(1−tμ−π(0))nk−1).

Proof: We assume that a neutral mutation accumulates in a population at the rate μ. Therefore, under the deterministic model, after t generations proportion of that mutation in a population is equal to tμ+π(0), where π(0) is the proportion of that mutation in the initial inoculum. Since a site where a particular neutral mutation can occur is sampled nk times and the probability of sampling this mutation in one trial equals tμ+π(0), the number of times it is observed in n replicates is a binomially distributed random variable, which we denote by W

W~ Bin (nk, tμ+π(0)).

Then, the probability of sampling that mutation more than once (denoted by γ) is given by

γ=P(W>1)=1−[P(W=0)+P(W=1)]=1−(1−tμ−π(0))nk−nk(tμ+π(0))(1−tμ−π(0))nk−1.

Assuming that there are r_new available neutral mutations, the number of neutral mutations observed more than once (denoted by V) has a binomial distribution

V~ Bin (rneu, γ).

Therefore, the expected number of neutral mutations observed more than once is given by

E(V)=rneuγ=rneu(1−(1−tμ−π(0))nk−nk(tμ+π(0))(1−tμ−π(0))nk−1).