PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of bmcgeneBioMed Centralsearchsubmit a manuscriptregisterthis articleBMC Genetics
 
BMC Genet. 2008; 9: 17.
Published online Feb 19, 2008. doi:  10.1186/1471-2156-9-17
PMCID: PMC2289835

A complete classification of epistatic two-locus models

Abstract

Background

The study of epistasis is of great importance in statistical genetics in fields such as linkage and association analysis and QTL mapping. In an effort to classify the types of epistasis in the case of two biallelic loci Li and Reich listed and described all models in the simplest case of 0/1 penetrance values. However, they left open the problem of finding a classification of two-locus models with continuous penetrance values.

Results

We provide a complete classification of biallelic two-locus models. In addition to solving the classification problem for dichotomous trait disease models, our results apply to any instance where real numbers are assigned to genotypes, and provide a complete framework for studying epistasis in QTL data. Our approach is geometric and we show that there are 387 distinct types of two-locus models, which can be reduced to 69 when symmetry between loci and alleles is accounted for. The model types are defined by 86 circuits, which are linear combinations of genotype values, each of which measures a fundamental unit of interaction.

Conclusion

The circuits provide information on epistasis beyond that contained in the additive × additive, additive × dominance, and dominance × dominance interaction terms. We discuss the connection between our classification and standard epistatic models and demonstrate its utility by analyzing a previously published dataset.

Background

The genetic dissection of complex traits is at the center of current research in human genetics. Complex traits are caused by multiple susceptibility genes and environmental factors, and mounting evidence from both human genetics and model organisms suggests that epistasis (gene × gene interaction) plays an important role [1,2]. Although the need to consider epistasis when mapping complex trait loci has been discussed by several authors [3-6], most statistical methods used in gene mapping, be it case-control association studies, quantitative trait loci (QTL) mapping, or linkage analysis, are based only on measures of marginal effects at individual loci and do not consider epistasis. Due to recent advances in genotyping technology many large case-control genome-wide association studies [7] have recently been completed, and there has been renewed interest in two-locus disease models and two-locus tests for association [8-11]. The application of two-locus models also arises in expression QTL mapping where thousands of gene expression traits are mapped with linkage analysis and it is imperative to study gene interactions [12].

Ideally, a test for epistasis between two loci A and B should test for biological interaction, or non-independence between the effects of locus A and locus B. Loci A and B are considered independent if the effect of the genotype at locus A does not depend on the genotype at locus B. This biologically motivated concept has been formalized in a variety of ways by different communities seeking simple, mathematically convenient definitions. In the statistical genetics literature the term epistasis is typically taken to mean that the effects at loci A and B are not additive (the "effect" of a locus is defined in terms of the statistical model used [13]). For a further discussion on epistasis see [14-16]. Fisher [17] considered a linear model for the contribution of different loci to a quantitative trait and used the term epistasy to describe a departure from additivity. In linkage analysis based on variance component models, a model without epistasis is a model in which all dominance variance components are zero. In case-control association studies of dichotomous traits it is common to use logistic regression, and additivity is measured in the log-odds of disease for a genotype [18]. A new test for epistasis was recently suggested in [10], it tests for departures from linkage disequilibrium (LD) in the cases, which is equivalent to testing for departure from additivity of the log penetrance values (i.e. departure from a multiplicative model for the two-locus penetrances). These tests all test for departure from additivity on a particular scale, but if an additive model is rejected they provide no information on the type of interaction present. Furthermore, it is not clear what the biological meaning of the interaction is.

With each of the nine two-locus genotypes we associate a genotype value. In the case of a dichotomous phenotype the genotype value can e.g. be the penetrance associated with the genotype, the logarithm of the penetrance, or the logarithm of the odds ratio. In the case of a quantitative trait a natural choice for the genotype value is the expected phenotype value of individuals with that genotype (sometimes called measured genotype). We will consider epistasis to be any deviation from additivity of the genotype values. This is consistent with the definition of epistasis given in [13], both for quantitative and dichotomous traits. In this paper we provide a framework within which one can study and classify the types of epistasis possible between two biallelic loci. Our results are based on recent work of [19] who provide a rigorous geometric approach to epistasis in the haploid case. We extend their results to the diploid case, and characterize all possible patterns of physical interactions among the 9 possible genotypes in the two locus case, showing that there are 387 classes of models that fall into 69 symmetry classes. We discuss the meaning of the different types of interaction and show how the interaction pattern can be effectively measured and visualized.

In genetic analysis it is common to test not only for departure from additivity, but also for whether the data fits a particular two-locus model (e.g. recessive or dominant). We discuss the models that are frequently used and show how they relate to the classification given here. In order to study a wider class of two-locus models [20] enumerated all two-locus, two-allele, two-phenotype disease models with penetrance values 0 or 1 for the nine possible phenotypes. There are 512 such models, which can be reduced to a list of 50 models after allowing for symmetry between alleles, loci and affection status. We classify models with continuously varying penetrances, overcoming the difficulty they highlight in their paper, and show that in fact their 50 models fall into 29 of the 69 symmetry classes.

We introduce the mathematical concepts used to derive the 387 classes of two-locus models and demonstrate on a real dataset how the shapes can be used to classify pairs of loci and identify pairs with similar genetic effects. Finally, we consider the two-locus models typically used in human genetics, the 50 models from [20], and some models with epistasis. We show that these models only represent a small fraction of all possible two-locus models.

Results and Discussion

Shapes of two-locus models

