Email updates

Keep up to date with the latest news and content from BMC Genetics and BioMed Central.

Open Access Methodology article

On coding genotypes for genetic markers with multiple alleles in genetic association study of quantitative traits

Tao Wang

Author affiliations

Division of Biostatistics, Institute for Health and Society, Medical College of Wisconsin, Milwaukee, WI 53226, USA

Citation and License

BMC Genetics 2011, 12:82  doi:10.1186/1471-2156-12-82

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2156/12/82


Received:31 May 2011
Accepted:21 September 2011
Published:21 September 2011

© 2011 Wang; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

In genetic association study of quantitative traits using Fmodels, how to code the marker genotypes and interpret the model parameters appropriately is important for constructing hypothesis tests and making statistical inferences. Currently, the coding of marker genotypes in building Fmodels has mainly focused on the biallelic case. A thorough work on the coding of marker genotypes and interpretation of model parameters for Fmodels is needed especially for genetic markers with multiple alleles.

Results

In this study, we will formulate Fgenetic models under various regression model frameworks and introduce three genotype coding schemes for genetic markers with multiple alleles. Starting from an allele-based modeling strategy, we first describe a regression framework to model the expected genotypic values at given markers. Then, as extension from the biallelic case, we introduce three coding schemes for constructing fully parameterized one-locus Fmodels and discuss the relationships between the model parameters and the expected genotypic values. Next, under a simplified modeling framework for the expected genotypic values, we consider several reduced one-locus Fmodels from the three coding schemes on the estimability and interpretation of their model parameters. Finally, we explore some extensions of the one-locus Fmodels to two loci. Several fully parameterized as well as reduced two-locus Fmodels are addressed.

Conclusions

The genotype coding schemes provide different ways to construct Fmodels for association testing of multi-allele genetic markers with quantitative traits. Which coding scheme should be applied depends on how convenient it can provide the statistical inferences on the parameters of our research interests. Based on these Fmodels, the standard regression model fitting tools can be used to estimate and test for various genetic effects through statistical contrasts with the adjustment for environmental factors.

Background

Genetic markers with multiple alleles are common phenomena in genetic studies. It is well known that the ABO blood types in human are determined by three alleles at a genetic locus on chromosome 9. Molecular markers such as microsatellites often have multiple alleles. The major histocompatibility complex (MHC), a highly polymorphic genome region that resides on the human chromosome 6, encompasses multiple genes that encode for many human leukocyte antigens (HLA) and play an important role in regulation of the immune responses. Depending on the resolution level of allele typing, each of the HLA-A, B, C, DR, DQ and DP gene loci could contain tens to hundreds of allele types. In addition, in the haplotype analysis of single-nucleotide polymorphisms (SNPs), various haplotypes from a set of SNPs can also be treated as different alleles from a 'super' marker locus that consists of the set of SNPs.

Presently, there are mainly three types of genetic models that are commonly used in the genetic analysis of quantitative traits. One is Fisher's analysis of variance (ANOVA) models that focus on a decomposition of the genotypic variance into genetic variance components contributed by various genetic effects at quantitative trait loci (QTL) [1-6]. Another is the Fmodels that concentrate on direct statistical modeling of the expected genotypic values at target genetic markers or QTL and the association testing of various genetic effects. The other one is the so-called functional genetic models that emphasize on modeling the functional effects of genes [7]. Both Fisher's and Fmodels can be referred to as statistical models, while the functional genetic models have fundamentally different objectives and estimation methods from the statistical models. A considerable amount of discussion has been made about the distinction between these different types of genetic models [8-11].

The Fmodels have been widely used in genetic association studies of quantitative traits. In building Fmodels, how to code genotypes at a marker (or QTL) and interpret the model parameters are fundamental issues for constructing appropriate testing hypotheses and making correct statistical inferences. While the Fisher's ANOVA models can be directly applicable to genetic markers with multiple alleles, the Fmodels by contrast have been mainly discussed in the biallelic case [1,9,12]. For haplotype analysis, Zaykin et al. in [13] proposed a simple coding which included only the additive effects of haplotypes but ignored their interactions. More recently, Yang et al. in [11] explored an extension of the biallelic Fmodels to multi-allele models with a focus on the definition of various genetic effects and their relationships with the average genetic effects defined in the Fisher's models. A thorough work on coding of marker genotypes and interpretation of model parameters for Fmodels has not been done in the past especially for genetic markers with multiple alleles.

In general, there are two different strategies in coding the marker or QTL genotypes. One is to treat each marker or QTL as a potential risk factor with its genotypes as the risk units. Then, similar to the strategy in handling categorical covariates in classical regression models, at each locus we can create one dummy variable per genotype and then include all but one (as the reference) of these dummy variables into a model. But this genotype coding is often limited by the available sample sizes especially when the number of alleles at the marker locus is large. Alternatively, as alleles are often supposed to be the basic genetic risk units that may contribute to disease phenotypes in genetic studies, we may want to treat alleles at each marker or QTL as the risk units and examine the effects of alleles. However, genetic data has some specialty that needs to be taken into account in order to build the allele-based models. In the genome of diploid species such as human being, alleles normally appear in pairs to form a genotype at each marker locus or QTL with one from the father and one from the mother, except for the sex chromosomes in males. That is, at each locus we have two within-locus risk factors that reside on a homologous pair of chromosomes. Unlike the classical two-way ANOVA model in which the two risk factors own different risk units, the paternal and maternal risk factors at a locus often share the same set of alleles. Besides, the parental origins (i.e., the phase) of the two alleles at each locus are quite often unknown. These features could sometimes complicate the allele-based coding of marker genotypes and generate confusion in interpretation of the model parameters.

In this study, we introduce three allele-based coding schemes for building Fmodels, namely allele, Fand allele-count codings. First, we formulate Fmodels under a general regression framework to model the expected genotypic values at given markers or QTL. Then, under a standard ANOVA model setting, we present several fully parameterized one-locus models using the three allele-based coding schemes. Some potential collinearity relationships among the coding variables of the marker genotypes are clarified. Strategies to avoid the redundant model parameters are also proposed. After that, we examine the definition of model parameters under a reduced one-locus model framework. The impact of a linear relationship among the coding variables of marker genotypes on the estimability of the model parameters is fully explored based on the linear model theory. Finally, we consider extension of the one-locus models to two-locus situation. Several fully parameterized as well as reduced two-locus models are addressed. A focus of this study is to establish the relationships between the model parameters and the expected genotypic values at given marker loci or QTL for various Fmodels from these three coding schemes under various different model frameworks, and explain how to estimate and test for various genetic effects through statistical contrasts. Relationships among different coding schemes and models are also illustrated through simulation.

Results

Fully parameterized one-locus models

In genetic studies, a quantitative trait Y is typically considered as a combination of a genetic component G and an environmental component E with perhaps the genetic by environmental interactions G × E, where G is the true genotypic value from a joint (unobservable) contribution of all the genetic factors to the quantitative trait Y. In practice, given a random sample of N individuals from a study population, let gi be the observed genotypes at certain target marker loci or QTL and zi be a vector of some environmental covariates that may contribute to the variation of the quantitative trait for individuals i = 1, ..., N. By ignoring the genetic by environmental interactions and assuming that the genotypic value G and environmental component E do not depend on the environmental covariates zi and gi, respectively, then the observed quantitative trait yi of an individual i can be expressed through a regression model as

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M1">View MathML</a>

(1)

