Theoretical Basis of the Beavis Effect

STATISTICAL THEORY

Marker analysis: The theory is developed using a backcross (BC) mating design, which provides the simplest genetic model in QTL mapping. The result is then extended to other progeny derived from common mating designs. Let y_j be the phenotypic value of a quantitative trait measured from individual j sampled from a BC population. The linear model describing y_j is yj=μ+xja+εj, (1) where μ is the population mean, a is the additive genetic effect at the locus of interest, and ϵ_j is the residual error, which is assumed to follow a normal distribution, N(0, σ₂). In the BC progeny, there are only two possible genotypes at the QTL, heterozygote and homozygote for the backcross parent allele. The independent variable x_j is an indicator variable defined as x_j = 0 for the homozygote and x_j = 1 for the heterozygote.

Let â be the estimated genetic effect and σa^2 be the variance of the estimate for a QTL. For simplicity, let us assume that the QTL is closely linked to a marker so that σa^2=σ2∑j=1n(xj−x‒)2≈σ2nσx2, (2) where n is the sample size (the number of BC progeny) and σx2=n−1∑j=1n(xj−x‒)2 (3) is the sampling variance of the genotype. Note that n may be replaced by n – 1 for a more precise result.

The variance of x depends on the type of progeny. In a BC population without segregation distortion, half of the individuals will be homozygous and half heterozygous, leading to σx2=14 in the limit of large n. We have implicitly made the assumption that there is no dominance effect. If there is dominance, it will be confounded with the additive effect; i.e., the genetic effect is actually a + d instead of a. Therefore, all subsequent analyses are based on the additive model.

If the positions of all QTL are known a priori (uncensored), the estimate of each QTL effect is approximately unbiased (depending on the method used) so that the distribution of â can be assumed to be approximately normal, i.e., N(a, σa^2 ). The normal distribution is mathematically convenient because it allows the standard statistical machinery for censored/truncated data to be used (Cohen 1991). Of course, if the assumption of normal distribution is violated, the theory developed here is not valid. The assumption of a normal distribution for â holds better if the residual error is normally, independently, and identically distributed. Some traits may be distributed in a discrete manner but analyzed as a quantitative trait, and some traits may have a skewed distribution with unequal variances. These traits may not have a normal residual error, and thus the theory developed in this study should be used with extra caution for such traits.

In the practice of QTL mapping, estimated genetic effects are reported only for significant QTL. Thus, the reported QTL consists of a censored sample so that the distribution of the estimated QTL effects conditional on statistical significance is a truncated normal distribution with a mean and variance different from those of the original normal distribution.

I now proceed to calculate this truncated distribution. Let us assume that a z-test statistic is used for the significance test so that the critical value for the z-test statistic is defined as Z_1–α/2, where α is a controlled type I error rate. A QTL is reported only if |z| > Z_1–α/2, where z = â/σ_â is the z-test statistic. In other words, all the QTL reported satisfy |â| > σ_âZ_1–α/2, i.e., â <–σ_âZ_1–α/2 or â > σ_âZ_1–α/2. The two-tailed test leads to the possibility that even if a is positive, it may be detected as a significantly negative effect due to sampling. Denote the truncated â by â_T so that a^T∼NT(aT,σa^T2) , where N_T represents a truncated normal distribution, a_T = E(â_T) is the expectation, and σa^T2=Var(a^T) is the variance of the truncated normal distribution. The truncated normal distribution is not symmetric. Let us define the two truncation points in the standardized (z-test statistic) scale as ξ1=−Z1−α∕2−aσa^=−Z1−α∕2−nσxaσ (4) and ξ2=Z1−α∕2−aσa^=Z1−α∕2−nσxaσ. (5) Further define ψi=ϕ(ξi)1+Φ(ξ1)−Φ(ξ2)fori=1,2, (6) where φ(ξ) and Φ(ξ) are the standardized normal probability density function (pdf) and cumulative distribution function (cdf), respectively.

Often, when using likelihood methods, the LOD score criterion may be used for QTL detection. The critical value of the LOD score can be converted into the critical value of the z-test statistic, using the following approximate relationship, Z1−α≈χ1,1−α2=4.605LOD=2.146LOD, (7) where χ1,1−α2 is the critical value of the chi-square distribution with 1 d.f. and is also the critical value for the likelihood-ratio test statistic. The relationship between the likelihood-ratio test statistic and LOD can be found in Lynch and Walsh (1998).

