METHOD

The principle of a genetic algorithm: The implementation of a genetic algorithm for QTL mapping is outlined in Figure 1. The genetic algorithm (ga) is initiated by randomly generating a population, each member of which constitutes a vector of coded parameter values. The parameter values in the vectors are linked together like genes in a linkage group. Each vector is named a ga-chromosome and each cell in the vector is named a ga-gene. We have used ga-chromosomes containing four ga-genes with real values between 0 and 1. The ga-genes are coded representatives for the chromosomal positions of the two QTL we are fitting (two ga-genes are used to represent the chromosome numbers and two are used to represent the respective locations on these chromosomes).

After initiation, the population member that constitutes the best solution is sought. This is done by selection of good individual ga-chromosomes on the basis of their evaluation values from an objective function (in our case the QTL mapping method). During the evaluation process, the real values of the ga-genes on each gachromosome are translated to putative QTL locations. The statistical support for QTL at the locations given by the ga-chromosome is evaluated and each ga-chromosome becomes associated with an evaluation value, named as its fitness. Ga-chromosomes can then be selected and used further in the algorithm on the basis of their fitness. The standard selection operator in PGA Pack is tournament selection, which operates like a tennis tournament, where the fitnesses of the ga-chromosomes are compared in pairs (matches) in several consecutive rounds. The winners of each match (the most fit ga-chromosomes) advance to the next round. The tournament continues until the desired number of gachromosomes are left. This scheme gives more variance in solutions as less fit ga-chromosomes can become selected, while the most fit solution is always selected.

The selected ga-chromosomes are used to generate new ga-chromosomes for the next generation by a process resembling genetic recombination. The ga-chromosomes are paired and the linkage between ga-genes is broken at one or several places along the ga-chromosome. New ga-chromosomes are generated by combining stretches of ga-genes (parameter values) from the old ga-chromosomes. The standard recombination operator in PGAPack is two-point crossover, which is defined as the process where the linkage is broken at two locations on the ga-chromosome before exchanging stretches of ga-genes. The standard probability of a crossover event occurring during generation of a new ga-gene is 0.85. A small amount of new variation is also introduced by mutation of the ga-genes, i.e., the value of the ga-gene is changed from its present value either using a mathematical function or assigning it to a new random value from the range of parameter values. The standard mutation frequency, defined as the probability of applying the mutation operator on a ga-gene, is 0.001. Propagation of the algorithm through sufficient generations is expected to lead to convergence to the optimal solution according to the objective function.

The implementation of the genetic algorithm: The genetic algorithm was used to find the location in the genome with the strongest support for a QTL effect and high robustness was our major objective when implementing the ga. The computational performance was compared to other search methods when robust ga settings had been found. The robustness of the method was controlled by evaluating the ability of the ga to repeatedly find the global extremes in the simulated data set with 18 QTL described below and real data from a wild boar × Large White intercross consisting of 191 individuals (Anderssonet al. 1994).

The ga was parameterized with the chromosome number and the chromosomal location as separate parameters to obtain maximal blending of parameters during the search. A real-valued genetic algorithm was chosen due to simple translation of ga-gene values to QTL locations. The real-valued ga using the standard settings in PGAPack (Levine 1996) had problems with premature convergence when mapping in populations where several pairs of QTL had effects of similar magnitude. To increase robustness, the ga settings were changed to maintain variation in the ga-population until termination. The mutation frequency was increased from 0.001 to 0.4 (0.4 was chosen, but any value from 0.2 to 0.5 works well) to introduce more new variation and a mutation-operator that chooses new ga-gene values in the whole parameter range was chosen to maximize the new variation generated by mutation. The ga was also set to not allow creation of new ga-chromosomes, that are duplicates of already existing ga-chromosomes. With these settings, the ga was reasonably robust, but convergence at local extremes was sometimes observed when the data contained several extremes of similar magnitude. The PGAPack manual indicated that practitioners often prefer to have recombination and mutation as exclusive operators. This means that if recombination has occurred during generation of a new ga-chromosome, mutation is not to occur, and if no recombination has occurred, mutation is to be performed with a certain probability for each ga-gene. Implementation of this procedure led to a slightly better performance.