where G(gi) = E(G|gi) is the expected genotypic value of G given the marker (or QTL) genotypes gi, β denotes the effects of the environmental covariates, and ei is the residual error of the model with E(ei) = 0. Similar to introducing dummy variables for the covariates zi which allow us to assess various environmental effects β in the model, it is convenient to further represent G(gi) as G(gi) = x(gi)α so that we can fit the regression model and assess the genetic effects α of the markers or QTL, where x(gi) is a coding function of the marker genotypes. When the marker locus is not associated with the phenotype, then G(gi) = E(G) is a constant which does not depend on gi. In the rest of the paper, we will focus on the interpretation of the marker effects α in terms of the expected genotypic values G(g) = E(G|g) according to different coding schemes. When certain genetic by environmental interactions are included in the model, the interpretation of α could be modified accordingly. It has to be pointed out that QTL are generally assumed to be unknown genomic regions that may contribute to the variation of the quantitative traits with their genotypes unobserved. But the results (i.e., the coding schemes and the relationships between the model parameters and the expected genotypic values) are held for QTL as well, although the expected genotypic values at a target QTL can no longer be directly estimated via fitting the regression models.

Now, consider one target marker locus with multiple alleles A1, ..., Am, m ≥ 2. In general, there are m possible homozygous genotypes AjAj, j = 1 ..., m, and m(m - 1)/2 possible heterozygous genotypes AjAk, j k. Let Gjk = E(G|g = AjAk) be the expected genotypic values, given the marker genotypes AjAk in a study population. Without knowing the parental origins of the alleles, we assume as usual that the parental origin of the alleles does not make a difference (i.e., no imprinting). We have then Gjk = Gkj for j, k = 1, ..., m, and there are totally m(m + 1)/2 possible distinctive expected genotypic values Gjk, j, k = 1, ..., m, which could be estimated through the means in the genotypic subgroups after adjustment for the environmental covariates. Here we assume no missing genotypes for the sampled individuals, and the random sample has its individuals carrying all possible genotypes. How to handle missing genotypes will be discussed in the discussion. To fully re-parameterize these expected genotypic values through a linear model, we then need totally m(m + 1)/2 parameters including the intercept in the model. By treating the paternal and maternal alleles as two independent risk factors and following the classical two-way ANOVA notation, we can represent the genotypic values Gjk as

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M2">View MathML</a>

(2)

where <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M3">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M4">View MathML</a> are the realized (but unobservable) additive effects of allele Aj and the allelic interaction between the two alleles Aj and Ak, respectively. The above model is different from the classical two-way ANOVA model in that here both the paternal and the maternal risk factors share the same set of alleles A1, ..., Am. As usual, with the unknown paternal origins of alleles at the locus, we assume the paternal and maternal alleles have the same genetic effect. More precisely, the paternal allele Aj and maternal allele Aj have the same additive allelic effects <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M5">View MathML</a> for j = 1, ..., m. Besides, the allelic interaction between a paternal allele Aj and a maternal allele Ak is the same as that between the paternal allele Ak and the maternal allele Aj; i.e., <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M6">View MathML</a>, for j, k = 1, ..., m. Still, with m additive allelic effects and m(m + 1)/2 allelic interactions plus the intercept, it is clear that model (2) is over-parameterized on modeling the m(m + 1)/2 expected genotypic values Gjk for j, k = 1, ..., m. As a result, the parameters μ*, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M7">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M8">View MathML</a> in model (2) are not all estimable in terms of the expected genotypic values Gjk (see [14,15]).

In order to avoid the inestimability issue, one way is to add constraints on the model parameters. However, those constraints, together with the symmetry property of <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M9">View MathML</a>, could make it difficult to fit the model using the standard software package such as SAS. Alternatively, we consider dropping certain redundant parameters in the model. Similar to the biallelic case [10], let us first introduce the following indicator variables to describe the transmission of alleles from parents to their offspring

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M10">View MathML</a>

and

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M11">View MathML</a>

for each allele type Aj, j = 1, ..., m. Then we define the following coding variables of the marker genotypes

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M12">View MathML</a>

for j, k = 1, ..., m, where <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M13">View MathML</a> denotes any other allele type except Aj. Note that z1j, z2jare not observable because we do not know exactly which allele is inherited from paternal or maternal gamete for the sampled individuals without their parental information. But this unknown phase problem does not affect the definitions of wj, vjk since wj only counts the number of allele Aj in the genotypes and the value of vjk is 1 when the genotype is AjAk and 0 otherwise regardless of where the two alleles come from. We refer to the above coding of marker genotypes as an allele coding scheme. Model (2) can then be re-written in a linear model form as

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M14">View MathML</a>

(3)

for i = 1, ..., N. As each individual always carries two alleles at a marker locus with one from the father and the other from the mother, we have <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M15">View MathML</a>, for any i = 1, ..., N. Therefore, given a particular j, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M16">View MathML</a>, which is a linear combination of the rest of {wk, k j}. For vjk, we also have <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M17">View MathML</a>, or <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M18">View MathML</a>. Hence, each of the vjk, k = 1, ..., m, is also a linear combination of the coding variables {wk, k j} and {vlk, l, k j}. To avoid the redundancy of parameters due to these collinearity relationships among the coding variables in model (3), without losing generality, we consider dropping wm and {vkm, k = 1, ..., m} in (3). Then

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M19">View MathML</a>

(4)

for i = 1, ..., N. Model (4) now provides a full re-parameterization of the m(m + 1)/2 expected genotypic values Gjk for j, k = 1, ..., m with its parameters αj can be referred to as the additive allelic effects and δjk the allelic interactions with respect to the reference allele Am. Given a random sample, we can then incorporate model (4) into (1) and fit the regression model (1) using the standard least-square approach. In terms of the expected genotypic values, it is easy to show that μ = Gmm, αj = Gjm - Gmm and δjk = (Gjk - Gkm) - (Gjm - Gmm), for j = 1, ..., m - 1 and k = j, ..., m - 1. Therefore, the additive allelic effect αj can be interpreted as the substitution effect of replacing allele Am by Aj when paired with another allele Am to form the genotypes. Meanwhile, the allelic interaction δjk is the difference between the substitution effect of replacing allele Am by Aj (or Ak) when paired with allele Ak (or Aj) and that when paired with allele Am. Or, in other words, δjk is the difference between the substitution effects of replacing allele Am by Aj (or Ak) with paired alleles Ak (or Aj) and Am. Note that dropping wj and {vkj, k = 1, ..., m} for a particular j m instead of wm and {vkm, k = 1, ..., m} can lead to similar interpretations of the model parameters with Aj being the reference allele. Using model (4), we can also estimate and test for various other genetic effects. For example, the so-called functional 'additive effects' <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M20">View MathML</a> and the 'dominance effects' <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M21">View MathML</a>, j k defined in [11] can be expressed as <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M22">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M23">View MathML</a>, j k, respectively, in terms of the above model parameters. So we can estimate <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M24">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M25">View MathML</a> using the fitted model parameters or test for the hypothesis of <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M26">View MathML</a> or <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M27">View MathML</a> through the general linear contrasts [15] using the standard software such as SAS. To test whether a particular allele Aj has an overall effect, the null hypothesis is H0 : αj = δjk = 0 for k = 1, ⋯, m - 1, which can be performed through either a general linear contrast (or likelihood ratio test) with the degrees of freedom being m for the test statistic. The association test for overall effects of the locus corresponds to the null hypothesis of H0 : αj = δjk = 0 for any j, k = 1, ⋯, m - 1, which has its degrees of freedom being m(m + 1)/2 - 1 for the test statistic. Currently, the so-called Fmodel has been widely used in genetic association studies. In the simple biallelic case with two alleles A and α, an Fmodel gives [16-19].

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M28">View MathML</a>

where GAA = E(G|AA), GAa = E(G|Aa) and Gaa = E(G|aa) are the three possible expected genotypic values at the marker. The parameters a, d are often referred to as the additive and dominance effects of the allele A over a, and in terms of the expected genotypic values we have a = (GAA - Gaa)/2 and d = GAa - (GAA + Gaa)/2. This Fmodel can also be written in a linear model form as [10]

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M29">View MathML</a>

where f, h are two coding variables of the marker genotypes that are defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M30">View MathML</a>

