Maximum likelihood (ML) estimators for scaled mutation parameters with a   strand symmetric mutation model in equilibrium

Claus Vogl; Lynette Caitlin Mikula

arXiv:1905.02015·q-bio.PE·May 7, 2019

Maximum likelihood (ML) estimators for scaled mutation parameters with a strand symmetric mutation model in equilibrium

Claus Vogl, Lynette Caitlin Mikula

PDF

Open Access

TL;DR

This paper derives maximum likelihood estimators for six mutation parameters in a strand-symmetric mutation model, extending beyond simpler parent-independent models, under equilibrium and small mutation rate assumptions.

Contribution

It introduces the first ML estimators for a complex mutation model that deviates from parent-independent mutation, broadening the scope of mutation parameter inference.

Findings

01

ML estimators derived for six mutation parameters

02

Applicable under mutation-drift equilibrium and small mutation rates

03

Extends mutation modeling beyond parent-independent assumptions

Abstract

With the multiallelic parent-independent mutation-drift model, the equilibrium proportions of alleles are known to be Dirichlet distributed. A special case is the biallelic model, in which the proportions are beta distributed. A sample taken from these models is then Dirichlet-multinomially or beta-binomially distributed, respectively. Maximum likelihood (ML) estimators for the mutation parameters of the biallelic parent-independent mutation model are available via an expectation maximization algorithm. Assuming small scaled mutation rates, the distribution of a sample of size $M$ can be expanded in a Taylor series of first order. Then the ML estimators for the two parameters in the biallelic model can be expressed using the site frequency spectrum. In this article, we go beyond parent-independent mutation and analyse a strand-symmetric mutation model with six scaled mutation parameters…

Equations130

{\Pr}(y{\,|\,}M,\alpha,\theta\to 0)=\begin{cases}(1-\alpha)\big{(}1-\alpha\theta H_{M-1}\big{)}+O(\theta^{2})&\text{for $y=0$,}\\ \alpha(1-\alpha)\theta\,\frac{M}{y(M-y)}+O(\theta^{2})&\text{for $1\leq y\leq(M-1)$,}\\ \alpha\big{(}1-(1-\alpha)\theta H_{M-1}\big{)}+O(\theta^{2})&\text{for $y=M$}\,.\end{cases}

{\Pr}(y{\,|\,}M,\alpha,\theta\to 0)=\begin{cases}(1-\alpha)\big{(}1-\alpha\theta H_{M-1}\big{)}+O(\theta^{2})&\text{for $y=0$,}\\ \alpha(1-\alpha)\theta\,\frac{M}{y(M-y)}+O(\theta^{2})&\text{for $1\leq y\leq(M-1)$,}\\ \alpha\big{(}1-(1-\alpha)\theta H_{M-1}\big{)}+O(\theta^{2})&\text{for $y=M$}\,.\end{cases}

\overset{α}{^} = \frac{L _{M} + L _{p}}{L}

\overset{α}{^} = \frac{L _{M} + L _{p}}{L}

\hat{θ} = \frac{1}{2 α ( 1 - α )} \frac{L _{p}}{H _{M - 1} L},

\hat{θ} = \frac{1}{2 α ( 1 - α )} \frac{L _{p}}{H _{M - 1} L},

Q = (- a_{12} - a_{13} - a_{14} a_{21} a_{31} a_{41} a_{12} - a_{21} - a_{23} - a_{24} a_{32} a_{42} a_{13} a_{23} - a_{31} - a_{32} - a_{34} a_{43} a_{14} a_{24} a_{34} - a_{41} - a_{42} - a_{43}) .

Q = (- a_{12} - a_{13} - a_{14} a_{21} a_{31} a_{41} a_{12} - a_{21} - a_{23} - a_{24} a_{32} a_{42} a_{13} a_{23} - a_{31} - a_{32} - a_{34} a_{43} a_{14} a_{24} a_{34} - a_{41} - a_{42} - a_{43}) .

Q = - a - c - e a b d a - a - c - e d b c e - b - d - f f e c f - d - b - f

Q = - a - c - e a b d a - a - c - e d b c e - b - d - f f e c f - d - b - f

v_{1} = (\frac{d + b}{2 ( b + c + d + e )} \frac{d + b}{2 ( b + c + d + e )} \frac{c + e}{2 ( b + c + d + e )} \frac{c + e}{2 ( b + c + d + e )}) .

v_{1} = (\frac{d + b}{2 ( b + c + d + e )} \frac{d + b}{2 ( b + c + d + e )} \frac{c + e}{2 ( b + c + d + e )} \frac{c + e}{2 ( b + c + d + e )}) .

v_{2} = (11 - 1 - 1) .

v_{2} = (11 - 1 - 1) .

Q = θ \cdot - α_{1} α_{21} α_{31} α_{41} α_{12} - α_{2} α_{32} α_{42} α_{13} α_{23} - α_{3} α_{43} α_{14} α_{24} α_{34} - α_{4},

Q = θ \cdot - α_{1} α_{21} α_{31} α_{41} α_{12} - α_{2} α_{32} α_{42} α_{13} α_{23} - α_{3} α_{43} α_{14} α_{24} α_{34} - α_{4},

f_{0} (θ) = Pr (y_{k} = M ∣ α, M, θ, no mutation) = Pr (k ∣ α) \cdot \frac{1}{1 + α _{k} θ} \cdot \frac{2}{2 + α _{k} θ} \dots \frac{M - 1}{( M - 1 ) + α _{k} θ} .

f_{0} (θ) = Pr (y_{k} = M ∣ α, M, θ, no mutation) = Pr (k ∣ α) \cdot \frac{1}{1 + α _{k} θ} \cdot \frac{2}{2 + α _{k} θ} \dots \frac{M - 1}{( M - 1 ) + α _{k} θ} .

Pr (y_{k} = M ∣ α, M, no mutation) = Pr (k ∣ α) - Pr (k ∣ α) \cdot m = 2 \sum M \frac{α _{k} θ}{m - 1} + O (θ^{2}) .

Pr (y_{k} = M ∣ α, M, no mutation) = Pr (k ∣ α) - Pr (k ∣ α) \cdot m = 2 \sum M \frac{α _{k} θ}{m - 1} + O (θ^{2}) .

f_{0} (θ = 0) = Pr (k ∣ α)

f_{0} (θ = 0) = Pr (k ∣ α)

\begin{split}f_{0}^{\prime}(\theta)&=\Pr(k\mid\vec{\alpha})\cdot\bigg{(}-\frac{\alpha_{k}}{(1+\alpha_{k}\theta)^{2}}\bigg{[}\frac{2}{2+\alpha_{k}\theta}\cdot\frac{3}{3+\alpha_{k}\theta}\cdots\frac{(M-1)}{(M-1)+\alpha_{k}\theta}\bigg{]}-\\ &\qquad-\frac{2\alpha_{k}}{(2+\alpha_{k}\theta)^{2}}\left[\frac{1}{1+\alpha_{k}\theta}\cdot\frac{3}{3+\alpha_{k}\theta}\cdots\frac{(M-1)}{(M-1)+\alpha_{k}\theta}\right]-\cdots\\ &\qquad\cdots-\frac{(M-1)\alpha_{k}}{((M-1)+\alpha_{k}\theta)^{2}}\bigg{[}\frac{1}{1+\alpha_{k}\theta}\cdot\frac{2}{2+\alpha_{k}\theta}\cdots\frac{(M-2)}{(M-2)+\alpha_{k}\theta}\bigg{]}\bigg{)}\,,\end{split}