Splitting of the ga-population into several smaller gapopulations has been suggested as a means to increase robustness without excessive increase in computational demand. During optimization, each ga-population is in turn allowed to converge at an extreme that is not identified by any previous populations. The detected extremes are then compared to identify the global extreme. The use of multiple populations is implemented by repetitive calls of PGAPack in a DO-loop structure. Three populations of size 50, 10 populations of size 20, and 20 populations of size 10 were compared and the best results were for 10 populations of size 20. With these settings, the ga-based search gave good performance in all our examined populations. In most cases, the ga found the global extreme in <5 populations, but occasionally all 10 populations were needed. A decrease in the number of populations should therefore be avoided unless the major limitation is computational performance.

Genetic algorithms are good at finding regions where the extremes are located but can have difficulty finding the precise location. One way to overcome this is to perform a local exhaustive search in the vicinity (±5 cM) of the extreme found by each ga-population. We noted a significant improvement in the location estimates after implementing this.

The objective function used was the residual sum of squared errors from a weighted least-squares approach to QTL mapping. The method is the extension of the method of Jansen (1992) to the two-loci linear model G = m + A₁ + A₂ + D₁ + D₂ + AA₁₂ + AD₁₂ + AD₂₁ + DD₂₂ as indicated by the author. The parameters of the model are explained below. Markers have not been used as cofactors and successive iterations in the EM algorithm have been removed to increase the computational efficiency during the evaluation procedure. The modifications needed for the single QTL mapping procedure described by Jansen and Stam (1994) when implementing the two-QTL model included duplication of each individual nine times (instead of three times, i.e., once for every possible two-QTL genotype) and the use of an expanded design matrix (X). The design matrix for the two-locus linear model has been described by Jana (1971). The weight for each observation was taken to be the product of the conditional probabilities of the single-QTL genotypes given the markers (Haley and Knott 1992) at each of the two fitted QTL. The estimates of the model parameters can be found as β^=(XTWX)−1XTWYσ2=(1∕N)(Y−Xβ^)TW(Y−Xβ^), where Y is the complete data vector, X is the design matrix for the complete data, W is the diagonal matrix of weights, β is the vector of the regression parameters, σ² is the normal variance, and N is the number of individuals (Jansen and Stam 1994).

The residual sums of squared errors can then be calculated as SSE=(Y−Xβ^)TW(Y−Xβ^).

The method described above can easily be extended to take account of background QTL in the analysis. Two extra ga-genes are added to the genetic algorithm and two extra columns are added to the X matrix for each background QTL. The extra ga-genes represent the chromosomal location for the QTL and the columns in the design matrix are to contain the QTL indicator variables a and d (Haley and Knott 1992) for a QTL at the location given by the ga-genes. The rest of the evaluation procedure is the same as before. We have evaluated the increase in computational demand for a simultaneous search for more than two QTL using this method, but have not investigated any other properties.

Simulations: To evaluate the robustness of the genetic algorithm and to show the advantages of mapping multiple QTL simultaneously, data were simulated for five genetic models (described below) with 100 replicates of 520 F₂ individuals from a cross between two inbred lines. The simulated genome consisted of 20 chromosomes, each 100 cM in length, carrying fully informative marker loci at the ends and at 10-cM intervals. Crossovers were generated using Haldane's mapping function without interference (Haldane 1919).

Four different combinations of interacting and additive QTL were evaluated. All contained two epistatically interacting QTL and a varying number of fully additive QTL, which were all fixed in the parental lines and simulated at random locations in the genome. Only one QTL was allowed on each chromosome. The number of additive QTL were in turn 16, 8, 4, and 0. The additive effects for the combination with 16 fully additive QTL were 0.4 for QTL 1–4, 0.3 for QTL 5–8, 0.2 for QTL 9–12, and 0.1 for QTL 13–16. These effects correspond to the genetic variances 0.08, 0.045, 0.020, and 0.005, respectively. The QTL effects for the models with 8 and 4 additive QTL were similar, but with 2 or 1 QTL in each of the four QTL effect classes. The genetic variance explained by the pair of interacting QTL is given in Table 2. The narrow-sense heritability was set to 0.5. The simulated effects of the QTL were chosen such that mapping of the QTL using an exhaustive enumerative search would detect the pair of interacting QTL in >80% of the simulated populations.

