In MoisesExpositoAlonso/popgensim: Wright Fisher simulations: C++ implementation with interactive Shiny app visualization

Genotypes and selection model

Given two loci, A and B, with alleles A and a, and B and b, there must be nine possible gammetes. If the A allele and B allele are selected with coefficients $s_{A}$ and $s_{B}$, and we assume no dominance and no epistasis, we can write the average fitness per genotype. For simplicity we use an haploid model, which would be similar to the diploid homozygotes case (also the case of highly self-fertilizing species):

\ \begin{center} \begin{tabular}{ c | c c} & B & b \ \hline A & $[(1+s_A)(1+s_B)]^{\epsilon} $&$ [(1+s_{A}) (1-s_{B})]^e $\
a & $[(1- s_{A})(1+ s_{B})]^e$&$ [(1-s_A)(1-s_B)]^{\epsilon} $
\end{tabular} \end{center} \

We could then expect that the observed fitness $y$ of one of $N$ number of haplotype of $p$ loci is expressed as:

$$ y = \prod_{i=1}^{p} (1+s_i \odot X_{i}) $$

The probabilistic model

We can assume the data $y_h$ follows a Gamma distribution (typically with $\alpha$ and $beta$ parameters or $\alpha = k$ and $\theta = 1/\beta$),

$$ Y \sim \Gamma ( \alpha , \beta ) $$

The first moment is:

$$ E[Y] = \alpha / \beta = 1 $$

The second moment is then:

$$ Var(y) = E[Y^2]-E[Y] = \alpha / \beta^2 = \sigma^{-1} = N^{-1} \sum_{h \in haps} \prod_{i=1}^{p} (1+s_i \odot X_{hi})^2- 1 $$ So,

$$ \sigma = \frac{N}{\sum_{h \in haps} \prod_{i=1}^{p} (1+s_i \odot X_{hi})^2} - 1 $$

This is the expected shape of the distribution given a configuration of genotypes and selection coefficients.

Thus, $\alpha = \beta$ so we use instead $\sigma$. This special type of $\Gamma$ distribution will always generate a distribution whose mean is in 1, as any relative fitness situation. What contains the information is the shape of the gamma distribution, hence our use of $\sigma$. Theory says that a sum of $h$ Gamma distributions $y_h \sim \Gamma(\alpha_h,\beta_h) $ produce a Gamma distribution with the same shape but with alpha equal to the sum of all alphas: $Y \sim \Gamma( \sum\alpha_h,\beta) $.

Therefore we can expect that depending on the ability of the genotype to survive, it might follow the shape of the population but have a different mean. Then:

$$\alpha_h = \prod_{i=1}^{p} (1+s_i \odot X_{hi}) $$ $$ \sigma=\alpha=\sum_{h \in haps} \alpha_h = \sum_{h \in haps} \prod_{i=1}^{p} (1+s_i \odot X_{hi})$$

The probability density function of Gamma in terms of $k=alpha$ and $\theta=1/ \beta$ is:

$$f(y) = \frac{1}{\Gamma(k) \theta^{k}} y^{k-1} e^{-\frac{y}{\theta}}$$ In log form this can be written as:

$$\ln f(y) = (k-1) \ln(y) {-\frac{y}{\theta}} - \ln \Gamma(k) - k \ln \theta $$ By substituting the observed $y$ fitness of every $h$ haplotype, the population distribution $\theta = 1/ \sigma$, and each haplotype's $\alpha_h = k_h$, we can optain the log-likelihood of all observations and sum them:

$$ - \ell (\vec{s} ) = - \sum_{h \in haps} \ln f(y|\sigma,\vec{\alpha})$$ So the vector $ \vec{s}=s(s_1, ..., s_n)$ are all selection coefficients that follow a multivariate $ \mathcal{N}(m,s)$

MoisesExpositoAlonso/popgensim documentation built on May 7, 2019, 8:57 p.m.