We refer to the above coding of the marker genotypes as the Fcoding. As a straightforward extension of the Fcoding scheme to multiple alleles, we can define the following coding variables

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M31">View MathML</a>

for each j = 1, ..., m. It is easy to see that fj, hj and the previous wj, vjk, j, k = 1, ..., m have the relationships: fj(g) = wj(g) - 1, hj(g) = wj(g) - 2vjj(g), and vjk(g) = hj(g)hk(g) as j k. Thus, for the same reason to avoid collinearity, we can exclude some redundant coding variables and write a fully parameterized one-locus model using the Fcoding as

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M32">View MathML</a>

(5)

for i = 1, ..., N. By having model (5) equivalent to (4), we can first build the relationships between the two model parameters and then establish the relationships between the parameters of model (5) and the expected genotypic values as following

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M33">View MathML</a>

Therefore, aj can be interpreted as a half of the difference between the two expected homozygous genotypic values Gjj and Gmm, which is the same as the additive effect <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M34">View MathML</a> defined in [11]. Besides, djj is the difference between the expected heterozygous genotypic value Gjm and the averaged expected homozygous genotypic value (Gjj + Gmm)/2, which is the same as the dominance effect <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M35">View MathML</a> defined in [11]. It is interesting to see that djk, j k, has the same interpretation as δjk in model (4), which is the difference between the substitution effects of replacing allele Am by Aj when paired with alleles Ak and Am. Note that djj can also be interpreted as the allelic interaction - the difference between the substitution effects of replacing allele Aj by Am when paired with another Aj and Am. In addition, based on model (5), the additive effects <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M36">View MathML</a> and the dominance effects <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M37">View MathML</a> proposed in [11] have the relationship with the model parameters: <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M38">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M39">View MathML</a>, j k. The overall effect of a particular allele Aj can be tested through the composite hypothesis of H0 : aj = djk = 0 for k = 1, ⋯, m - 1, and the overall effects of the locus can be tested via the null hypothesis of H0 : aj = djk = 0 for any j, k = 1, ⋯, m - 1.

In addition to the allele and Fcodings, another way of coding the marker genotypes which occasionally appears in practice is to count the number of alleles in marker genotypes for each specific allele Aj. As each individual can have 0, 1 or 2 copies of an allele Aj, by taking the genotypic group with 0 copy of allele Aj as the baseline, we can introduce the following two indicator (or dummy) variables for the genotypic groups with 1 and 2 copies of the allele Aj, respectively.

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M40">View MathML</a>

for each j = 1, ..., m - 1. These coding variables of marker genotypes have relationships h1j(g) = hj(g) = wj(g) - 2vjj(g) and h2j(g) = vjj(g) with previous ones. We refer to this coding of marker genotypes as the allele-count coding. Similar to models (4) and (5), by excluding some redundant coding variables, the allele-count coding leads to another fully parameterized one-locus model as

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M41">View MathML</a>

(6)

for i = 1, ..., N. Similarly, by having model (6) equivalent to (4), we can establish the following relationships

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M42">View MathML</a>

Therefore, πj in model (6) can still be interpreted as the substitution effect of replacing allele Am by Aj when paired with allele Am, or the difference between the genotypic values of the genotype group AjAm with one copy of Aj versus the genotype group AmAm (baseline). ηjj is the difference between the expected genotypic value Gjj in the homozygous genotypic group AjAj with two copies of Aj and Gmm in the baseline group AmAm. Besides, ηjk in model (6) has the same interpretation as δjk (or djk) before. From model (6), the general additive effects <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M43">View MathML</a> and the dominance effects <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M44">View MathML</a>, j k, which can be tested either separately or jointly. The overall effect of a particular allele Aj can be tested through the composite hypothesis of H0 : πj = ηjk = 0 for k = 1, ⋯, m - 1. The overall effects of the locus can also be tested via the null hypothesis of H0 : πj = ηjk = 0 for any j, k = 1, ⋯, m - 1.

Each of the three models (4), (5) and (6) provides a full re-parameterization of the m(m + 1)/2 expected genotypic values under the same model framework (3). The relationships between their model parameters and the expected genotypic values are summarized in Table 1. It is interesting to see from Table 1 that the null hypothesis of αj = δjj = 0 is equivalent to either aj = djj = 0 or πj = ηjj = 0, which implies Gjj = Gjm = Gmm. So the three models above should provide the same test statistics for testing αj = δjj = 0, aj = djj = 0 or πj = ηjj = 0.

Table 1. Parameterization of fully parameterized one-locus models (4), (5), (6).

For a biallelic locus with alleles A (or A1) and a (or A2), we have m = 2 with three possible genotypic values GAA = E(G|AA), GAa = E(G|Aa) and Gaa = E(G|aa). If we adopt the allele coding, then w2(g) = 2 - w1(g), v12(g) = w1(g) - v11(g), and v22(g) = 1 - w1(g) + v11(g). For the Fcoding, we have f2(g) = -f1(g) and h2(g) = h1(g). So we can further drop d2 in model (5). For the allele-count coding, we have h12(g) = h11(g) and h22(g) = 1 - h11(g) - h21(g). The interpretation of model parameters for these three biallelic QTL models are summarized in Table 2, which is a special case of Table 1.

Table 2. Parameterization of one-locus models (4), (5), (6) when m = 2.

For a locus with three alleles A1, A2 (i.e., m = 3), we have six possibly distinctive expected genotypic values G11, G22, G33, G12, G13 and G23. Each of the three fully parameterized models (4), (5) and (6) can provide a full re-parameterization of the six expected genotypic values. In a matrix form, from the allele coding model (4), we have

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M54">View MathML</a>

From the Fcoding model (5), we have

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M55">View MathML</a>

And the allele-count coding model (6) gives

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M56">View MathML</a>

By multiplying the design matrices on the left side of the equations, we can show that the model parameters and the expected genotypic values have the relationships as summarized in Table 3, which is consistent with that in Table 1.

Table 3. Parameterization of one-locus models (4), (5), (6) when m = 3.

Reduced one-locus models

Due to limited available sample sizes in practice, it may not always be feasible to use the fully parameterized models. Quite often, one may want to check the main effects of alleles first before including all possible allelic interactions. Here we consider the case of including possible interactions between Aj and itself for the homozygous genotypes Aj Aj, j = 1, ..., m, but ignore other interactions between different alleles Aj and Ak (j k). Then we obtain a reduced case of model (2) as below

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M60">View MathML</a>

(7)

for j, k = 1, ..., m. Similarly, using the allele coding, we can present this model in a linear model form as

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M61">View MathML</a>

(8)

for i = 1, ..., N, where vj(g) = vjj(g) for j = 1, ..., m, with vjj(g) defined as before.

Model (8) contains only one redundant parameter in the α*'s due to the fact that <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M62">View MathML</a> for i = 1, ..., N. In this case, as shown in Appendix A, the parameters <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M63">View MathML</a> in model (8) are estimable but the parameters μ* and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M64">View MathML</a> are not estimable. To overcome the redundant parameter problem, we can drop wm from model (8) and consider

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M65">View MathML</a>

(9)

for i = 1, ..., N. Note that <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M66">View MathML</a>, which cannot be completely determined by {wj, vj, j = 1, ..., m - 1}. Therefore, dropping {δjk, j, k = 1, ..., m - 1, j < k} from model (4) does not directly lead to an equivalent model of (9) as the latter contains vm. In fact, as further dropping vm in (9), it will lead to a more restricted model structure for the expected genotypic values with the similar interpretation of its model parameters as presented in model (4). It is also interesting to see that the haplotype coding proposed in [13] is a special case of model (9) when we further ignore all the allelic interactions and drop all the {vj, j = 1, ..., m} in the model.