To evaluate the sensitivity of the search methods to differences in detection power of the statistical mapping methods, lower QTL effects were simulated for the recessive epistatic model by using three different heritabilities, 0.5, 0.1, and 0.01. The different heritabilities were obtained by increasing the residual variance relative to the QTL effects. The detection powers for these heritabilities were 80, 70, and 50%, using the exhaustive enumerative search.

Each simulated population was analyzed using two different methods. Method 1 was simultaneous mapping of two interacting QTL using a genetic algorithm. Method 2 was mapping of two interacting QTL by first performing a one-dimensional search, using the method described by Haley and Knott (1992), to identify the most significant single QTL and then subsequently searching for a second QTL showing epistatic interactions with the previously detected QTL. This latter procedure is from here on called a conditional QTL search. The second search is in one dimension and uses the two-QTL epistatic model, where the location of the first QTL is fixed at the location of the most significant single QTL. Both methods used the previously described weighted least-squares method. Detection of the epistatic QTL was declared upon identification of two QTL on the chromosomes where the simulated major QTL were located, i.e., no significance thresholds have been used at this stage. The comparison between the ga-based and the conditional search is not fair, since the conditional search is model dependent in the sense that it relies on the contribution of epistasis to additive and dominance variances for detection of the first QTL. Thus, it does not guarantee giving the right answer when searching for two QTL, while the ga does. This comparison has been included anyway to show how the more theoretically correct genetic algorithm compares to the conditional search that is frequently used in practice.

Epistatic model: The epistatic model used in the regression was the two-locus linear model G = m + A₁ + A₂ + D₁ + D₂ + AA₁₂ + AD₁₂ + AD₂₁ + DD₂₂, where A₁ is the additive effect at locus 1, A₂ is the additive effect at locus 2, D₁ is the dominance effect at locus 1, D₂ is the dominance effect at locus 2, AA₁₂ is the interaction between A₁ and A₂, AD₁₂ is the interaction between A₁ and D₂, AD₂₁ is the interaction between A₂ and D₁, and DD₁₂ is the interaction between D₁ and D₂ (Jana 1971). The epistatic interactions used for the simulation study were the classical models of epistasis: complementary, dominant, duplicate, recessive, and inhibitory epistasis (Jana 1971). Complementary epistasis is observed when a defect in either of two genes gives the same mutant phenotype, giving an expected Mendelian segregation ratio of 9:7 (Table 1). In this case, functional copies of both genes must be present to produce the dominant phenotype. Duplicate epistasis is observed when a defect in two genes gives a mutant phenotype and the expected segregation ratio is 15:1. In this case, a functional copy of only one of the two genes must be present to produce the dominant phenotype. Dominant, recessive, and inhibitory epistasis occurs when one gene blocks the phenotypic expression of a second gene. For dominant epistasis, the dominant allele at the first locus is also dominant over the alleles at the second locus. The phenotypic effects of the second locus are therefore expressed only when the individual is recessive homozygote at the first locus. This gives an expected segregation ratio of 12:3:1. Recessive epistasis occurs when the recessive homozygote at one locus is dominant over the alleles at the other locus. The phenotypic effects for the second locus are therefore expressed only when the individual is dominant homozygous or heterozygous at the first locus. The expected segregation ratio is in this case 9:3:4. Inhibitory epistasis works in the same way as dominant epistasis and is the special case when the two genes have equal-sized effects with opposite signs. The expected segregation ratio is here 13:3. The relationships among the genetic parameters for these five genetic models are given in Table 1. The translation of the genetic parameters to the genotypic effects of the two interacting QTL are given in Figure 2. The genotypic effects and genotypic variances for the two interacting QTL for the five epistatic models are given in Table 2.

Computational efficiency: The computational efficiency was compared by calculating the number of statistical evaluations performed during QTL mapping. Three search methods were compared: an exhaustive enumerative search, a genetic algorithm-based search, and a conditional search. Comparisons were made for mapping two, three, and four QTL in a genome size of 2000 cM using 1-cM resolution.