Introduction

Epistasis plays a fundamental role in the genetic control and evolution of complex traits. However, epistatic effects are hard to detect (Cheverud and Routman, 1995) because an epistatic genetic model potentially contains a large number of model effects. One choice is to use a model selection technique to eliminate the spurious effects so that the number of effects is reduced to a manageable level. Several model selection methods have been developed recently. A maximum likelihood (ML) based stepwise regression method has been developed by Kao et al (1999) for model selection. A one-dimensional search method, also based on ML, was developed by Jannink and Jansen (2001) and Boer et al (2002). More recently, a Bayesian method based on the stochastic search variable selection (SSVS) has been applied to QTL mapping (Oh et al, 2003; Yi et al, 2003). These various selection methods are still open to discussion because the criteria of variable inclusion and exclusion are somewhat subjective (Balding et al, 2002; Broman and Speed, 2002; Sillanpaa and Corander, 2002; Kadane and Lazar, 2004).

Ridge regression represents another class of methods for handling oversaturated models. Whittaker et al (2000) adopted the original ridge regression idea of Hoerl and Kennard (1970) to shrink marker effects proportionally in the context of marker-assisted selection. Gianola et al (2003) claimed that the mixed model analysis of genetic effects is the same as the ridge regression analysis, except that the ridge factors vary across model effects and can be estimated from the data. Xu (2003) found that ridge regression works only if the number of model effects is in the same order as the number of observations. Xu (2003) modified the ridge regression by allowing the ridge factor to vary across different model effects. The difference between Xu (2003) and Gianola et al (2003) is that Xu's (2003) method can estimate the QTL variance using only a single regression coefficient whereas the method of Gianola et al (2003) estimates the QTL variance using a batch of regression coefficients. The modified ridge regression methods turn out to be equivalent to the Bayesian analysis with different model effects taking different prior distributions. The model-selection-free method of Xu (2003) has successfully detected multiple QTL with main effects. Extension of the method to an epistatic effects model has not been explored, although it is straightforward. One concern about the extension is the intensive computing time because the Markov Chain Monte Carlo (MCMC) algorithm requires repeatedly sampling a huge number of model effects. If we incorporate the idea of estimating the parameters of the prior distribution from the data, a kind of empirical Bayesian analysis, the method becomes a penalized ML method (Boer et al, 2002).

In this study, we develop such a penalized likelihood method, with the penalty being a function of the parameters. The method can handle an oversaturated model, with the number of model effects many times larger than the number of observations. The method allows spurious effects to be shrunk towards zero, while QTL with large effects is subject to virtually no shrinkage. Therefore, model selection is no longer needed, and all problems associated with model selection, for example, lack of full exploration of the parameter space of models and intensive computing time, are not of concern.

We use the backcross (BC) design as an example to demonstrate the new method. We also fix the QTL positions at markers so that estimation of QTL positions is irrelevant and we can concentrate on evaluating the performance of the method based on the estimated effects. The method is essentially the multiple marker analysis of Xu (2003) incorporating epistatic effects from the ML perspective.

Theory and methods

Epistatic effect model

Let yi (i=1,…,n) be the phenotypic value of the ith individual in a BC mapping population of size n. The epistatic effect model is

where a0 is the population mean, Zij is a dummy variable indicating the genotype of the jth marker for individual i, aj is the effect of marker j (j=1, …, p), p is the total number of markers on the entire genome, ars is the epistatic effect between markers r and s (r=1, …, p–1; s=r+1, …, p), and ɛj is the residual error with a N(0, σ2) distribution. In a BC population, an individual can take one of two genotypes, heterozygote and homozygote. The dummy variable is defined as Zij=1 for heterozygote and Zij=−1 for homozygote.

Methods of estimating the main effects and interaction effects are the same. For the sake of clarity of notation, we redefine the design matrix and the regression coefficients as follows. Let b0=a0, bj=aj (j=1, …, p), and

where q=p(p+1)/2. Similarly, we define xij=zij (j=1, …, p) and

Model (1) is now rewritten as

We now have a simple model that includes both the main and the interaction effects.

Penalized likelihood function