A two-locus disease model on two biallelic loci is specified by the genotype values of the 9 two-locus genotypes. We consider two loci with genotypes aa, Aa, and AA, and bb, Bb, and BB, respectively, where A and B are the susceptibility alleles. The genotype values, fij, i, j = 0, 1, 2, are represented by a 3 × 3 table, where i and j refer to the number of disease alleles at loci A and B, respectively:

equation M1

In the case of a dichotomous trait, fij can e.g. be a penetrance, the probability that an individual with genotype ij will get the disease. For a quantitative trait, fij can e.g. be the expected phenotypic value for an individual with genotype ij.

In an additive model, the genotype values can be written as a sum of the effect at each locus, fij = αi + βj, where αi is the effect associated with having i disease alleles at the first locus, and βj is the effect associated with having j disease alleles at the second locus. An epistatic model is any non-additive two-locus model. To study epistasis we consider the interaction space, which is the space of all two-locus models modulo the space spanned by all additive two-locus models. The interaction space is spanned by a set of linear forms in the {fij} called circuits. There is a circuit for each set of 3 collinear points, and for each set of four points in the plane such that no three of the points are collinear, resulting in a total of 86 circuits. The coefficients in the linear form are such that the sum of the points in the circuit, when scaled by these coefficients, is zero. For example, the circuit arising from the points f00, f01, f20, and f12 is

-3f00 + 4f01 + f20 - 2f12,

since

-3·(0, 0) + 4·(0, 1) + (2, 0) - 2·(1, 2) = 0.

Every circuit with four points can be seen as a contrast between two pairs of genotype values and measures a specific deviation from additivity. For example, the above circuit is positive if 4f01 + f20 ≥ 3f00 + 2f12 and negative otherwise. For some circuits this contrast has a simple interpretation, e.g. the circuit arising from f00, f01 and f02 is f00 - 2f01 + f02. It compares the genotype value f01 (for genotype aa/Bb) to the average of the genotypic values f00 and f02 (for genotypes aa/bb and aa/BB), i.e. it measures deviation from additivity at locus B in individuals with genotype aa at locus A.

To more easily interpret the meaning of the circuits we perform a change of coordinates. In quantitative trait genetics the phenotypic value is often decomposed into additive (fa and fb) and dominance (δa and δb) main effects at loci A and B respectively, and four epistatic effects, additive × additive (IAA), additive × dominance (IAD), dominance × additive (IDA), and dominance × dominance (IDD). We will use the same notation here to decompose the genotype values into main and epistatic effects. We write the two-locus model as

equation M2

where

equation M3
fa = -f00 + f02 - f20 + f22
fb = -f00 - f02 + f20 + f22
δa = -f00 + 2f01 - f02 - f20 + 2f21 - f22
δb = -f00 + 2f10 - f20 - f02 + 2f12 - f22
IAA = f00 - f02 - f20 + f22
IAD = f00 - 2f01 + f02 - f20 + 2f21 - f22
IDA = f00 - 2f10 + f20 - f02 + 2f12 - f22
IDD = f00 - 2f01 + f02 - 2f10 + 4f11 - 2f12 + f20 - 2f21 + f22

Note that with this choice the additive effect is scaled so that the contribution is -fa, 0, and fa for genotypes aa, Aa, and AA respectively, and similarly for the second locus. This is a simple linear transformation of the genotype values which can be used both for penetrances and phenotypic means. The space of all two-locus models has dimension 9 and the interaction space has dimension 6. A natural choice of a basis for the interaction space is given by the interaction coordinates (δa, δb, IAA, IAD, IDA, IDD) where δa and δb measure within-locus interaction and IAA, IAD, IDA, and IDD measure between-loci interaction.

A full list of the 86 circuits in the new coordinates is given in Appendix B. Although the circuits are fully specified by the six interaction coordinates they do contain important information on the type of interaction present. The circuits measure interesting contrasts and the new parameterization allows us to interpret them. For example, circuit c30 = 2δa - 2δb measures the difference between the dominance effects, circuit c1 = -2δa + 2IAD measures the difference between the dominance effect at the first locus and the additive × dominance interaction, etc. The sign of a circuit specifies whether the type of epistasis measured by the circuit is positive or negative, and its magnitude measures the degree of interaction. The circuits contain detailed information on the interaction in a model and to fully describe the pattern of interaction we can consider the sign pattern of all 86 circuits, however, this leads to a very large number of categories. For a more useful classification of all two-locus models according to the type of interaction present we consider the triangulation induced by the penetrances. The connection between a triangulation and the circuits will be discussed further below.

The mathematical definition of a triangulation is given in Appendix A but an informal description is provided here. We represent the 9 genotypes by 9 points in the plane on a 3 × 3 grid and the genotypic values by heights above these points. If the values come from an additive model it is possible to fit a plane through the height points. For any non-additive model we consider the surface given by the upper faces of the convex hull of the heights. Intuitively this is the surface formed if we were to drape a piece of stiff cloth on top of the heights and consider its shape. Any departure from additivity in the model becomes apparent in this surface. The triangulation, or shape, of a model is obtained by projecting these upper faces (the "creases" in the surface) onto the xy-plane.