By definition, a reduced model can be derived from its original model by adding certain restrictions on the model parameters. Typically, the model parameters in a reduced model could be interpreted similarly as that in its original model when these restrictions are simple enough (e.g., by setting a subset of them being zero). When the restrictions on the original model parameters are complicated, however, the interpretation of the reduced model parameters could be different from that presented in the original model. For model (9), we can establish the relationship between its model parameters and the expected genotypic values using a classical matrix approach, as shown in Appendix B. An alternative way of building this relationship is to simply treat model (9) as a reduced form of model (8) by adding a restriction <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M67">View MathML</a> and taking μ = μ*, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M68">View MathML</a> for j = 1, ..., m - 1, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M69">View MathML</a> for j = 1, ..., m. Note that adding the restriction <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M70">View MathML</a> on (8) does not change the modeling structure of the expected genotypic values because <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M71">View MathML</a> is a redundant parameter given the others. Therefore,

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M72">View MathML</a>

Comparing with the parameters in model (4), we can see that the interpretation of the parameters in model (9) have changed slightly. The intercept μ now becomes <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M73">View MathML</a> instead of Gmm, the αj is the substitution effect of replacing allele Am by Aj when paired with any allele Ak (k j, m) instead of just Am, while the δj is the difference between the substitution effect of replacing any allele Ak by Aj when paired with Aj itself and that when paired with another allele Al (l j, k). If both αj and δj are zero for a particular j < m, then Gjj = Gjm = μ and Gjk = Gkm for any k j, m.

Under the same model framework (8), the Fcoding leads to the following model

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M74">View MathML</a>

(10)

for i = 1, ..., N. By applying the relationship fj(g) = wj(g) - 1 and hj(g) = wj(g) - 2vj(g) for j = 1, ..., m, we can show that for models (10) and (8) to be equivalent their model parameters have the relationship

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M75">View MathML</a>

In other words, model (10) leads to a restriction <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M76">View MathML</a> on the parameters in model (8) which makes <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M77">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M78">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M79">View MathML</a>, j = 1, ..., m - 1. Thus,

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M80">View MathML</a>

Now dj becomes a half of the difference between the substitution effect of replacing any allele Ak by Aj when paired with another Aj and that when paired with an allele Al (l j, k), which can no longer be referred to as a dominance effect.

With the allele-count coding, we can actually construct two equivalent models in this case

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M81">View MathML</a>

(11)

and

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M82">View MathML</a>

(12)

for i = 1, ..., N. Similarly, we can show that model (11) can be treated as a reduced model by adding the restriction <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M83">View MathML</a> on parameters in model (8) with the following relationships

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M84">View MathML</a>

On the other hand, model (12) can be treated as a reduced model by adding the restriction <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M85">View MathML</a> on parameters in model (10) with the following relationships

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M86">View MathML</a>

While the effect ηjj in model (6) is the difference between the two expected homozygous genotypic values Gjj and Gmm, the effect ηj in model (11) becomes the sum of the substitution effects of replacing allele Am by Aj when paired with Aj itself and when paired with another allele Ak (k j, m. It is also interesting to see that the definition of parameters in models (11) and (12) are quite different. A null hypothesis of <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M87">View MathML</a> for a particular j < m in model (12) implies that Gjj = Gmm and Gjm - Gmm = Gjk - Gkm for any k j, m, while the null hypothesis of H0 : πj = ηj = 0 for a j < m in model (11) implies that Gjj = Gjm and Gjk = Gkm for any k j, m, which has nothing to do with Gmm.

Under the same model framework (8), each of the above four models (9), (10), (11) and (12) contains 2m non-redundant parameters (including the intercept) to model the m(m + 1)/2 expected genotypic values. When m > 3, we have m(m + 1)/2 > 2m. Therefore, the model framework (7) enforces certain constraints on the m(m + 1)/2 genotypic values. If m = 3, then each of the four models actually provides a full re-parameterization of the six expected genotypic values G11, G22, G33, G12, G13 and G23. The relationships between the four model parameters and the expected genotypic values are summarized in Table 4.

Table 4. Parameterization of one-locus models (9), (10), (11), (12) when m ≥ 3.

Comparing Table 4 with Table 1, we can see that the definition of model parameters depends not only on the coding schemes of marker genotypes but also on the underlying framework for the structure of the expected genotypic values. From Table 4, it is also interesting to see that the null hypothesis of H0 : αj = δj = 0 (j < m) in model (9) is equivalent to πj = ηj = 0 in model (11), which implies <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M96">View MathML</a> in model (8) with restriction <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M97">View MathML</a>, or Gjk = Gkm for any k = 1, ..., m. On the other hand, the null hypothesis of H0 : aj = dj = 0 (j < m) in model (10) is equivalent to <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M98">View MathML</a> in model (12), which implies <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M99">View MathML</a> in model (8) with a restriction <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M100">View MathML</a>, or Gjj = Gmm and Gjj - Gjm = Gjk - Gkm for any k m. In general, the two null hypotheses of αj = δj = 0 and aj = dj = 0 may not always be equivalent. For example, when m = 3, similar to the three-allele models discussed in the previous section, we can show that the four model parameters and the expected genotypic values have the relationships as shown in Table 5, which is a special case of Table 4. We can see from Table 5 that α1 = δ1 = 0 is equivalent to π1 = η1 = 0 which implies G12 = G23 and G11 = G13; while a1 = d1 = 0 is equivalent to <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M101">View MathML</a> which implies G11 = G33 and G12 + G13 = G11 + G23. So, depending on the underlying true setting of the expected genotypic values, the null hypotheses of α1 = δ1 = 0 in model (9) could be different from that of a1 = d1 = 0 in model (10).

Table 5. Parameterization of one-locus models (9), (10), (11), (12) when m = 3.

Extension to two-locus models

In this section, we further explore some extensions of the previous one-locus models to two-locus models. Consider two marker loci with alleles <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M110">View MathML</a> at locus 1 and alleles <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M111">View MathML</a> at locus 2, respectively. Without distinguishing the parental origins of the alleles, there are totally m1m2(m1 + 1)(m2 + 1)/4 possible distinctive expected genotypic values: Gjkrs = E(G|A1jA1kA2rA2s) for j, k = 1, ..., m1, j k; and r, s = 1, ..., m2, r s. Using the allele coding, we introduce the following coding variables

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M112">View MathML</a>

j, k = 1, ..., m1, for marker genotypes at locus 1 and

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M113">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M114">View MathML</a>

r, s = 1, ..., m2, for marker genotypes at locus 2, where <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M115','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M115">View MathML</a> (or <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M116','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M116">View MathML</a>) denotes any other allele type except A1j(or A2r) at locus 1 (or 2). A fully parameterized two-locus model for Gjkrs can then be presented as

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M117">View MathML</a>

(13)

for i = 1, ..., N. Similar to the one-locus models, we can establish the relationship between the model parameters and the expected genotypic values as shown in (C.1) of Appendix C. A nice property of this allele coding model is that a higher order effect is simply the deviation of its corresponding expected genotypic value from an approximation of the other lower order effects. Here the corresponding expected genotypic value of a marker effect is determined by the position of alleles that differ from the two reference alleles <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M118">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M119">View MathML</a>. So, starting from the lowest order parameter μ, it seems straightforward to build the relationships between the model parameters and the expected genotypic values starting from the low-order effect parameters up to the high-order effect parameters.

For the Fcoding, we can define the following coding variables for the genotypes at the two marker loci separately.

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M120">View MathML</a>

for j = 1, ..., m1, and

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M121','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M121">View MathML</a>

for r = 1, ..., m2. A fully parameterized two-locus model using this Fcoding is then

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M122','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M122">View MathML</a>

(14)

