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:
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
since
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
where
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 . 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 , 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 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.
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 ). We classify all two-locus models according to which of the 387 (or 69) triangulations they belong to.
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,
, the pair
and
are the "winners", i.e.
. 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
and
are the "losers". For example, in Figure , there is a line between (1, 0) and (0, 1) and
f01 +
f10 ≥
f00 +
f11, and also 2
f01 + 2
f10 ≥ 3
f00 +
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 . 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,
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 , 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 . 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 1Classical 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)(β (more ...) |
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 .
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 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 . 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 or any rotation thereof. Thus from the shape we can often infer what type of interaction is the strongest in the data.
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 . 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 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 may still be close. Consider the last row in Figure . 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 ). In the last two observations in the row,
dom × dom is the strongest interaction, whereas the first three also show strong
add × dom interaction.
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
where the coefficients of the model are the coordinates
,
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:
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 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.
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 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 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.
1 and Debbie S Yuster2,3
This article has been 
