Fast and Accurate Detection of Multiple Quantitative Trait Loci

3.1. The original DIRECT algorithm

In the DIRECT scheme, the search space is successively divided into progressively smaller boxes. The search effort is focused in the most promising regions, and the subdivision of less promising boxes is postponed. For a traditional exhaustive search over a d-dimensional hypercube, the objective function is evaluated in a brute-force fashion at every point of a fine lattice covering the search space, using no further information on the objective function. If DIRECT is run to completion using a minimum resolution criterion matching the step length in the exhaustive search lattice, it will eventually explore exactly the same points as the corresponding exhaustive search. The efficiency of the DIRECT algorithm comes from exploiting the property that the most promising regions are explored first. A heuristic termination criterion is then needed to stop the search well before it devolves into an exhaustive search. As the heuristic criterion has no firm mathematical foundation, this implies that there is no well-established guarantee that the result from DIRECT is equal to the result from an exhaustive search over the same space. This is a typical result for global optimization schemes aimed at solving general problems. To be able to provide guaranteed accuracy without exploring the full search space, more information on the properties of the objective functions must be provided and used in the optimization procedure.

The original DIRECT algorithm initially creates a single search box covering the full search space. The objective function, that is, in our case the RSS, is evaluated in the center of this box. The box is then split into three equally sized boxes along the majoring dimension. This trinary split results in the centroid of the original box coinciding with the centroid of one of the new boxes. Therefore, only two additional function evaluations are required for the three resulting boxes. DIRECT then continues iteratively splitting the boxes. At the end of each DIRECT iteration, the convex hull is determined among the remaining boxes, in a space of box radii vs. objective function values. This hull determines which boxes to split in the following iteration. By computing the hull, the RSS value for the sequence of boxes will be monotonously decreasing if tracing along the box radii. The idea is that a promising box is characterized by either being large, so there is a high possibility that exploring it further might recover a new optimum, or that the value at the centroid is close to the current optimum and hence a new optimum might be found, even if the radius is small. This qualitative argument can also be represented as different assumptions of the value of a Lipschitz constant K, that is, on the maximal slope of the objective function (Jones et al., 1993). Figure 1 illustrates a few iterations of DIRECT in a simple one-dimensional space. The hull is “peeled off” when those boxes are split, making new boxes available in the next iteration. This process is repeated until a suitable termination condition has been reached.

OPEN IN VIEWER

FIG. 1.

To the left, three boxes in the optimization space are illustrated. As all boxes are of equal size, only one will be selected for splitting in the convex hull, resulting in the boxes on the right. If the splitting continues, two boxes would be split, one from each size, as the smallest function value in the smaller box size is slightly lower than the smallest value in the larger one. Dashed lines indicate possible minimum function values found at each distance from the box centroid, assuming a strict Lipschitz bound of K = 0.04.

3.2. Adaption of DIRECT to solve QTL mapping problems

In Ljungberg et al. (2004), the original DIRECT algorithm was adapted for solving QTL search problems. A heuristic termination criterion, as described above, was used to stop the search, and the resulting scheme was proven to be orders of magnitudes faster than the corresponding exhaustive search. This made QTL searches for at least d = 3, 4, and 5 possible using a reasonable computational effort. However, the results in Ljungberg et al. (2004) also show that using DIRECT for the permutation tests can result in a bias in the significance thresholds compared to those from an exhaustive search. The reason for this is that the termination criterion used in Ljungberg et al. (2004) results in premature termination for the very flat optimization landscape present in most permuted cases. The end result is that a putative 95% threshold instead would give, for example, 94.8% significance. It should also be noted that no bias was detected in the DIRECT searches for nonpermuted phenotype data, probably due to the more structured nature of these optimization landscapes. A main result of the analysis presented in later sections is that it is possible to terminate our new PruneDIRECT process at a much earlier point than the corresponding exhaustive search while still being able to guarantee, up to some defined residual probability, that the same global optimum is found as for the exhaustive search.

We argue that the successful results for DIRECT presented in Ljungberg et al. (2004) originate from the fact that the probability of co-inheritance of two genomic loci on a single chromosome is related to their physical distance from each other. The standard unit for genetic distances, mapping distance, is even directly defined from such probabilities. In Ljungberg et al. (2004), this fact is not explicitly used in the algorithm, but since the end result is that the slope of the objective function is limited, the performance of DIRECT ends up being quite good. It is also clear that this notion of co-inheritance cannot be extended across chromosome boundaries. Separate chromosomes segregate independently and can be considered to have an infinite mapping distance. In the search space for the QTL mapping problems, all genome locations are lined up as a single line, one chromosome after another. There is a clear disjunction at the chromosome borders, corresponding to a discontinuity of the RSS at these boundaries. In Ljungberg et al. (2004), this situation is handled by introducing what is essentially several DIRECT searches governed by a common priority queue. In this version of the algorithm, each chromosome combination is considered to correspond to a separate search space, called a chromosome combination box (cc-box). In the first step of the DIRECT scheme, the RSS is evaluated at the centroids of all these cc-boxes.