for i = 1, ..., N. Still, using the relationships w1j= 1 + f1j, w2r= 1 + f2r, v1jj= (1 + f1j- h1j), v2rr= (1 + f2r- h2r), v1jk= h1jh1kfor j k, and v2rs= h2rh2sfor r s between the Fcoding variables and the allele coding variables, we can establish the relationships between the model parameters and the expected genotypic values as shown in (C.2) of Appendix C. We can easily verify that the biallelic two-locus effects <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M123','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M123">View MathML</a> in [9] is a special case of our results with m1 = m2 = 2. It is also interesting to see that the interpretation of model parameters in terms of the expected genotypic values becomes much more complicated than that in the previous allele coding model. When m1, m2 > 2, the low-order within-locus main effect a1jis a weighted combination of the differences <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M124">View MathML</a>, where r = 1, ..., m2 refer to various homozygous genotypes A2rA2rat locus 2. The within-locus effect d1jjis a weighted combination of the allelic interactions <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M125">View MathML</a>, r = 1, ..., m2, at locus 1 with reference A2rA2rat locus 2. Even the intercept τ of the model becomes a complex function of various homozygous genotypic values.

Applying the allele-count coding, we can define

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M126','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M126">View MathML</a>

for j = 1, ..., m1, and

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M127">View MathML</a>

for r = 1, ..., m2. Another fully parameterized two-locus model for Gjkrs can be written as

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M128">View MathML</a>

(15)

for i = 1, ..., N. In this case, the allele-count coding variables and the allele coding variables have the relationships <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M129">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M130">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M131','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M131">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M132','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M132">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M133','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M133">View MathML</a> for j k, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M134','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M134">View MathML</a> for r s. Through the equivalence of the two models (13) and (15), we can also construct relationships between the parameters in model (15) and the expected genotypic values as shown in (C.3) of Appendix C. We can see that the interpretation of parameters in the allele-count coding model (15) are as simple as that in the allele coding model (13) with the same intercept being <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M135','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M135">View MathML</a>. Besides, it seems that some parameters such as (η1jjη2rr), (η1jkη2rs) and (η1jkη2rr) have simpler relationships than the corresponding ones in the allele coding model (13).

Finally, let us consider some reduced cases of the two-locus models. By ignoring locus-by-locus interactions (i.e., epistases), we have the following simplified two-locus model framework

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M136','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M136">View MathML</a>

(16)

for j, k = 1, ..., m1 and r, s = 1, ..., m2. If we further ignore the within-locus allelic interactions between different alleles, then another reduced two-locus model framework is

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M137','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M137">View MathML</a>

(17)

Similar to the one-locus models, under each of the two reduced model frameworks we can construct the two-locus models from the three coding schemes. The relationships between the model parameters and the expected genotypic values under framework (14) are summarized in Table 6, which can be treated as an extension of Table 1 to the two-locus case. The relationships between the model parameters and the expected genotypic values under framework (17) are also summarized in Table 7, which is a straightforward extension of Table 4. Further dropping <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M138','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M138">View MathML</a> for j = 1, ..., m1 and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M139','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M139">View MathML</a> for r = 1, ..., m2 in (15) will lead to an additive model framework, which has its model parameters interpretable similar to that in Table 6. From Tables 6 and 7, we can see that both the allele and allele-count coding models have their lower-order main effects keep similar interpretation as to that in the previous fully parameterized case with epistases, while the Fcoding models have the definition of their lower-order main effects vary depending on whether there are epistases involved in the models.

Table 6. Parameterization of two-locus models under model framework (16).

Table 7. Parameterization of two-locus models under model framework (17) when m1, m2 ≥ 3.

As pointed out in [9], the genetic effects of a marker may have different interpretation depending upon whether the marker is fitted in a one-locus model or a two-locus model. From the linear model theory, the genetic effects of a marker in a one-locus model are defined based on the expected genotypic values of certain genotypes at this particular marker locus with genotypes at the other marker loci being averaged out based on the joint genotype distribution. For instance, marker 1 in the two-locus setting above has its effects defined in a one-locus model based on the one-locus genotypic values <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M151','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M151">View MathML</a>, which could depend on the LDs of alleles between the two loci. When the same marker is fitted in a two-locus model, its effects are usually functions of the expected genotypic values with their joint genotypes taking certain reference alleles or genotypes at the other marker loci. So, in general, even without locus-by-locus interactions, a single marker's effects could be different from the one defined in a multi-locus model when the alleles at different loci are in linkage disequilibrium (LD). Consider a 2-locus haploid model with alleles A, a at locus 1 and B, b at locus 2. If we ignore the locus-by-locus interaction, it is easy to show that the additive allelic effects are α1 = GAB - GaB = GAb - Gab and α2 = GAB - GAb = GaB - Gab at locus 1 and 2, respectively. In a one-locus model at locus 1, however, we can show that the locus has its additive allelic effect <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M152','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M152">View MathML</a>, where D = PAB - pApB is the LD between the two loci.

Simulation Examples

We use some numerical examples to illustrate properties of the models we have discussed. First, we consider the same example discussed in [11] of a three-allele locus with allele frequencies p1 = 0.2 for A1, p2 = 0.3 for A2, and p3 = 0.5 for A3. The six genotypic values are G11 = 10, G12 = 30, G22 = 50, G13 = 36, G23 = 46 and G33 = 42. We adopt a similar strategy to specify the genotype frequencies as: <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M153','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M153">View MathML</a> for j = 1, 2, 3 and Pjk = 2pjpk + D for j k, where D is a measure of departure from Hardy-Weinberg equilibrium (HWE) for the three alleles at the locus and D- D D+ with

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M154','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M154">View MathML</a>

and

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M155','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M155">View MathML</a>

We consider two cases: i) D = 0 for HWE, and ii) D = 0.02 for Hardy-Weinberg disequilibrium (HWD). The phenotypic value of an individual is simulated as a sum of its true genotypic value and an environmental noise from N(0, σ2), where the σ2 is chosen to be either 0 or σ2 = 288 with the latter one corresponds to a 20% heritability level when D = 0. For each of the four configurations, we simulate 10,000 random samples with 1000 individuals each. For each random sample, we fit the three fully parameterized one-locus models (4), (5) and (6) under model framework (2) using the least square approach and estimate the model parameters as well as the six genotypic values. The means and standard deviations (SD) of the least square estimates (LSE) of the model parameters and the six genotypic values from the 10,000 random samples in fitting these three models are summarized in Table 8.

Table 8. Means (SD) of LSE for three one-locus models (4), (5) and (6) when m = 3.

As each of the three models provides a re-parameterization of the six genotypic values, for each random sample the three models always give exactly the same estimates of the six genotypic values and the residual variance as we expected, even though their model parameters are defined in different ways. As a result, under each configuration, the three models have the same means and SD for the LSE of the six genotypic values and the residual variance. Without environmental variation, each model can accurately estimate its model parameters and the six genotypic values for each random sample regardless of whether there is HWE or HWD. When there is environmental variation on the phenotypes, it is known that the least square estimators of the model parameters are unbiased under either HWE or HWD. However, the HWD may affect the variance of the least square estimators of the model parameters and the six genotypic values. Note that the genotypic frequencies are P11 = 0.04, P22 = 0.09, P33 = 0.25, P12 = 0.12, P13 = 0.20 and P23 = 0.30 under HWE, while with D = 0.02 the genotypic frequencies become P11 = 0.02, P22 = 0.07, P33 = 0.23, P12 = 0.14, P13 = 0.22 and P23 = 0.32. So, under HWD, we tend to have more individuals carrying genotypes A1A2, A1A3, A2A3 but less individuals carrying genotypes A1A1, A2A2, A3A3 in the random samples than that under HWE. Without knowing the accurate genotypic values, more individuals with certain genotypes in a random sample can then provide better estimates of the corresponding genotypic values. This explains why under HWD the estimates of G11, G22 and G33 have larger SD (or variances) than that under the HWE, and the estimates of G12, G13 and G23 under HWD have smaller variances than that under the HWE.