The penalized likelihood is similar to the posterior distribution of the parameters, with the prior distribution of the parameters serving as the penalty. The difference between the penalized likelihood method and the Bayesian method is that the parameters in the prior distributions are estimated simultaneously along with the parameters of interest. Let θ={b0, b1, …, bq, σ2} be the vector of parameters of interest. The log likelihood function is

where and is the normal density with mean βi and variance σ2.

We now introduce a factor to penalize the large number of model effects. This penalty should be a function of the parameters. The prior density of the parameters in the Bayesian framework is an ideal choice for the penalty factor. Let us introduce the following prior density for each of the parameters. Parameters b0 and σ2 are always included in the model and thus their inclusion should not be penalized. We introduce a normal prior for each of the regression coefficients,

In classical Bayesian regression analysis, μj and are hyperparameters. In the oversaturated model, however, the choice of the parameters in the prior distribution is very important. Therefore, we will estimate these hyperparameters from the data. Our experience shows that a prior distribution should also be assigned to μj and the normal prior given below is necessary

where η>0 serves as a prior sample size for accessing μj. Let be the hyperparameters that are subject to estimation. The logarithm of the prior density is used as the penalty and it has the following form;

The penalized log likelihood is defined as

Parameter estimation

The parameters are estimated by maximizing ψ(θ, ξ) with respect to θ and ξ simultaneously. The solutions are called the penalized maximum likelihood estimates (PMLE) of the parameters. The PMLE of ξ are not of direct interest, but provided estimates of nuisance parameters. We now describe an iterative algorithm to solve the PMLE of the parameters.

The PMLE of the intercept is found by setting

and solving for b0, which is

Setting

and solving for bj, we get

The residual variance is estimated by setting

and solving for σ2, which is

Although the nuisance parameters ξ are not of direct interest, they must be estimated from the data. The PMLE of ξ are all obtained by setting , where ξj is the jth component of the nuisance parameters. The PMLE of the nuisance parameters are

and

Summary of iterations

  1. 1

    Set η>0 and provide initial values for θ and ξ;

  2. 2

    Update b0 using Eq. (9);

  3. 3

    Update bj (j=1, …, q) using Eq. (11);

  4. 4

    Update σ2 using Eq. (13);

  5. 5

    Update μj (j=1, …, q) using Eq. (14);

  6. 6

    Update (j=1, …, q) using Eq. (15);

  7. 7

    Repeat step 2 to step 6 until a certain criterion of convergence is satisfied.

The initial values are suggested for θ, where ȳ and are sample mean and variance of the phenotypic values. The initial values are suggested for ξ.

Statistical test

The proposed penalized likelihood method is intended to include all markers in a single model. Theoretically, no statistical tests are necessary because all marker effects are estimated and none of them are missing. However, investigators may only want to report markers with relatively large effects. Fortunately, the penalized likelihood method often provides extremely small estimates for markers not closely linked to QTL. The signals (estimated effects) of these null markers are almost negligible compared to those of the markers that are linked to QTL (Xu, 2003). It is not hard to pick up ‘significant’ markers just by visually scanning on the genome. The majority of the markers should have no effects and their estimated effects should be close to zero. When we plot the estimated marker effects against the genome location, the effects of the null markers should serve as background. Our simulation experiments show that the background noise is indeed very close to zero, making the signals of markers linked to QTL very clear. How large an estimated marker effect is large enough to warrant a spot in the final list of markers associated with the phenotype? An objective statistical test may be helpful. Unfortunately, a usual likelihood ratio test cannot be performed with the penalized likelihood method because of the overparameterization. Therefore, we propose the following two-stage selection process to screen the markers. All markers with are deemed to have passed the first round of selection. If even if it was significant statistically, it would not be interesting biologically. In a BC population, a QTL of this size would explain less than 10−10% of the phenotypic variance. This selection criterion is already quite stringent because very few spurious markers will survive this selection owing to the enforced stringent penalty (shrinkage). In the second stage of the selection, we are more careful on choosing the criterion of selection. We now modify our epistatic model so that only effects that have passed the first round of selection are included in the model because the dimensionality of such a model is quite small compared to the original oversaturated model. Owing to the dimension of the modified model being small, we can use a regular (unpenalized) ML method to reanalyze the data and perform a likelihood ratio test for each QTL. The estimated QTL effects from the penalized likelihood using the oversaturated model are almost identical to the effects estimated from the likelihood analysis using the modified model that includes only the QTL surviving the first round of selection (see results of simulation).