3.3. The PruneDIRECT algorithm

The original DIRECT algorithm is based on an assumption of Lipschitz continuity, where the value of the Lipschitz constant is unknown. This means that we know there exists a constant K such that no partial first-order derivative of the RSS will exceed K in any position but that the value on this bound of the slope of the objective function is not known.

Our improvement of the DIRECT procedure for QTL mapping is based on the fact that a Lipschitz bound for a QTL problem can be computed based on the relation between co-inheritance and physical distance mentioned above. Parts of the search space can then be permanently excluded from further subdivisions at each iteration in the DIRECT procedure, resulting in an efficient algorithm where a heuristic termination criterion does not have to be used. The resulting scheme has a well-defined error bound and is equivalent to the corresponding exhaustive search.

If a bound on the Lipschitz constant K is known, it is possible to compute upper and lower bounds for the objective function within any box in the search space, given knowledge of the value at the centroid. When using the original DIRECT algorithm, K is unknown, but it is still possible to impose a partial ordering of boxes. If a box A is both smaller and has a larger value of RSS than another box B, then no value of K can result in a smaller minimum bound within A than within B. This partial ordering gives rise to the choice of a convex hull of boxes to split in the next iteration.

However, DIRECT can also be considered in the more general context of Lipschitz optimization schemes. Here, the more traditional schemes (Shubert, 1972; Pinter, 1986; Galerpin, 1985; Mladineo, 1987) assume that a bound of the Lipschitz constant is known. By introducing an upper bound on K in DIRECT, the resulting algorithm is brought close to classic Lipschitz optimization. A known value of the objective function in small boxes will also enable exclusion of larger boxes, hereby introducing a pruning of the search tree, indicating the choice of name for our new algorithm. In contrast to the peeling effect of successive convex hulls in DIRECT, which will eventually result in every existing box getting split, box volumes pruned due to a Lipschitz criterion are permanently irrelevant for further evaluation of the RSS and can be removed from all data structures.

The natural termination condition for the PruneDIRECT algorithm is given by the finite resolution criterion corresponding to the lattice in the underlying exhaustive search. This results in an error bound of Kh on the value of the objective function, where h is the step length in the lattice. Thus, a bound on the Lipschitz constant K leads to the new PruneDIRECT algorithm having improved performance compared to exhaustive search (by excluding parts of the search space) and a well-defined error bound. Since no heuristic termination criterion is used anymore, a main new result is also that PruneDIRECT can be guaranteed to be equivalent to the corresponding exhaustive search, effectively removing any bias of the results caused by the optimization procedure.

For multidimensional QTL mapping that explores the full search space, it is possible to directly use the permutation test methodology for single-QTL models. However, performing hundreds or thousands of multidimensional QTL searches for permuted data in order to get a significance threshold can be a very computationally demanding task, even when an efficient global optimization scheme is employed. Here, it should be noted that since the purpose of doing the random permutations is to determine a significance threshold empirically, it is not necessary to locate the location of the best model fit for every set of randomized data. Rather, it is enough to answer the yes/no question if it is possible to find a location in the permuted dataset with a residual variance below that of the putative set of QTL. If the optimum value used to determine pruning is replaced by the value of the QTL candidate tested, rather than the optimum found in the current permute search, significant decreases in computational effort are possible, even compared to using the same algorithm to solve the full QTL search problem for each permuted dataset. The significance levels derived are identical. We use this approach to derive a very fast scheme for performing permutation tests, described in more detail later.

4.1. Infinite size population

Consider an infinite size population. The total variance is the sum of genetic variance and environmental variance. If there are d QTL, then all genetic variance is explainable by them. We start by analyzing d = 1 and then show how these results can be generalized to higher dimensions. Assume there is a QTL at position x₀. The total variance is the “mean squared error.” We split the total variance into a sum of undiscovered genetic variance and the discovered genetic variance at position x + x₀: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}V_{total} =& V_{genetic} + V_{enviromental} = V_{g - discovered} ( x + x_0 ) + V_{g - undiscovered} ( x + x_0 ) + V_{e} , \\& \hbox{which can be abbreviated into\ } V_{t} = V_{gd} ( x + x_0 ) + V_{gu} ( x + x_0 ) + V_{e}.\end{split} \tag{8}\end{align*} \end{document}

Then define V_r(x + x₀) = V_gu(x + x₀) + V_e as the unexplainable variance at position x + x₀. The goal is to express V_r(x + x₀) as a function of the recombination frequency and then find a bound for it. For simplicity, we start by letting V_e = 0 and later we add V_e back to the calculation.