A visual representation of a two-locus model is given in Figure Figure1.1. The data comes from an example that will be discussed further later. The genotype values, relative to the value of aa/BB, are listed in Panel (a). Panel (b) shows the classical visualization of this table, where each line corresponds to one row in the table. In Panel (c) there is a bar-chart of the data, and the corresponding shape is shown in Panel (d). There is clearly epistatic interaction in the model in Figure Figure1,1, as the genotypes aa/bb, aa/Bb, Aa/bb, and AA/BB have much higher means than the remaining 5 genotypes. The shape shows the four planes of the upper convex hull of the heights. It includes a plane through the genotypes Aa/bb, aa/Bb, and AA/BB, which is given by the middle triangle in the picture, and three planes corresponding to the outer three triangles. Although the classical visualization in Panel (b) of Figure Figure11 contains complete information on the relative genotype values it is hard to grasp what types of interactions occur just by glancing at the figure. The bar-chart is a very good visual representation of the 9 values, however, any comparison between two different datasets based on bar-charts would be not only tedious, but hard to define. Some information is lost by considering only the shape of the model, but since it summarizes the epistasis that is present, the shape enables us to easily compare and classify different models.

Figure 1
Example of epistasis. Example of epistasis in QTL data. The data is on chicken growth [24]. (a) The phenotypic means of the two-locus genotypes, (b) a wiggle plot of the data, where each line corresponds to a row in the table, (c) bar plot of the data, ...

We used TOPCOM [21] to compute all possible triangulations, or shapes, and found that there are 387, however, many are equivalent when we account for symmetry. By symmetry we mean i) the interchange of locus 1 and locus 2, or ii) the interchange of two alleles at one or both loci. These same symmetry conditions were used in [20]. After accounting for symmetry, there are 69 shapes (see Figure Figure2).2). We classify all two-locus models according to which of the 387 (or 69) triangulations they belong to.

Figure 2
Symmetry classes. Shapes. The 69 symmetry classes of the shapes of two-locus models.

A sign pattern for the circuits specifies a model shape, but the converse is not true. Thus considering the shape of a model, rather than the sign pattern of the 86 circuits, gives a coarser model classification, but it provides a very useful description of the type of epistasis in the model. A shape contains information about the signs of some of the 86 circuits. Every group of points in a circuit can be triangulated in exactly two ways [22] corresponding to the type of epistasis. If a model shape has a line connecting the points (i1, j1) and (i2, j2) then for some circuit, equation M4, the pair equation M5 and equation M6 are the "winners", i.e. equation M7. Similarly, if there is no line connecting the points (i1, j1) and (i2, j2), and it is not possible to add one without crossing an existing line segment, then there is some circuit such that equation M8 and equation M9 are the "losers". For example, in Figure Figure1,1, there is a line between (1, 0) and (0, 1) and f01 + f10 f00 + f11, and also 2f01 + 2f10 ≥ 3f00 + f22.

Note that the model shape gives information about the types of interaction present in the model, but does not reveal the magnitude of the interaction (for that we need the actual value of the circuits). For generic models we always get a triangulation of the 3 × 3 grid, but for some models the resulting shape is not a triangulation but a subdivision, where not all cells in the shape are 3-sided (this happens e.g. when many of the genotype values are identical). These coarse subdivisions are not counted in our 387 models, however each coarse subdivision is refined by two or more of our models. The model shape provides information that is complementary to that given by the values of the interaction coordinates. Looking at a specific triangulation or subdivision tells us which way some (but not all) of the circuits are triangulated, thus giving information about interaction for that particular model, in particular the triangulation allows us to identify the dominating interactions. Consider e.g. the example in Figure Figure1.1. Although there is additive × additive interaction present (IAA = 210.9), and the circuit c17 = 4IAA is clearly positive, the corresponding line between (0, 0) and (2, 2) is not included. This is because this interaction is dominated by other types of interaction. The two circuits with the largest values are c71 and c72. The first of these contrasts f01 and f22 with f02 and f10 and thus the line between (0, 1) and (2, 2) is included in the triangulation, the second contrasts f10 and f22 with f02 and f21 and thus the line between (1, 0) and (2, 2) is included.

It is useful to have a notion of when two shapes are "close" or "similar". We say that two shapes are adjacent if one can move from one to the other by changing the sign of one of the 86 circuits (note that most sign changes do not result in a move between shapes). Out of 387 shapes, 350 are adjacent to 6 other shapes, 16 are adjacent to 7 other shapes, and 21 are adjacent to 8 other shapes. We define the distance between two shapes as the minimum number of circuit changes that are necessary to get from one to the other. In the set of 387 shapes the maximum distance between two shapes is 9, and around 70% of all pairs of shapes are distance 4 to 6 apart. Two-locus models which fall into adjacent model shapes share many of the same two-locus interactions, and in general the shorter the distance between two shapes, the more similar the genetic effects. For a further discussion on the shapes of genetic models see [23].

Each shape divides the 3 × 3 grid into 2 to 8 triangles (the numbers in each category are 2, 11, 38, 68, 96, 108 and 64 out of 387). Each shape corresponds to a subspace of 9-dimensional space and the volume of this subspace measures how much of the parameter space the shape inhabits. We obtained an estimate of this by generating 1,000,000 random vectors of length 9 and calculating the shape that each of them falls into. The model shape for a 9-vector is conserved under shifting and scaling so it suffices to consider vectors in [0, 1]9 and each of the 9 numbers were drawn uniformly at random from the interval [0, 1]. The fraction of observations that fell into shapes which divide the grid into 2, 3, 4, triangles, etc., was 6.4%, 17.2%, 28.3%, 24.9%, 15.0%, 6.1% and 2.0%, very different from the fraction of shapes in each category, which is 0.5%, 2.8%, 9.8%, 17.6%, 24.9%, 29.9% and 16.5%. Two-locus models where one, or a few, genotype values are larger than the remaining values induce shapes which contain fewer triangles. However, if the genotype values show only slight deviations from falling on a plane (i.e. δa, δb, IAA, IAD, IDA, and IDD are small), the surface is not dominated by a few genotypes and the resulting shape will be more subdivided.

