A revised Fisher model on analysis of quantitative trait loci with multiple alleles

Zeng et al. (2005) proposed a general two-allele (G2A) model to model bi-allelic quantitative trait loci (QTL). Comparing with the classical Fisher model, the G2A model can avoid using redundant parameters and be fitted directly using standard least square (LS) approach. In this study, we further extend the G2A model to general multi-allele (GMA) model. First, we propose a one-locus GMA model for phase known genotypes based on modeling the inheritance of paternal and maternal alleles. Next, we develop a one-locus GMA model for phase unknown genotypes by treating it as a special case of the phase known one-locus GMA model. Thirdly, we extend the one-locus GMA models to multiple loci. We discuss how the genetic variance components can be analyzed using these GMA models in equilibrium as well as disequilibrium populations. Finally, we apply the GMA model to a published experimental data set.

1. Introduction

Currently there are two types of statistical genetic models that are commonly used in genetic analysis of quantitative traits. One is the F_∞ type models that concentrate on direct modeling of the expected genotypic values at targeted quantitative trait loci (QTL) or genetic markers and association testing for various allelic effects and interactions (Fisher, 1918; Cheverud, 2000; Hansen and Wagner, 2001; Wang, 2011). Another is Fisher's analysis of variance (ANOVA) models that focus on assessing variations contributed by some genetic components (i.e., grouped allelic effects or allelic interactions) at targeted QTL or genetic markers (Fisher, 1918; Cockerham, 1954, 1963; Kempthorne, 1969; Weir and Cockerham, 1977; Wang and Zeng, 2006). A considerable amount of discussion has been made about the distinction between these two different types of genetic models (Zeng et al., 2005; Álvarez-Castro and Carlborg, 2007; Yang and Álvarez-Castro, 2008; Wang and Zeng, 2009). More recently, a comprehensive review of various genetic models was also given in Álvarez-Castro (2012).

In genetic association studies, we are often interested in direct comparison of the expected genotypic values at certain QTL or marker loci. The F_∞ models are appealing in this setting due to their simplicity in interpretation of their model parameters, which are often referred as the fixed genetic effects such as the additive and dominance effects or the allelic effects and allelic interactions in terms of the expected genotypic values. By applying the F_∞ models, we can compare the expected genotypic values via hypothesis tests on various fixed genetic effects. However, as pointed out in Wang and Zeng (2009), the p-value based association tests on these fixed genetic effects could highly depend upon the regression model, the distribution assumption and the available sample size. Besides, a statistically significant genetic effect with a small enough p-value may not necessarily imply a clinically important finding. On the other hand, the Fisher type models allow us to assess the genetic variations contributed by certain genetic components to the overall variation of a quantitative trait. By definition, these genetic variance components do not depend on the sample size, and they can provide additional information on better understanding the genetic etiology and assessing for clinical importance. Nonetheless, both the F_∞ and the Fisher type models form basis in the analysis of quantitative traits. They provide different perspectives in assessing the genetic effects of QTL or genetic markers.

The basic genetic model on assessing the genetic variance components was first proposed by Fisher (1918). Cockerham (1954, 1963) extended Fisher's one-locus model to two bi-allelic QTL with a particular focus on epistatic variance components. Kempthorne (1954, 1957) further extended the two-locus model to multiple alleles. Wang and Zeng (2006) also explored the Fisher type multi-allele two-locus model on partition of the genotypic variance in the presence of linkage disequilibrium (LD). With their model parameters referred as the average effect of the gene substitution (see Falconer and Mackay, 1996) or the average allelic effects and interactions (Wang and Zeng, 2009), these classical Fisher type models often contain constraints on their model parameters due to an over-parameterization of the expected genotypic values. These constraints could make it difficult to fit this type of models using standard least square (LS) regression approach. To avoid this problem, under a regression model framework, Zeng et al. (2005) proposed a general two-allele (G2A) model for variance component analysis on bi-allelic QTL. The G2A model retains the nice property of the classical Fisher's model on orthogonal partition of the genotypic variance in equilibrium populations. Meanwhile, without using redundant parameters, the G2A model can be fitted directly using the standard LS approach. Wang and Zeng (2009) further explored the origin of the G2A model and clarified its theoretical basis on pertaining the orthogonal partition of the genotypic variance in equilibrium populations. Álvarez-Castro and Carlborg (2007) also proposed a different way for re-parameterization of the expected genotypic values, which the authors referred as the natural and orthogonal interactions (NOIA) model. It appears that the NOIA model first defines the genetic effects in terms of the expected genotypic values, and then derives the design matrix for the genetic effects via an inverse of the matrix. Though the matrix formulation is helpful for facilitating the transformation between different parameterizations, making inverse of the matrix to construct the design matrix is difficult to implement analytically. For multiple QTL, how to define various genetic effects in terms of the expected genotypic values could also be a challenge especially when locus-by-locus interactions are involved. Besides, the NOIA model assumes that the design matrix is of full rank, which makes it unsuitable for reduced re-parameterization of the genotypic values. Currently, no NOIA model has been proposed for reduced re-parameterization of the expected genotypic values, or for multiple QTL in the presence of LD.

In this study, we further extend the G2A model to QTL with multiple alleles and multiple loci. In bi-allelic case, only one additive effect and one dominance effect are needed at each locus, and the locus-by-locus interactions can be easily included for constructing a full re-parameterization of the genotypic values. For one QTL with multiple alleles, how to define the dominance effects for various allelic interactions is not straightforward especially when phases of its genotypes are unknown. The extension to multiple loci is also cumbersome by the much more complex structure of locus-by-locus interactions. How to present the model and define various genetic variance components are not trivial tasks. To construct one-locus general multi-allele (GMA) model, we overcome the phase problem by appropriately merging the paternal and maternal allelic effects and allelic interactions in the phase-known situation. Typically, with phase unknown genotypes at a locus, we may have to assume that the paternal and maternal alleles have the same frequencies and contribute the same genetic effects so that we could merge them without distinguishing their parental origins. With phase-known genotypes, we can further break down the additive variance component into paternal and maternal variance components. For multiple QTL with multiple alleles, we develop concise expressions for constructing multi-locus GMA models and defining various genetic components. Explicit formulas for calculating various genetic variance components in equilibrium population are also derived.

The structure of this manuscript is organized as the following. First, we consider one multi-allele QTL with phase known genotypes. Following the same strategy as adopted in Wang and Zeng (2009), we start by introducing indicator variables that describe the inheritance of paternal and maternal alleles. Then we make mean corrections on these indicator variables and build one-locus GMA model based on the mean-corrected index variables. Next, we consider one QTL with phase unknown genotypes. We construct a one-locus GMA model for unphased genotypes by appropriately combining the paternal and maternal allelic effects in the phase known one-locus GMA model. We derive formulas on partitioning the genotypic variance into the additive and dominance variance components under Hardy-Weinberg equilibrium (HWE) as well as in Hardy-Weinberg disequilibrium (HWD) for both the phase-known and phase-unknown one-locus GMA models. Thirdly, we extend the one-locus GMA models to multiple loci with either phase known or phase unknown genotypes. Based on these multi-locus GMA models, we describe how the various genetic components can be defined. An orthogonal partition on the genotypic variance in an equilibrium population is also presented in both the phase known and unknown cases. In addition, we discuss how to construct the reduced multi-locus GMA models in practice. The difference in using F_∞ models and GMA models to build reduced models for the expected genotypic values is explored. Finally, we apply the GMA model to a published experimental data set.

2. Methods and Results

The variation of a quantitative trait Y is usually assumed to be contributed by both genetic and environmental effects. Let G denote the true (unobservable) genotypic value from the joint contribution of all the genetic factors to the quantitative trait Y. Given genotypes at targeted QTL or genetic marker loci, we focus on assessing the variations (i.e., variances) contributed by the allelic effects and interactions of the QTL or marker loci to the total genotypic variance V_G.

2.1. One-locus Model for Phase Known

First, let us consider a single QTL with multiple alleles A₁, …, A_m (m ≥ 2). Suppose we know the parental origins of the alleles for observed genotypes. Then we can distinguish the paternal and maternal allelic effects separately. With m alleles, there are in total m² possible phased genotypes: (A_i, A_j), i, j = 1, …, m, where A_i and A_j represent the paternal and maternal allele, respectively. Let Gⁱ_j = E[G|g = (A_i, A_j)] denote the expected genotypic value for a phased genotype (A_i, A_j) in a study population. Then there are totally m² possible expected genotypic values Gⁱ_j, i, j = 1, …, m. To model these expected genotypic values, Fisher's classical one-locus model is given by

$\begin{array}{cc}{G}_j^i=\mu +{\alpha }^i+{\alpha }_j+{\delta }_j^i, i,j=1,\dots ,m,& (1)\end{array}$

where αⁱ (or α_j) is the so-called average additive or main allelic effect of a paternal (or maternal) allele A_i (or A_j), and δⁱ_j is the average allelic interaction between a paternal allele A_i and a maternal allele A_j. The above model is a typical two-way ANOVA model with the paternal and maternal gametes being treated as two independent risk factors. Although the paternal and maternal gametes often share (but do not have to) the same set of alleles A₁, …, A_m at the QTL, it allows the paternal and maternal gametes to have different allele frequencies and allelic effects at the QTL.

Let pⁱ be the frequency of allele A_i on paternal gametes $\left({\displaystyle {\sum }_{i=1}^m{p}^i}=1\right)$, and p_j be the frequency of allele A_j on the maternal gametes $\left({\displaystyle {\sum }_{j=1}^m{p}_j}=1\right)$. One nice feature of the Fisher model above is that it can assess the additive and dominance variance components V_A and V_D, which are defined as variations contributed by the additive allelic effects from the paternal and maternal alleles and the allelic interactions between the paternal or maternal alleles, respectively. For example, under HWE, it is well known that the total genotypic variance V_G = ∑_i,jpⁱp_j(Gⁱ_j − μ)² has an orthogonal partition V_G = V_A + V_D, where

$V_A={\displaystyle \sum _i{p}^i}{({\alpha }^i)}^2+{\displaystyle \sum _j{p}_j}{({\alpha }_j)}^2 , V_D={\displaystyle \sum _{i,j}{p}^i}{p}_j{({\delta }_j^i)}^2.$

Note that there are in total m² + 2m + 1 = (m + 1)² parameters being involved in Fisher model (1) including the intercept μ, which is more than the total number m² of the expected genotypic values Gⁱ_j, i, j = 1, …, m. As a result, not all the model parameters are estimable. To avoid this problem, some constraints on the model parameters are often required. It is usually assumed that all the genetic effects are averaged to zero over any index; i.e.,