Let s be the total number of QTL effects that have passed the first round of selection and be the parameters that are subject to the ML analysis for significance test. To test the null hypothesis that H0:bj=0, that is, the jth surviving QTL (passed the first round of selection) is not true, we use the following likelihood ratio test statistic,

(Lander and Botstein, 1989), where is the vector of parameters that excludes bj. As pointed out by Kao et al (1999), the choice of critical value for claiming a significant QTL becomes complicated for multiple QTL tests. For simplicity, we use the usual LODj≥3 as the criterion, where LODj=LRj/(2ln 10)≈LRj/4.61. Application of the permutation test (Churchill and Doerge, 1994) is discussed later.

Simulation studies

We conducted three simulation experiments to evaluate the performance of the method. In the first experiment, we simulated a single genome of 200 cM long with 21 evenly spaced markers, with equal marker distance of 10 cM. We put four main QTL effects and four pairwise interaction effects, all of which overlap with markers. The positions and effects of the simulated QTL are given in Table 1 along with the simulated residual variance. We simulated a BC population with sample size of n=200 for one case and n=500 for the other case. Each case was replicated 100 times to evaluate the accuracy, the precision, and the statistical power for each estimated QTL effect. The total number of QTL effects included in the model is 21(21+1)/2=231. We used the two-stage screening process to select markers and further tested the selected marker using the likelihood ratio tests. For each simulated QTL, we counted the samples in which the LOD statistic had passed 3. The ratio of the number of such samples to the total number of replicates (100 in this case) represented the empirical power for this QTL. When the sample size was small, we noticed that a marker with a simulated QTL effect was not always significant, but a significant LOD occurred in a nearby marker. This reflected the uncertainty of the estimated QTL position. In this case, the simulated QTL was also counted as significant (detected). This is why there is an average estimate of QTL position shown in Table 1. The table shows that the larger sample size does have a higher power than the smaller sample size. QTL with small effects tend to be associated with lower powers. The method can detect the smallest QTL (explaining 2.5% of the phenotypic variance) with 63% power even when n=200.

Table 1 Effect of sample sizes on the results of epistatic QTL analysis (100 replicates)
Full size table

The prior value η=5 was used in the simulation experiment. We also tried η=10 and 30, which had virtually no effect on the result. The convergence criterion was chosen as where θ(t) is the vector of parameter values at the iteration. The convergence was usually very fast, taking only 40–70 iterations to the convergence criterion.

In the second simulation experiment, we doubled the chromosome size (400 cM long) but simulated the same number of QTL with the same positions and effects as those given in the first experiment. The total number of markers was 41, with the total number of marker effects (including all pairwise interactions) being 41(41+1)=861. The sample size was now 300 in the second simulation experiment. We also simulated the chromosome length 200 cM (the same as that in the first experiment) with n=300 for comparison. Again, all QTL resided at marker positions. The objective of this experiment was to evaluate the performance of the new method on a more saturated model. Our prediction was that the method would still have a satisfactory performance, even though the number of model effects was almost three times as large as the sample size. The results are given in Table 2 and consistent with the above prediction.

Table 2 Effect of the number of variables on the results of epistatic QTL analysis (100 replicates)
Full size table