We will further discuss how the shapes can be used to characterize the type of interaction in a dataset using an example on QTL mapping in chicken.

Two-locus models

In this section we study a number of model classes that are often used in genetic analysis, and the shapes that they induce. We show that each of the model classes restricts the analysis to a small subset of all possible two-locus models. Furthermore, because these models are very specific, they limit the types of interaction that can be modeled and only represent a small fraction of the 69 shapes.

A two-locus penetrance model can be defined by specifying single locus penetrance factors, (α0, α1, α2) and (β0, β1, β2), and combining them in one of three ways,

equation M10

The penetrance factors are typically chosen from a recessive (0, 0, α), dominant (0, α, α) or additive (0, α/2, α) model. For an additive two-locus model with additive penetrance factors, the interaction coordinates δa, δb, IAA, IAD, IDA, IDD are all zero and the circuits all vanish. For all additive two-locus models IAA = IAD = IDA = IDD = 0, but δa and δb depend on the penetrance factors. The heterogeneous model is often viewed as an approximation to the additive model because if the same penetrance factors are used, the models give very similar penetrances. However, in terms of the type of interaction that can be modeled, the multiplicative and heterogeneous models are very similar. In Table Table1,1, we list the values of the interaction coordinates for some common multiplicative models. If we consider the corresponding heterogeneous models, the single locus dominance terms, δa and δb, have the same value, as listed in the table, and the interaction terms, IAA, IAD, IDA, and IDD, all have the same absolute value but opposite signs. The shapes induced by these models are shown in Figure Figure3.3. Note that only 8 shapes can be induced and 6 of the 8 shapes are not generic models (they are subdivisions rather than triangulations of the 3 × 3 grid).

Table 1
Classical two-locus models. The table lists the values of the interaction coordinates for multiplicative two-locus models. The parameters are γ = (α2 - α0)(β2 - β0), η1 = (α0 + α2)(β ...
Figure 3
Shapes of two-locus models. Two-locus models. The model shapes for multiplicative and heterogeneous two-locus models.

In [20] a classification of all two-locus disease models with 0/1 penetrance values is given. Although this classification is useful to generate data under various scenarios and to study general properties of two-locus models, it cannot be used to classify observed data. This class of models is much larger than the class of disease models discussed above, yet they only cover a small part of all two-locus models. The 50 models represent only 29 unique subdivisions, and only 10 out of those 29 are among the 69 model shapes, see Figure Figure44.

Figure 4
Subdivisions. 0/1 models. The subdivisions for the 50 Li and Reich 0/1 penetrance models.

In population genetics and in the study of quantitative traits, two-locus models are classified according to the type of epistatic effects. Four commonly studied patterns of epistasis are additive × additive, additive × dominance, dominance × additive and dominance × dominance. In an additive × additive model two double homozygotes, aa/bb and AA/BB, have higher phenotypic mean (or fitness) than expected, but the other two, aa/BB and AA/bb, have lower phenotypic mean than expected. A numeric representation of the four types is given in Figure Figure55 and the corresponding shapes are also shown. If these epistatic effects are added to a fully additive two-locus model, the resulting shape will be the one shown in Figure Figure5.5. However, the epistasis observed in real data is seldom purely of one type and although e.g. dominance × dominance epistasis is present in the data, the resulting shape can be different. If the dominance terms, δa and δb, are non-zero, the resulting shape will be the dominance × dominance shape, with the possible addition of one or both of the horizontal and vertical lines through the middle of the shape (depending on the magnitude of the dominance terms). A model with both additive × dominance and dominance × additive interaction can fall into one of three shapes. If either the additive × dominance or the dominance × additive interaction is much stronger than the other, the corresponding shape will dominate. If the magnitude of both types of interaction is similar, the resulting shape will be the shape shown in Figure Figure11 or any rotation thereof. Thus from the shape we can often infer what type of interaction is the strongest in the data.

Figure 5
Shapes of two-locus models. Epistastic models. The tables list the genotype values associated with four epistatic models and below each table is the shape induced by a model with purely add × add, add × dom, dom × add or dom × ...

Classification of epistatic effects

Model organisms such as yeast, mouse or chicken are frequently used in genetic analysis, and several recent studies have shown that epistatic effects contribute greatly to observed genetic variability. When pairs of interacting loci have been found, using either QTL mapping, linkage analysis, or association analysis, it is of interest to describe the epistasis in the data. If many pairs of interacting loci have been found, it is of interest to identify pairs with similar genetic effects. This classification can be based on finding, for each pair, the model which best fits the data, out of the classical two-locus models. However, many datasets do not fall into any one of these classes (e.g. more than one type of epistasis can be present in the data). Another option is to base the classification on visual inspection, but that can be inaccurate and very time consuming, especially since in most applications the two alleles at a locus are interchangeable, so one would have to consider many rotations of the 3 × 3 data matrix.

We propose classifying observations according to the shape that they induce, and measuring the similarity of the genetic effects observed in two different datasets by the minimum distance between their induced shapes. This allows us to quickly and automatically identify observations with similar genetic effects. Here we only consider the shape of a model for classification but a more robust classification, outside the scope of this paper, could be obtained by testing which circuits are non-zero and considering the shape induced after circuits which are not significant have been set to zero. This would help in reducing mis-classification due to measurement error in the data and in particular this would reveal whether the data comes from an additive (or near additive) model.

In a study of growth traits in chickens, [24] measured various growth and body weight variables on 546 chickens from an F2 cross between two lines, a commercial broiler sire line and a White Leghorn line. The alleles at each locus are labeled with B and L, according to which line they came from. A method for simultaneous mapping of interacting QTLs was used to do a genome-wide analysis of five growth traits which identified 21 QTL pairs with a significant genetic effect. Some of the 21 QTL pairs were associated with more than one growth trait, resulting in 30 combinations of traits and QTL pairs. For each trait and QTL pair the phenotypic means of each of the nine two-locus genotypes were estimated using linear regression (see Table 2 in [24]). They noted that the standard models for epistasis do not adequately describe the types of interaction present in their data, and classified the QTL pairs into groups with similar genetic effect by visual inspection. They identified 4 general classes of models in this dataset, and classified 16 out of the 21 QTL pairs into one of these classes (when a QTL pair was associated with more than one trait the observations from both traits were considered to be in the same class). The classes are H) some of the homozygote/heterozygote combinations are lower than expected, B) the phenotype value associated with the genotype BB/BB is lower than expected, A) the data fits an additive model, by visual inspection, L) there is a set of genotypes with a high value, a set with a low value associated with it, and the value associated with the genotype LL/LL is between the two, and U) the 5 QTL pairs which did not fit into any of the four classes were left unclassified.