$\begin{array}{cc}{\displaystyle \sum _i{p}^i}{\alpha }^i=0, {\displaystyle \sum _j{p}_j}{\alpha }_j=0, {\displaystyle \sum _i{p}^i}{\delta }_j^i=0, {\displaystyle \sum _j{p}_j}{\delta }_j^i=0.& (2)\end{array}$

With these constraints, Fisher (1918) showed that the least square estimates (LSE) of the model parameters are given by

$\begin{array}{l}\mu =E(G), {\alpha }_i=G_{.}^i-G_{.}^{ ·} , {\alpha }_j=G_j^{ ·}-G_{.}^{ ·} , {\delta }_j^i=G_j^i-G_{.}^i\\ -G_j^{ ·}+G_{.}^{ ·}\end{array}$

where G^·_. = E(G), Gⁱ_. = E[G|g = (A_i, −)] and G^·_j = E[G|g = (−, A_j)]. In simple cases, we could estimate the model parameters using these formulas. In general, however, those “irregular” constraints in (2) make it difficult to fit model (1) using the standard LS approach via commonly used software such as SAS (SAS Institute INC, Raleigh, NC), especially when we need to adjust for certain environmental covariates.

Here we propose a way to get rid of the redundant parameters. Let us first introduce the following indicator variables that describe the transmission inheritance of the paternal and maternal alleles.