Finally, we simulated a single large chromosome 1800 cM long, covered by 121 evenly spaced markers with a 15 cM per marker interval. The total number of QTL effects included in the model was 121(121+1)/2=7381. We increased the sample size to 600. The number of model effects was about 12 times as large as the number of observations. The simulated parameters (positions and effects of QTL) are given in Table 3 for the main effects and Table 4 for the epistatic effects. We simulated nine main-effect QTL and 13 interacting QTL effects. The sizes of QTL (measured by the proportions of phenotypic variance explained by QTL) varied from 0.5 to 20%. Residual variance was set at 10. The data were analyzed with two methods: a Bayesian method and the penalized likelihood method. The Bayesian method was implemented via the MCMC algorithm and it is a simple extension of Xu (2003) by incorporating epistatic effects into the oversaturated model. The initial values and prior distribution of the parameters for the Bayesian analysis were the same as those given by Xu (2003). The length of the Markov chain was of 20 000 iterations, excluding 4000 iterations for the burn-in period. The chain was trimmed by keeping one observation in every 20 iterations so that the posterior sample size for the post-MCMC analysis was 1000. The penalized likelihood method took about 58 iterations to converge and the total computing time with our (SAS, 1999) program run at PC Dell Optiplex GX 400 was about 14 h. The Bayesian analysis, however, took about 3 weeks in the same computer. The estimated main QTL effects are plotted in Figure 1 and the interaction effects are plotted (3D plot) in Figure 2. For the penalized likelihood method, the number of estimated effects (including both the main and interaction effects) that had passed the first round of selection was 23 (nine main effects and 14 interaction effects). In the second round of selection, all the nine main effects and 11 of the interaction effects surpassed the critical value of the test statistic (Table 5). One interaction effect was partitioned into two significant effects (Table 5) and three interaction effects were not detected (Table 4). The penalized likelihood analysis generated results that are quite consistent with the Bayesian analysis, except that the former missed a few QTL effects. This observation was expected because we gain the fast speed with the cost of failing to detect a few very small QTL effects (all explaining <1% of the phenotypic variance). Among the three undetected epistatic QTL effects, two are the interactions between two closely linked markers, that is, marker pairs 20–21 and 101–105, and one is an effect explaining only 0.5% of the phenotypic variance.

Table 3 Simulated and estimated QTL positions and effects from a single dataset of a large genome
Full size table
Table 4 Simulated and estimated positions and effects of interacting QTL from a single data set of a large genome
Full size table
Figure 1
figure 1

Plot of the estimated main QTL effects against the genome location using the penalized likelihood method.

Full size image
Figure 2
figure 2

Plot (3-D) of the interaction effects against the genome location. The left-hand side of the figure shows the true effects (red) and the right-hand side of the figure shows the estimated effects (blue or green). The blue and green colored prisms represent positive and negative estimates of the interaction effects, respectively.

Full size image
Table 5 Likelihood ratio test for the significance of the estimated QTL effects
Full size table

If we compare results of Tables 3 and 4 (penalized likelihood estimates, last columns) with the results of Table 5 (the reduced likelihood estimates, columns 2 and 5) methods), we can see that the penalized estimates and the reduced likelihood estimates are almost identical (differ only by <10−6). Therefore, the two-stage analysis did not change the original estimates of the genetic effects other than facilitate a way to test the significances of the estimated genetic effects.

Discussion

We used a BC design as an example to demonstrate the penalized likelihood analysis. The method can be directly applied to double haploids (DH) and recombinant inbred lines (RIL). Extension to more complicated designs is possible. For example, to analyze data for an F2 design, we need to partition the genetic effect of a single locus into an additive effect (a) and a dominance effect (d). Correspondingly, the epistatic effect can be partitioned into additive-by-additive (aa), additive-by-dominance (ad), dominance-by-additive (da), and dominance-by-dominance (dd). The dimensionality of the model will be doubled for the main effect QTL (a and d) and quadrupled for the epistatic effect QTL (aa, ad, da, and dd), but the method remains the same. The model appears to be

where Zij={1, 0, −1} and Wij={0, 1, 0} for the three genotypes {QQ, Qq, qq}. One can adopt the notation, as given in equation (2), for the above model and thereafter use the same method to estimate all the model effects. For a four-way cross design, the genetic effect of a single locus can be partitioned into three terms, allelic effect from the male parent (am), allelic effect from the female parent (af), and the dominance effect (δ). Correspondingly, the epistatic effect can be partitioned into 3 × 3=9 terms. The incidence variables (coefficients of the main QTL effects) given by Xu et al (2003) for a four-way cross design may be adopted here. The incidence variables for the epistatis effects simply take the products of the corresponding incidence variables of the main effects. No additional theory and methods are involved other than that the dimension of the model should be expanded. Extension of the method to pedigree data is not obvious and deserves further investigation.