According to the standard statistical machinery of truncated normal distributions E(a^T)=a+σnσx(ψ2−ψ1) (8) (Cohen 1991), so that the bias is b=σnσx(ψ2−ψ1). (9) The variance of the truncated sample is Var(a^T)=σ2nσx2[(1+ξ2ψ2−ξ1ψ1)−(ψ2−ψ1)2]. (10) This truncated variance is used later when we discuss the bias in the QTL variance estimate.

The phenotypic variance may be partitioned into the true QTL variance (σG2 ) and the residual variance (σ²). The genetic variance of a QTL in a BC population is σG2=a2∕4 . The corresponding variance in an F₂ population is σG2=a2∕2 . In general, the genetic variance can be expressed as σG2=γa2 , where γ is a constant depending on the filial relationship among segregating progeny; e.g., γ = ½ for BC and γ = ½ for F₂. Therefore, a commonly used approach to estimating the genetic variance is σG2=γa^T2 . However, this estimation is also biased because E(σ^G2)=E(γa^T2)=γ[E2(a^T)+V(a^T)]=γ[(a+σnσx(ψ2−ψ1))2+σ2nσx2[(1+ξ2ψ2−ξ1ψ1)−(ψ2−ψ1)2]]=γ[a2+σ2nσx2[1+ψ2(ξ2+2nσxaσ)−ψ1(ξ1+2nσxaσ)]]=σG2+γσ2nσx2[1+(ξ2′ψ2−ξ1′ψ1)], (11) where ξi′=ξi+nσxaσfori=1,2. (12) Therefore, the bias in the estimated variance is B=E(σ^G2)−σG2=γσ2nσx2[1+(ξ2′ψ2−ξ1′ψ1)]. (13) Note that there are two sources of bias: the contribution of environmental variance to the estimate of a² given by [nσx2]−1γσ2 plus the Beavis effect given by [nσx2]−1γσ2(ξ2′ψ2−ξ1′ψ1) . Thus, this estimator is biased even in the absence of the Beavis effect. The bias due to the first source of error has been investigated by Luo et al. (2003), who also provided a method to correct the bias.

The one-tailed test is a special case of the two-tailed test in that we simply set ξ₁ = – ∞ and ψ₁ = 0. For the two-tailed test, we can further simplify B into B=γ(σ2∕nσx2)[1+(ξ2ψ1−ξ1ψ2)] due to ξ1′=−ξ2 and ξ2′=−ξ1 . We did not present this simplified version of B as a general formula because it would not allow us to extend the results simply to the one-tailed test.

The proportion of phenotypic variance explained by the QTL is defined as hQ2=σG2σG2+σ2=γa2γa2+σ2=γa2σP2, (14) where σP2=γa2+σ2 is the total phenotypic variance among the test progeny. The hQ2 notation is adopted from classical quantitative genetics and is known as the heritability. When used as a proportion of variance contributed by an individual QTL, it is no longer called the heritability. The typical estimator for hQ2 is h^Q2=γa^T2∕σ^P2 . It is hard to find the expectation of a ratio. However, if we assume that the estimation error of the denominator is negligible, we can take E(h^Q2)=E(γa^T2)∕σP2 . This assumption is justified for many quantitative traits because the phenotypic variance is usually accurately estimated. On the basis of this assumption, we get E(h^Q2)≈γa2+BσP2=hQ2+BσP2. (15) Letting σ2=(1−hQ2)σP2 , the bias of the genetic variance can be rewritten as B=γ(1−hQ2)σP2nσx2[1+(ξ2′ψ2−ξ1′ψ1)], (16) which leads to E(h^Q2)≈hQ2+γnσx2[1+(ξ2′ψ2−ξ1′ψ1)](1−hQ2). (17) Therefore, the proportion of the phenotypic variance associated with the QTL also is biased. Note that the term (ξ2′ψ2−ξ1′ψ1) in Equation 17 involves the genetic effect, which can be expressed using (14) as a=σhQ2γ(1−hQ2). (18) If we substitute Equation 18 into Equations 4 and 5, σ cancels out from (ξ2′ψ2−ξ1′ψ1) . Therefore, the bias given in Equation 17 is only a function of hQ2 and the parameters associated with the sampled progeny.