Assume that the phenotype values for the two QTL genotypes are 0 and 1. Due to the symmetric structure of individuals with phenotype value 0 and individuals with phenotype value 1, V_r(x + x₀) for each class is equal to the total residual variance. Hence: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}V_r ( x + x_0 ) = V_r ( x + x_0 \ \hbox{for class 1} ) =& \hbox{Prob ( no recombination ) } \times ( 1 - \mu ( x + x_0 \ \hbox{for class 1} ) ) ^2 \\& + \hbox{Prob ( recombination ) } \ \ \times ( 0 - \mu ( x + x_0 \ \hbox{for class} 1 ) ) ^2\end{split} \tag{9}\end{align*} \end{document}

where: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} \mu ( x + x_0 \ \hbox{for class 1} ) & = 1 \times \hbox{Prob(No recombination)} + 0 \times \hbox{Prob(recombination)} \\ & = p ( x + x_0 ) .\end{split} \tag{10}\end{align*} \end{document}

Hence, from substituting (10) in (9) and denoting the recombination frequency at position x + x₀ by p(x + x₀) we get: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}V_r ( x + x_0 ) & = ( 1 - p ( x + x_0 ) ) ( 1 - p ( x + x_0 ) ) ^2 + p ( x + x_0 ) ( 0 - p ( x + x_0 ) ) ^2 \\ & = p ( x + x_0 ) - p ( x + x_0 ) ^2.\end{split} \tag{11}\end{align*} \end{document}

This result for V_r should be related to the total phenotypic variance, which is the genetic variance at position x₀ since we defined V_e = 0. We know that all the variance at a QTL point is the discovered variance, hence: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}V_t = V_g ( x_0 ) = V_{gd} ( x_0 ) + v_{gu} ( x_0 ) = [ 0.5 \times ( 1 - 0.5 ) ^2 + 0.5 \times ( 0 - 0.5 ) ^2 ] + 0 = 0.25. \tag{12}\end{align*} \end{document}

At last, we need to have a recombination map, relating genetic distance to the recombination frequency. We use Haldane's mapping function (Haldane, 1919): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}p ( x + x_0 ) = 0.5 + 0.5e^ { - \frac { 2x } { 100 } } , \tag { 13 } \end{align*} \end{document}

where x is the genetic distance from a fixed point x₀ measured in centimorgan. Inserting p(x + x₀) from (13) into equation (11) we get: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}V_r ( x + x_0 ) = & ( 0.5 + 0.5e^ { - \frac { 2x } { 100 } } ) ( 1 - ( 0.5 + 0.5e^ { - \frac {2x} {100}} ) ) \qquad\qquad(14) \\& = 0.25 - 0.25e^ {- \frac {4x} {100}} \qquad\qquad\qquad\qquad\qquad(15) \\ & = V_t - V_g ( x_0 ) e^ { - \frac { 4x } { 100 } } \ ( { \rm from } 12 ) . \qquad\qquad\qquad\qquad(16) \end{align*} \end{document}

Here, V_e can be added to the final calculation since variance is additive. Alternatively, it is safe to assume that V_e is already included in calculating V_t. Now that we have the formula for calculating the explainable variance as a function of genetic distance, we can introduce a transformation of the objective function in the DIRECT optimization. Ignoring a scaling factor, the RSS is equivalent to the residual variance V_r. Adding a constant will not affect the location of minima, so we can instead consider \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$g ( x ) = V_r ( x ) - V_t = - V_g e^ { - \frac { \mid 4x \mid } { 100 } } $$ \end{document} based on (16), assuming x₀ = 0. Furthermore, the function g(x) is always negative, so f (x) = −ln(−g(x)) will always be defined and share the locations of minima with g(x). Hence, we can introduce a transformed objective function; call it LogVar objective function: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}f ( x ) = - \ln ( - V_r ( x ) + V_t ) = - \ln \left( V_g e^ { - \frac { \mid 4x \mid } { 100 } } \right). \tag { 17 } \end{align*} \end{document}

Here, it is easy to verify the derivative of the transformed objective function is bounded: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { df ( x ) } { dx } = \frac { d } { dx } - \ln \left( V_g e^ { - \frac { \mid 4x \mid } { 100 } } \right) = \pm 0.04. \tag { 18 } \end{align*} \end{document}

The use of |x| in the definition of g(x) is related to the fact that x is defined as the distance from the QTL, while a position in the chromosome can naturally be both upstream and downstream from this position. It is possible to shift the function by introducing the true QTL position y, resulting in f (x) = ln V_g− 0.04|x − y|.