\begin{split}f_{0}^{\prime}(\theta)&=\Pr(k\mid\vec{\alpha})\cdot\bigg{(}-\frac{\alpha_{k}}{(1+\alpha_{k}\theta)^{2}}\bigg{[}\frac{2}{2+\alpha_{k}\theta}\cdot\frac{3}{3+\alpha_{k}\theta}\cdots\frac{(M-1)}{(M-1)+\alpha_{k}\theta}\bigg{]}-\\ &\qquad-\frac{2\alpha_{k}}{(2+\alpha_{k}\theta)^{2}}\left[\frac{1}{1+\alpha_{k}\theta}\cdot\frac{3}{3+\alpha_{k}\theta}\cdots\frac{(M-1)}{(M-1)+\alpha_{k}\theta}\right]-\cdots\\ &\qquad\cdots-\frac{(M-1)\alpha_{k}}{((M-1)+\alpha_{k}\theta)^{2}}\bigg{[}\frac{1}{1+\alpha_{k}\theta}\cdot\frac{2}{2+\alpha_{k}\theta}\cdots\frac{(M-2)}{(M-2)+\alpha_{k}\theta}\bigg{]}\bigg{)}\,,\end{split}

f_{0}^{'} (θ = 0) = - Pr (k ∣ α) \cdot m = 2 \sum M \frac{α _{k}}{m - 1} .

f_{0}^{'} (θ = 0) = - Pr (k ∣ α) \cdot m = 2 \sum M \frac{α _{k}}{m - 1} .

f_{1} (θ) = Pr (y_{k} = M - y_{l} ∣ α, M, θ, one mutation at m = y_{k} + 1) = Pr (k ∣ α) \cdot \frac{1}{1 + α _{k} θ} \cdot \frac{2}{2 + α _{k} θ} \dots \frac{y _{k} - 1}{y _{k} - 1 + α _{k} θ} \cdot \frac{α _{k l} θ}{m - 1 + α _{k} θ} \times \frac{1}{y _{k} + \frac{y _{k}}{y _{k} + 1} α _{k} θ + \frac{1}{y _{k} + 1} α _{l} θ} \dots \frac{y _{l} - 1}{M - 1 + \frac{y _{k}}{M} α _{k} θ + \frac{y _{l}}{M - y _{k}} α _{l} θ}

f_{1} (θ) = Pr (y_{k} = M - y_{l} ∣ α, M, θ, one mutation at m = y_{k} + 1) = Pr (k ∣ α) \cdot \frac{1}{1 + α _{k} θ} \cdot \frac{2}{2 + α _{k} θ} \dots \frac{y _{k} - 1}{y _{k} - 1 + α _{k} θ} \cdot \frac{α _{k l} θ}{m - 1 + α _{k} θ} \times \frac{1}{y _{k} + \frac{y _{k}}{y _{k} + 1} α _{k} θ + \frac{1}{y _{k} + 1} α _{l} θ} \dots \frac{y _{l} - 1}{M - 1 + \frac{y _{k}}{M} α _{k} θ + \frac{y _{l}}{M - y _{k}} α _{l} θ}

Pr (y_{k} = M - y_{l} ∣ α, M, one mutation) = Pr (k ∣ α) \cdot \frac{α _{k l}}{m - 1} \frac{( y _{k} - 1 )! ( y _{m} - 1 )!}{( M - 1 )!} + O (θ^{2})

Pr (y_{k} = M - y_{l} ∣ α, M, one mutation) = Pr (k ∣ α) \cdot \frac{α _{k l}}{m - 1} \frac{( y _{k} - 1 )! ( y _{m} - 1 )!}{( M - 1 )!} + O (θ^{2})

g_{y} (θ) = Pr (k ∣ α) \cdot \frac{1}{1 + α _{k} θ} \frac{2}{2 + α _{k} θ} \dots \frac{y _{k} - 1}{y _{k} - 1 + α _{k} θ} \times \frac{1}{y _{k} + \frac{y _{k}}{y _{k} + 1} α _{k} θ + \frac{1}{y _{k} + 1} α _{l} θ} \dots \frac{y _{l} - 1}{M - 1 + \frac{y _{k}}{M} α _{k} θ + \frac{y _{l}}{M - y _{k}} α _{l} θ},

g_{y} (θ) = Pr (k ∣ α) \cdot \frac{1}{1 + α _{k} θ} \frac{2}{2 + α _{k} θ} \dots \frac{y _{k} - 1}{y _{k} - 1 + α _{k} θ} \times \frac{1}{y _{k} + \frac{y _{k}}{y _{k} + 1} α _{k} θ + \frac{1}{y _{k} + 1} α _{l} θ} \dots \frac{y _{l} - 1}{M - 1 + \frac{y _{k}}{M} α _{k} θ + \frac{y _{l}}{M - y _{k}} α _{l} θ},

g_{y} (θ = 0) = Pr (k ∣ α) \frac{( y _{k} - 1 )! ( y _{l} - 1 )!}{( M - 1 )!} .

g_{y} (θ = 0) = Pr (k ∣ α) \frac{( y _{k} - 1 )! ( y _{l} - 1 )!}{( M - 1 )!} .

\begin{split}f_{1}(\theta)&=\frac{\alpha_{kl}\theta}{(m-1)+\alpha_{k}\theta}g_{y}(\theta)\\ f_{1}^{\prime}(\theta)&=\bigg{(}\frac{\alpha_{kl}\theta}{(m-1)+\alpha_{k}\theta}\bigg{)}^{\prime}g_{m}(\theta)+\frac{\alpha_{kl}\theta}{(m-1)+\alpha_{k}\theta}g_{m}^{\prime}(\theta)\,.\end{split}

\begin{split}f_{1}(\theta)&=\frac{\alpha_{kl}\theta}{(m-1)+\alpha_{k}\theta}g_{y}(\theta)\\ f_{1}^{\prime}(\theta)&=\bigg{(}\frac{\alpha_{kl}\theta}{(m-1)+\alpha_{k}\theta}\bigg{)}^{\prime}g_{m}(\theta)+\frac{\alpha_{kl}\theta}{(m-1)+\alpha_{k}\theta}g_{m}^{\prime}(\theta)\,.\end{split}

f_{1}^{'} (θ = 0) = \frac{α _{k l}}{m - 1} g_{m} (θ = 0) + 0 = Pr (k ∣ α) \cdot \frac{α _{k l}}{m - 1} \frac{( y _{k} - 1 )! ( y _{m} - 1 )!}{( M - 1 )!} .

f_{1}^{'} (θ = 0) = \frac{α _{k l}}{m - 1} g_{m} (θ = 0) + 0 = Pr (k ∣ α) \cdot \frac{α _{k l}}{m - 1} \frac{( y _{k} - 1 )! ( y _{m} - 1 )!}{( M - 1 )!} .

Pr (X_{1} = 1, \dots, X_{k} = 1, X_{k + 1} = 0, \dots, X_{n} = 0) = \int_{0}^{1} ν^{k} (1 - ν^{n - k}) d F (ν) = E [U^{k} (1 - U^{n - k})]

Pr (X_{1} = 1, \dots, X_{k} = 1, X_{k + 1} = 0, \dots, X_{n} = 0) = \int_{0}^{1} ν^{k} (1 - ν^{n - k}) d F (ν) = E [U^{k} (1 - U^{n - k})]

Pr (U \leq ν) = F (ν)

Pr (U \leq ν) = F (ν)

Pr (S_{n} = i \sum n X_{i} = k) = \int_{0}^{1} (k n) ν^{k} (1 - ν^{n - k}) d F (ν) .

Pr (S_{n} = i \sum n X_{i} = k) = \int_{0}^{1} (k n) ν^{k} (1 - ν^{n - k}) d F (ν) .

(1, m - 1) (1001),

(1, m - 1) (1001),

U \sim B (1, m - 1) .

U \sim B (1, m - 1) .

Pr (S_{M - m} ∣ M, m, 1) = \int_{0}^{1} (y _{l} - 1 M - m) ν^{y_{l} - 1} (1 - ν)^{M - y_{l} - m} (1 - ν)^{m - 1} ν^{0} \frac{Γ ( m )}{Γ ( m - 1 ) Γ ( 1 )} d ν = (y _{l} - 1 M - m) \frac{Γ ( M - y _{l} - m + m - 1 ) Γ ( y _{l} )}{Γ ( M - y _{l} - m + m - 1 + y _{l} )} \frac{Γ ( m )}{Γ ( m - 1 ) Γ ( 1 )} = (y _{l} - 1 M - m) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} \frac{Γ ( m )}{Γ ( m - 1 ) Γ ( 1 )} = (y _{l} - 1 M - m) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} \frac{( m - 1 )!}{( m - 2 )! ( 1 - 1 )!} = (y _{l} - 1 M - m) (m - 1) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} .

Pr (S_{M - m} ∣ M, m, 1) = \int_{0}^{1} (y _{l} - 1 M - m) ν^{y_{l} - 1} (1 - ν)^{M - y_{l} - m} (1 - ν)^{m - 1} ν^{0} \frac{Γ ( m )}{Γ ( m - 1 ) Γ ( 1 )} d ν = (y _{l} - 1 M - m) \frac{Γ ( M - y _{l} - m + m - 1 ) Γ ( y _{l} )}{Γ ( M - y _{l} - m + m - 1 + y _{l} )} \frac{Γ ( m )}{Γ ( m - 1 ) Γ ( 1 )} = (y _{l} - 1 M - m) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} \frac{Γ ( m )}{Γ ( m - 1 ) Γ ( 1 )} = (y _{l} - 1 M - m) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} \frac{( m - 1 )!}{( m - 2 )! ( 1 - 1 )!} = (y _{l} - 1 M - m) (m - 1) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} .