We computed the shapes of the 30 observations and found that 23 of the 387 shapes occurred, or 16 out of 69 up to symmetry. The data are shown in Figure Figure6.6. For each observation we show a bar-chart of the phenotype means and the corresponding shape. The point in the upper left corner of the shape corresponds to the genotype BB/BB, and the point in the lower right corner corresponds to LL/LL. Although in most applications one would consider the two alleles at a locus to be interchangeable we do not here, since they come from different chicken lines. To group together observations with similar genetic effects we clustered the shapes based on the pair-wise distances between them, using complete linkage hierarchical clustering. There are four main clusters in the resulting dendrogram (not shown). Under each panel in Figure Figure66 we list which cluster it falls into, and in parentheses we list which group it belongs to according to [24]. The observations are ordered based on the hierarchical clustering with observations in the same cluster listed together and observations within each cluster listed according to the distance between them, as far as possible. For four observations we switched the order of the first and second locus, compared to the order in [24], in order to minimize the distance to the closest observation. Within a cluster, the distance between the shapes in side-by-side panels is typically one but occasionally two. Many of the observed shapes are adjacent to more than one other shape, so two shapes that are not adjacent in Figure Figure66 may still be close. Consider the last row in Figure Figure6.6. In all five panels the values of the genotypes BB/BL, BL/BB and LL/LL dominate the shape, resulting in a central triangular plane. The value at BB/BB varies considerably but does not affect the shape. The shape that each of the observations fall into is, however, affected by the values of BL/LL and LL/BL. When they are relatively high an additional partition is added in the shape. Recall from the previous section that this shape is observed when there is both additive × dominance and dominance × additive interaction in the data. The shapes in the second-to-last row indicate strong dominance × dominance interaction (compare to the shape given in Figure Figure5).5). In the last two observations in the row, dom × dom is the strongest interaction, whereas the first three also show strong add × dom interaction.

Figure 6
Example. Example of observed epistasis. A visual representation of the 30 trait/QTL pairs. The phenotype average for each genotype is given by the heights of the bars, the corresponding shape is also given, and the trait (A-E) and QTL pair (1–21) ...

The visual classification corresponds very well to the classification based on shapes. All observations labeled H fall into clusters 1 and 2 (which are close to each other in the dendrogram) and all observations labeled B fall into clusters 3 and 4. The observations in group A (additive model) fall into two different clusters. An additive model has no shape (one can fit a plane through the points) but due to measurement error in real data this will not be the case. Note that 3 of the 5 observations in group A induce shapes which are very subdivided, as can be expected when there are no genotypes with very high values which dominate the shape. The observations in group U, which were previously unclassified, have now been grouped with the observations they are closest to. Two QTL pairs (4 and 6) were grouped together in category L. The two observations on QTL pair 6 are in cluster 4 and the observation on QTL pair 4 in cluster 1.

The power to detect epistasis depends on the model shape

In the previous sections we have discussed how to visualize and classify interaction, however, the first step in a two-locus analysis is typically to identify pairs of loci with statistically significant interaction. We now ask whether (and how) the power to detect interaction depends on the true model shape. To fully answer that question it is necessary to perform a thorough simulation study which is outside the scope of this paper, but we have performed a preliminary analysis with the goal of comparing the relative power to detect interaction under different model shapes. We considered three different situations: QTL mapping, association analysis using logistic regression, and association analysis using an LD based measure for interaction. We consider the power to detect interacting loci as a function of only the true model shape although the power will also depend on the minor allele frequency at the two loci, the sample size in the study, the number of genotyped markers, and the prevalence of disease/phenotypic mean in the population. In two-locus QTL mapping, the phenotype is typically modeled as a function of the genotype using a linear model. If y is the phenotype, the model is

equation M11

where the coefficients of the model are the coordinates equation M12, fa, fb, δa, δb, IAA, IAD, IDA and IDD, and ε is Gaussian. The x* are dummy variables; xA takes the values -1, 0, and 1 for individuals with genotypes aa, Aa, and AA, respectively, and xAa takes the value 1 for individuals with genotype Aa. The variables xB and xBb are defined similarly. To test for epistasis, the fit of this model is compared to an additive model where IAA = IAD = IDA = IDD = 0. The test statistic for a likelihood ratio test is minus twice the difference between the log-likelihood of the additive and the full model. This is equivalent to testing if the circuits c7 = c8 = c9 = c10 = 0.