The result above can be further generalized, maintaining the bound on K. First, we have assumed that all genetic variance was attributable to a single locus (at y). We can now assume a single-locus model for analysis, but that the true QTL, with respective components of V_g are represented as a vector \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bf y} = \left\{ y_1 , y_2 , \ldots , y_n \right\} $$ \end{document} , resulting in the following expression for the residual variance (assuming all y_i being unlinked): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}V_r ( x ) = V_t - \sum_ { i = 1 } ^n V_ { g_i } e^ { - \frac { 4 \mid x - y_i \mid } { 100 } } . \tag { 19 } \end{align*} \end{document}

If all QTL are indeed unlinked, the positions y_i relative to any point considered in a single-QTL model will be ±∞, except for at most one y_j. Since linkage is transitive, the observation position x can at most be linked to a single QTL. Thus, (19) reduces to: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}V_r ( x ) = V_t - V_ { g_j } e^ { - \frac { 4 \mid x - y_j \mid } { 100 } } , \tag { 20 } \end{align*} \end{document}

and the derivative bound on (20) follows from the result in (18).

The next extension is to make the search landscape itself multidimensional, replacing the scalar x with a vector \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bf x} = \left\{ x_1 , x_2 \ldots , x_m \right\} $$ \end{document} . Assume that each x_i is defined with a point of reference in linkage with the corresponding QTL z_i. Modeled locations not in linkage with any true QTL will result in no explainable variance, and therefore do not need to be considered. As the numbering of both vectors is essentially arbitrary, all other cases are also symmetrical to this one. Using the earlier result, each x_i will only have a term for the corresponding z_i, as all other mapping distances |x_i − z_j| would be +∞. This results in: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}V_r ( { \rm x } ) = V_t - \sum_ { i = 1 } ^n V_ { g_i } e^ { - \frac { 4 \mid x_i - z_i \mid } { 100 } } . \tag { 21 } \end{align*} \end{document}

We then have: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}f ( { \rm x} ) = - \ln ( - V_r ( { \rm x} ) + V_t ) , \tag{22}\end{align*} \end{document}

where the logarithm does not affect the location of minima.

Equation (21) can be rewritten as: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}V_r ( { \rm x } ) = V_t - V_ { g_j } e^ { - \frac { 4 \mid x_i - z_i \mid } { 100 } } - \underbrace { \sum_ { i = 1 , i \ne j } ^n V_ { g_i } e^ { - \frac { 4 \mid x_i - z_i \mid } { 100 } } } _ { C } , \tag { 23 } \end{align*} \end{document}

for any arbitrary \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$j \in {N} , j \leq n$$ \end{document} . Based on (23), (22) becomes: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}f ( { \rm x } ) = - \ln \left( V_ { g_j } e^ { - \frac { 4 \mid x_j - y_j \mid } { 100 } } + C \right). \tag { 24 } \end{align*} \end{document}

From basic calculus, we know that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mid \frac { d } { dx } \ln ( e^ { kx } + C ) \mid = \mid \frac { k e^ { kx } } { e^ { kx } + C } \mid \leq \frac { d } { dx } \ln ( e^ { kx } ) = \frac { k e^kx } { e^kx } = k$$ \end{document} for any \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$C \in {R} , C \geq 0$$ \end{document} . Hence, all partial derivatives \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { \partial f } { \partial x_j } $$ \end{document} are confined by the bound given for the single-dimensional first derivative presented earlier. Thus, the LogVar objective function f (x) as defined above will again have a well-defined Lipschitz bound.

4.2. Explainable variance as a function of genetic distance (finite size population)

The bound derived for an infinite size population does not apply directly to the derivative of the actual residual variance in experimental data, but to the expected value of the derivative, corresponding to the relation between the mapping distance and the expected number of crossover events. Depending on what recombinations are actually present (i.e., in which individuals, with accompanying phenotype values), the actual residual variance can, and will, be different from the relationship predicted. This is an effect of sampling, which should decrease with an increasing size in population and vanish at a theoretical infinite population size. However, experimental populations tend to be rather small, and thus the infinite-size approximation cannot be used directly.

In this section, we derive an approximation to the distribution of the residual sum of errors in the vicinity of any point in the search space for a backcross. We revisit the linear model for a single QTL. Assuming there is a QTL at point x₀ = 0, the model is presented as: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}y_i = \mu ( x ) + a ( x ) Z_{i} ( x ) + e_i ( x ) , { \rm where :\,} Z_{i} ( x ) = \begin{cases}1 , \quad { \rm Qq \ \ at\ \ } x; \\ 0 , \quad { \rm qq \ \ at\ \ } x.\end{cases} \tag{25}\end{align*} \end{document}