Pr (y_{l} ∣ M, α, k, θ \to 0) = m = 2 \sum M - y_{l} + 1 (y _{l} - 1 M - m) (m - 1) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} \frac{α _{k l} θ}{m - 1} + O (θ^{2}) = α_{k l} θ m = 2 \sum M - y_{l} + 1 (y _{l} - 1 M - m) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} + O (θ^{2}) = \frac{α _{k l} θ}{y _{l}} + O (θ^{2}) .

Pr (y_{l} ∣ M, α, k, θ \to 0) = m = 2 \sum M - y_{l} + 1 (y _{l} - 1 M - m) (m - 1) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} \frac{α _{k l} θ}{m - 1} + O (θ^{2}) = α_{k l} θ m = 2 \sum M - y_{l} + 1 (y _{l} - 1 M - m) \frac{Γ ( M - y _{l} ) Γ ( y _{l} )}{Γ ( M )} + O (θ^{2}) = \frac{α _{k l} θ}{y _{l}} + O (θ^{2}) .

(x N) = (x - 1 N - 1) + (x N - 1) \Leftrightarrow (x - 1 N - 1) = (x N) - (x N - 1),

(x N) = (x - 1 N - 1) + (x N - 1) \Leftrightarrow (x - 1 N - 1) = (x N) - (x N - 1),

m = 2 \sum M - y_{l} + 1 (y _{l} - 1 M - m) = (y _{l} - 1 M - 2) + (y _{l} - 1 M - 3) + \dots + (y _{l} - 1 y _{1}) + (y _{l} - 1 y _{l} - 1) = (y _{l} M - 1) - (y _{l} M - 2) + (y _{l} M - 2) - (y _{l} M - 3) + \dots + (y _{l} y _{l} + 1) - (y _{l} y _{l}) + (y _{l} y _{l}) - (y _{l} y _{l} - 1) = (y _{l} M - 1),

m = 2 \sum M - y_{l} + 1 (y _{l} - 1 M - m) = (y _{l} - 1 M - 2) + (y _{l} - 1 M - 3) + \dots + (y _{l} - 1 y _{1}) + (y _{l} - 1 y _{l} - 1) = (y _{l} M - 1) - (y _{l} M - 2) + (y _{l} M - 2) - (y _{l} M - 3) + \dots + (y _{l} y _{l} + 1) - (y _{l} y _{l}) + (y _{l} y _{l}) - (y _{l} y _{l} - 1) = (y _{l} M - 1),

⎩ ⎨ ⎧ Pr (k ∣ α) - Pr (k ∣ α) \cdot \sum_{k \neq = l}^{K} α_{k} θ H_{M - 1} + O (θ^{2}) Pr (k ∣ α) \frac{α _{k l} θ}{y _{l}} + Pr (l ∣ α) \frac{α _{l k} θ}{y _{k}} + O (θ^{2}) y_{k} = M, y_{l \neq = k} = 0 1 \leq y_{k}, y_{l \neq = k} \leq M - 1 and y_{l} + y_{k} = M .

⎩ ⎨ ⎧ Pr (k ∣ α) - Pr (k ∣ α) \cdot \sum_{k \neq = l}^{K} α_{k} θ H_{M - 1} + O (θ^{2}) Pr (k ∣ α) \frac{α _{k l} θ}{y _{l}} + Pr (l ∣ α) \frac{α _{l k} θ}{y _{k}} + O (θ^{2}) y_{k} = M, y_{l \neq = k} = 0 1 \leq y_{k}, y_{l \neq = k} \leq M - 1 and y_{l} + y_{k} = M .