As another example, let us consider the statistical modeling of two-locus genotypic values Gjkrs, where the first locus have three alleles A1, A2, A3 and the second locus have two alleles B1, B2. Assume that the alleles at locus 1 have the same allele frequencies as that in the previous example; i.e., p1 = 0.2 for A1, p2 = 0.3 for A2, and p3 = 0.5 for A3, while the two alleles at locus 2 have frequencies q1 = 0.2 for B1 and q2 = 0.8 for B2. The two-locus genotypic values G2 = (Gjkrs), j, k = 1, 2, 3; r, s = 1, 2 are given by

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M156','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M156">View MathML</a>

which are modified values from the previous one-locus model in a way that the Gjk11 = Gjk, Gjk12 = Gjk + e1jkand Gjk22 = Gjk - e2jkwith e1jkand e2jkbeing some small positive fluctuations according to the genotypes B1B2 and B2B2 at locus 2. We assume Hardy-Weinberg equilibria at both loci and specify their haplotype frequencies as: h11 = p1q1 - D1, h12 = p1q2 + D1, h21 = p2q1 - D2, h22 = p2q2 - D2, h31 = p3q1 + (D1 - D2), h32 = p3q2 - (D1 - D2), where D1 (and D2) are the linkage disequilibria (LD) between alleles A1 and B2 (and A2 and B1) at the two loci. We consider two scenarios: i) D1 = D2 = 0 for linkage equilibrium (LE); and ii) D1 = 0, D2 = 0.03 for LD. The phenotypic value of an individual is still simulated as a sum of its genotypic value and an environmental noise from N(0, σ2), where the σ2 was chosen to be either 0 or σ2 = 286 with the latter one corresponds to a 20% heritability level when D1 = D2 = 0. For each of the four configurations, we simulate 10,000 random samples with 1000 individuals each. For each random sample, we consider fitting models under three model frameworks: i) one-locus models (4), (5) and (6) at locus 1 under model framework (2); ii) two-locus models without epistases from the three coding schemes under model framework (14); iii) fully parameterized two-locus models (13), (14) and (15) with epistases. Still, for each random sample, the three allele coding models under the same model framework give exactly the same estimates of the 18 genotypic values as we expected (results not shown here). As the result, under each model framework, the three models have the same means and SD for the LSE of the 18 genotypic values and the residual variance, although the means and SD for the LSE of their model parameters are different. To compare the LSE of model parameters for models from the same coding under different model frameworks, we summarize in Table 9 the means and SD of the LSE of the model parameters from the 10,000 random samples in fitting the three allele-coding models: the one-locus model (4), the two-locus model under model framework (14), and the two-locus model under model framework (13). Models from the other two coding schemes behave similarly.

Table 9. Means (SD) of LSE for three allele-coding models regarding the two-locus genotypic values

As we mentioned before, the one-locus models are actually modeling the expected genotypic values given the genotypes at locus 1. When D1 = D2 = 0, we can show that the expected genotypic values at locus 1 are G11 = 10.03, G22 = 50.03, G33 = 41.68, G12 = 29.90, G13 = 35.87 and G23 = 45.71, which correspond to μ = 41.68, α11 = -5.81, α12 = 4.03, δ111 = -20.03, δ122 = 0.29 and δ112 = -10 as the true parameters in the allele coding one-locus model. When D1 = 0, D2 = 0.03, the expected genotypic values at locus 1 become G11 = 10.03, G22 = 50.08, G33 = 41.55, G12 = 29.97, G13 = 35.81 and G23 = 45.77, which correspond to μ = 41.55, α11 = -5.74, α12 = 4.21, δ111 = -20.04, δ122 = 0.09 and δ112 = -10.06 as the true parameters in the allele coding one-locus model. In both cases, the least square estimators of the one-locus model parameters are unbiased estimators of the true parameters. Note that, unlike the one-locus model in the previous example, the LSE of the model parameters are no longer exactly the same as the true values even when no environmental noises are involved. The reason is that the expected genotypic values at locus 1 depend on not only the genotypic values but also the joint genotype frequencies in the sample, which may change slightly from sample to sample due to the sampling variation.

For the two-locus model without epistases, it cannot provide unbiased estimators for all the genotypic values because of the model mis-specification. However, the LSE of its parameters associated with locus 1 are similar to the ones in the one-locus model at locus 1. In fact, as we know from the linear model theory, the true values of its parameters associated with locus 1 are the same as the ones defined in the one-locus model at locus 1 when the two loci are in LE. Under LD, the least square estimators of its model parameters associated with locus 1 could be biased, and the biasness depends on the LD setting.

The two-locus model with epistases gives a full re-parameterization of the 18 genotypic values. Therefore, when no environmental noises are involved, the LSE of its model parameters are exactly the same as their true values for each random sample regardless of the LD between the two loci. It has to be pointed out that this phenomenon holds only when the random sample contains all the 18 possible genotypes. In our simulation setting, the frequencies for certain genotypes such as A1A1B1B1, A1A3B1B1 and A2A2B1B1 are pretty small. As the result, we occasionally (about 22-23% of the 1000 random samples) may obtain a random sample that has no individuals carrying certain genotypes. In this case, the design matrix in the fully parameterized model becomes singular and the LSE of the model parameters are no longer unique. To keep our illustration of the model properties simple, we excluded those random samples in fitting the two-locus model with epistases (reduced models are less likely to have singular design matrices). Other techniques such as ridge regression could be applied to handle those skewed random samples. In the presence of environmental noises, it is also noted that the LSE for some of its model parameters such as δ111, (δ111α21) and (δ111δ211) have much larger SD than the LSE of other parameters. This is due to the low frequencies of genotypes A1A1B1B1, A1A3B1B1 and A2A2B1B1. As a random sample has few individuals carrying these genotypes, it has reduced accuracy in estimation of their corresponding true genotypic values to which the model parameters δ111, (δ111α21) and (δ111δ211) are related.

Discussion

In this study, we introduced three genotype coding schemes to build Fmodels for multi-allele markers. The relationship between the model parameters and the expected genotypic values were established in some fully parameterized as well as reduced one-locus and two-locus Fmodels. Our results showed that the relationships between the model parameters and the expected genotypic values could become more intricate in the multi-allele case than that in the biallelic case, even though the extension of the coding schemes from biallelic to multiple alleles appears straightforward. We built the relationships between different model parameters mainly through their coding variables of marker genotypes, which simplified the tedious derivation process comparing with the classical matrix approach. The Fmodels we proposed can be used directly for association testing of multi-allele markers and their possible interactions with quantitative traits using random unrelated samples. These Fmodels could also be applied to test for the risk haplotypes and their interactions when incorporated with the likelihood approach (e.g., [20]), or analyze family data by combining them with the likelihood to account for the transmission probability of alleles from parents to their offspring. Although our discussion focused on genetic modeling of quantitative traits, the results can be extended to other phenotypic traits such as binary outcomes in case-control studies using logistic regression models or time-to-event data using the Cox proportional hazard models.

Throughout the paper, we assumed that all the possible genotypes are available from the sampled individuals. If certain genotypes are not observable, then the expected genotypic values on these genotypes will not be estimable by themselves, which could change the interpretation of the model parameters as well. The models we have presented can also be modified to handle the situation when some individuals have missing genotypes at certain marker loci. When the missing genotypes at a marker locus have both alleles missing at the same time, we can simply introduce an indicator variable to code for the missing genotype at the marker. The regression coefficient of this indicator variable for this missing genotype can usually be interpreted as the difference between the expected genotypic value with missing genotype at the marker locus and the intercept of the model, while the other regression coefficients would keep the same interpretation as before.