Our goal is to calculate the distribution of the residual sum of squares (RSS) at point x. At first, we introduce some definitions: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}n = \hbox {total number of individuals} \tag{26}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}n_0 ( x ) = \hbox{number of 0's at} \ x \tag{27}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}n_1 ( x ) = \hbox {number of 1's at} \ x \tag{28}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}m_{10} ( x ) = \hbox{number of 0's at} \ x \ \hbox {that are 1 at zero} \tag{29}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}m_{01} ( x ) = \hbox{number of 1's at} \ x \ \hbox {that were 0 at zero} \tag{30}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}M_{10} ( x ) = \hbox {set of indices of individuals that recombined from 1 at zero to 0 at} \ x \tag{31}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}M_{01} ( x ) = \hbox{set of indices of individuals that recombined from 0 at zero to 1 at} \ x \tag{32}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}p = \hbox{recombination frequency} \tag{33}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}S_0 = S_1 = \hbox{within - class environmental variance} = 1 \tag{34}\end{align*} \end{document}

The last definitions are simplifying assumptions relating to the normal residual assumption underlying linear regression. Based on the above definitions we get: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}n_0 ( x ) = n_0 ( 0 ) + m_{10} ( x ) - m_{01} ( x ) \tag{35}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}n_1 ( x ) = n_1 ( 0 ) + m_{01} ( x ) - m_{10} ( x ) \tag{36}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\overline { Z } ( x ) = \frac { n_1 ( x ) } { n } \tag { 37 } \end{align*} \end{document}

For now, consider m₀₁(x) and m₁₀(x) to be fixed and known. We also assume that n₀(0) = n₁(0) = n/2 for simplicity. We introduce the estimates for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\hat \mu ( x ) , \ \hat a ( x )$$ \end{document} , and RSS_x. Set \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$a_1 = \sum ( Z_i ( x ) - \overline{Z} ( x ) ) y_i$$ \end{document} (proportional to the covariance between Z and y) and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$a_2 = \sum ( Z_i ( x ) - \overline{Z} ( x ) ) ^2$$ \end{document} (proportional to the variance of Z), then: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\hat a ( x ) = \frac { a_1 } { a_2 } \tag { 38 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\hat \mu ( x ) = \overline{y} - \hat{a} ( x ) \overline{Z} ( x ) \tag{39}\end{align*} \end{document}

After some calculation we get the values for a₁ and a₂: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}a_1 = \sum_{i \in N \backslash N_r ( x ) } ( Z_i ( x ) - \overline{Z} ( x ) ) y_i + \sum_{i \in N_r ( x ) } ( Z_i ( x ) - \overline{Z} ( x ) ) y_i \tag{40}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}= \sum_{i \in N_r ( x ) } ( 1 - 2Z_i ( 0 ) ) y_i + \sum_{i \in N} ( Z_i ( 0 ) - \overline{Z} ( x ) ) y_i \tag{41}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$N = \{ 1 , 2 , \ldots , n \} $$ \end{document} and N_r(x) = {indices for recombined individuals at point x}

As Equation (41) shows, we split a₁ into two sums, one of recombinants and the other nonrecombinant, which can be summarized as: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}a_1 = a_{11} + a_{12} \tag{42}\end{align*} \end{document}

where: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}a_{11} = \sum_{i \in N_r ( x ) } ( 1 - 2Z_i ( 0 ) ) y_i = \sum_{i \in M_{01}}y_i - \sum_{i \in M_{10}}y_i \tag{43}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}a_{12} = \sum_{i \in N} ( Z_i ( 0 ) - \overline{Z} ( x ) ) y_i \tag{44}\end{align*} \end{document}

a₁₁ is a random weighted sum of phonotypes of recombinants and is the only stochastic variable involved in \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\hat{a} ( x )$$ \end{document} . The stochastic behavior of a₁₁ comes from the fact that N_r(x) is a random variable related to the recombination process. Generally speaking, one can say that a₁₁ captures all the stochasticity of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\hat{a} ( x )$$ \end{document} . We approximate the value of a₁₁ by assuming y_i's in a₁₁ to be normally distributed: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\hbox{Assume ( only for recombinant individuals ) }: \ y_{i} \sim N ( \mu ( 0 ) + a ( 0 ) Z_i ( 0 ) , 1 ) \tag{45}\end{align*} \end{document}

then: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}a_{11} \sim N ( \mu_{a_{11}} , \sigma^2_{a_{11}} ) \tag{46}\end{align*} \end{document}

where from (43): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\mu_{a_{11}} = \mu ( 0 ) ( m_{01} ( x ) - m_{10} ( x ) ) - a ( 0 ) m_{10} ( x ) , \tag{47}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\sigma^2_ { a_ { 11 } } = S_0m_ { 01 } ( x ) \left( 1 - \frac { m_ { 01 } ( x ) - 1 } { n_0 ( x ) - 1 } \right) + S_1m_ { 10 } ( x ) \left( 1 - \frac { m_ { 10 } ( x ) - 1 } { n_1 ( x ) - 1 } \right) \tag { 48 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}S_0 = 1 , S_1 = 1 \tag{49}\end{align*} \end{document}