Pr (y ∣ M, α, θ) = ⎩ ⎨ ⎧ \frac{b + d}{2 ( b + d + c + e )} - \frac{b + d}{2 ( b + d + c + e )} (a + c + e) H_{M - 1} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} - \frac{b + d}{2 ( b + d + c + e )} (a + c + e) H_{M - 1} + O (θ^{2}) \frac{c + e}{2 ( b + d + c + e )} - \frac{c + e}{2 ( b + d + c + e )} (b + d + f) H_{M - 1} + O (θ^{2}) \frac{c + e}{2 ( b + d + c + e )} - \frac{c + e}{2 ( b + d + c + e )} (b + d + f) H_{M - 1} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{a}{y _{A}} + \frac{b + d}{2 ( b + d + c + e )} \frac{a}{y _{T}} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{e}{y _{A}} + \frac{c + e}{2 ( b + d + c + e )} \frac{d}{y _{C}} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{c}{y _{A}} + \frac{c + e}{2 ( b + d + c + e )} \frac{b}{y _{G}} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{c}{y _{T}} + \frac{c + e}{2 ( b + d + c + e )} \frac{b}{y _{C}} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{e}{y _{T}} + \frac{c + e}{2 ( b + d + c + e )} \frac{d}{y _{G}} + O (θ^{2}) \frac{c + e}{2 ( b + d + c + e )} \frac{f}{y _{G}} + \frac{c + e}{2 ( b + d + c + e )} \frac{f}{y _{C}} + O (θ^{2}) y_{A} = M, y_{T, G, C} = 0 y_{T} = M, y_{A, G, C} = 0 y_{G} = M, y_{A, T, C} = 0 y_{C} = M, y_{A, T, G} = 0 1 \leq y_{A}, y_{T} \leq M - 1, y_{A} + y_{T} = M 1 \leq y_{A}, y_{C} \leq M - 1, y_{A} + y_{C} = M 1 \leq y_{A}, y_{G} \leq M - 1, y_{A} + y_{G} = M 1 \leq y_{T}, y_{C} \leq M - 1, y_{T} + y_{C} = M 1 \leq y_{T}, y_{G} \leq M - 1, y_{T} + y_{G} = M 1 \leq y_{G}, y_{C} \leq M - 1, y_{G} + y_{C} = M .

Pr (y ∣ M, α, θ) = ⎩ ⎨ ⎧ \frac{b + d}{2 ( b + d + c + e )} - \frac{b + d}{2 ( b + d + c + e )} (a + c + e) H_{M - 1} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} - \frac{b + d}{2 ( b + d + c + e )} (a + c + e) H_{M - 1} + O (θ^{2}) \frac{c + e}{2 ( b + d + c + e )} - \frac{c + e}{2 ( b + d + c + e )} (b + d + f) H_{M - 1} + O (θ^{2}) \frac{c + e}{2 ( b + d + c + e )} - \frac{c + e}{2 ( b + d + c + e )} (b + d + f) H_{M - 1} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{a}{y _{A}} + \frac{b + d}{2 ( b + d + c + e )} \frac{a}{y _{T}} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{e}{y _{A}} + \frac{c + e}{2 ( b + d + c + e )} \frac{d}{y _{C}} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{c}{y _{A}} + \frac{c + e}{2 ( b + d + c + e )} \frac{b}{y _{G}} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{c}{y _{T}} + \frac{c + e}{2 ( b + d + c + e )} \frac{b}{y _{C}} + O (θ^{2}) \frac{b + d}{2 ( b + d + c + e )} \frac{e}{y _{T}} + \frac{c + e}{2 ( b + d + c + e )} \frac{d}{y _{G}} + O (θ^{2}) \frac{c + e}{2 ( b + d + c + e )} \frac{f}{y _{G}} + \frac{c + e}{2 ( b + d + c + e )} \frac{f}{y _{C}} + O (θ^{2}) y_{A} = M, y_{T, G, C} = 0 y_{T} = M, y_{A, G, C} = 0 y_{G} = M, y_{A, T, C} = 0 y_{C} = M, y_{A, T, G} = 0 1 \leq y_{A}, y_{T} \leq M - 1, y_{A} + y_{T} = M 1 \leq y_{A}, y_{C} \leq M - 1, y_{A} + y_{C} = M 1 \leq y_{A}, y_{G} \leq M - 1, y_{A} + y_{G} = M 1 \leq y_{T}, y_{C} \leq M - 1, y_{T} + y_{C} = M 1 \leq y_{T}, y_{G} \leq M - 1, y_{T} + y_{G} = M 1 \leq y_{G}, y_{C} \leq M - 1, y_{G} + y_{C} = M .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsEvolution and Genetic Dynamics · Genetic Associations and Epidemiology · Genetic diversity and population structure

Full text

Maximum likelihood (ML) estimators for scaled mutation parameters with a strand symmetric mutation model in equilibrium

Claus Vogl

[email protected]

Lynette Caitlin Mikula

[email protected]

Department of Biomedical Sciences, Vetmeduni Vienna, Veterinärplatz 1, A-1210 Wien, Austria

Centre for Biological Diversity, School of Biology, University of St. Andrews, St Andrews KY16 9TH, UK

Abstract

With the multiallelic parent-independent mutation-drift model, the equilibrium proportions of alleles are known to be Dirichlet distributed. A special case is the biallelic model, in which the proportions are beta distributed. A sample taken from these models is then Dirichlet-multinomially or beta-binomially distributed, respectively. Maximum likelihood (ML) estimators for the mutation parameters of the biallelic parent-independent mutation model are available via an expectation maximization algorithm. Assuming small scaled mutation rates, the distribution of a sample of size $M$ can be expanded in a Taylor series of first order. Then the ML estimators for the two parameters in the biallelic model can be expressed using the site frequency spectrum. In this article, we go beyond parent-independent mutation and analyse a strand-symmetric mutation model with six scaled mutation parameters that deviates from parent independent mutation and, generally, from detailed balance. We derive ML estimators for these six parameters assuming mutation-drift equilibrium and small scaled mutation rates. This is the first time that ML estimators are provided for a mutation model more complex than parent-independent mutation.

keywords:

strand-symmetric mutation, mutation-drift model , scaled mutation parameters , maximum likelihood inference , expectation-maximization algorithm.

††journal: Theoretical Population Biology

1 Introduction

With the parent-independent mutation-drift model, theoretical results for the estimation of mutation parameters are available: The proportions of alleles of a biallelic model in equilibrium are beta distributed [18]. Data for the inference of parameters are usually in the form of site frequency spectrum (also called allele frequency spectrum) data taken from one population. Such data are a vector of allele frequencies for samples of $M$ haploid individuals at a total of $L$ sites (or loci) whereby there are $L_{0},L_{1},\dots,L_{M}$ sites with $y=0,1,\dots,M$ alleles of the focal type. Given a sample of size $M$ and a binomial sampling distribution conditional on the allele proportion, one obtains a beta-binomial distribution for the allele frequencies. Maximum likelihood (ML) estimators can be constructed from the beta-binomial distribution via an expectation maximization algorithm, where a root of a polynomial of the order of the sample size needs to be evaluated at each iteration [16]. With a multiallelic parent-independent mutation model, a Dirichlet-multinomial distribution for the allele frequencies follows analogously.