It has to be pointed out that the relationships between the model parameters and the expected genotypic values are based on the assumption that the models can correctly specify the structure of the expected genotypic values. When a fully parameterized model is applied, the definition of its model parameters do not depend on the allele frequencies, HWD among alleles within a locus, or LD structure between alleles at different loci. In fitting a reduced model, however, a simplified model may not be totally correct in modeling all the expected genotypic values. In this case, depending on how accurate the simplified model is on approximating the expected genotypic values, the allele frequencies, HWD and LD structure between marker alleles could affect the definition and LSE of its model parameters. In the presence of environmental variation on the phenotypic values, regardless of whether a fully parameterized or reduced model is applied, the allele frequencies, HWD or LD between marker alleles may affect the LSE of the model parameters and the power in detection of the associated marker alleles as shown in our simulation studies.

All the models we have discussed so far are Fmodels. Statistically, these Fmodels are fixed-effect models which focus on modeling the expected genotypic values directly. On the other hand, the Fisher's ANOVA models, which target on evaluation of the variations contributed by various allelic effects and interactions, can be treated as random-effect models (see [21]) in which the expected genotypic values come from a discrete random variable G(g) = E(G|g) with its limited genotypes g being randomly sampled from a study population. Both the Fand the Fisher type models form basis in the analysis of quantitative traits and they provide different perspectives in assessing the genetic effects of QTL and markers. For biallelic markers, we proposed in [10] a 'mean corrected' Fisher (mc-Fisher) model for decomposition of the genotypic variances. In the multi-allele marker case, we can also construct similar mc-Fisher models by applying mean corrections on all the indicator variables of the paternal and maternal alleles in the allele coding Fmodels. For example, based on the allele coding model (4), we can construct its corresponding mc-Fisher model by replacing the coding variables wj and vjk with <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M157','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M157">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M158','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M158">View MathML</a>, respectively; where pj is the allele frequency of Aj. Then the genetic additive and dominant variance components VA and VD of G(g), which are defined as variations contributed by the additive allelic effect and allelic interactions respectively, can be estimated from <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M159','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M159">View MathML</a>'s and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M160','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M160">View MathML</a>'s separately. As pointed out in [10], the mc-Fisher model can provide an orthogonal partition of V(G) into the sum of VA and VD under Hardy-Weinberg equilibrium, and it can be fitted through the standard least-square regression approach. Similar to the Fmodels, the definition of the model parameters in such a mc-Fisher model also depend on the choice of the reference allele 'Am'. But the estimates of the additive and dominant variance components VA and VD do not depend on such a choice. In addition, when a fully parameterized model is applied, the mc-Fisher model is equivalent to its original Fmodel in modeling the expected genotypic values. Therefore, both models have the same residual variance and the F-statistics in testing for the overall effect of the marker locus. When reduced models are applied, the mc-Fisher model could become inequivalent to its original Fmodel especially when allelic interactions are involved.

Of the three coding schemes that we have discussed, the Fcoding is perhaps the most widely used in current genetic association studies of quantitative traits. From what we have shown, the three coding schemes can essentially lead to equivalent models and have the same power in detection of various genetic effects. In practice, just like the various existing coding schemes such as 'Reference', 'GLM' and 'Effect' that are commonly used in the analysis of categorical covariates [22], we usually only need to adopt one specific coding scheme in building the regression models. Which coding scheme should be applied depends on how convenient it can provide the statistical inferences on the parameters of our research interests. In general, the allele coding models can provide direct estimates of certain substitution effects of alleles and allelic interactions and, in the two-locus case, allele coding models are perhaps the easiest among the three codings in building the relationships between their model parameters and the expected genotypic values. Besides, they are generically linked to the genetic variance components as we have shown above. On the other hand, the allele-count coding models are attractive in that it often leads to simple comparisons among the three genotypic groups with 0, 1 or 2 copies of a particular allele. In the two-locus case, the allele-count coding models also have the definition of their model parameters remain as simple as (if not simpler than) that in the allele coding models even in the presence of epistases. Meanwhile, both the allele and allele-count coding show an advantage that their lower-order main effects in the models can keep the same interpretation regardless of whether there are epistases involved in the model or not. In contrast, the Fcoding models may have the definition of their lower-order main effects vary depending on the absence or presence of epistases in the models. Even though the one-locus Fcoding model parameters are closely related to the additive and dominance effects, the two-locus Fcoding model parameters including the lower-order main effects have more complicated interpretations than that in the allele or allele-count coding models especially when epistases are involved.

The coding of marker genotypes are not limited to the three allele-based coding schemes that we have discussed. Application of a coding scheme could also be subject to the number of individuals available in each genotype group. For example, under the model framework (7), the allele coding scheme typically creates wj(g), vj(g) for each allele type Aj, j = 1, ..., m. When the group of a homozygous genotype AjAj includes very few individuals for a particular allele Aj, we may want to combine this genotypic group with another genotype such as the one carrying one copy of the allele Aj. Then we can replace the original wj(g) and vj(g) by an allele presence-absence coding variable dj(g) for this specific allele Aj while keeping two coding variables wk(g), vk(g) for other alleles Ak, which leads to a mixed use of the allele coding and this allele presence-absence coding variable. In certain situations, the genotype-based coding could also be very useful as it can provide direct tests on pair-wise comparisons of certain genotypic values. Comparing with the genotype-based coding, the allele coding has the advantage of further dissecting the genetic effects into the allelic effects and allelic interactions, which allow us to specify reduced models with varying degrees of interactions among the main allelic effects - a useful tool in the model building procedures. Given a fixed coding, the likelihood ratio test can be applied to compare a full model with its reduced models. Statistical model selection tools such as AIC and BIC criteria, which provide a balance between the goodness of model fitting to the data and the complexity of the models in terms of the number of parameters, could also be used to compare some non-nested reduced models or frameworks. The current study focuses on establishing the theoretical relationships between the model parameters and the expected genotypic values according to different coding schemes under various model frameworks. A power comparison of some reduced models from different coding schemes under various scenarios with respect to the allele frequencies and possible HWD or LDs between marker alleles is beyond the scope of this study and might be worth of further exploration.

Conclusions

In summary, we introduced three allele-based coding schemes to construct Fmodels for association testing of multi-allele genetic markers with quantitative traits. Depending upon whether certain allelic effects or comparisons between genotypic groups are of the main research interest, investigators may adopt one of the three allele-based codings (i.e., allele, For allele-count), or perhaps a genotype-based coding in building an Fmodel. Based on the Fmodel from a given coding scheme, standard regression model fitting tools can then be applied to estimate or test for various genetic effects. Understanding the definition of model parameters from different coding schemes under various model frameworks are crucial for constructing appropriate testing hypothesis and making the correct statistical inferences in the genetic association studies.

Authors' contributions

TW planned the study, conducted the derivation and wrote the manuscript.

Appendices

A. Estimability of parameters in model (8)

Let G = (G(g1), ...,G(gN))T denote a vector of the expected genotypic values of all the individuals in the sample, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M161','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M161">View MathML</a> be a vector of all the model parameters. We can rewrite model (8) in a matrix form as G = * +e, where e = (e1, ..., eN) and the design matrix X is

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M162','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M162">View MathML</a>

(18)

with Wj = (wj(g1), ..., wj(gN))T and Vj = (vj(g1), ..., vj(gN))T for j = 1, ..., m. As every individual carries two and only two alleles at the locus, we have <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M163','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M163">View MathML</a>, which means that the first (m+1) column vectors 1N, W1, W2, ..., Wm of the design matrix X are linearly dependent. So, rank(X) ≤ 2m; i.e., X is not a full column rank matrix.

From (7), we have <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M164','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M164">View MathML</a> for j k, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M165','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M165">View MathML</a>. If we write G0 = (G12, G13, ..., G1m, G23, ..., GL-1, L, G11, ..., Gmm)T, then this model gives

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M166','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M166">View MathML</a>

where s = m(m - 1)/2, and

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M167','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M167">View MathML</a>