We are using population correction coefficients \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$( 1 - \frac { m_ { 01 } ( x ) - 1 } { n_0 ( x ) - 1 } )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$( 1 - \frac { m_ { 10 } ( x ) - 1 } { n_1 ( x ) - 1 } )$$ \end{document} , since we are sampling from a finite size population. If, for example, all individuals would recombine, all randomness in a₁₁ would also disappear, which this correction reflects.

We also need to simplify a₂. After some calculations, one gets \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}a_2 = n_0 ( x ) \overline{Z} ( x ) ^2 + n_1 ( x ) ( 1 - \overline{Z} ( x ) ) ^2 \tag{50}\end{align*} \end{document}

Now that we have a₁ and a₂, we should calculate RSS. If a dataset is divided into disjoint categories, it is known that one can write RSS_x as: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}RSS_x = \sum ( y_i - \overline{y} ) ^2 - n_0 ( x ) ( \hat \mu_0 ( x ) - \overline{y} ) ^2 - n_1 ( x ) ( \hat \mu_1 ( x ) - \overline{y} ) ^2 \tag{51}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\hat \mu_0 ( x )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\hat \mu_1 ( x )$$ \end{document} are the means of each category. Using this formula we get: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}RSS_x & = RSS_{total} - n_0 ( x ) ( \hat \mu_0 ( x ) - \overline{y} ) ^2 - n_1 ( x ) ( \hat \mu_1 ( x ) + \hat a ( x ) - \overline{y} ) ^2 \\ & = RSS_t - a_2 \hat a ( x ) ^2 & ( 52 ) \end{align*} \end{document}

Given the above description, we calculate the cumulative distribution function (CDF) of RSS_x: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} F_ { RSS_ { x } } ( rss_x ) = & P_ { a_ { 11 } } \left( \hat a ( x ) ^2 \leq \frac { rss_x - RSS_t } { - a_ { 2 } } \right) \\ = &1 - P_ { a_ { 11 } } \left( a_ { 11 } \leq \sqrt { ( RSS_t - rss_x ) a_2 } - a_ { 12 } \right) \\ &+ P_ { a_ { 11 } } \left( a_ { 11 } \leq - \sqrt { ( RSS_t - rss_x ) a_2 } - a_ { 12 } \right) & ( 53 ) \end{align*} \end{document}

Set \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$R = \sqrt{ ( RSS_t - rss_x ) a_2}$$ \end{document} and then normalize the CDF to get the standard normal: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm Set :\ } z = \frac { a_ { 11 } - \mu_ { a_ { 11 } } } { \sigma_ { a_ { 11 } } } , { \rm then: } \tag { 54 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}F_ { RSS_x } ( rss_x \mid m_ { 01 } ( x ) , m_ { 10 } ( x ) , y , Z ( 0 ) , p ) = 1 - \Phi \left( \frac { R - a_ { 12 } - \mu_ { a_ { 11 } } } { \sigma_ { a_ { 11 } } } \right) + \Phi \left( \frac { - R - a_ { 12 } - \mu_ { a_ { 11 } } } { \sigma_ { a_ { 11 } } } \right) \tag { 55 } \end{align*} \end{document}

where Φ is the CDF of standard normal distribution. Given this result, we remove the constrains on m₀₁(x) and m₁₀(x) by summing over all values they can take. Since n₀(0) = n₁(0) = n/2, we have \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$m_{01} ( x ) , m_{01} ( x ) \in \{ 0 , 1 , \ldots , n / 2 \} $$ \end{document} . Hence, one gets: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} F_{RSS_x} &( rss_x \mid y , Z ( 0 ) , p ) = \\ & \sum_{m_{01} ( x ) } [ B ( m_{01} ( x ) , n_0 ( 0 ) , p ) \sum_{m_{10} ( x ) }B ( m_{10} ( x ) , n_1 ( 0 ) , p ) F_{RSS_x} ( rss_x \mid m_{01} ( x ) , m_{10} ( x ) , y , Z ( 0 ) , p ) ] \end{split} \tag{56}\end{align*} \end{document}

Where B is the binomial probability mass function, and the sums are over all values of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$m_{01} ( x ) , m_{01} ( x ) \in \{ 0 , 1 , \ldots , n_{0 / 1} ( 0 ) \} $$ \end{document} .

In this section, we derived the CDF of RSS_x(x + x₀), the residual sum of squares at point x in the vicinity of x₀, conditional on known m₀₁ and m₁₀, by approximating the phenotypes of recombinants with a normal distribution. Then, we removed the condition on m₀₁ and m₁₀ by computing a probability-weighted sum over all values they can take. This sum forms an approximation of the CDF for RSS_x for fixed y_i's. This CDF is then a two-leveled binomially weighted sum of mixture normals. One can find the value of rss_x for which this CDF is arbitrarily close to 1, for example, the 1−ε quantile. We thus have an upper bound for the value of rss_x at location x + x₀, which holds with probability 1−ε.