In case-control association studies the penetrances of the genotypes are not observed, only the counts of cases and controls that have each genotype. To study the shape of a two-locus disease model we can fit a full two-locus model using logistic regression and obtain the fitted log odds-ratio for each genotype, which can then be used to obtain an estimate of the penetrances. The model is:

equation M13

where the fij are penetrances and the dummy variables x* are defined as above. By using logistic regression the log-odds scale is chosen as the scale of interest, and additivity on that scale corresponds to no interaction. A likelihood ratio test for epistasis compares the fit of the full model to an additive model where IAA = IAD = IDA = IDD = 0. This test is equivalent to testing c7 = c8 = c9 = c10 = 0 where the circuits are obtained by replacing fij with log(fij/(1 - fij)).

Recently [10] proposed a new test to detect unlinked interacting disease loci. They use an LD based interaction measure, I = h00h11 - h01h10, where hij is defined as the penetrance of a haplotype hij (h00 is the haplotype ab, h01 is aB, etc.). The haplotype penetrance depends on the two locus penetrances as well as the allele frequencies. It is easy to show that the interaction measure, I, vanishes if c7 = c8 = c9 = c10 = 0 when the circuits are calculated using the log penetrance values. In other words, this interaction measure tests for multiplicative penetrances.

We generated 50, 000 random vectors of length 9. For the QTL analysis we fixed the population mean of the phenotype, fixed the allele frequencies of A and B, and then normalized each random vector to give the desired population mean. For each vector we generated 10 datasets, each with sample size 300, and fit both the full model and an additive model (note that for all of the models the 6 interaction coordinates are non-zero so the tests all have the same degrees of freedom). We used the average likelihood ratio statistic as an indicator of the power to detect interaction for that particular model. For each random model we then recorded which of the 387 model shapes it fell into and for each shape looked at maximum of the likelihood ratio statistic. In the first panel of Figure Figure77 we show the maximum for each shape. These maxima are highly variable between shapes, indicating that some types of interactions are easier to detect than others. We also observed that there is a strong association between large values of the likelihood ratio test statistic and the number of polygons a shape divides the square into.

Figure 7
Power to detect association. Power to detect association. The plots show the maximum value of the likelihood ratio test statistic observed for randomly generated data from each of the 387 shapes.

We also generated case-control data from an association study. The random vectors were normalized so that they all give the same population prevalence of disease. In the middle panel of Figure Figure77 we have plotted the maximum of the likelihood ratio test statistic as a function of the shape induced by the penetrances. The test measures deviation from additivity on a log-odds scale, so the difference between the different shapes is relatively small. When the shape is calculated based on the log-odds the results are the same as before. Finally, in Panel 3 of Figure Figure77 we plot the maximum absolute value of the interaction measure I. This test measures deviation from additivity on a log scale, yet the results seem to be more similar to the QTL mapping case.

Conclusion

The multitude of terms used to describe gene interactions are a testament not only to the importance of describing and classifying gene interaction, but also to the fact that even in a two-locus model the types of interactions that can and do occur are diverse and difficult to classify. Most examples of gene interactions that are observed in real data do not fall into any one of the categories typically used to describe interactions. Our approach overcomes this limitation and provides a complete classification of all two-locus models with continuous genotypic values into 69 (or 387) classes. The shape of a two-locus model reveals information about the types of gene interaction present and provides a visual representation of epistasis. By comparing an observed shape to the shapes of standard epistatic models we see which type of interaction is strongest in the data. Moreover, the values of the individual circuits listed in Appendix B provide a complete description of the epistasis in a two-locus system. The observed shape can differ from the true underlying model shape due to noise in the data. Rather than assign an observation to a shape based on the observed genotype values, one could test which circuits are significantly different from zero and use only those circuits to obtain the shape.

Two-locus models are frequently used to generate simulated datasets that form the basis for studies of the power of single-locus and two-locus methods. These can then be used e.g. to choose between exhaustive two-locus searches or two-stage two-locus analyses. There are many examples, both for linkage analysis and association analysis, where the results and ensuing recommendations depend on the models, and types of gene interactions, that are considered [8,9]. With our complete classification it is possible to generate data from each model class (while varying parameters such as population prevalence and allele frequencies) and subsequently a more thorough analysis than previously possible can be performed.

As observed in [24] "there are no striking similarities with a Mendelian pattern of digenic epistasis" in the QTL example and we found many types of nontrivial interaction, including models which cannot easily be described using existing models. The fact that our classification is purely mathematical lends it strength, since we can describe all possible models and categorize them according to the relative genotypic values. It can easily be extended to three or more loci. It remains to be seen whether all of the 69 types occur in nature. Our results also provide a formalism for identifying types of epistasis that may play a role in determining genetic variability in populations [25], but we do not address these implications in this paper.

Authors' contributions

IBH and DSY both analysed the data, wrote the software and prepared the manuscript. DSY studied the triangulations of the two-locus model and IBH worked on the connections to biology. IBH also guided the project. Both authors read and approved the manuscript.

Appendix A: Polyhedral subdivisions

Our classification is based on the theory of regular polyhedral subdivisions.

Definition 1 A polyhedral subdivision of a point set A is a decomposition of conv(A), the convex hull of A, into a finite number of bounded polyhedra, such that the union of these polyhedra is conv(A), and the intersection of any two polyhedra is a common face of each (possibly the empty face).

A polyhedral subdivision where all the polyhedra are simplices is called a triangulation. We can construct a regular polyhedral subdivision of a point set A using the following construction: Assign to every point ai in A a 'height', hi. Then lift each point in A to its specified height by forming the new point set

equation M14