In the limit of small scaled mutation rates, it is convenient to reparametrize the two parameters of the beta distribution for equilibrium allele proportions with the mutation bias $\alpha$ and the overall mutation rate $\theta$ . Expanding the beta-binomial distribution in a Taylor series up to first order of the scaled mutation rate $\theta$ then results in a relatively simple equilibrium sample distribution [16]:

[TABLE]

whereby $H_{m}$ denotes the $m$ th harmonic number. Note that since $H_{M-1}$ approaches infinity logarithmically, the sample size must stay below a limit of approximately $M<e^{\frac{1}{\max(1-\alpha,\alpha)\theta}}$ . Furthermore, simulations have shown that the mutation rates for the first order Taylor series expansion must be $2\alpha(1-\alpha)\theta\leq 0.025$ [17]; this upper boundary for the scaled mutation rate was further substantiated by Schrempf and Hobolth [14]. It is higher than the scaled mutation rate of most eukaryotes [12]. With this approximation, simple ML estimators given site-frequency spectrum data are obtained:

[TABLE]

and

[TABLE]

respectively [16].

Vogl and Clemente [17] proposed a Moran model, in which mutations are restricted to monomorphic states. In this so-called boundary-mutation Moran model, the same equilibrium sample distribution is obtained as with the first order Taylor series expansion in formula (1). Schrempf and Hobolth [14] were able to explicitly derive the approximate neutral multiallelic stationary distribution of allele frequencies as the equilibrium state of a discrete boundary mutation Moran model with a completely general mutation model. Starting from a Wright-Fisher model and passing to the diffusion limit, Burden and Griffiths [1] showed that a Taylor series expansion in the scaled mutation rate can be used to obtain the stationary distribution for mutation models more complex than the parent-independent mutation model. This is remarkable because the equilibrium distribution of such models is not known. Below, we provide an alternative proof.

Our main results, however, are ML estimators for the scaled mutation parameters of a strand-symmetric mutation model. The two strands of the DNA double helix are held together in anti-parallel orientation by hydrogen bonds. Assuming symmetry between the two strands, the twelve mutation rate parameters between the four bases adenine, thymine, cytosine, and guanine ( $A$ , $T$ , $C$ , and $G$ ) are reduced to six because the transition rates of complementary bases are identical, e.g., the transition rate from $A$ to $C$ is identical to that from $T$ to $G$ . For these six parameters, we provide estimators using expectation-maximization (EM) algorithms that are guaranteed to converge to the global maximum.

2 The strand-symmetric mutation-drift model

At each site, we assume $K=4$ alleles corresponding to the four bases. The evolutionary process is described by a continuous time Markov chain with instantaneous transition rate matrix $\mathbf{Q}$ of $dim(\mathbf{Q})=4$ x $4$ and corresponding forward transition probability matrix $\mathbf{P}$ . Under the most general mutation model, the off-diagonal entries of $\mathbf{Q}$ are the strictly positive transition rates among the four bases $A$ , $T$ , $G$ , and $C$ :

[TABLE]

Thus the general transition rate matrix has twelve parameters. Finding the unique stationary state for a general transition rate matrix $\mathbf{Q}$ requires scaling the individual entries to ensure they are $\leq 1$ . Applying the Perron-Frobenius theorem to the stochastic matrix given by the sum of this scaled transition rate matrix and the identity matrix yields a dominant eigenvalue of zero. All further eigenvalues have a negative real component. Provided the process is irreducible, the eigenvector associated with the zero-valued eigenvalue corresponds to the unique stationary distribution of the process.

With the strand-symmetric mutation model, similar conclusions can be drawn without invoking the Perron-Frobenius theorem. The number of parameters of $\mathbf{Q}$ is reduced to six because the transition rates between bases are identical on complementary DNA strands. The system is thereby simplified and $\mathbf{Q}$ can be written as in Lobry [11] (note that our $\mathbf{Q}$ is the transpose of Lobry’s):

[TABLE]

. Again $a,..,f$ are strictly positive and the rows of $\mathbf{Q}$ sum to zero. The leading eigenvalue of $\mathbf{Q}$ is $\lambda_{1}=0$ with the associated normed eigenvector

[TABLE]

This vector corresponds to the stationary distribution [11]. Clearly, the probabilities of the four bases at stationarity follow Chargaff’s second parity rule: For each of the two DNA strands, the expectations of the proportions of $A$ and $T$ are identical, as are those of $C$ and $G$ . Chargaff’s second parity rule has been shown to hold for all types of double-stranded DNA except organellar DNA [13].

The second eigenvalue [11] is $\lambda_{2}=-(b+d+c+e)$ with associated eigenvector

[TABLE]

This eigenvector contrasts the sum of $A$ and $T$ with the sum of $C$ and $G$ .

The next two eigenvalues ( $\lambda_{3},\lambda_{4}$ ) may have complex components that introduce a probability flow into the system [11]. However, it must be noted that if all transition rates are similar in magnitude, this potentially imaginary term fluctuates around zero. The stationary distribution given by the leading eigenvector is stable despite the complex eigenvalues because the real components of $\lambda_{2,r;3,r;4,r}$ are strictly negative [11]. Provided the individual transition rates are indeed similar in magnitude, the real components of $\lambda_{3,r}$ and $\lambda_{4,r}$ are equal to $-a-f+\frac{\lambda_{2}}{2}$ and are of the same order of magnitude as $\lambda_{2}$ (The sign in front of the fraction is incorrect in Lobry’s listing of eigenvalues; his further arguments indicate that this is merely a typo). Thus the evolutionary process converges exponentially towards a stationary state at a rate that depends on $\lambda_{2,3,4}$ .

On the basis of Lobry’s work, we can be confident of a unique stationary distribution given the transition rates among alleles in the strand-symmetric mutation model except in degenerate and thus biologically meaningless cases.

3 The stationary distribution of the strand-symmetric mutation-drift model with small scaled mutation rates

In this section, we derive the stationary distribution for a sample of size $M$ taken from the general multiallelic boundary-mutation Moran model in the limit of small scaled overall mutation rates. It is assumed that the mutation matrix gives rise to a unique stationary distribution. This is a similar approach to that of Burden and Griffiths [1] and recovers the full stationary distribution for the proportions of alleles (their formulas 7-9 in Theorem 1) [see also 14].

We reparametrize the general mutation matrix as follows:

[TABLE]

where $\alpha_{k}=\sum_{l=1,l\neq k}^{K}\alpha_{kl}$ .