4.3. Applications

Based on the derivation in Section 3.3, the PruneDIRECT algorithm, and the Lipschitz bound for the LogVar objective function in Equation (17), a method can be devised for accelerating multidimensional QTL searches and especially permutation testing for significance of an existing candidate QTL.

As indicated in the previous section, Equation (56) can be used to compute a quantile for the distribution of the residual sum of squares for a finite-size population, that is, F_RSSx, at a specific distance from a hypothetical QTL explaining all genetic variance. We define a bound for the LogVar objective function based on the quantile of F_RSSx instead of a Lipschitz bound. All boxes are compared against the LogVar transform for an F_RSSx distribution based on the currently found optimum. Pruning of a box is possible if the 1 − ε quantile for the LogVar distribution at distance x lies below the value evaluated in the centroid of a box with Manhattan radius x. If that condition is fulfilled, then the probability that there would be a new optimum, surpassing the current one, within the box is less than ε. The value for ε will then need to be chosen so that the aggregate probability of missing a minimum over all splits in a QTL search (and any permutation testing) is limited to an acceptable level.

In practice, the set of possible box radii appearing in the DIRECT search for a specific dataset is limited by the structure of the marker map. For a specific minimum value, the set of quantile limits can therefore be reused after first calculated, reducing the load of computing the sum of Gaussian-binomial products. Furthermore, the set of all binomial distribution coefficients for a specific distance can be computed in almost linear time with respect to the total number of individuals (see Appendix B). If, instead, each coefficient would be determined individually when computing the sum, the complexity per term would be on the order of O(log n), with a rather larger constant term, even when using an efficient implementation such as the ones in the Boost Project (2012) and R Development Core Team (2011).

Finally, the binomial distribution is rather thin-tailed. Therefore, it is not necessary to sum over all m₀₁, m₁₀, but rather only the subset from, for example, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$np / 2 - 8 \sqrt{np ( 1 - p ) / 2}$$ \end{document} to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$np / 2 + 8 \sqrt{np ( 1 - p ) / 2}$$ \end{document} , still capturing almost all of the probability mass. This reduces the growth rate of the number of terms in the sum from O(n²) to O(n). The resulting bound is also only slightly more conservative, as the quantiles will get shifted upward. Details for calculating the quantiles can be found in Appendix A.

In all, these adaptations and implementation aspects make it possible to use the defined bound for finite-population online within the Prune-DIRECT algorithm, while avoiding it representing a major part of the total computational load. Using LogVar as the objective function rather than the RSS directly is also beneficial even when the bound is not as simple as the Lipschitz bound, since the expectation will still be a matter of (close to) linear behavior with a bounded slope. This fact makes the splitting strategy used within DIRECT more efficient.

For intercrosses and other configurations with more than two states, we are still using the finite-size distribution determined for the backcross as an approximation, with an added condition. When more parameters are added, the computed RSS can fluctuate more. If the RSS and then LogVar value computed in a box is higher than what would be determined by average in any single linear regression with the same number of degrees of freedom, that is, the mean of the χ² distribution, then the LogVar transform of the χ² mean, rather than the actually computed LogVar value, is used in determining the pruning condition. At this level, small disturbances due to the specific phenotype values overtake the general Gaussian assumptions.

5.1. Simulations

The simulations were performed on the Tintin cluster at the UPPMAX computational resource center, running as single threads on nodes with AMD Opteron 6220 CPUs. The code is parallel, but since a very large number of replicates were used, many jobs corresponding to individual replicates were executed simultaneously, where each job executed the serial version of the code.

For each run, first the non permuted QTL model of specified size was fit, allowing all levels of interaction (each genotype–phenotype mean was a free parameter). Then, permutations were created. For exhaustive search, all possible candidate loci sets were explored in a 1 cM lattice. For PruneDIRECT, the minimum resolution was that same lattice, but in addition the bound was used to avoid splitting of some boxes, if it was definite that the values within that box could not improve upon the minimum found in the original main run. The ε value used in PruneDIRECT was 10⁻⁹.

Table 1 presents the specific number of replicates, the simulated broad-sense heritability h², and the number of permutations done for each replicate. Table 2 presents average, minimum and median number of function evaluations for complete sets of main run plus permutations. Timings show that objective function evaluations exceed 90% of the time used in the PruneDIRECT version. The maximum number of function evaluations can exceed the number of an exhaustive search, as the DIRECT implementation we started out from did not natively implement a discretized search grid and thus points coinciding in the discrete lattice could be evaluated multiple times.