We now extend the results to other types of progeny. For simplicity, we assume that dominance is absent. In an F₂ population, there are three possible genotypes whose genotypic values are defined as a for the homozygote with the “high allele,” 0 for the heterozygote, and –a for the homozygote with the “low allele.” The linear model given in Equation 1 applies here in the F₂ population except that the x variable is now defined as x_j = 1, 0, –1 for the three genotypes, respectively. Without segregation distortion, the variance of x in an F₂ population is σx2=12 . Therefore, the results derived above apply here with σx2=14 in BC progeny replaced by σx2=12 in F₂ progeny.

The method also can be applied to double-haploid (DH) and recombinant inbred line (RIL) populations. In both DH and RIL, the heterozygous genotype is absent. The x variable is defined as x_j = 1 for one homozygote and x_j =–1 for the other homozygote. The two types of homozygote have an equal frequency, and thus σx2=1 in both DH and RIL. The results derived above apply by substituting σx2=14 in the BC by σx2=1 in the DH and RIL.

Note that σx2 and γ are identical for all types of progeny if a marker provides a fully informative genotype for the QTL. Different notation is used because there can be differences when the QTL is not tightly linked to a marker.

Interval mapping: In interval mapping, a QTL may be identified at an intermediate position between markers by inferring the genotype of the QTL from flanking marker information. This will affect both the parameter estimates and the statistical tests of inference. We need to substitute the variance of x by the variance of the estimated x, denoted by σx^2 . This variance depends on the relative position of the QTL within the interval. The explicit expression is tedious, but numerical values may be computed conveniently. For example, in a BC population, there are four possible marker classes, and the estimated x is denoted by x̂_j = E(x_j|m_1j, m_2j) = p_j(1), where m_1j and m_2j indicate the flanking marker genotypes and p_j(1) = Pr(x_j = 1|m_1j, m_2j). Let Pr(m_1j, m_2j) be the joint probability of the flanking marker genotype. These values are given in Table 1, where r₁ and r₂ are the recombination fractions between the QTL and the two markers, respectively, and r is the recombination between the two markers. From this table, we can calculate σx^2 using σx^2=E(x^2)−E2(x^)=∑j=14Pr(m1j,m2j)x^j2−[∑j=14Pr(m1j,m2j)x^j]2. (19) Note that when σx^2 is used in place of σx2 , the assumption of a normal distribution for â may be less valid.

The method used to derive σx^2 implicitly assumes that the least-squares method (Haley and Knott 1992) is used because the expectation of the QTL genotype indicator is calculated conditional on markers only. In maximum-likelihood analysis, the expectation is actually calculated conditional on both markers and the phenotypic values. However, it is extremely difficult to derive the variance using expectation conditional on both markers and phenotypes. Therefore, the way we calculate σx^2 is considered as an approximation.

For a DH population, we define x̂_j in a slightly different way (see Table 2) but still use (19) to calculate σx^2 . For RIL, we use the same table (Table 2), but replace r₁, r₂, and r by c₁, c₂, and c, where c = 2r/(1 + 2r) (Lynch and Walsh 1998).

In an F₂ population, there are nine flanking marker genotypic classes. The definitions of x̂_j's and their probabilities for the nine classes are given in Table 3. The same formula (Equation 19) is used to calculate σx^2 except that the summation is taken over nine categories. For example, if the QTL position is in the middle of a 20-cM interval, we have r₁ = r₂ = 0.0906 and r = 0.1648, leading to σx^2=0.4013 . Figure 1 shows the relationship between σx^2 and the position of the QTL in the 20-cM interval flanked by two fully informative markers. Note, as the QTL position approaches a marker, σx^2 approaches ½.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

—Relationship between and the position of a QTL flanked by two fully informative markers in an interval of 20 cM in an F₂ population in the limit of large n.

NUMERICAL EVALUATION

Bias was evaluated numerically by considering the following factors: the sample size, the genetic effect measured by hQ2 , and the LOD score criterion. In evaluating the bias of QTL effects, the residual variance was set at unity, i.e., σ² = 1.0. If the actual residual variance is not unity, one can always standardize the genetic effect using a/σ in place of a.