Assuming a unique stationary distribution for $\mathbf{Q}$ and recalling the duality between the Moran model and Kingman’s coalescent [6, chapt. 2], we can use a sampling algorithm proposed by Stephens and Donnelly [15] to build a genealogical realization of a Moran model of sample size $M$ forwards in time. This sampling algorithm is essentially an urn sampling process. Hoppe [9] first drew the connection between Polya-like urn processes and Kingman’s coalescent; Donnelly and Kurtz [5] proposed an algorithm by which one can obtain the realization of a sample from a particle Fleming-Viot process at stationarity forwards in time (The Moran model falls within the attraction of the Fleming-Viot processes). The Stephens-Donnelly algorithm is the continuous version of such particle sampling algorithms:

Start with a sample of size $m=1$ . Randomly select an allele of type $k$ from $1\leq k\leq K$ with probability $\Pr(k{\,|\,}\vec{\alpha})$ taken from the stationary distribution characterized by $\mathbf{Q}$ . Split immediately into two lineages of this type. 2. 2.

At this point there are $m$ lineages in the ancestry. Wait for an exponentially distributed time with rate $m(m-1+\theta)$ and then select an ancestral lineage at random. Split it with probability $m(m-1+\theta)$ , otherwise introduce a mutation. 3. 3.

Hit a predetermined stopping criterion at sample size $M+1$ , then go back to the time when there were $M$ ancestral lineages.

4 Sampling Paths

4.1 Taylor Expansion

Given a sample of alleles drawn with the Stephens-Donnelly algorithm, one can explicitly calculate the probabilities of the possible sampling paths that lead to the observed sample.

A monomorphic sample of size $M$ is either the result of a sampling path that consists purely of splitting ancestries or of one that has an even number of reversible mutations (which we also take to include the unlikely case of four mutations through all variants back to the original allele). In the former case, the sampling path can be written as follows:

[TABLE]

With $\theta\ll 1$ and $\alpha_{i}$ , $\alpha_{ij}$ of order one or smaller, a Taylor expansion of the full sampling path at $\theta=0$ is to first order:

[TABLE]

Here the full derivation: We have

[TABLE]

and

[TABLE]

such that

[TABLE]

Note that $m$ can be interpreted as the sample size at which a mutation could potentially occur.

A sampling path with an even number of mutations that restores the monomorphic condition is unlikely with low mutation rates. Indeed it is easy to show that, in this case, probabilities of such sampling paths are at least of second order in $\theta$ and thus do not contribute to the first order approximation.

Consider now the sampling paths that create polymorphic samples of size $M$ : The simplest of these includes one mutation, e.g., when the sample size $m$ is equal to $y_{k}+1$ :

[TABLE]

A Taylor expansion of the probability of this sampling path at $\theta=0$ yields:

[TABLE]

The derivation is similar to those above: Set

[TABLE]

such that

[TABLE]

Then we have

[TABLE]

Setting $\theta=0$ , we get

[TABLE]

The occurrence of two or three mutations that lead to a polymorphic sample with two, three, or four segregating alleles is a theoretical possibility. However, the expansion of the sampling path probabilities would again approximate to zero for the same reason as in the case of (multiple) reversible mutations in a monomorphic sampling path.

Note that the mutation may in principle occur at any sample size between $2\leq m_{k}\leq(y_{k}+1)$ . Hence, while the series expansion of the monomorphic sampling paths accounts for the only possible sampling scenario, the expansion of a polymorphic sampling path represents one of several feasible branching structures along the sampling algorithm that result in the given configuration of alleles. We thus need to sum over these possibilities.

4.2 Sum of Ordered Probabilities

In the following subsection, we show that the sampling distribution of a polymorphic sampling path, conditional on the mutation occurring at sample size $m$ , is beta-binomial. The argument is similar to that in Burden and Griffiths [1]. These authors use a direct result equivalent to the deFinetti density of a Polya urn model. We start from the properties of the sampling path and build up the deFinetti representation.

According to deFinetti’s theorem [3], there exists for every infinite sequence of Bernoulli random variables a probability distribution $F$ on $\left[0,1\right]$ so that

[TABLE]

for a random variable $U$ , whereby

[TABLE]

and the $(X_{i})_{i}\leq 1$ are i.i.d. when conditioned on a $U$ that fulfills $\Pr(X_{i}=k)=U$ . The following also holds:

[TABLE]

The joint distribution of the polymorphic sampling path, which is a finite number of draws of alleles, does not change if the order of alleles within the path is changed. We start from an initial allele of type $k$ , then there is a mutational event to allele $l$ at a random time point. At this point, $m-1$ alleles of type $k$ are already in the sample. Each subsequent draw can increase the number of sampled alleles of either type by one. This is exactly the specification of the following Polya urn process:

[TABLE]

which in our case runs until we have $y_{l}-1$ alleles of type $l$ and $y_{k}-m=M-y_{l}-m$ alleles of type $k$ . By noting that the probabilities of specific sampling paths are sequences of beta distributed moments and recalling Hausdorff’s moment problem [7, 8], it follows that

[TABLE]

Therefore, the probability of a process that yields $y$ alleles of type $l$ in a sample size of $M-m$ (denoted by $S_{M-m}$ ) is given by the following:

[TABLE]

For a given proportion $y_{l}$ of mutant alleles, we need to sum over all possible mutation points, i.e., $2\leq m\leq(M-y_{l}+1)$ . Given that the first allele sampled at sample size $m=1$ is of type $k$ , we have:

[TABLE]

Line three follows from line two by repeatedly applying the identity:

[TABLE]

from which it follows that

[TABLE]

since $\binom{y_{l}-1}{y_{l}}=0$ . Note that equation (19) also appears in Burden and Griffiths [1], but is solved differently.

This provides an alternative proof of Theorem 1 in Burden and Griffiths [1] and concludes the characterization of the polymorphic sampling paths.

5 Stationary Distribution

On slightly reparametrizing the results of the previous two subsections, the probability of obtaining the overall configuration $y$ of alleles can be concisely formulated: Given a fixed sample size $M<e^{\frac{1}{\alpha_{k}\theta}}$ , mutation biases of maximal order $1$ , and a small scaled overall mutation rate $\theta=N\mu<0.1$ , the stationary allelic configuration $Y$ to the first order of $\theta$ is the distribution $\Pr(y\mid M,\vec{\alpha},\theta)$ :

[TABLE]

All other possibilities have probabilities of at most $\mathcal{O}(\theta^{2})$ and thus do not contribute to the distribution. This stationary distribution is equivalent to that derived by Burden and Griffiths [1], as expected.

6 Maximum likelihood (ML) estimators of the scaled mutation rate parameters

Starting from the general stationary distribution and recalling that we are certain of a unique stationary distribution of $\mathbf{Q}$ in the case of strand symmetry (i.e., we can build a sample of size $M$ using the Stephens-Donnelly algorithm if we assume strand symmetry), we now aim to infer the six parameters of the strand-symmetric boundary-mutation Moran model given site frequency spectrum data.