Out of 500 simulated 2D QTL, the loci recovered in 25 of them resulted in more than 10 permuted datasets surpassing the simulated QTL, that is, less than 99% significant. These runs were also more time-consuming, since the efficiency from the PruneDIRECT algorithm arises from the optimum being rare. If the QTL is not significant, then more boxes will be similar by random and thus fewer boxes can be pruned. If these are removed from the computed number of function evaluations, the average number of evaluations decreases to 7529 thousand. Acceleration compared to exhaustive search is only possible when the result from the main run rises above the noise floor defined by the quantiles and χ² distributions. The possible accelerations increase rapidly with stronger QTL signals. Finally, we tested simulating a single four-locus network with a total heritability of 0.30. Finding this took 9.0 million function evaluations, which compares favorably to the 963 millions just finding that network, that is, omitting any permutation testing, would require using an exhaustive search.

We also created 10,000 random QTL searches in a backcross templated on this intercross dataset (by fixing allele origin for either parent). Full concordance was achieved in classical one-dimensional QTL searches between the minimum found through exhaustive search and using PruneDIRECT.

A. Quantile calculations

We want to find the value of rss_x such that: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}F_{RSS_x} ( rss_x \mid y , Z ( 0 ) , p ) = 1 - \epsilon \tag{57}\end{align*} \end{document}

One knows that B(x,n,p) for x outside the interval [μ − 8σ, μ + 8σ] is almost 0, where μ and σ are the mean and standard deviations of the binomial distribution considered. By using this, one gets: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} F_{RSS_x} &( rss_x \mid y , Z ( 0 ) , p ) = & ( 58 ) \\ & \sum_{m_{01}^{ \prime} ( x ) } [ B ( m_{01}^{ \prime} ( x ) , n_0 ( 0 ) , p ) \sum_{m_{10}^{ \prime} ( x ) } B ( m_{10}^{ \prime} ( x ) , n_1 ( 0 ) , p ) F_{RSS_x} ( rss_x \mid m_{01}^{ \prime} ( x ) , m_{10}^{ \prime} ( x ) , y , Z ( 0 ) , p ) ] + \delta\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm where :\ } & \delta \approx 0 , & ( 59 ) \\ & m_ { 01 } ^ { \prime } ( x ) , m_ { 10 } ^ { \prime } ( x ) \in \left[ \frac { np } { 2 } - 8 \sqrt { \frac { np ( 1 - p ) } { 2 } } , \frac { np } { 2 } + 8 \sqrt { \frac { np ( 1 - p ) } { 2 } } \,\right] & ( 60 ) \end{align*} \end{document}

An alternate region for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$m_{01}^{ \prime} ( x ) , m_{10}^{ \prime} ( x )$$ \end{document} is to choose \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$m_{10}^{ \prime} ( x ) \in [ np / 2 - 8 \sqrt{np ( 1 - p ) / 2} , np / 2 + 8 \sqrt{np ( 1 - p ) / 2} ]$$ \end{document} and then \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$m_{01}^{ \prime} ( x )$$ \end{document} in a circle around \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$m_{01}^{ \prime} ( x )$$ \end{document} with radius \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$8 \sqrt{np ( 1 - p ) / 2}$$ \end{document} . After all, one should find rss_x such that: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} \sum_{m_{01}^{ \prime} ( x ) } [ B ( m_{01}^{ \prime} ( x ) , n_0 ( 0 ) , p ) \sum_{m_{10}^{ \prime} ( x ) } B ( m_{10}^{ \prime} ( x ) , n_1 ( 0 ) , p ) F_{RSS_x} ( rss_x \mid m_{01}^{ \prime} ( x ) , m_{10}^{ \prime} ( x ) , y , Z ( 0 ) , p ) ] = 1 - \epsilon\end{split} \tag{61}\end{align*} \end{document}

This equation can be solved efficiently numerically by, for example, the bisection method. In our implementation, we only solve it down to a resolution level of 0.04, overestimating the location of the LogVar quantile by at most the equivalent of 1 cM in the inifite-size population case. Furthermore, since we know the quantile location to be monotonous in terms of increasing x, our precalculation scheme for different box radii can select the low and high initial bisection bounds based on previously calculated solutions for lower and higher x, respectively. In all, this can reduce the amortized number of bisection steps for computing the sum to a handful per unique x. For different x, the φ values going into the sum are also unchanged and can thus be stored precalculated, since evaluating the normal CDF is relatively time-consuming.

B. Calculating the full set of binomial coefficients

Most library implementations supporting the binomial distribution compute single values individually. Even if one is using, for example, R (R Development Core Team, 2011) and provides a vector as x when computing B(x,n,p), each step will be performed independently. The computational complexity for these steps is nontrivial and on the order of O(log(x) + log(n)).

If one knows that all binomial coefficients will be used in a sum, it is instead far more efficient to compute them all as part of the same process. The total computational complexity then becomes linear, with only a few multiplications and divisions per iteration, augmented by a few exponentializations to renormalize the carry values in order to avoid falling outside the dynamic range of conventional floating point implementations.