The numerical evaluation was conducted only in BC populations because the general trends are similar for all types of progeny (data not shown). In addition, only a one-tailed test was evaluated. The one-tailed test is a special case of the two-tailed test with ξ₁ = – ∞ and ψ₁ = 0. The functional relationships between the size of the detected QTL and the true sizes are shown in Figure 2. The diagonal lines in the first column of Figure 2 represent the case where a_T = a and those in the second column represent the case where h^Q2=hQ2 , which holds only when the sample size is infinitely large. With finite sample sizes, the curves deviate from this straight line and the deviation increases as the sample size decreases. The deviation also increases as the LOD criterion increases. The deviation is negligible when hQ2=0.5 (corresponding to a/σ= 0.2), even if the sample size is as small as n = 50. Assume that the commonly used LOD criterion is 3 and the sample size is 200; the bias is within 7% of the true value as long as hQ2≥0.10 . The bias becomes more severe, however, for small detectable QTL. This observation is consistent with that observed by Beavis (1998). It is interesting to note that when hQ2=0 and n = 50, the bias in h^Q2 is as high as 0.18 for LOD = 2 and 0.33 for LOD = 5.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

—Functional relationships of the detected QTL parameters to the true parametric values under various sample sizes and LOD score critical values in BC mapping. The diagonal lines represent the relationship between the parameter of detected QTL and the true parametric value when the sample size is infinitely large (unbiased relationship). Note that as the LOD criterion increases, the bias associated with detecting a very small QTL (a → 0) becomes worse but the probability of detecting such a QTL becomes much smaller (not shown). Note that in the Beavis experiment, LOD = 2.5 and n = 100, 500, 1000.

The functional relationships between h^Q2 and the sample size under various LOD criteria and hQ2 are shown in the right column of Figure 3. The corresponding plots for the effects are given in the left column. When n ≥ 200 and hQ2 ≥ 0.10, all the biases are negligible (within 10% of the true value), regardless of the LOD criterion. For small hQ2 , a large sample size, even as large as n = 500, is not sufficient to eliminate the bias, again consistent with results of the Beavis experiment (Beavis 1998).

If the average estimated h^Q2 is 0.14 among all the experiments surveyed where the average sample size is ∼100 and the average LOD criterion is ∼3, the true hQ2 may actually be zero (found from the second panel from the top of the right column of Figure 2). If h^Q2 is 0.15 (just a slight increase), however, the true hQ2 is ∼0.08. If h^Q2=0.25 , the true hQ2 is about the same as h^Q2 ; i.e., very small bias is expected. From these graphs or using Equation 17, one can find the true hQ2 retrospectively from h^Q2 for all other settings.

We now use the parameter values of the Beavis experiment to compare the biases with those observed by Beavis (1998). Table 4 shows the original data reported by Beavis as well as the predicted biases for both the QTL effects and their variances. The average QTL position is in the middle of a 20-cM interval. However, for an interval this short, the position of a detected QTL rarely coincided with the true position. We observed that the detected position within the interval varied almost uniformly across the interval. Therefore, we choose an average σx^2 = (0.5 + 0.4)/2 = 0.45 as the estimated variance of x, where 0.5 is the variance when the QTL is completely linked with a marker and 0.4 is the variance when the QTL is exactly in the middle of the interval. The agreements between the observed and the predicted biases, judged by the differences between the two, are quite good except when the sample size is small (n = 100). However, the percentage differences, (observed-predicted)/predicted, show an opposite trend with >8% discrepancies in the variance observed only when n = 1000. The percentage difference may be misleading, however, as the predicted effects become smaller and more sensitive to error with increasing sample size and because, ultimately, one is interested in the absolute error made in inferring the effects of QTL. The large absolute deviations of the predicted biases from the observed values for small sample size may be explained as follows. The critical value of LOD = 2.5 was used by Beavis for a test statistic involving both the additive and dominance effects. Our prediction, however, was based on a test statistic for additive effects only. For small samples, some QTL were detected due to large estimated dominance effects, even though the dominance was absent in the simulations (see Table 10.5 of Beavis 1998). To verify this, we redid the simulations for all the cases with n = 100 and used the pure additive model (dominance was included in neither the simulation nor the model). We used our own interval-mapping program known as QTL-BY-SAS to detect QTL (Xu and Xu 2003). Here, we replicated the experiment 500 times, instead of 200 times. The results are listed in Table 5, which shows very good agreement between the predicted and observed biases. Even though the LOD test statistic of the Beavis experiment involved 2 d.f. (additive and dominance), our model still predicted the biases quite accurately when the sample size was relatively large (n ≥ 500). This was because for large samples, the LOD score test statistic was determined primarily by the additive effect. The 2-d.f. LOD score and 1-d.f. LOD score were virtually identical.