$\begin{array}{l}{z}_{P_i}(g)=\{\begin{array}{l}1,\text{the}\text{paternal}\text{allele}\text{is}{A}_i\hfill \\ 0,\text{the}\text{paternal}\text{allele}\text{is}\text{not}{A}_i ,\hfill \end{array}\\ z_{M_j}(g)=\{\begin{array}{l}1,\text{the}\text{maternal}\text{allele}\text{is}{A}_j\hfill \\ 0,\text{ the}\text{maternal}\text{allele}\text{is}\text{not}{A}_j ,\hfill \end{array}\end{array}$

for i, j = 1, 2, …, m at the QTL. Next, using the same strategy as adopted in Wang and Zeng (2006), we further make mean corrections on these indicator variables z_Pi, z_Mj and define the following mean-corrected index variables

$\begin{array}{l}{x}_{P_i}(g)=z_{P_i}(g)-E[z_{P_i}(g)]\\ =\{\begin{array}{l}1-p^i,\text{the}\text{paternal}\text{allele}\text{is}{A}_i\hfill \\ -p^i,\text{the}\text{paternal}\text{allele}\text{is}\text{not}{A}_i ,\hfill \end{array}\\ x_{M_j}(g)=z_{M_j}(g)-E[z_{M_j}(g)]\\ =\{\begin{array}{l}1-p_j,\text{the}\text{paternal}\text{allele}\text{is}{A}_j\hfill \\ -p_j,\text{the}\text{paternal}\text{allele}\text{is}\text{not}{A}_j ,\hfill \end{array}\end{array}$

Then we can re-write the Fisher model (1) as

$\begin{array}{cc}\begin{array}{l}E(G|g)=\mu +{\displaystyle \sum _{i=1}^m{\alpha }^i}{x}_{P_i}(g)+{\displaystyle \sum _{j=1}^m{\alpha }_j}{x}_{M_j}(g)\\ +{\displaystyle \sum _{i,j}{\delta }_j^i}{x}_{P_i}(g)x_{M_j}(g),\end{array}& (3)\end{array}$

where E(G|g) denotes the expected genotypic value given a genotype g at the targeted QTL. As shown in Wang and Zeng (2006), the above model is simply a different representation of the Fisher model (1) with the same constraints being applied on the model parameters. However, as the HWD measurements can be quantified as covariances between those mean-corrected index variables x_Pi's and x_Mj's, this model expression can facilitate derivation on examining the genetic variance components under HWD.

Now, we further exclude the redundant parameters in model (3). For a diploid subject such as human being, his or her genotype at a locus on a pair of homologous chromosomes consists of two alleles with one from the father and the other one from the mother. Therefore, we always have ${\displaystyle {\sum }_{i=1}^m{z}_{P_i}}(g)\equiv 1$ and ${\displaystyle {\sum }_{j=1}^m{z}_{M_j}}(g)\equiv 1$. Thus, $x_{P_m}=-{\displaystyle {\sum }_{i=1}^{m-1}{x}_{P_i}}$ and $x_{M_m}=-{\displaystyle {\sum }_{j=1}^{m-1}{x}_{M_j}}$. So we can simply replace x_Pm by $(-{\displaystyle {\sum }_{i=1}^{m-1}{x}_{P_i}})$, and x_Mm by $(-{\displaystyle {\sum }_{j=1}^{m-1}{x}_{M_j}})$ in model (3), which leads to the following revised Fisher model

$\begin{array}{cc}\begin{array}{l}E(G|g)=\mu +{\displaystyle \sum _{i=1}^{m-1}{\alpha }^{\ast i}}{x}_{P_i}(g)+{\displaystyle \sum _{j=1}^{m-1}{\alpha }_{\ast j}}{x}_{M_j}(g)\\ +{\displaystyle \sum _{i=1}^{m-1}{\displaystyle \sum _{j=1}^{m-1}{\delta }_{\ast j}^{\ast i}}}{x}_{P_i}(g)x_{M_j}(g).\end{array}& (4)\end{array}$

Model (4) provides a full re-parameterization of the m² expected genotypic values Gⁱ_j, i, j = 1, …, m, without using redundant parameters. We refer it as a revised one-locus Fisher model or a general multi-allele (GMA) model. We also refer its model parameters α*ⁱ (or α_*j) as the average (additive) allelic effect of the paternal (or maternal) allele A_i (or A_j), and δ*ⁱ_*j the average allelic interaction between the paternal allele A_i and maternal allele A_j, with respect to the baseline allele A_m. Here, we choose allele A_m as the baseline allele but it could be any other allele as well. It is easy to show that, in terms of the parameters in the original Fisher model (1), we have α*ⁱ = αⁱ − α^m, α_*j = α_j − α_m and δ*ⁱ_*j = δⁱ_j − δ^m_j − δⁱ_m + δ^m_m, for i, j = 1, …, m.

Model (4) retains most of the nice features of the original Fisher model (1) on partition of the genotypic variance. Based on this model, we have the genetic additive variance components $V_{A_P}=\text{Var}({\displaystyle {\sum }_{i=1}^{m-1}{\alpha }^{\ast i}}{x}_{P_i})$ and $V_{A_M}=\text{Var}({\displaystyle {\sum }_{j=1}^{m-1}{\alpha }_{\ast j}}{x}_{M_j})$, which denote variations contributed by the additive effects of the paternal and maternal alleles, respectively. The genetic dominant variance component $V_D=\text{Var}({\displaystyle {\sum }_{i=1}^{m-1}{\displaystyle {\sum }_{i=1}^{m-1}{\delta }_{\ast j}^{\ast i}}}{x}_{P_i}{x}_{M_j})$, which accounts for a variation contributed by the interactions between paternal and maternal alleles in addition to the additive variance components. Under HWE, the paternal index variables {x_Pi, i = 1, …, m} are independent of the maternal index variables {x_Mj, j = 1, …, m}. Therefore, we have μ = E(G) and an orthogonal partition on the variance of the expected genotypic values: V(E(G|g)) = V_AP + V_AM + V_D, where

$\begin{array}{l}{V}_{A_P}={\displaystyle \sum _{i=1}^{m-1}{p}_i}{({\alpha }^{\ast i})}^2-{\left({\displaystyle \sum _{i=1}^{m-1}{p}^i}{\alpha }^{\ast i}\right)}^2,\\ V_{A_M}={\displaystyle \sum _{j=1}^{m-1}{p}_j}{({\alpha }_{\ast j})}^2-{\left({\displaystyle \sum _{j=1}^{m-1}{p}_j}{\alpha }_{\ast j}\right)}^2 \\ V_D={\displaystyle \sum _{i,j=1}^{m-1}{p}^i}{p}_j{({\delta }_{\ast j}^{\ast i})}^2-{\displaystyle \sum _{i=1}^{m-1}{p}^i}{\left({\displaystyle \sum _{j=1}^{m-1}{p}_j}{\delta }_{\ast j}^{\ast i}\right)}^2\\ -{\displaystyle \sum _{j=1}^{m-1}{p}_j}{\left({\displaystyle \sum _{i=1}^{m-1}{p}^i}{\delta }_{\ast j}^{\ast i}\right)}^2+{\left({\displaystyle \sum _{i,j=1}^{m-1}{p}^i}{p}_j{\delta }_{\ast j}^{\ast i}\right)}^2.\end{array}$

In HWD, the disequilibrium measurements can be captured by the covariances between the index variables x_Pi's and x_Mj's. In this case, we no longer have an orthogonal partition on the variance of the genotypic values for the additive and dominance variance components (see Appendix A). It should be pointed out that although the model parameters defined in model (4) depend upon the choice of the reference allele A_m, the additive and dominant variance components V_AP, V_AM, V_D as well as their covariance components do not depend on such a choice.

Note that the partition of the total genotypic variance V_G based on the Fisher model (1) assumes that the genotypic values are completely determined by the QTL. In general, beyond the targeted QTL, the genotypic value G could also be affected by other unobserved QTL. In this case, the total genotypic variance V_G = V(E(G|g)) + E(V(G|g)), where V(G|g) is a variation that cannot be explained by the targeted QTL. So the Fisher model (1), or its revised model (4), actually provides a partition on the variance V(E(G|g)) of the expected genotypic values at the targeted QTL instead of the total genotypic variance V_G. In practice, given a random sample from the study population, if we ignore the genetic by environmental interaction, a quantitative trait can typically be expressed in a linear regression model form as

$y_k=\beta z_k+E(G|g_k)+{ϵ}_k, k=1,\dots ,n,$

where, for k = 1, …, n, g_k is the observed QTL genotype of subject k, z_k is a vector of subject k from some adjusted environmental covariates with fixed effects β, and ϵ_k is a model residual contributed by other environmental and genetic factors that cannot be captured by z_k and g_k at the targeted QTL. Assuming that ϵ_k, k = 1, …, n, are independent and identically distributed (i.i.d) with a variance σ²_ϵ, we have V_y = V(E(G|g)) + σ²_ϵ. To further partition V(E(G|g)) based on the QTL genotypes, we can first calculate the allele frequencies and compute the values of the index variables {x_Pi(g_k), x_Mj(g_k)} for each subject k = 1, …, n, in the sample. Then, by treating these index variables as fixed covariates, we can incorporate the GMA model (4) into the regression model above and fit the model using the standard LS approach. Based on the LSE ${\widehat{\alpha }}^{\ast i}$, ${\widehat{\alpha }}_{\ast j}$ and ${\widehat{\delta }}_{\ast j}^{\ast i}$ of the model parameters, we can compute the additive and dominance genetic components $A_P(k)={\displaystyle {\sum }_{i=1}^{m-1}{\widehat{\alpha }}^{\ast i}}{x}_{P_i}(g_k),A_M(k)={\displaystyle {\sum }_{j=1}^{m-1}{\widehat{\alpha }}_{\ast j}}{x}_{M_j}(g_k)$ and $D(k)={\displaystyle {\sum }_{i=1}^{m-1}{\displaystyle {\sum }_{j=1}^{m-1}{\widehat{\delta }}_{\ast j}^{\ast i}}}{x}_{P_i}(g_k)x_{M_j}(g_k)$, for k = 1, …, n. Finally, the genetic variance components V_AP, V_AM and V_D can be estimated as the sample variances of {A_P(k), k = 1, …, n}, {A_M(k), k = 1, …, n}, and {D(k), k = 1, …, n}, respectively. Meanwhile, the genetic covariance components Cov(A_P, A_M), Cov(A_P, D), and Cov(A_M, D) can be estimated through the sample covariances

$\begin{array}{l}\widehat{\text{Cov}(A_P,A_M)}=\frac{1}{n}{\displaystyle \sum _{k=1}^n(A_P(}k)-{\overline{A}}_P)(A_M(k)-{\overline{A}}_M),\\ \widehat{\text{Cov}(A_P,D)}=\frac{1}{n}{\displaystyle \sum _{k=1}^n(A_P(}k)-{\overline{A}}_P)(D(k)-\overline{D}),\\ \widehat{\text{Cov}(A_M,D)}=\frac{1}{n}{\displaystyle \sum _{k=1}^n(A_M(}k)-{\overline{A}}_M)(D(k)-\overline{D}),\end{array}$

respectively, where ${\overline{A}}_P={\displaystyle {\sum }_{k=1}^n{A}_P}(k)/n,{\overline{A}}_M={\displaystyle {\sum }_{k=1}^n{A}_M}(k)/n$ and $\overline{D}={\displaystyle {\sum }_{k=1}^{n}D}(k)/n$.

2.2. One-locus Model for Phase Unknown

In this subsection, we consider modeling a multi-allele QTL with phase unknown genotypes—a more common situation in practice. As we cannot distinguish the parental origins of alleles in QTL genotypes, as usual, we assume that the paternal and maternal gametes share the same set of alleles with the same allele frequencies. Let A₁, …, A_m (m ≥ 2) denote the alleles at a target QTL or marker locus with allele frequencies p_i, i = 1, …, m. Ignoring the phases, there are m possible homozygous genotypes A_iA_i, i = 1 …, m, and m(m − 1)/2 possible heterozygous genotypes A_iA_j, i ≠ j. We also assume that these alleles contribute the same genetic effects regardless of their parental origins, which implies that the expected genotypic values G_ij = E(G|g = A_iA_j) satisfy the symmetric property: G_ij = G_ji, for i ≠ j. So there are totally m(m + 1)/2 possible distinctive expected genotypic values G_ij, i, j = 1, …, m. In this case, by treating the paternal and maternal gametes as two independent risk factors, the Fisher's ANOVA model for the expected genotypic values G_ij can be written as

$\begin{array}{cc}{G}_{ij}=\mu +{\alpha }_i+{\alpha }_j+{\delta }_{ij}, i,j=1,\dots ,m,& (5)\end{array}$

where α_i is the average (additive) allelic effect of the paternal or maternal allele A_i (i = 1, …, m), and δ_ij is the average allelic interaction between two alleles A_i and A_j (i, j, = 1, …, m). As pointed out in Wang (2011), the above model is different from the classical two-way ANOVA model in that the paternal and the maternal gametes share the same set of alleles and have the same allelic effects at the locus. To avoid the inestimability of model parameters in model (5) due to over-parameterization, the following constraints are often added on the model parameters

${\displaystyle \sum _i{p}_i}{\alpha }_i=0, {\displaystyle \sum _i{p}_i}{\delta }_{ij}=0.$

From the symmetric property of G_ij, we also assume δ_jk's are symmetric. Similarly, the above constraints together with the symmetry property of δ_jk make it difficult to fit model (5) using the standard LS approach.

Note that model (5) can be treated as a special case of model (1) or (3) with αⁱ = α_i and δⁱ_j = δ^j_i for i, j = 1, …, m. To construct a similar GMA model for model (5), we can combine the term α*ⁱx_Pi with α_*ix_Mi, and δ*ⁱ_*jx_Pix_Mj with δ*^j_*ix_Pjx_Mi in model (4). By denoting α*ⁱ = α_*i as α*_i for i = 1, …, m − 1, and δ*^j_*i = δ*ⁱ_*j as δ*_ij for i ≤ j, we obtain the following GMA model

$\begin{array}{cc}\begin{array}{l}E(G|g)=\mu +{\displaystyle \sum _{i=1}^{m-1}{\alpha }_i^{\ast }}{w}_i(g)+{\displaystyle \sum _{i=1}^{m-1}{\delta }_{ii}^{\ast }}{v}_{ii}(g)\\ +{\displaystyle \sum _{j=2}^{m-1}{\displaystyle \sum _{i<j}{\delta }_{ij}^{\ast }}}{v}_{ij}(g),\end{array}& (6)\end{array}$

where, for i = 1, …, m − 1,

$\begin{array}{l}{w}_i(g)=x_{P_i}+x_{M_i}=\{\begin{array}{l}2(1-p_i),\text{if}g=A_i{A}_i\hfill \\ 1-2p_i,\text{if}g=A_i{A}_i^c\hfill \\ -2p_i,\text{if}g=A_i^c{A}_i^c ,\hfill \end{array}\\ v_{ii}(g)=x_{P_i}{x}_{M_i}=\{\begin{array}{l}{(1-p_i)}^2,i\text{f}g=A_i{A}_i\hfill \\ -p_i(1-p_i),\text{ if}g=A_i{A}_i^c\hfill \\ p_i^2,\text{ if}g=A_i^c{A}_i^c ,\hfill \end{array}\end{array}$

and for i < j

$v_{ij}(g)=x_{P_i}{x}_{M_j}+x_{P_j}{x}_{M_i}=\{\begin{array}{l}(1-p_i)(1-p_j)+p_i{p}_j,\text{ if}g=A_i{A}_j\hfill \\ -2p_j(1-p_i),\text{ if}g=A_i{A}_i\hfill \\ -p_j(1-2p_i),\text{ if}g=A_i{A}_{ij}^c\hfill \\ -2p_i(1-p_j),\text{ if}g=A_j{A}_j\hfill \\ -p_i(1-2p_j),\text{ if}g=A_j{A}_{ij}^c\hfill \\ 2p_i{p}_j,\text{ if}g=A_{ij}^c{A}_{ij}^c .\hfill \end{array}$

Here A^c_ij denotes an allele which is different from A_i and A_j. Note that the combined index variables w_i(g), v_ii(g), and v_ij(g) above are well defined on unphased genotypes, although x_Pi, x_Mj are not. We refer them as genotype coding variables, and the model parameter α*_i as the average allelic effect of allele A_i (i = 1, …, m), and δ*_ij as the average allelic interaction between two alleles A_i and A_j (i, j, = 1, …, m), with respect to the baseline allele A_m. In terms of the parameters in the original model (5), we can show that α*_i = α_i − α_m, for i = 1, …, m − 1; and δ*_ij = δ_ij − δ_im − δ_jm + δ_mm, for i ≤ j, i, j = 1, …, m − 1.

Model (6) is an extension of the one-locus G2A model proposed in Zeng et al. (2005) to one QTL with multiple alleles. Note that no v_ij(g)'s for i < j are needed in the bi-allelic case. The combined index variables v_ii(g), i = 1, …, m − 1, are also slightly different from the ones defined in Zeng et al. (2005) by a scalar of (−2) folds. Here we drop the scalar (−2) so that the coefficient δ*_ii can keep the same interpretation as δ*ⁱ_*i in model (4). In addition, the combined index variables v_ij (i < j) defined above are not exactly the same as the ones we suggested in the discussion section of Wang (2011). Based on the previously defined inheritance indicator variables, we define w*_i = z_Pi + z_Mi and v*_ij(g) = z_Piz_Mj + z_Pjz_Mi for i < j. Then here we have v_ij = v*_ij − (p_iw*_j + p_jw*_i) + 2p_ip_j, for i < j.

Still, model (6) retains the nice feature of the classical Fisher's model on partition of the genotypic variance. The additive variance component $V_A=V\left({\displaystyle {\sum }_{i=1}^{m-1}{\alpha }_i^{\ast }}{w}_i(g)\right)$, which is contributed by the additive effects of both paternal and maternal alleles. The dominant variance component $V_D=V\left({\displaystyle {\sum }_{i=1}^{m-1}{\delta }_{ii}^{\ast }}{v}_{ii}(g)+{\displaystyle {\sum }_{i<j}{\delta }_{ij}^{\ast }}{v}_{ij}(g)\right)$, which is contributed by all the interactions between paternal and maternal alleles. Under HWE, we have μ = E(G) and an orthogonal partition on the variance of the expected genotypic values V(E(G|g)) = V_A + V_D, where

$\begin{array}{l}{V}_A=2{\displaystyle \sum _{i=1}^{m-1}{p}_i}{({\alpha }_i^{\ast })}^2-2{\left({\displaystyle \sum _{i=1}^{m-1}{p}_i}{\alpha }_i^{\ast }\right)}^2,\\ V_D={\displaystyle \sum _{i,j=1}^{m-1}{p}_i}{p}_j{({\delta }_{ij}^{\ast })}^2-2{\displaystyle \sum _{j=1}^{m-1}{p}_j}\text{}{\left({\displaystyle \sum _{i=1}^{m-1}{p}_i}{\delta }_{ij}^{\ast }\text{}\right)}^2+{\left({\displaystyle \sum _{i,j=1}^{m-1}{p}_i}{p}_j{\delta }_{ij}^{\ast }\right)}^2\text{}\text{}.\end{array}$

Here we define δ*_ji = δ*_ij, for i < j. In HWD, the desired orthogonal partition on the variance of the expected genotypic values is no longer held. But the model can still allow us to capture the disequilibria via covariances between those index variables x_Pi's and x_Mj's (see Appendix B).

As an example, let us consider a QTL with 3 alleles A₁, A₂, and A₃. By taking A₃ as the baseline allele, model (6) leads to

$\begin{array}{cc}\begin{array}{l}E(G|g)=\mu +{\alpha }_1^{\ast }{w}_1(g)+{\alpha }_2^{\ast }{w}_2(g)+{\delta }_{11}^{\ast }{v}_{11}(g)+{\delta }_{22}^{\ast }{v}_{22}(g)\\ +{\delta }_{12}^{\ast }{v}_{12}(g).\end{array}& (7)\end{array}$

Or, in a matrix form, we have

$\begin{array}{l}\text{}\left(\text{}\begin{array}{l}{G}_{11}\hfill \\ G_{12}\hfill \\ G_{22}\hfill \\ G_{13}\hfill \\ G_{23}\hfill \\ G_{33}\hfill \end{array}\text{}\right)\text{}=\left(\begin{array}{llllll}1\hfill & 2(1-p_1)\hfill & -2p_2\hfill & {(1-p_1)}^2\hfill & p_2^2\hfill & -2p_2(1-p_1)\hfill \\ 1\hfill & 1-2p_1\hfill & 1-2p_2\hfill & -p_1(1-p_1)\hfill & -p_2(1-p_2)\hfill & 1-p_1-p_2+2p_1{p}_2\hfill \\ 1\hfill & -2p_1\hfill & 2(1-p_2)\hfill & p_1^2\hfill & {(1-p_2)}^2\hfill & -2p_1(1-p_2)\hfill \\ 1\hfill & 1-2p_1\hfill & -2p_2\hfill & -p_1(1-p_1)\hfill & p_2^2\hfill & -p_2(1-2p_1)\hfill \\ 1\hfill & -2p_1\hfill & 1-2p_2\hfill & p_1^2\hfill & -p_2(1-p_2)\hfill & -p_1(1-2p_2)\hfill \\ 1\hfill & -2p_1\hfill & -2p_2\hfill & p_1^2\hfill & p_2^2\hfill & 2p_1{p}_2\hfill \end{array}\right)\\ \left(\begin{array}{r}\hfill \mu \\ \hfill {\alpha }_1^{\ast }\\ \hfill {\alpha }_2^{\ast }\\ \hfill {\delta }_{11}^{\ast }\\ \hfill {\delta }_{22}^{\ast }\\ \hfill {\delta }_{12}^{\ast }\end{array}\right)\end{array}$

If we choose A₁ (or A₂) instead of A₃ as the baseline allele, we can obtain different re-parameterizations of the six expected genotypic values. But they all give the same partition on the variance of the expected genotypic values. Álvarez-Castro and Yang (2011) also presented a similar re-parameterization of the six expected genotypic values from their NOIA model. It appears that their one-locus 3-allele NOIA model and the GMA model (7) above are likely equivalent on partitioning the variance of the six expected genotypic values, as we will see later in the Example Subsection 2.4.

Similar to the phase-known case, we can estimate the additive, dominance variance components and the covariance Cov(A, D) from a random sample when the QTL genotypes are available. First, we can compute the values of the combined index variables w_i(g_k) and v_ij(g_k) at the target QTL for each subject k = 1, …, n. Next, we can incorporate model (6) into a regression model with other possible adjusted covariates. By treating those combined index variables as regular fixed covariates, we can fit the regression model using the standard LS approach. Then, based on the fitted model, we compute the additive and dominance components $A(k)={\displaystyle {\sum }_{i=1}^{m-1}{\widehat{\alpha }}_i^{\ast }}{w}_i(g_k)$ and $D(k)={\displaystyle {\sum }_{i=1}^{m-1}{\widehat{\delta }}_{ii}^{\ast }}{v}_{ii}(g_k)+{\displaystyle {\sum }_{i<j}{\widehat{\delta }}_{ij}^{\ast }}{v}_{ij}(g_k)$ for each subject k = 1, …, n, where ${\widehat{\alpha }}_i^{\ast },{\widehat{\delta }}_{ii}^{\ast }$ and ${\widehat{\delta }}_{ij}^{\ast }$ are LSE of the model parameters. Finally, we can estimate V_A and V_D as the sample variances of {A(k), k = 1, …, n} and {D(k), k = 1, …, n}, respectively. Meanwhile, an estimate of Cov(A, D) is given by the sample covariance between {A(k), k = 1, …, n} and {D(k), k = 1, …, n}.

2.3. Multi-locus Models

The one-locus GMA models can be extended to multiple loci. Typically, for each locus k, we can create the mean-corrected index variables x^(k)_Pi, x^(k)_Mj (or combined index variables w_k,i, v_k,ij for phase-unknown genotypes) in the same way as in the one-locus case. Then we can build a multi-locus GMA model by including the within locus additive and dominance effects as well as the locus-by-locus interactions (i.e., epistases). The proposed multi-locus GMA model allows LD among multiple loci even though the LD may affect the partition on the variance of the expected genotypic values.

Consider L loci with alleles Ak1, …, Akmk at the k-th locus for k = 1, …, L. With phase known, there are totally m21··· m2L possible expected genotypic values: Gs1···sLt1···tL = E[G|g = (A1s1···ALsL, A1t1···ALtL)], where the joint genotype (A1s1···ALsL, A1t1···ALtL) is formed by the union of a paternal gamete A1s1···ALsL and a maternal gamete A1t1···ALtL with sk, tk ∈ {1, …, mk} for k = 1, …, L. Let pPksk and pMktk be the frequencies of the paternal allele Aksk and maternal allele Aktk, respectively. At each locus k = 1, …, L, we define the mean-corrected index variables x(k)Psk for paternal alleles Aksk (sk = 1, …, mk) and x(k)Mtk for maternal alleles Aktk (tk = 1, …, mk) in the same way as in the one-locus case. Then, when we choose A1m1, …, ALmL as the baseline alleles, a fully parameterized L-locus GMA model can be expressed as

$\begin{array}{cc}\begin{array}{l}\text{}\text{}\text{}E(G|g)={\displaystyle \sum _{i_1=0}^1{\displaystyle \sum _{j_1=0}^1\cdots }}{\displaystyle \sum _{i_L=0}^1{\displaystyle \sum }_{j_L=0}^1}{\displaystyle \sum _{\{s_1,\text{if} i_1=1\}}{\displaystyle \sum _{\{t_1,\text{if} j_1=1\}}\text{}}}\text{}\cdots {\displaystyle \sum _{\{s_L,\text{if} i_L=1\}}\text{}}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\\ {\displaystyle \sum _{\{t_L,\text{if} j_L=1\}}\text{}}\text{}{\alpha }_{t_1^{\ast }\cdots t_L^{\ast }}^{s_1^{\ast }\cdots s_L^{\ast }}·{\left(x_{P{s}_1}^{(1)}\right)}^{i_1}{\left(x_{M{t}_1}^{(1)}\right)}^{j_1}\cdots {\left(x_{P{s}_L}^{(L)}\right)}^{i_L}{\left(x_{M{t}_L}^{(L)}\right)}^{j_L}\text{}\text{}\text{} \end{array}& (8)\end{array}$

where the summation of sk (or tk) is from 1 up to (mk − 1) for k = 1, …, L; sk (or tk) refers to a particular paternal allele Aksk (or maternal allele Aktk) at the k-th locus; and ik (or jk) is an indicator variable for the presence or absence of a paternal (or maternal) allele at locus k which is involved in a term. The coefficient αs*1···s*Lt*1···t*L in each term represents an average allelic effect of a single paternal or maternal allele, or an allelic interaction from a set of paternal and maternal alleles that are involved in this term with respect to the baseline alleles A1m1, …, ALmL. The superscripts (or subscripts) in αs*1···s*Lt*1···t*L are defined as s*k = ik · sk (or t*k = jk · tk), indicating which paternal (or maternal) allele at locus k is involved in this term. If the term does not have any paternal (or maternal) allele at locus k being involved, then we have s*k = 0 (or t*k = 0). Note that $n{a}_P={\displaystyle {\sum }_{k=1}^L{i}_k}$ (or $n{a}_M={\displaystyle {\sum }_{k=1}^L{j}_k}$) specifies the total number of paternal (or maternal) alleles being involved in a term. We refer the total number of both paternal and maternal alleles that are involved in a term na = naP + naM as the order of this term.

The multi-locus GMA model (8) provides a full re-parameterization of the m21···m2L expected genotypic values Gs1···sLt1···tL with phase-known genotypes without using redundant parameters. Note that E(G|g) can also be partitioned into a sum of the following genetic components

$\begin{array}{l}{C}_{j_1\cdots j_L}^{i_1\cdots i_L}={\displaystyle \sum _{\{s_1,\text{if} i_1=1\}}}{\displaystyle \sum _{\{t_1,\text{if} j_1=1\}}}\mathrm{...}{\displaystyle \sum _{\{s_L,\text{if} i_L=1\}}}{\displaystyle \sum _{\{t_L,\text{if} j_L=1\}}}\\ {\alpha }_{t_1^{\ast }\cdots t_L^{\ast }}^{s_1^{\ast }\cdots s_L^{\ast }}·{\left(x_{P{s}_1}^{(1)}\right)}^{i_1}{\left(x_{M{t}_1}^{(1)}\right)}^{j_1}\cdots {\left(x_{P{s}_L}^{(L)}\right)}^{i_L}{\left(x_{M{t}_L}^{(L)}\right)}^{j_L},\end{array}$

where i₁, j₁, …, i_L, j_L = 0, 1 are indicator variables, which specify the parental origins of the contributing alleles. In other words, each genetic component consists of the terms with the same order and having their alleles all coming from the same subset of loci with the same parental origins but allowing varied allelic types. The classical genetic variance components can then be defined as variances of these genetic components. Note that the component C^0···0_0···0 is a constant, which corresponds to an intercept without any alleles being involved. Therefore, for a L-locus GMA model with phase known genotypes, there are in total (2^2L−1) genetic variance components $V\left(C_{j_1\cdots j_L}^{i_1\cdots i_L}\right)$, where i₁, j₁, …, i_L, j_L = 0, 1 but not all being zeros. Among them, there are L paternal and L maternal variance components of order 1 with each having a single paternal or maternal allele being involved, L dominance variance components of order 2 with each having two alleles within the same locus being involved, L(L − 1)/2 paternal by paternal (or maternal by maternal) variance components of order 2 with two paternal (or maternal) alleles coming from different loci, and L(L − 1) paternal by maternal variance components of order 2 with one paternal and one maternal allele coming from two different loci. For examples, the within-locus paternal, maternal and dominance variance components at a locus k can be written as

$\begin{array}{l}{V}_{A_P^{(k)}}=V\left({\displaystyle \sum _{s_k=1}^{m_k-1}{\alpha }_{0\cdots 0\cdots 0}^{0\cdots s_k\cdots 0}}·x_{P{s}_k}^{(k)}\right)=V\left(C_{0\cdots 0\cdots 0}^{0\cdots 1\cdots 0}\right),\\ V_{A_M^{(k)}}=V\left({\displaystyle \sum _{t_k=1}^{m_k-1}{\alpha }_{0\cdots t_k\cdots 0}^{0\cdots 0\cdots 0}}·x_{M{t}_k}^{(k)}\right)=V\left(C_{0\cdots 1\cdots 0}^{0\cdots 0\cdots 0}\right),\\ V_{D^{(k)}}=V\left({\displaystyle \sum _{s_k=1}^{m_k-1}{\displaystyle \sum _{t_k=1}^{m_k-1}{\alpha }_{0\cdots t_k\cdots 0}^{0\cdots s_k\cdots 0}}}·x_{P{s}_k}^{(k)}{x}_{M{t}_k}^{(k)}\right)=V\left(C_{0\cdots 1\cdots 0}^{0\cdots 1\cdots 0}\right),\end{array}$

respectively (k = 1, …, L). For a pair of loci j ≠ k, the two-locus paternal by paternal variance component is

$\begin{array}{l}{V}_{A_P^{(j)}\times A_P^{(k)}}=V\left({\displaystyle \sum _{s_j=1}^{m_j-1}{\displaystyle \sum _{s_k=1}^{m_k-1}{\alpha }_{0\cdots 0\cdots 0\cdots 0}^{0\cdots s_j\cdots s_k\cdots 0}}}·x_{P{s}_j}^{(j)}{x}_{P{s}_k}^{(k)}\right)\\ =V\left(C_{0\cdots 0\cdots 0\cdots 0}^{0\cdots 1\cdots 1\cdots 0}\right),\end{array}$

and the two-locus paternal by maternal variance component is

$\begin{array}{l}{V}_{A_P^{(j)}\times A_M^{(k)}}=V\left({\displaystyle \sum _{s_j=1}^{m_j-1}{\displaystyle \sum _{t_k=1}^{m_k-1}{\alpha }_{0\cdots 0\cdots t_k\cdots 0}^{0\cdots s_j\cdots 0\cdots 0}}}·x_{P{s}_j}^{(j)}{x}_{M{t}_k}^{(k)}\right)\\ =V\left(C_{0\cdots 0\cdots 1\cdots 0}^{0\cdots 1\cdots 0\cdots 0}\right).\end{array}$

The variance component of epistases with the highest order is

$\begin{array}{l}\text{}\text{}\text{}\text{}{V}_{D^{(1)}\times \cdots \times D^{(L)}}\\ =V\left({\displaystyle \sum _{s_1=1}^{m_1-1}{\displaystyle \sum _{t_1=1}^{m_1-1}\cdots }}{\displaystyle \sum _{s_L=1}^{m_L-1}{\displaystyle \sum _{t_L=1}^{m_L-1}{\alpha }_{t_1\cdots t_L}^{s_1\cdots s_L}}}{x}_{P{s}_1}^{(1)}{x}_{M{t}_1}^{(1)}\cdots x_{P{s}_L}^{(L)}{x}_{M{t}_L}^{(L)}\right)\\ =V(C_{1\cdots 1}^{1\cdots 1}),\end{array}$

which has 2L alleles being involved with two alleles per locus.

Based on these genetic components, we can partition the variance of the expected genotypic values into the variances and covariances of the genetic components. Still, the coefficients αs*1···s*Lt*1···t*L are defined based on the baseline alleles A1m1, …, ALmL. But the variances and covariances of these genetic components do not depend on the choice of the baseline alleles. Under HWE, the within-locus paternal index variables {x(k)Psk, sk = 1, …, mk} are independent of the maternal index variables {x(k)Mtk, tk = 1, …, mk} for each k = 1, …, L. With LE, the index variables {x(j)Psj, x(j)Mtj, sj, tj = 1, …, mj} at a locus j are also independent of the index variables {x(k)Psk, x(k)Mtk, sk, tk = 1, …, mk} at a different locus k (k ≠ j). To achieve an orthogonal partition on the variance of the expected genotypic values, we need to further assume that all the paternal index variables {x(k)Psk, sk = 1, …, mk, k = 1, …, L} are independent of the maternal index variables {x(k)Mtk, tk = 1, …, mk, k = 1, …, L} across all the loci; i.e., the so-called gametic equilibrium (see Wang and Zeng, 2006). Under both the gametic and linkage equilibria, each genetic component has its mean $E\left(C_{j_1\cdots j_L}^{i_1\cdots i_L}\right)=0$, and the covariances between different genetic components are zeros. Thus, we have E(G) = α0···00···0 and an orthogonal partition on the variance of the expected genotypic values is given by

$\begin{array}{cc}V(E(G|g))={\displaystyle \sum _{i_1=0}^1{\displaystyle \sum _{j_1=0}^1\cdots }}{\displaystyle \sum _{i_L=0}^1{\displaystyle \sum _{j_L=0}^{1}V}}\left(C_{j_1\cdots j_L}^{i_1\cdots i_L}\right),& (9)\end{array}$

where

$\begin{array}{l}\text{}\text{}\text{}\text{}V\left(C_{j_1\cdots j_L}^{i_1\cdots i_L}\right)={\displaystyle \sum _{\{s_1,{s^{\prime }}_1,\text{if} i_1=1\}}{(p_{1s_1}^P{1}_{\{s_1={s^{\prime }}_1\}}-p_{1s_1}^P{p}_{1{s^{\prime }}_1}^P)}^{i_1}}\\ {\displaystyle \sum _{\{t_1,{t^{\prime }}_1,\text{if} j_1=1\}}{(p_{1t_1}^M{1}_{\{t_1={t^{\prime }}_1\}}-p_{1t_1}^M{p}_{1{t^{\prime }}_1}^M)}^{j_1}}\cdots \\ {\displaystyle \sum _{\{s_L,{s^{\prime }}_L,\text{if} i_L=1\}}{(p_{L{s}_L}^P{1}_{\{s_L={s^{\prime }}_L\}}-p_{L{s}_L}^P{p}_{L{s^{\prime }}_L}^P)}^{i_L}}\\ {\displaystyle \sum _{\{t_L,{t^{\prime }}_L,\text{if} j_L=1\}}{(p_{L{t}_L}^M{1}_{\{t_L={t^{\prime }}_L\}}-p_{L{t}_L}^M{p}_{L{t^{\prime }}_L}^M)}^{j_L}}·{\alpha }_{t_1^{\ast }\cdots t_L^{\ast }}^{s_1^{\ast }\cdots s_L^{\ast }}{\alpha }_{{t^{\prime }}_1^{\ast }\cdots {t^{\prime }}_L^{\ast }}^{{s^{\prime }}_1^{\ast }\cdots {s^{\prime }}_L^{\ast }}.\end{array}$

Similarly, we can construct multi-locus GMA models for QTL with phase unknown genotypes. Without distinguishing the parental origin of the alleles, there are totally ∏^L_{k = 1} m_k(m_k + 1)/2^L possible expected genotypic values: G_{s1t1···sLtL} = E[G(g)|g = (A_1s1A_1t1, …, A_LsLA_LtL)], for s_k, t_k = 1, …, m_k and k = 1, …, L. We assume that the paternal and maternal allele frequencies are the same and denoted by p_ksk. We define the combined index variables w_{k, i}(g) and v_k,ij at each locus k in the same way as the one-locus case for k = 1, …, L. Then, when we choose A_1m1, …, A_LmL as the baseline alleles, a fully parameterized L-locus GMA model for QTL with phase unknown genotypes can be expressed as

$\begin{array}{cc}\begin{array}{l}E(G|g)={\displaystyle \sum _{i_1=0}^2\cdots }{\displaystyle \sum _{i_L=0}^2{\displaystyle \sum _{\{s_1,\text{if} i_1=1\}}\cdots }}\\ {\displaystyle \sum _{\{s_L,\text{if}{i}_L=1\}}}{\displaystyle \sum _{\{s_1\le t_1,\text{if} i_1=2\}}}\mathrm{...}{\displaystyle \sum _{\{s_L\le t_L,\text{if} i_L=2\}}}\\ {\alpha }_{s_1^{\ast }\cdots s_L^{\ast }}·\left(w_{1,s_1}^{1_{\{i_1=1\}}}·v_{1,s_1{t}_1}^{1_{\{i_1=2\}}}\right)\cdots \left(w_{L,s_L}^{1_{\{i_L=1\}}}·v_{L,s_L{t}_L}^{1_{\{i_L=2\}}}\right),\end{array}& (10)\end{array}$

where 1_{ik=j} is the Kronecker function which equals 1 when i_k = j and 0 otherwise, for k = 1, …, L and j = 1 or 2; the summation of s_k (or t_k) is from 1 up to m_k − 1 with s_k (or t_k) referring to allele A_ksk (or A_ktk) at a locus k (k = 1, …, L); and i_k specifies how many alleles at a locus k are involved in this term. If i_k = 1, this term has only one allele A_ksk (either paternal or maternal) being involved via w_k,sk and we set s_k* = s_k. If i_k = 2, both the paternal and maternal alleles A_ksk and A_ktk are involved via v_k,sktk in this term and we set s*_k = s_kt_k regardless of the order of s_k and t_k. When i_k = 0, this term does not have any alleles at locus k being involved and we have $\left(w_{k,s_k}^{1_{\{i_k=1\}}}·v_{k,s_k{t}_k}^{1_{\{i_k=2\}}}\right)=1$. The coefficient α_s*1···s*L represents the average allelic effect of a single allele, or an allelic interaction from all the alleles that are involved in this term, with respect to the baseline alleles A_1m1, …, A_LmL. We still refer the total number of alleles being involved in each term; i.e., $na={\displaystyle {\sum }_{k=1}^L{i}_k}$, as the order of this term.

Based on the above model, we can define the genetic components as

$\begin{array}{l}{C}_{i_1\cdots i_L}={\displaystyle \sum _{\{s_1,\text{if} i_1=1\}}}\mathrm{...}{\displaystyle \sum _{\{s_L,\text{if} i_L=1\}}}{\displaystyle \sum _{\{s_1\le t_1,\text{if} i_1=2\}}}\mathrm{...}{\displaystyle \sum _{\{s_L\le t_L,\text{if} i_L=2\}}}\\ {\alpha }_{s_1^{\ast }\cdots s_L^{\ast }}·\left(w_{1,s_1}^{1_{\{i_1=1\}}}·v_{1,s_1{t}_1}^{1_{\{i_1=2\}}}\right)\cdots \left(w_{L,s_L}^{1_{\{i_L=1\}}}·v_{L,s_L{t}_L}^{1_{\{i_L=2\}}}\right),\end{array}$

for i₁, …, i_L = 0, 1, 2, with each component being a summation of the terms with their alleles all coming from the same subset of loci and having the same number of alleles from each locus. Excluding the one with i₁ = ··· = i_L = 0, which corresponds to an intercept, there are in total 3^L − 1 genetic variance components V(C_i1···iL) for i₁, …, i_L = 0, 1, 2. Among them, there are L additive variance components with each having a single paternal or maternal allele being involved, L dominance variance components with each having both paternal and maternal alleles within the same locus being involved, L(L − 1)/2 additive by additive variance components with two alleles coming from different loci, etc. For examples, the within-locus additive and dominance variance components at a locus k can be written as

$\begin{array}{l}{V}_{A^{(k)}}=V\left({\displaystyle \sum _{s_k=1}^{m_k-1}{\alpha }_{0\cdots s_k\cdots 0}}·w_{k,s_k}\right)=V(C_{0\cdots 1\cdots 0}),\\ V_{D^{(k)}}=V\left({\displaystyle \sum _{s_k\le t_k}{\alpha }_{0\cdots s_k{t}_k\cdots 0}}·v_{k,s_k{t}_k}\right)=V(C_{0\cdots 2\cdots 0}),\end{array}$

respectively, for k = 1, …, L. For a pair of loci j ≠ k, the two-locus additive by additive interaction is

$\begin{array}{l}{V}_{A^{(j)}\times A^{(k)}}=V\left({\displaystyle \sum _{s_j=1}^{m_j-1}{\displaystyle \sum _{s_k=1}^{m_k-1}{\alpha }_{0\cdots s_j\cdots s_k\cdots 0}}}·w_{j,s_j}{w}_{k,s_k}\right)\\ =V(C_{0\cdots 1\cdots 1\cdots 0}).\end{array}$

The variance component of epistases with the highest order of 2L is given by

$\begin{array}{l}\text{}\text{}\text{}\text{}{V}_{D^{(1)}\times \cdots \times D^{(L)}}\\ =V\left({\displaystyle \sum _{s_1\le t_1}\cdots }{\displaystyle \sum _{s_L\le t_L}{\alpha }_{s_1{t}_1,\cdots ,s_L{t}_L}}{v}_{1,s_1{t}_1}\cdots v_{L,s_L{t}_L}\right)\\ =V(C_{2\cdots 2}).\end{array}$

Under both the gametic and linkage equilibria, we have E(G) = α_0···0 and an orthogonal partition on the variance of the expected genotypic values

$\begin{array}{cc}V(E(G|g))={\displaystyle \sum _{i_1=0}^2\cdots }{\displaystyle \sum _{i_L=0}^{2}V}(C_{i_1\cdots i_L}),& (11)\end{array}$

where

$\begin{array}{l}\text{}\text{}\text{}\text{}V(C_{i_1\cdots i_L})={\displaystyle \sum _{\{s_1,{s^{\prime }}_1,\text{if} i_1=1\}}{\left\{2\left(p_{1s_1}{1}_{\{s_1={s^{\prime }}_1\}}-p_{1s_1}{p}_{1s_1}\right)\right\}}^{1_{\{i_1=1\}}}}\cdots \\ {\displaystyle \sum _{\{s_L,{s^{\prime }}_L,\text{if} i_L=1\}}{\left\{2\left(p_{L{s}_L}{1}_{\{s_L={s^{\prime }}_L\}}-p_{L{s}_L}{p}_{L{s^{\prime }}_L}\right)\right\}}^{1_{\{i_L=1\}}}}\\ {\displaystyle \sum _{\{s_1,{s^{\prime }}_1,t_1,{t^{\prime }}_1,\text{if} i_1=2\}}\{\left(p_{1s_1}{1}_{\{s_1={s^{\prime }}_1\}}-p_{1s_1}{p}_{1{s^{\prime }}_1}\right)}\\ {\left(p_{1t_1}{1}_{\{t_1={t^{\prime }}_1\}}-p_{1t_1}{p}_{1{t^{\prime }}_1}\right)\}}^{1_{\{i_1=2\}}}\cdots \\ {\displaystyle \sum _{\{s_L,{s^{\prime }}_L,t_L,{t^{\prime }}_L,\text{if} i_L=2\}}\{\left(p_{L{s}_L}{1}_{\{s_L={s^{\prime }}_L\}}-p_{L{s}_L}{p}_{L{s^{\prime }}_L}\right)}\\ {\left(p_{L{t}_L}{1}_{\{t_L={t^{\prime }}_L\}}-p_{L{t}_L}{p}_{L{t^{\prime }}_L}\right)\}}^{1_{\{i_L=2\}}}{\alpha }_{s_1^{\ast }\cdots s_L^{\ast }}{\alpha }_{s_1^{\ast }{}^{\prime }\cdots s_L^{\ast }{}^{\prime }}.\end{array}$

Here, when i_k = 2, we have s*_k = s_kt_k and s*′_k = s′_kt′_k. We also define α_s*1···s*L (or α_{s*′1···s*′L}) to be the same if we switch the order of s_k and t_k in s*_k (or s_k' and t_k' in s*′_k) for k = 1, …, L. In the presence of disequilibria, the desired orthogonal partition may no longer hold. However, regardless of the equilibrium, the coefficients α_s*1···s*L in model (10) are defined based on the baseline alleles A_1m1, …, A_LmL, while the variances and covariances of the genetic components do not depend on the choice of these baseline alleles.

In practice, we do not have to rely on the derived formula to estimate the genetic variance or covariance components. Similar to the one-locus case, given the observed QTL genotypes for a random sample from a study population, we can always incorporate model (8) or (10) into a regression model with other possible adjusted covariates and fit the model using standard LS approach. Then we can estimate various genetic variance components as well as the covariances among different genetic components based on the fitted model. A good fit of a fully parameterized GMA model often requires that the expected genotypic values for all possible joint genotypes of the QTL are estimable from the study sample. If certain genotypes are not observable or rarely present in subjects from the study sample, a situation which likely happens when the number of alleles or the number of QTL is large with moderate or small sample size, the design matrix for the genetic effects could become singular which implies that some genetic variance components cannot be estimated reliably. But we do not have to use fully parameterized GMA models to model the expected genotypic values. In this case, we may want to build a reduced GMA model that can provide a good approximation to the expected genotypic values overall and meanwhile has a less complicated model structure. The fact that two terms within the same genetic component are unavoidably correlated suggests that we should perhaps treat each genetic component as a whole and keep or drop its terms all at once in building a GMA model. As genetic components of lower orders tend to have bigger impact on the expected genotypic values than the higher order ones, one way to construct a reduced GMA model is perhaps to go through a stepwise forward selection procedure by hierarchically adding the lowest order genetic component that can achieve a nominal significance level (e.g., 5%) but has not yet been selected in the model into the model one at a time. Here, the classical likelihood ratio statistic can be used to assess each genetic component for entering into or dropping from the model.

It has been known that the model building procedure is often sensitive to potential confounding among the selected variables. A GMA model uses the mean-corrected index variables x^(k)_Psk, x^(k)_Mtk (or their combined index variables w_k,sk, v_k,sktk in the phase unknown case) as the basic units in constructing all its model terms. At least in an equilibrium population, these mean-corrected index variables can reduce the confounding among different genetic components and make them uncorrelated. This orthogonality implies that each genetic component can be assessed independently regardless of which components are presented in the model. On the other hand, an F_∞ model can be thought of directly using the inheritance indicator variables z^(k)_Psk and z^(k)_Mtk (or their merged genotype coding variables w*_k,sk(g) = z^(k)_Psk + z^(k)_Msk and v*_k,sktk = z^(k)_Pskz^(k)_Mtk + z^(k)_Ptkz^(k)_Msk in phase unknown case) as the basic units in constructing its terms (Wang and Zeng, 2009; Wang, 2011). Even in an equilibrium population, an F_∞ model could have its low-order terms being highly confounded with other high-order terms when they contain shared alleles (see Zeng et al., 2005). As the result, in building a predictive model based on F_∞ models, the stepwise forward selection procedure could be problematic because failing to include a significant higher order term (or component) in a reduced F_∞ model could make the assessment of some low-order terms (or components) unreliable. On selecting significant QTL from a set of loci without having the locus-by-locus interactions being involved, the choice of using GMA or F_∞ model in building a reduced model for the expected genotypic values should not matter much because using mean-corrected index variables mainly affects the intercept in this case. However, when we consider including epistases for a given set of QTL, the GMA model can appropriately use the orthogonal property among different genetic components to dissect the confounding at least in equilibrium populations, while it appears that the F_∞ models cannot make full use of the equilibrium information. When disequilibria are present, as Hardy-Weinberg equilibria are expected to be held in most of the human genomic regions and LD mainly present for closely linked loci, we would expect that most of the genetic components in a GMA model are likely uncorrelated. Therefore, in most cases, using the GMA model could still be preferable to using F_∞ model in building reduced models for expected genotypic values especially when epistases are involved.

2.4. Example

As an example, we apply the GMA model to a published experimental data set on the polymorphism at the human acid phosphatase locus (ACP1). The analysis of this data set was first conducted by Greene et al. (2000), and recently re-analyzed as an example in Álvarez-Castro and Yang (2011). The ACP1 gene involves 3 alleles A, B, C. Two phenotypic traits are considered, which measure the ACP1 enzyme activity (y^ac) and inhibition (yⁱⁿ). The estimates of the expected genotypic values and the genotype frequencies are summarized in Table 1.

TABLE 1

Table 1. Expected genotypic values and observed genotype frequencies.

From the genotype frequencies, we first estimate the allele frequencies as p_A = 0.3534, p_B = 0.5818, and p_C = 0.0647. Then, for each trait, we fit a separate GMA model (7) to its expected genotypic values by taking C as the baseline allele and creating the combined index variables w₁(g), w₂(g), v₁₁(g), v₂₂(g), v₁₂(g). For ACP1 enzyme activity, we obtain LSE of the model parameters as μ = 167.735, α*₁ = −59.260, α*₂ = −26.254, δ*₁₁ = −4.800, δ*₂₂ = 3.700 and δ*₁₂ = −2.000. For ACP1 enzyme inhibition, we have μ = 39.386, α*₁ = −16.149, α*₂ = −19.714, δ*₁₁ = −0.200, δ*₂₂ = 4.200, and δ*₁₂ = 2.100.

Next, for each trait separately, we calculate A(g_k) = α*₁w₁(g_k) + α*₂w₂(g_k) and D(g_k) = δ*₁₁v₁₁(g_k) + δ*₂₂v₂₂(g_k) + δ*₁₂v₁₂(g_k), k = 1, …, 6, for the six ACP1 genotypes. Based on the genotype frequencies, we then calculate the genetic variance components. For ACP1 enzyme activity, we obtain V_A = 658.868, V_D = 0.973 and the covariance Cov(A, D) = −0.356. The total variance of the expected genotypic values is V(E(G^ac|g)) = 659.129, and V_A/V(E(G^ac|g)) = 99.96%. For ACP1 enzyme inhibition, we have V_A = 44.920, V_D = 0.146 and the covariance Cov(A, D) = −0.217. The total variance of the expected genotypic values is V(E(Gⁱⁿ|g)) = 44.632, and V_A/V(E(Gⁱⁿ|g)) = 100.65%. Note that the partition V(E(G|g)) = V_A + V_D + 2Cov(A, D), which is not orthogonal for both traits due to a slight deviation from HWE at ACP1 locus as we have observed from the genotype frequencies in Table 1. It is also interesting to see that for ACP1 enzyme inhibition, V_A/V(E(Gⁱⁿ|g)) is bigger than 100% due to HWD and the fact that V_D + 2Cov(A, D) < 0.

Using the same allele frequencies but assuming HWE, we would have slightly different genotype frequencies. Note that the LSE of the model parameters keep the same and do not depend on the genotype frequencies because they are completely determined by the allele frequencies and the six expected genotypic values when a fully-parameterized GMA is used. But the total variance of the expected genotypic values and its variance components will be different. For ACP1 enzyme activity, we obtain V_A = 660.588, V_D = 0.971, Cov(A, D) = 0.006, V(E(G^ac|g)) = 661.573 and V_A/V(E(G^ac|g)) = 99.85%. For ACP1 enzyme inhibition, we have V_A = 46.422, V_D = 0.152, Cov(A, D) = 0.001, V(E(Gⁱⁿ|g)) = 46.573 and V_A/V(E(Gⁱⁿ|g)) = 99.68%. The minor deviations from orthogonal partitions are likely due to some numerical round-off and the fact that the summation of the three allele frequencies is 0.9999 instead of exactly 1.

For this three-allele example, we also applied the NOIA model using the formulas (10) and (11) provided in Álvarez-Castro and Yang (2011), and in both HWE and HWD cases we obtained exactly the same results as above based on our own calculation, which however appear to be slightly different from what had been reported in Álvarez-Castro and Yang (2011). For example, based on the observed genotype frequencies, an estimate of V_D = 1.15 for the trait of ACP1 enzyme inhibition was reported in Álvarez-Castro and Yang (2011), which is noticeably larger than V_D = 0.146 that we obtain above.

3. Discussion

In the analysis of genetic variance components, a separation of the variations contributed by the additive allelic effects and allelic interactions is complicated by the fact that the observed genotypes are often phase-unknown. In this study, by appropriately merging the paternal and maternal allelic effects and allelic interactions in the phase-known situation, we propose a way to construct one-locus and multi-locus GMA models on analysis of genetic variance components for QTL with multiple alleles. In the same way as building a G2A model, we construct a GMA model by first specifying its design matrix for the genetic effects via some mean-corrected index variables. As these mean-corrected index variables are well defined based on the observed genotypes and allele frequencies, they can be treated as regular covariates for coding QTL genotypes. These one-locus or multi-locus GMA models can then be incorporated into standard regression models with other possible adjusted covariates and fitted using standard LS approach. Based on the fitted models, we can further estimate the genetic variance and covariance components through the sample variances and covariances of various genetic components. As we have pointed out, these GMA models can be applied to equilibrium populations as well as populations in Hardy-Weinberg and/or linkage disequilibria. By using the full set or a low-order subset of the index variables (or genetic components), the GMA model allows us to make either full or reduced re-parameterization of the genotypic values. When some loci have phase known genotypes while other loci have phase unknown genotypes (a possible hypothetical situation), a mixed GMA model could also be constructed by adopting the same modeling strategy.

Sometimes we may want to perform hypothesis tests on the existence of certain genetic variance and covariance components. Note that the GMA models have allele frequencies being involved in their design matrices. As allele frequencies often need to be estimated from the genotype data, they could contribute another source of variation in the LSE ($\widehat{\beta }$) of the model parameters as well as the genetic variance components. When the allele frequencies can be accurately estimated, we could simply treat them as fixed constants. When the residuals in a regression model do not depend on the allele frequency estimates, based on the linear model theory (see Ravishanker and Dey, 2002), $\widehat{\beta }$ are known to be unbiased with its covariance matrix cov($\widehat{\beta }$) = σ²E[(X′X)⁻¹], where X is the design matrix of a GMA model. In this case, we can assess the existence of variance components by performing traditional hypothesis tests on $\widehat{\beta }$. In general, we could also assess the existence of these genetic variance and covariance components through a bootstrap procedure. By repeatedly drawing random samples of the same size from the observed random sample with replacement, we can estimate the genetic variance and covariance components for each bootstrap sample and meanwhile assess the variances in estimates of these genetic variance and covariance components and test for their existence.

In genetic studies, QTL with missing genotypes is a common phenomenon. GMA model can be used to fit QTL with missing genotypes. Rather than excluding patients with missing QTL genotypes, we could treat “missing” as an allele although this strategy may induce potential bias as we assume that all the missing alleles have the same genetic effect. GMA models could also be applied incorporation with various imputation methods. In recent years, there has been a great deal of interest in developing methodologies for QTL mapping using recombinant intercrosses from multiple inbred lines. In this case, the putative QTL often have their locations and genotypes unknown. But the allele frequencies of QTL could probably be inferred from the study design and the QTL genotypes might be imputable from their neighboring genetic markers. How to apply GMA models to this type of experimental crosses for QTL mapping could be a research topic for further exploration.

In summary, the analysis of genetic variance components for multi-allele QTL has been challenging due to complex allelic interactions and locus-by-locus interactions. In this study, we thoroughly explored the architecture of one-locus and multi-locus GMA models with either phase known or unknown genotypes. Particularly, we described in detail the architecture of the multi-locus GMA model, and how the model terms can be grouped into various genetic components. Under equilibria populations, we also derived formulas for orthogonal partition of the genetic variance components, which could be useful for analytical assessment of the variance components. Comparing to the classical Fisher model, the GMA models can estimate the genetic variance and covariance components more conveniently via standard LS approach for either one or multiple QTL with multiple alleles, in equilibrium as well as disequilibrium populations.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The author would like to acknowledge Dr. Zhao-Bang Zeng at Bioinformatics Research Center, North Carolina State University, for his thoughtful comments and suggestions on an earlier version of the manuscript.

References

Álvarez-Castro, J. M. (2012). Current applications of models of genetic effects with interactions across the genome. Curr. Genomics 13, 163–175. doi: 10.2174/138920212799860689

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Álvarez-Castro, J. M., and Carlborg, Ö. (2007). A unified model for functional and statistical epistasis and its application in quantitative trait loci analysis. Genetics 176, 1151–1167. doi: 10.1534/genetics.106.067348

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Álvarez-Castro, J. M., and Yang, R.-C. (2011). Multiallelic models of genetic effects and variance decomposition in non-equilibrium populations. Genetica 139, 1119–1134. doi: 10.1007/s10709-011-9614-9

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheverud, J. M. (2000). “Detecting epistasis among quantitative trait loci,” in Epistasis and the Evolutionary Process, eds J. B. Wolf, E. D. Brodie, and M. J. Wade (New York, NY: Oxford University Press), 58–81.

Cockerham, C. C. (1954). An extension of the concept of partitioning hereditary variance for analysis of covariances among relatives when epistasis is present. Genetics 39, 859–882.

Pubmed Abstract | Pubmed Full Text

Cockerham, C. C. (1963). “Estimation of genetic variances,” in Statistical Genetics and Plant Breeding, Vol. 982, eds W. D. Henson and H. F. Robinson (Washington, DC: National Academy of Science - National Research Council), 53–94.

Falconer, D. S., and Mackay, T. F. C. (1996). Introduction to Quantitative Genetics, 4th Edn. Harlow: Longman Group Ltd.

Fisher, R. A. (1918). The correlation between relatives on the supposition of mendelian inheritance. Trans. R. Soc. Edinburgh 52, 399–433. doi: 10.1017/S0080456800012163

CrossRef Full Text

Greene, L. S., Bottini, N., Borgiani, P., and Gloria-Bottini, F. (2000). Acid phosphatase locus 1 (acp1): possible relationship of allelic variation to body size and human population adaptation to thermal stressa theoretical perspective. Am. J. Hum. Biol. 12, 688–701. doi: 10.1002/1520-6300(200009/10)12:5<688::AID-AJHB14>3.0.CO;2-C

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hansen, T. F., and Wagner, G. P. (2001). Modeling genetic architecture: a multilinear theory of gene interaction. Theor. Popul. Biol. 59, 61–86. doi: 10.1006/tpbi.2000.1508

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kempthorne, O. (1954). The correlation between relatives in a random mating population. Proc. R. Soc. Lond. B Biol. Sci. 143, 103–113. doi: 10.1098/rspb.1954.0056

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kempthorne, O. (1957). An Introduction to Genetic Statistics. New York, NY: John Wiley and Sons, Inc.

Kempthorne, O. (1969). An Introduction to Genetic Statistics. New Haven, CT: Iowa State University Press Ames.

Ravishanker, N., and Dey, D. K. (2002). A First Course in Linear Model Theory. Boca Raton, FL: Chapman and Hall, CRC.

Wang, T. (2011). On coding genotypes for genetic markers with multiple alleles in genetic association study of quantitative traits. BMC Genet. 12:82. doi: 10.1186/1471-2156-12-82

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wang, T., and Zeng, Z. B. (2006). Models and partition of variance for quantitative trait loci with epistasis and linkage disequilibrium. BMC Genet. 7:9. doi: 10.1186/1471-2156-7-9

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wang, T., and Zeng, Z. B. (2009). Contribution of genetic effects to genetic variance components with epistasis and linkage disequilibrium. BMC Genet. 10:52. doi: 10.1186/1471-2156-10-52

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Weir, B. S. (1996). Genetic Data Analysis II: Methods for Discrete Population Genetic Data. Sunderland, MA: Sinauer Associates, Inc.

Weir, B. S., and Cockerham, C. C. (1977). “Two-locus theory in quantitative genetics,” in Proceedings of the International Conference on Quantitative Genetics, eds E. Pollak, O. Kempthorne, and T. B. Bailey Jr. (Ames, IA: Iowa State University Press), 247–269.

Yang, R.-C., and Álvarez-Castro, J. M. (2008). Functional and statistical genetic effects with multiple alleles. Curr. Top. Genet. 3, 49–62.

Zeng, Z.-B., Wang, T., and Zou, W. (2005). Modeling quantitative trait loci and interpretation of models. Genetics 169, 1711–1725. doi: 10.1534/genetics.104.035857

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Appendix A. One-locus Model for Phase Known – in Hardy-Weinberg Disequilibrium

In HWD, we can represent the genotype frequencies as P_(Ai,Aj) = pⁱp_j + Dⁱ_j, where Dⁱ_j measures the departure from HWE with constraints ${\displaystyle {\sum }_{i=1}^m{D}_j^i}=0$ for j = 1, …, m, and ${\displaystyle {\sum }_{j=1}^m{D}_j^i}=0$ for i = 1, …, m. Note that Dⁱ_j = P_(Ai,Aj) − pⁱp_j = Cov(x_Pi, x_Mj), we can show that $E(G)=\mu +{\displaystyle {\sum }_{i=1}^{m-1}{\displaystyle {\sum }_{j=1}^{m-1}{\delta }_{\ast j}^{\ast i}}}{D}_j^i$, and the variance of the expected genotypic values becomes

$V(E(G|g))=V_{A_P}+V_{A_M}+V_D+2\text{Cov}(A_P,A_M)+2\text{Cov}(A_P,D)+2\text{Cov}(A_M,D)$

where $A_P={\displaystyle {\sum }_{i=1}^{m-1}{\alpha }^{\ast i}}{x}_{P_i}(g),A_M={\displaystyle {\sum }_{j=1}^{m-1}{\alpha }_{\ast j}}{x}_{M_j}(g)$ and $D={\displaystyle {\sum }_{i=1}^{m-1}{\displaystyle {\sum }_{j=1}^{m-1}{\delta }_{\ast j}^{\ast i}}}{x}_{P_i}(g)x_{M_j}(g)$. The formulas on calculating the additive variance components V_AP and V_AM are the same as in HWE case. But the formula for dominance variance component V_D becomes

$\begin{array}{l}{V}_D={\displaystyle \sum _{i,j=1}^{m-1}{p}^i}{p}_j{({\delta }_{\ast j}^{\ast i})}^2-{\displaystyle \sum _{i=1}^{m-1}{p}^i}{\left({\displaystyle \sum _{j=1}^{m-1}{p}_j}{\delta }_{\ast j}^{\ast i}\right)}^2-{\displaystyle \sum _{j=1}^{m-1}{p}_j}{\left({\displaystyle \sum _{i=1}^{m-1}{p}^i}{\delta }_{\ast j}^{\ast i}\right)}^2+{\left({\displaystyle \sum _{i,j=1}^{m-1}{p}^i}{p}_j{\delta }_{\ast j}^{\ast i}\right)}^2\\ +{\displaystyle \sum _{i,j=1}^{m-1}{({\delta }_{\ast j}^{\ast i})}^2}{D}_{ij}-2{\displaystyle \sum _{i=1}^{m-1}\left({\displaystyle \sum _{j=1}^{m-1}{p}_j}{\delta }_{\ast j}^{\ast i}\right)}\left({\displaystyle \sum _{j=1}^{m-1}{D}_j^i}{\delta }_{\ast j}^{\ast i}\right)-2{\displaystyle \sum _{j=1}^{m-1}\left({\displaystyle \sum _{i=1}^{m-1}{p}^i}{\delta }_{\ast j}^{\ast i}\right)}\left({\displaystyle \sum _{i=1}^{m-1}{D}_j^i}{\delta }_{\ast j}^{\ast i}\right)\\ +2\left({\displaystyle \sum _{i,j=1}^{m-1}{p}^i}{p}_j{\delta }_{\ast j}^{\ast i}\right)\left({\displaystyle \sum _{i,j=1}^{m-1}{D}_j^i}{\delta }_{\ast j}^{\ast i}\right)+2{\displaystyle \sum _{i,s=1}^{m-1}{p}^i}\left({\displaystyle \sum _{j=1}^{m-1}{p}_j}{\delta }_{\ast j}^{\ast s}\right)\left({\displaystyle \sum _{j=1}^{m-1}{D}_j^s}{\delta }_{\ast j}^{\ast i}\right)-{\left({\displaystyle \sum _{i,j=1}^{m-1}{\delta }_{\ast j}^{\ast i}}{D}_j^i\right)}^2\end{array}$

For the covariances, we have

$\begin{array}{l}\text{Cov}(A_P,A_M)={\displaystyle \sum _{i,j=1}^{m-1}{D}_j^i}{\alpha }^{\ast i}{\alpha }_{\ast j}\\ \text{Cov}(A_P,D)={\displaystyle \sum _{i,j=1}^{m-1}{D}_j^i}{\alpha }^{\ast i}{\delta }_{\ast j}^{\ast i}-\left({\displaystyle \sum _{i=1}^{m-1}{p}^i}{\alpha }^{\ast i}\right)\left({\displaystyle \sum _{j,s=1}^{m-1}{D}_j^s}{\delta }_{\ast j}^{\ast s}\right)-{\displaystyle \sum _{j=1}^{m-1}\left({\displaystyle \sum _{i=1}^{m-1}{\alpha }^{\ast i}}{D}_j^i\right)}\left({\displaystyle \sum _{s=1}^{m-1}{p}^s}{\delta }_{\ast j}^{\ast s}\right)\\ \text{Cov}(A_M,D)={\displaystyle \sum _{i,j=1}^{m-1}{D}_j^i}{\alpha }_{\ast j}{\delta }_{\ast j}^{\ast i}-\left({\displaystyle \sum _{j=1}^{m-1}{p}_j}{\alpha }_{\ast j}\right)\left({\displaystyle \sum _{i,t=1}^{m-1}{D}_t^i}{\delta }_{\ast t}^{\ast i}\right)-{\displaystyle \sum _{i=1}^{m-1}\left({\displaystyle \sum _{j=1}^{m-1}{\alpha }_{\ast j}}{D}_j^i\right)}\left({\displaystyle \sum _{t=1}^{m-1}{p}_t}{\delta }_{\ast t}^{\ast i}\right)\end{array}$

As the paternal and maternal alleles are correlated in HWD, we will likely have non-zero covariances among A_P, A_M and D.

Appendix B. One-locus Model for Phase Unknown – in Hardy-Weinberg Disequilibrium

In HWD, we can represent the genotype frequencies as P_AiAi = p²_i + D_ii for j = 1, …, m, and P_AiAj = 2p_ip_j + 2D_ij for i ≠ j. Since p_i = P_AiAi + ∑_j≠iP_AiAj/2, we have ${\displaystyle {\sum }_{i=1}^m{D}_{ij}}=0$ for j = 1, …, m; ${\displaystyle {\sum }_{j=1}^m{D}_{ij}}=0$ for i = 1, …, m; and D_ij = D_ji. As the disequilibrium measures can be represented as D_ij = Cov(x_Pi, x_Mj) for i, j = 1, …, m, we can show that $E(G)=\mu +{\displaystyle {\sum }_{i=1}^{m-1}{\delta }_{ii}^{\ast }}{D}_{ii}+2{\displaystyle {\sum }_{j=2}^{m-1}{\displaystyle {\sum }_{i<j}{\delta }_{ij}^{\ast }}}{D}_{ij}$, and the variance partition of the expected genotypic values V(E(G|g)) = V_A + V_D + 2Cov(A, D), where

$\begin{array}{l}{V}_A=2{\displaystyle \sum _{i=1}^{m-1}{p}_i}{({\alpha }_i^{\ast })}^2-2{\left({\displaystyle \sum _{i=1}^{m-1}{p}_i}{\alpha }_i^{\ast }\right)}^2+2\left({\displaystyle \sum _{i=1}^{m-1}{D}_{ii}}{({\alpha }_i^{\ast })}^2+2{\displaystyle \sum _{i<j}{\alpha }_i^{\ast }}{\alpha }_j^{\ast }{D}_{ij}\right)\\ V_D={\displaystyle \sum _{i,j=1}^{m-1}{p}_i}{p}_j{({\delta }_{ij}^{\ast })}^2-2{\displaystyle \sum _{j=1}^{m-1}{p}_j}{\left({\displaystyle \sum _{i=1}^{m-1}{p}_i}{\delta }_{ij}^{\ast }\right)}^2+{\left({\displaystyle \sum _{i,j=1}^{m-1}{p}_i}{p}_j{\delta }_{ij}^{\ast }\right)}^2+{\displaystyle \sum _{i,j=1}^{m-1}{D}_{ij}}{({\delta }_{ij}^{\ast })}^2-4{\displaystyle \sum _{j=1}^{m-1}\left({\displaystyle \sum _{i=1}^{m-1}{p}_i}{\delta }_{ij}^{\ast }\right)}\left({\displaystyle \sum _{i=1}^{m-1}{D}_{ij}}{\delta }_{ij}^{\ast }\right)\\ +2\left({\displaystyle \sum _{i,j=1}^{m-1}{p}_i}{p}_j{\delta }_{ij}^{\ast }\right)\left({\displaystyle \sum _{i,j=1}^{m-1}{D}_{ij}}{\delta }_{ij}^{\ast }\right)+2{\displaystyle \sum _{i,s=1}^{m-1}{p}_i}\left({\displaystyle \sum _{j=1}^{m-1}{p}_j}{\delta }_{sj}^{\ast }\right)\left({\displaystyle \sum _{j=1}^{m-1}{D}_{sj}}{\delta }_{ij}^{\ast }\right)-{\left({\displaystyle \sum _{i,j=1}^{m-1}{D}_{ij}}{\delta }_{ij}^{\ast }\right)}^2\\ \text{Cov}(A,D)=2{\displaystyle \sum _{i,j=1}^{m-1}{D}_{ij}}{\alpha }_i^{\ast }{\delta }_{ij}^{\ast }-2\left({\displaystyle \sum _{i=1}^{m-1}{p}_i}{\alpha }_i^{\ast }\right)\left({\displaystyle \sum _{s,j=1}^{m-1}{D}_{sj}}{\delta }_{sj}^{\ast }\right)-2{\displaystyle \sum _{j=1}^{m-1}\left({\displaystyle \sum _{i=1}^{m-1}{\alpha }_i^{\ast }}{D}_{ij}\right)}\left({\displaystyle \sum _{s=1}^{m-1}{p}_s}{\delta }_{sj}^{\ast }\right)\end{array}$

In HWD, the genotype frequencies can also be parameterized as P_AiAi = p²_i + p_i(1 − p_i)f, and P_AiAj = 2p_ip_j(1 − f) for i ≠ j (see Weir, 1996). This is a special case of the above parameterization with D_ii = p_i(1 − p_i)f and D_ij = −p_ip_jf, for i ≠ j.