Take conv(equation M15), and consider its "upper faces", that is, the faces whose outward-pointing normal vector has its last coordinate positive. Project each upper face onto conv(A), by dropping the final coordinate of each point. In this manner, we obtain a polyhedral subdivision of A. Note that some points of A may not be used in this subdivision.

Remark 2 In the construction of an induced subdivision there is some ambiguity as to the whether to project with the lower or upper faces of conv(equation M16). Both conventions are commonplace. We chose to use the upper faces in order to stay consistent with literature on induced subdivisions and gene epistasis [19].

If the set of heights {hi} is sufficiently generic, then the subdivision induced by the heights will be a triangulation. We will only consider regular subdivisions and triangulations, thus we will use the term "subdivision" to mean "regular polyhedral subdivision", and "triangulation" to mean "regular triangulation". For more on polyhedral geometry see the book [26].

Appendix B: Circuits

c1 = -2δa + 2IAD
c2 = -2δa - 2IDD
c3 = -2δa - 2IAD
c4 = -2δb + 2IDA
c5 = -2δb - 2IDD
c6 = -2δb - 2IDA
c7 = IAA + IAD + IDA + IDD
c8 = IAA - IAD + IDA - IDD
c9 = IAA - IAD - IDA + IDD
c10 = IAA + IAD - IDA - IDD
c11 = -2δa - 2δb + 2IAA - 2IDD
c12 = -2δa - 2δb - 2IAA - 2IDD
c13 = 2IAA + 2IDA
c14 = 2IAA - 2IDA
c15 = 2IAA + 2IAD
c16 = 2IAA - 2IAD
c17 = 4IAA
c18 = -2δa + IAA + IDA - IDD + IAD
c19 = 2δa + IAA + IDA + IDD - IAD
c20 = -2δa + IAA - IDA - IDD - IAD
c21 = 2δa + IAA - IDA + IDD + IAD
c22 = -2δb + IAA + IDA - IDD + IAD
c23 = -2δb - IAA + IDA - IDD - IAD
c24 = -2δb + IAA - IDA - IDD - IAD
c25 = -2δb - IAA - IDA - IDD + IAD
c26 = -2δa + 2IAA
c27 = 2δa + 2IAA
c28 = -2δb + 2IAA
c29 = -2δb - 2IAA
c30 = 2δa - 2δb
c31 = 2δa + 2IAA + 2IDA + 2IDD
c32 = -2δa + 2IAA + 2IDA - 2IDD
c33 = 2δa + 2IAA - 2IDA + 2IAD
c34 = -2δa + 2IAA - 2IDA - 2IAD
c35 = -2δb - 2IAA - 2IDD + 2IAD
c36 = 2δb + 2IAA + 2IDD + 2IAD
c37 = -2δa + 2IAA + 2IDA + 2IAD
c38 = 2δa + 2IAA + 2IDA - 2IAD
c39 = -2δb + 2IAA - 2IDD + 2IAD
c40 = 2δb + 2IAA + 2IDA - 2IAD
c41 = -2δb + 2IAA - 2IDD - 2IAD
c42 = -2δb - 2IAA + 2IDA - 2IAD
c43 = 2δa + 2IAA - 2IDA + 2IDD
c44 = -2δb + 2IAA + 2IDA + 2IAD
c45 = -2δb + 2IAA - 2IDA - 2IAD
c46 = -2δa + 2IAA - 2IDA - 2IDD
c47 = -2δa + 4IAA + 2IAD
c48 = -2δb + 4IAA - 2IDA
c49 = -2δa + 4IAA - 2IAD
c50 = 2δa + 4IAA - 2IAD
c51 = -2δb - 4IAA + 2IDA
c52 = 2δa + 4IAA + 2IAD
c53 = 2δb + 4IAA + 2IDA
c54 = -2δb + 4IAA + 2IDA
c55 = -2δa + 2δb + 2IDA + 2IAD
c56 = 2δa - 2δb + 2IDA + 2IAD
c57 = 2δa - 2δb + 2IDA - 2IAD
c58 = 2δa - 2δb - 2IDA + 2IAD
c59 = 2δa + 2δb + 4IAA + 2IDA - 2IAD
c60 = -2δa - 2δb + 4IAA - 2IDA - 2IAD
c61 = -2δa - 2δb + 4IAA + 2IDA + 2IAD
c62 = -2δa - 2δb - 4IAA + 2IDA - 2IAD
c63 = 2δa - 4δb - 2IAA - 2IDA + 2IAD
c64 = 4δa - 2δb + 2IAA + 2IDA - 2IAD
c65 = -4δa + 2δb + 2IAA + 2IDA + 2IAD
c66 = 2δa - 4δb - 2IAA + 2IDA - 2IAD
c67 = 2δa - 4δb + 2IAA - 2IDA - 2IAD
c68 = 4δa - 2δb + 2IAA - 2IDA + 2IAD
c69 = 2δa - 4δb + 2IAA + 2IDA + 2IAD
c70 = 4δa - 2δb - 2IAA + 2IDA + 2IAD
c71 = 4δa - 2δb + 4IAA + 2IDA - 4IAD
c72 = 4δa - 2δb - 4IAA + 2IDA + 4IAD
c73 = 4δa - 2δb + 4IAA - 2IDA + 4IAD
c74 = 2δa - 4δb - 4IAA + 4IDA - 2IAD
c75 = -2δa + 4δb + 4IAA + 4IDA - 2IAD
c76 = -4δa + 2δb + 4IAA + 2IDA + 4IAD
c77 = 2δa - 4δb + 4IAA + 4IDA + 2IAD
c78 = 2δa - 4δb + 4IAA - 4IDA - 2IAD
c79 = -4δa - 2δb - 2IDA - 4IDD
c80 = -2δa - 4δb - 4IDD + 2IAD
c81 = -4δa - 2δb + 2IDA - 4IDD
c82 = -2δa - 4δb - 4IDD - 2IAD
c83 = -2δa - 2δb + IAA + IDA - 3IDD + IAD
c84 = -2δa - 2δb - IAA + IDA - 3IDD - IAD
c85 = -2δa - 2δb + IAA - IDA - 3IDD - IAD
c86 = -2δa - 2δb - IAA - IDA - 3IDD + IAD