Kao and Zeng (2002) proposed the use of Cockerham's (1954) model to define the epistatic and other model effects. Cockerham's model in its original form only applies to two loci. For interactions involving multiple loci, substantial additional work may be required to construct a Cockerham's model. The key of the Cockerham's model is the orthogonality of the model effects. Cockerham (1954) defined model effects by orthogonal linear contrasts of the original genotypic effects. The linear contrasts may be called the statistical parameters, while the genotypic effects may be called the genetic parameters. The transformation from genetic parameters into statistical parameters has several desirable properties, including additivity of the variance components and ease of parameter estimation. However, the additive effects (defined as the linear contrasts) also contain dominance and epistatic effects defined in the original genotypic scale if the population is in linkage disequilibrium. This, unfortunately, has caused some problem in interpreting the biological meanings of the statistical parameters. Nonetheless, we might use the term Cockerham's model if orthogonal linear contrasts are estimated as the model effects. In that case, Cockerham's model in the context of multiple QTL should be described using the following orthogonal transformations. Let us rewrite model (2) as

which can be further expressed in matrix notation by

where ψi=[ψi1, …, ψiq] and β=[β1, …, βq]T. Model (18) may be interpreted as the Cockerham's model if the following constraints are enforced,

where 0q is a 1 × q vector of zeros and Iq × q is an identity matrix with dimension q. Such a ψi can be found using

where L is a generalized inverse of the Choleskey decomposition (upper triangular matrix) of

ie, L=∑−1/2 and LTL=∑−1. By taking such an orthogonal transformation, we can see that and β=L−1b. Substituting ψi, β0 and β into equation (18), we get

and thus we have recovered model (2). One can directly deal with the Cockerham's model using our penalized likelihood method for estimating β, or estimate the original genetic parameters b̂ with our method and then convert them into Cockerham's statistical parameters using . With the Cockerham's orthogonal transformation, we can see that the total phenotypic variance can be partitioned into independent variance components, as shown below,

The Bayesian method developed by Xu (2003) is a model-selection-free method for multiple QTL mapping. The method is very simple so that it can be easily extended to mapping epistatic effects with very little additional effort. However, because the method is implemented via the MCMC, computing time then becomes a major concern for that extension. Within each iteration, a huge number of parameters need to be updated (sampled) and the posterior sample size should be of the order of tens of thousands. Although the computing times spent on updating parameters in each iterations are almost the same for the proposed penalized method and the Bayesian method, the former only takes a few (<70 usually) iterations to converge, while the latter takes several orders of magnitude of iterations to converge to a stationary distribution. The time-saving factor was the major motivation for developing the current method. In addition, the method facilitates an approximate hypothesis test for each putative QTL effect, while the Bayesian method of Xu (2003) does not.

The penalized likelihood method bears all the shrinkage property of the method of Xu (2003). This explains why both methods can handle an extremely oversaturated linear model. However, the shrinkage factor is even more selective in the penalized method. This is reflected by the additional prior distribution assigned to μj, which is the prior mean for bj and is equal to zero in the Bayesian method. Here, in the penalized analysis, the prior for μj is with η>0, which facilitates a mechanism to estimate μj from the data. It is different from Gianola et al (2003). With an estimated μj away from zero, the shrinkage for large QTL effects is not as stringent as when μj=0. However, when the QTL effects are indeed very small, the estimated μj is closer to zero than bj is, and thus the shrinkage becomes more stringent for smaller bj than for larger bj. This more selectively different shrinkage is one of the major differences of the penalized likelihood method from the Bayesian method of Xu (2003). The subjective parameter η>0 does not have a major influence on the result as long as . We tested a wide range of values for η>0, for example, 10, 30, 50, and so on; the results were all comparable to that when η=5. We understand that if this is of no difference from setting a prior for the regression coefficient, and the result is not as good as that when . Our experience shows that 1≤η≤30 works well for all the simulated data we have examined.