Recalling the parametrization for $\mathbf{Q}$ in formula (5) as visualized in Fig 1, we set $\theta=-\lambda_{2}=b+d+c+e$ . Then $\alpha_{1}=\frac{a+c+e}{b+d+c+e}$ , $\alpha_{12}=\frac{a}{b+d+c+e}$ , $\alpha_{13}=\frac{c}{b+d+c+e}$ , $\alpha_{14}=\frac{e}{b+d+c+e}$ , etc.

The full stationary distribution for the strand-symmetric model is then:

[TABLE]

In the next subsections, we will derive ML estimators for the parameter vector $(a,b,c,d,e,f)$ in the distribution (21) by linear transformation.

6.1 Biallelic Model

Grouping the bases $A$ and $T$ as well as $G$ and $C$ together results in a biallelic model [16]: If one interprets the allelic state $k$ as representing the grouped state $A-T$ including those sites polymorphic for $A$ / $T$ and the allelic state $l$ as representing the $G-C$ alleles and the polymorphic $G/C$ sites, $\alpha_{kl}$ and $\alpha_{lk}$ become equal to $\alpha_{k}$ and $\alpha_{l}$ respectively (see Figure 2).

In particular, this biallelic model enables us to infer $\Pr(k=AT\mid\vec{\alpha})=\frac{b+d}{b+d+c+e}$ and $\Pr(l=GC\mid\vec{\alpha})=\frac{e+c}{b+d+c+e}$ . Using the frequencies of $A$ and $T$ , one can then infer $a$ using a symmetric version of the biallelic model and similarly infer $f$ from the frequencies of $G$ and $C$ . The remaining transition rates can be estimated via the expectation-maximization algorithm.

Following Vogl [16], we will first find the ML estimators for $\Pr(k=AT\mid\vec{\alpha})=\frac{b+d}{b+d+c+e}=\beta$ and $\Pr(l=GC\mid\vec{\alpha})=\frac{e+c}{b+d+c+e}=1-\beta$ using $\alpha_{13}+\alpha_{14}=1-\beta$ and starting from the following biallelic version of 20:

[TABLE]

The estimators are to be expressed via site frequency data: We assume knowledge of the allele frequency counts at $L<\infty$ loci indexed as $1\leq l\leq L$ and further assume that these counts are based on a sample of size $M$ at each locus (the constant sample size merely simplifies the notation and calculations and is not a theoretical necessity). $L_{0}$ is the number of loci at which there are no occurrences of the state $A-T$ (i.e., the number of sites monomorphic for $G$ or $C$ , or $G/C$ polymorphic), and $L_{M}$ is the number of loci either monomorphic for $A$ or $T$ , or $A/T$ polymorphic. The remaining $L_{y}$ with $L_{p}=\sum_{y=1}^{M-1}L_{y}$ are the observed polymorphic counts between the $A-T$ and $G-C$ states. We define

[TABLE]

to efficiently handle products of random variables.

The likelihood of the site frequency spectrum then becomes:

[TABLE]

whereby $H_{m}$ again denotes the $m$ th harmonic number.

Then the log-likelihood can be written as proportional to:

[TABLE]

With asymmetric mutation rates between the two classes, we can take the derivative of the log-likelihood by $\beta$ :

[TABLE]

It follows that

[TABLE]

and therefore

[TABLE]

and similarly

[TABLE]

Taking the derivative of the log-likelihood by $2\gamma H_{M-1}$ and substituting we get:

[TABLE]

It follows that

[TABLE]

This corresponds to the estimator in formula (36) in [16]. Analogously, we get:

[TABLE]

Similarly, we have

[TABLE]

respectively yielding

[TABLE]

which corresponds to the estimator in formula (37) in [16], and

[TABLE]

Thus we provide ML estimators of the parameters of the combined states $A-T$ and $G-C$ in the biallelic model resulting from combining $A-T$ and $C-G$ under stationarity.

6.2 Symmetric Biallelic Model

Next we recover estimators for the transition rates $a$ and $f$ from symmetric biallelic systems [see 2]. Let us focus on $a$ (see Figure 3): Recalling Chargaff’s second parity rule and its correspondence with the stationary distribution of base frequencies, we can write

[TABLE]

Furthermore, we know that

[TABLE]

from which it follows that

[TABLE]

Subsequently the distribution (20) can be rewritten for a symmetric biallelic A-T model using the previously estimated parameters:

[TABLE]

This yields the following log-likelihood of the site frequency spectrum:

[TABLE]

Taking the derivative by $a$ , one can calculate the ML estimator for $a$ as follows:

[TABLE]

From earlier we have

[TABLE]

such that

[TABLE]

and furthermore

[TABLE]

As a result the ML estimator for $a$ becomes:

[TABLE]

Analogously, we get

[TABLE]

6.3 Expectation-Maximization Algorithm

The remaining transition rates $(b,d,e,c)$ cannot be disentangled and expressed individually through reparametrization in the same way as $a$ and $f$ . Instead, we will use the expectation-maximization (EM) algorithm [4] to obtain estimators of these transition rates as sequential updates that cycle through the parameter pairs $d,e$ and $b,c$ , respectively.

Let us focus on a biallelic A-C system (see Figure 4) and rewrite 20 to obtain the appropriate stationary distribution:

[TABLE]

Let $L_{i}$ denote the number of polymorphic samples in the frequency spectrum, whereby $y_{A}=i$ and $y_{C}=M-i$ . In order to distinguish between mutations of different directions, we introduce the auxiliary variable $z_{i}$ ( $0\leq z_{i}\leq L_{i}$ ) that counts the mutations from $A$ to $C$ . Conversely, $(L_{i}-z_{i})$ counts the mutations from $C$ to $A$ . This yields the complete data log-likelihood:

[TABLE]

The expectation of $z_{i}$ corresponds to the mean of a binomial distribution with sample size $L_{i}$ and $p_{i}=\frac{\frac{e}{i}}{\frac{e}{i}+\frac{d}{M-i}}$ :

[TABLE]

The expectation step of the EM-algorithm constitutes taking the expectation of the full log-likelihood and noting that only the part $\operatorname{Q}(d,e{\,|\,}d^{t},e^{t})=\operatorname{E}(\Pr(y,z_{i}{\,|\,}L_{0},\dots,L_{i},\dots,L_{M},\dots))$ of the function needs to be maximized:

[TABLE]

To calculate the parameter updates according to the maximization step, we must take the appropriate derivatives of $\operatorname{Q}(d,e{\,|\,}d^{t},e^{t})$ . For the calculation of $e^{t+1}$ , we take the derivative by $e$ and set it to zero:

[TABLE]

where we substituted $(1-\beta)\theta$ for $c+e$ to obtain the third line.

Substituting the ML estimators for $a$ and $(1-\beta)\theta$ in the numerator of the first factor yields:

[TABLE]

Furthermore, we have:

[TABLE]

Therefore, we get:

[TABLE]

Similarly, one obtains $d^{t+1}$ by taking the derivative of $Q(d,e{\,|\,}d^{t},e^{t})$ by $d$ and setting it to zero.