Acknowledgements

The authors thank Bernd Sturmfels and Lior Pachter for helpful discussions and the two anonymous reviewers whose comments greatly improved this paper. I.B.H. was supported by grant 512066 (LSHG-CT-2004) from the European Union FP6 programme.

References

  • Clark A, Wang L. Epistasis in measured genotypes: Drosophila P-element insertions. Genetics. 1997;147:157–63. [PubMed]
  • Storey J, Akey J, Kruglyak L. Multiple locus linkage analysis of genomewide expression in yeast. PLOS Biology. 2005;3 [PMC free article] [PubMed]
  • Carlborg O, Haley C. Epistasis: too often neglected in complex trait studies? Nat Rev Genet. 2004;5 [PubMed]
  • Wade M. Epistasis, complex traits, and mapping genes. Genetica. 2001;112–113:59–69. doi: 10.1023/A:1013316611768. [Cross Ref]
  • Hoh J, Ott J. Mathematical multi-locus approaches to localizing complex human trait genes. Nature Rev Genet. 2003;4:701–709. doi: 10.1038/nrg1155. [PubMed] [Cross Ref]
  • Moore J. A global view of epistasis. Nature Genetics. 2005;37
  • Hirschhorn J, Lohmueller K, Byrne E, Hirschhorn K. A comprehensive review of genetic association studies. Genet Med. 2002;4:45–61. [PubMed]
  • Marchini J, Donnelly P, Cardon L. Genome-wide strategies for detecting multiple loci that influence complex diseases. Nature Genetics. 2005;37:413–417. doi: 10.1038/ng1537. [PubMed] [Cross Ref]
  • Evans D, Marchini J, Morris A, Cardon L. Two-stage two-locus models in genome-wide association. PLOS Genetics. 2006;2:1424–1432. doi: 10.1371/journal.pgen.0020157. [Cross Ref]
  • Zhao J, Jin L, Xiong M. Test for interaction between two unlinked loci. Am J of Hum Genet. 2006;79:831–845. doi: 10.1086/508571. [PubMed] [Cross Ref]
  • Pérez-Enciso M. Multiple association analysis via simulated annealing (MASSA) Bioinformatics. 2006;22:573–580. doi: 10.1093/bioinformatics/btk050. [PubMed] [Cross Ref]
  • Brem R, Kruglyak L. The landscape of genetic complexity acreoss 5,700 gene expression traits in yeast. PNAS. 2005;102:1572–1577. doi: 10.1073/pnas.0408709102. [PubMed] [Cross Ref]
  • Cordell H. Epistasis: what it means, what it doesn't mean, and statistical methods to detect it in humans. Hum Mol Gen. 2002;11:2463–2468. doi: 10.1093/hmg/11.20.2463. [PubMed] [Cross Ref]
  • Wade M. A gene's eye view of epistasis, selection and speciation. J Evol Biol. 2002;15:337–346. doi: 10.1046/j.1420-9101.2002.00413.x. [Cross Ref]
  • Phillips P. The language of gene interaction. Genetics. 1998;149:1167–1171. [PubMed]
  • Cheverud J, Routman J. Epistasis and its contribution to genetic variance components. Genetics. 1995;139:1455–61. [PubMed]
  • Fisher R. The correlations between relatives on the supposition of Mendelian inheritance. Trans R Soc Edinburgh. 1918;52:399–433.
  • Gauderman W. Sample size requirements for association studies of gene-gene interaction. Int J of Epidemiology. 2002;155:478–484. doi: 10.1093/aje/155.5.478. [Cross Ref]
  • Beerenwinkel N, Pachter L, Sturmfels B. Epistasis and shapes of fitness landscapes. Statistica Sinica. 2007;17:1317–1342.
  • Li W, Reich J. A complete enumeration and classification of two-locus disease models. Hum Hered. 2000;50:334–349. doi: 10.1159/000022939. [PubMed] [Cross Ref]
  • Rambau J. TOPCOM: Triangulations of point configurations and oriented matroids. Proc Int Congress of Mathematical Software, ICMS. 2002.
  • Ziegler G. Lectures on Polytopes. Vol. 152. Springer, New York, NY; 1995.
  • Beerenwinkel N, Pachter L, Sturmfels B, Elena S, Lenski R. Analysis of epistatic interactions and fitness landscapes using a new geometric approach. BMC Evolutionary Biology. 2007;7 [PMC free article] [PubMed]
  • Carlborg O, Hocking P, Burt D, Haley C. Simultaneous mapping of epistatic QTL in chickens reveals clusters of QTL pairs with similar genetic effects on growth. Genet Res. 2004;83:197–209. doi: 10.1017/S0016672304006779. [PubMed] [Cross Ref]
  • Turelli M, Barton H. Will population bottlenecks and multilocus epistasis increase additive genetic variance? Evolution. 2006;60:1763–1776. [PubMed]
  • DeLoera JA, Rambau J, Santos F. Triangulations: Applications, Structures, Algorithms Springer to appear.

Articles from BMC Genetics are provided here courtesy of BioMed Central