The proposed penalized likelihood method also differs from that of Boer et al (2002), who defined the epistatic effect of each marker as a variance component of the current marker with all other markers (background). Our method actually pinpoints the occurrence of the interaction between specific markers. Is it possible to include higher-order interactions in the epistatic model? Theoretically, it is possible because our model does not have restriction on the number of model effects. Practically, however, there is a concern with the lengthy computing time. Based on the assumption of higher-order interactions being less important than lower-order interactions (Falconer, 1989), the pairwise epistatic genetic model should suffice. The proposed penalized likelihood method also differs from that of Kao and Zeng (2002). Although the latter makes use of the orthogonal property of Cockerham's model, variable selection is still needed.

The penalized likelihood method is similar to the BIC or other information-criterion-based methods for parameter estimation (Akaike, 1973; Schwarz, 1978; Jansen, 1994; Broman and Speed, 2002; Sillanpaa and Corander, 2002). However, these methods were designed mainly for model selection. We adopted a similar idea, but allowed the penalty to be a function of the parameters. This modification appears to be trivial, but it has played an important role in improving the performance of the method. All variables are included in the model and no variable selection is performed, and thus we take no risk of missing any important QTL effects.

We used LOD≥3 as the criterion to declare statistical significance. We found no false positive QTL. LOD=3 is somewhat arbitrary. We then took a permutation approach (Churchill and Doerge, 1994) to find the empirical critical values. These empirical 95% critical values tend to be small (all less than 1.0, data not shown). Based on the empirical critical value, we did find a few false-positive QTL. For the large-genome simulation experiment (the third simulation experiment), three false-positive QTL were found. Considering a total number of 7381 effects included in the model, we think that the false positive rate is extremely low. In the first and second simulation experiments (small genome simulations), we only found eight false-positive QTL. Therefore, the probability of false positive is very low. Our experience indicates that the penalized likelihood analysis usually generates such clear signals of QTL effects that a statistical test may not even be required.

The method was validated using simulated data. When applied to real data for estimating epistatic effects, most of the favorable properties will remain. However, some minor modifications are required before the method is applied to real data. (1) In reality, markers may not be evenly distributed along the genome. Although our method does not depend on the uniformity of marker distribution, very tightly linked markers may cause poor estimates of the marker effects due to high degree of multicollinearity. Therefore, it is recommended to use only one marker from a cluster of markers. (2) Missing marker genotypes may occur in real data analysis. Multiple imputations for the missing marker genotypes (Sen and Churchill, 2001) may be adopted here to simulate the missing genotypes. This requires multiple analyses of the data, each for one imputed data set. In all, 10–20 imputed data sets may suffice (Sen and Churchill, 2001). Alternatively, we may replace the indicator variables for the missing marker genotypes by their conditional expectations calculated with the multipoint method (Rao and Xu, 1998). (3) When two markers are far away, it is possible to insert a virtual marker in the middle of the interval bracketed by the markers. The genotypes of the virtual markers are missing across all individuals. Imputations of the virtual marker genotypes are required to complete the data analysis. (4) In real-data analysis, the expected outcome will be slightly different from what was observed in the simulation study in that the background noise of the plots (see Figure 2) will be larger in the real data analysis than in the simulation study. This is not a deficit of the method; rather, it is due to the polygenic nature of quantitative traits. The high ‘noise’ may not be the true noise but caused by the polygenic effects. Excluding these background polygenic effects from the model, as done in any model selection approach, may be detrimental because the polygenic effects, collectively, may significantly contribute to the residual variance.

Finally, we have paid all our attention to developing the penalized likelihood method for handling oversaturated model. We only evaluated marker effects so that we could focus on the performance of the penalized likelihood method on an oversaturated model rather than digressing to estimating QTL positions. Broman and Speed (2002) took the same approach when they investigated various model selection algorithms on QTL mapping. They fixed QTL at marker positions so that they could concentrate on the main issue of model selection rather than addressing estimation of QTL positions. The natural next step would be to develop a true QTL-mapping method incorporating the penalized likelihood method, in which we would allow QTL to move away from markers positions. This has been carried out in the Bayesian shrinkage analysis of Wang et al (2005) for main effect QTL. Extension to QTL with epistatic effects, making use of this penalized likelihood framework, is underway and will be reported in a subsequent paper.