The overall iteration scheme is then given by:

[TABLE]

Considering a biallelic A-G system and using the EM-algorithm analogously, the parameter updates for $c$ and $b$ can also be determined:

[TABLE]

For both pairs of parameters, cyclical calculation of estimators guarantees convergence towards a local maximum of the marginal likelihood by properties of the EM-algorithm. Furthermore: $\Pr(y{\,|\,}L_{0},\dots,L_{i},\dots,L_{M},\dots)$ is the distribution of counts of each type of segregating allele. The configurations of the types of alleles themselves are constructed via Polya-like urn processes. As such, the marginal likelihood $\Pr(y{\,|\,}L_{0},\dots,L_{i},\dots,L_{M},\dots)$ takes the form of a Dirichlet-multinomial distribution, which is known to be unimodal [10]. Therefore the estimators determined via the expectation-maximization algorithm converge towards the global optimum of the marginal distribution.

7 Summary

An alternative derivation of the distribution of a sample of size $M$ is derived given a general mutation model under the assumptions of small scaled mutation rates and mutational equilibrium. Further assuming a standard four allele DNA model with strand-symmetric mutation on complementary DNA strands and available site frequency spectrum data, ML estimators for the six scaled mutation parameters are determined. This is the first time such estimators are provided for a mutation model more complex than the parent-independent mutation model.

Acknowledgments

This article is a condensed version of LCM’s master’s thesis supervised by CV. CV wants to thank Juraj Bergman and Sandra Peer for discussions. CV’s research is supported by the Austrian Science Fund (FWF): DK W1225-B20; LCM’s by the School of Biology at the University of St.Andrews.

References

Burden and Griffiths [2018]

Burden, C. and Griffiths, R. (2018).

The stationary distribution of a sample from the Wright-Fisher diffusion model with general small mutation rates.

Journal of Mathematical Biology, 78, 1211–1224.

Burden and Tang [2016]

Burden, C. and Tang, Y. (2016).

An approximate stationary solution for multi-allele neutral diffusion with low mutation rates.

Theoretical Population Biology, 112, 22–32.

de Finetti [1931]

de Finetti, B. (1931).

Funzione caratteristica di un fenomeno aleatorio.

Atti della R. Academia Nazionale dei Lincei, Serie 6. Memorie, Classe di Scienze Fisiche, Mathematice e Naturale, 4, 251–299.

Dempster et al. [1977]

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).

Maximum Likelihood from Incomplete Data via the EM Algorithm.

Journal of the Royal Statistical Society. Series B, 39(1), 1–38.

Donnelly and Kurtz [1996]

Donnelly, P. and Kurtz, T. (1996).

A Countable Representation of the Fleming-Viot Measure-Valued Diffusion.

The Annals of Probability, 24(2), 698–742.

Etheridge [2012]

Etheridge, A. (2012).

Some Mathematical Models from Population Genetics: Lecture Notes in Mathematics.

Springer Verlag. Berlin, Heidelberg.

Hausdorff [1921a]

Hausdorff, F. (1921a).

Summationsmethoden und Momentenfolgen 1.

Mathematische Zeitschrift, 9, 74–109.

Hausdorff [1921b]

Hausdorff, F. (1921b).

Summationsmethoden und Momentenfolgen 2.

Mathematische Zeitschrift, 9, 280–299.

Hoppe [1987]

Hoppe, F. M. (1987).

The Sampling Theory of Neutral Alleles and an Urn Model in Population Genetics.

Journal of Mathematical Biology, 25, 123–159.

Levin and Reeds [1977]

Levin, B. and Reeds, J. (1977).

Compound Multinomial Likelihood Functions are Unimodal: Proof of a Conjecture of I.J. Good.

The Annals of Statistics, 5(1), 79–87.

Lobry [1995]

Lobry, J. (1995).

Properties of a General Model of DNA Evolution Under No-Strand-Bias Conditions.

Journal of Molecular Evolution, 40, 326–330.

Lynch et al. [2016]

Lynch, M., Ackerman, M., Gout, J., Long, H., Sung, W., Thomas, W., and Foster, P. (2016).

Genetic drift, selection and the evolution of the mutation rate.

Nature, 17, 704–714.

Mitchell and Bridge [2006]

Mitchell, D. and Bridge, R. (2006).

A test of Chargaff’s second rule.

Biochemical and Biophysical Research Communications, 340(1), 90–94.

Schrempf and Hobolth [2017]

Schrempf, D. and Hobolth, A. (2017).

An alternative derivation of the stationary distribution of the multivariate neutral Wright-Fisher model for low mutation rates with a view to mutation rate estimation from site frequency data.

Theoretical Population Biology, 114, 88–94.

Stephens and Donnelly [2000]

Stephens, M. and Donnelly, P. (2000).

Inference in Molecular Population Genetics.

Journal of the Royal Statistical Society, Series B, page Discussion Paper.

Vogl [2014]

Vogl, C. (2014).

Estimating the scaled mutation rate and mutation bias with site frequency data.

Theoretical Population Biology, 98, 19–27.

Vogl and Clemente [2012]

Vogl, C. and Clemente, F. (2012).

The allele-frequency spectrum in a decoupled Moran model with mutation, drift, and directional selection, assuming small mutation rates.

Theoretical Population Biology, 81(3), 197–209.

Wright [1931]

Wright, S. (1931).

Evolution in Mendelian populations.

Genetics, 16, 97–159.

Bibliography18

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Burden and Griffiths [2018] Burden, C. and Griffiths, R. (2018). The stationary distribution of a sample from the Wright-Fisher diffusion model with general small mutation rates. Journal of Mathematical Biology , 78 , 1211–1224.
2Burden and Tang [2016] Burden, C. and Tang, Y. (2016). An approximate stationary solution for multi-allele neutral diffusion with low mutation rates. Theoretical Population Biology , 112 , 22–32.
3de Finetti [1931] de Finetti, B. (1931). Funzione caratteristica di un fenomeno aleatorio. Atti della R. Academia Nazionale dei Lincei, Serie 6. Memorie, Classe di Scienze Fisiche, Mathematice e Naturale , 4 , 251–299.
4Dempster et al. [1977] Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society. Series B , 39(1) , 1–38.
5Donnelly and Kurtz [1996] Donnelly, P. and Kurtz, T. (1996). A Countable Representation of the Fleming-Viot Measure-Valued Diffusion. The Annals of Probability , 24(2) , 698–742.
6Etheridge [2012] Etheridge, A. (2012). Some Mathematical Models from Population Genetics: Lecture Notes in Mathematics . Springer Verlag. Berlin, Heidelberg.
7Hausdorff [1921 a] Hausdorff, F. (1921 a). Summationsmethoden und Momentenfolgen 1. Mathematische Zeitschrift , 9 , 74–109.
8Hausdorff [1921 b] Hausdorff, F. (1921 b). Summationsmethoden und Momentenfolgen 2. Mathematische Zeitschrift , 9 , 280–299.