Assume that the genotypes of the sampled individuals cover all possible genotypes AjAk for j, k = 1, ..., m. Then the design matrix X includes all the row vectors of X0, which implies that rank(X) ≥ rank(X0). It is clear that rank(X0) = m + rank([1s XA]), and it can be shown that rank([1s XA]) = rank(XA) = m when m ≥ 3. Therefore, rank(X) = 2m as m ≥ 3. Note that when m = 2, we have s = 1 and rank(X) = 3.

From the linear models theory, we know that for a vector λ = (λ0, ..., λ2m)T R2m+1, a linear function λTβ* of β* is estimable if and only if <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M168','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M168">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M169','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M169">View MathML</a> is the null space of the design matrix X. It is also known that <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M170','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M170">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M171','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M171">View MathML</a> is a linear space generated by the row vectors of X. Hence, we have <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M172','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M172">View MathML</a>. Note that <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M173','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M173">View MathML</a> due to the linear dependency among the column vectors 1N, W1, W2, ..., Wm in the design matrix X. Therefore, for a vector λ = (λ0, λ1, ..., λ2m)T R2m+1, the linear function λTβ* is estimable if and only if λ c, or equivalently, <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M174','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M174">View MathML</a>. As a result, we know that in model (8) the functions of model parameters <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M175','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M175">View MathML</a> for j k, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M176','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M176">View MathML</a> for j = 1, ..., m are estimable, and the parameters <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M177','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M177">View MathML</a> as j k, l and k l (or in abbreviation, j k, ≠ l) for j = 1, ..., m are also estimable. But the parameters μ* and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M178','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M178">View MathML</a> themselves are not estimable.

B. Estimability of parameters in model (9)

For model (9), we have its design matrix

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M179','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M179">View MathML</a>

where Wj = (wj(g1), ..., wj(gN))T and Vj = (vj(g1), ..., vj(gN))T for j = 1, ..., m. It can be shown that the W and the design matrix X defined in (16) for model (8) have the following relationship W = XT or X = WST, where

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M180','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M180">View MathML</a>

and

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M181','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M181">View MathML</a>

with d = (2, -1, ..., -1)' ∈ Rm. Let β = (μ, α1, ..., αm-1, δ1, ..., δm). Therefore, as (8) and (9) are two equivalent models, we have G = * = WSTβ* = , which yields

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M182','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M182">View MathML</a>

From this relationship, we have <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M183','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M183">View MathML</a>, j = 1, ..., m, which are estimable as shown in Appendix A. Besides, the intercept <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M184','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M184">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M185','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M185">View MathML</a>, k j, j = 1, ..., m - 1, are also estimable.

C. Relationships for fully parameterized two-locus models

(C.1) Relationships between parameters of the fully parameterized two-locus model (13) and the expected genotypic values are

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M186','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M186">View MathML</a>

for j, k = 1, ..., m1 - 1; r, s = 1, ..., m2 - 1 and j k, r s.

(C.2) Relationships between parameters of the fully parameterized two-locus model (14) and the expected genotypic values are

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M187','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M187">View MathML</a>

for j, k = 1, ..., m1 - 1 and r, s = 1, ..., m2 - 1, where the relationships between the parameters of model (14) and model (13) are built based on the equivalency between the two models. The relationships between the parameters of model (14) and the expected genotypic values can then be derived by replacing the parameters of model (13) with the expected genotypic values from the previous established results in (C.1).

(C.3) Relationships between parameters of the fully parameterized two-locus model (15) and the expected genotypic values are

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/82/mathml/M188','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/82/mathml/M188">View MathML</a>

for j, k = 1, ..., m1 - 1 and r, s = 1, ..., m2 - 1, where the relationships between the parameters of model (15) and model (13) are built based on the equivalency between the two models. The relationships between the parameters of model (15) and the expected genotypic values are then derived by replacing the parameters of model (13) with the expected genotypic values from the previous established results in (C.1).

References

  1. Fisher RA: The correlation between relatives on the supposition of Mendelian inheritance.

    Trans Roy Soc Edinburgh 1918, 52:399-433. OpenURL

  2. Cockerham CC: An extension of the concept of partitioning hereditary variance for analysis of covariances among relatives when epistasis is present.

    Genetics 1954, 39:859-882. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  3. Cockerham CC: Estimation of genetic variances. In Statistical Genetics and Plant Breeding Natl Acad Sci Natl Res. Edited by Henson WD, Robinson HF. Council publ No. 982. Washington, D.C.; 1963:53-94. OpenURL

  4. Weir BS, Cockerham C: Two-locus theory in quantitative genetics. In Proceedings of the international conference on quantitative genetics. Edited by Pollack EBT Kempthorne O. Iowa State University Press; 1977:247-269. OpenURL

  5. Kempthorne O: An introduction to Genetic Statistics. New Haven: Iowa State University Press, Ames; 1969. OpenURL

  6. Wang T, Zeng ZB: Models and partition of variance for quantitative trait loci with epistasis and linkage disequilibrium.

    BMC Genetics 2006, 7:Article 9. OpenURL

  7. Hansen TF, Wagner GP: Modeling genetic architecture: a multilinear theory of gene interaction.

    Theor Popul Biol 2001, 59:61-86. PubMed Abstract | Publisher Full Text OpenURL

  8. Alvarez-Castro JM, Carlborg O: A unified model for functional and statistical epistasis and its application in quantitative trait Loci analysis.

    Genetics 2007, 176(2):1151-1167. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  9. Zeng ZB, Wang T, Zou W: Modeling quantitative trait Loci and interpretation of models.

    Genetics 2005, 169(3):1711-1725. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  10. Wang T, Zeng ZB: Contribution of genetic effects to genetic variance components with epistasis and linkage disequilibrium.

    BMC Genetics 2009, 10:Article 52. OpenURL

  11. Yang RC, Alvarez-Castro JM: Functional and statistical genetic effects with multiple alleles.

    Current Topics in Genetics 2008, 3:49-62. OpenURL

  12. Lynch M, Walsh B: Genetics and Analysis of Quantitative Traits. Sunderland, MA: Sinauer; 1998. OpenURL

  13. Zaykin DV, Westfall PH, Young SS, Karnoub MA, Wagner MJ, Ehm MG: Testing association of statistically inferred haplotypes with discrete and continuous traits in samples of unrelated individuals.

    Hum Hered 2002, 53(2):79-91. PubMed Abstract | Publisher Full Text OpenURL

  14. Searle SR: Linear Models. John Wiley & Sons Inc., New York, NY; 1971. OpenURL

  15. Ravishanker N, Dey DK: A First Course in Linear Model Theory. Chapman & Hall, CRC, Boca Raton, Florida; 2002. OpenURL

  16. Van Der Veen JH: Tests of non-allelic interaction and linkage for quantitative characters in generations derived from two diploid pure lines.

    Genetica 1959, 30:201-232. PubMed Abstract | Publisher Full Text OpenURL

  17. Mather K, Jinks JL: Biometrical Genetics. 3rd edition. Landon: Chapman and Hall; 1982. OpenURL

  18. Falconer DS, Mackay TFC: Introduction to Quantitative Genetics. fourth edition. Harlow, UK: Longman; 1996. OpenURL

  19. Hayman BI, Mather KM: The description of genetic interactions in continuous variation.

    Biometrics 1955, 11:69-82. Publisher Full Text OpenURL

  20. Liu T, Johnson JA, Casella G, Wu R: Sequencing complex diseases with HapMap.

    Genetics 2004, 168:503-511. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  21. Searle SR, Casella G, McCulloch CE: Variance Components. John Wiley & Sons, NIC, Hoboken, NJ; 1992. OpenURL

  22. Stokes ME, Davis CS, Koch GG: Categorical Data Analysis using the SAS System. 2nd edition. SAS Institute Inc., Cary, NC; 2001. OpenURL