A Bayesian approach for constructing genetic maps when markers are miscoded

The advent of molecular markers has created opportunities for a better understanding of quantitative inheritance and for developing novel strategies for genetic improvement of agricultural species, using information on quantitative trait loci (QTL). A QTL analysis relies on accurate genetic marker maps. At present, most statistical methods used for map construction ignore the fact that molecular data may be read with error. Often, however, there is ambiguity about some marker genotypes. A Bayesian MCMC approach for inferences about a genetic marker map when random miscoding of genotypes occurs is presented, and simulated and real data sets are analyzed. The results suggest that unless there is strong reason to believe that genotypes are ascertained without error, the proposed approach provides more reliable inference on the genetic map.


INTRODUCTION
The advent of molecular markers has created opportunities for a better understanding of quantitative inheritance and for developing novel strategies for genetic improvement in agriculture. For example, the location and the effects of quantitative trait loci (QTL) can be inferred by combining information from marker genotypes and phenotypic scores of individuals in a population in linkage disequilibrium, such as in experiments with line crosses, e.g., using backcross or F2 progenies. A QTL analysis relies on the availability of accurate estimates of the genetic marker map, which includes information on the order and on genetic distances between marker loci order. Genetic maps are inferred from recombination events between markers, which are genotyped for each individual. Several statistical methods have been suggested for map construction. Lathrop et al. [14], Ott [17] and Smith and Stephens [21] discussed maximum likelihood procedures for marker map inferences, and George et al. [9] presented a Bayesian approach for ordering gene markers. Jones [10] reviewed a variety of statistical methods for gene mapping. At present, most statistical methods used for map construction ignore the possibility that molecular (marker) data may be read with error. Often, however, there is ambiguity about genotypes and, if ignored, this can adversely affect inferences [3,15]. The problem of miscoded genotypes has received the attention of some investigators. Most of their research, however, has focused on error detection and data cleaning [4,11,15]. The objective of our work is to discuss possible biases in marker map estimates when miscoding of genotypes is ignored and to suggest a robust approach for more realistic inferences about marker positions and their distances. The approach simultaneously estimates the genotyping error rate and corrects for possible miscoded genotypes, while making inferences on the order and distances between genetic markers.
The plan of the paper is as follows. In Section 2, the problem of miscoding genotypes is discussed, as well as the systematic bias that this imposes on genetic map estimation. In Section 3, a Bayesian approach for inferences about a genetic map, when miscoding is ignored, is reviewed. In Section 4, the methodology is extended to handle situations with miscoded genotypes, when these occur at random. Simulated and real data are analyzed in Sections 5 and 6, respectively, and the results are discussed. Concluding remarks are presented in Section 7.

THE PROBLEM CAUSED BY MISCODED GENOTYPES
First consider the estimation of the genetic distance between two marker loci having a recombination rate θ. In simple situations, e.g., with double haploid or backcross designs, each individual has one of two possible genotypes (say 0 or 1) at each marker locus. Inferences about genetic distance between loci are based on recombination events, which are observed by genotyping individuals. If marker genotypes could be read without error, the probability of observing a recombination event in a randomly drawn individual would be θ. However, it will be supposed that there is ambiguity in the assignment of genotypes to individuals. For example, a genotype 0 may be coded as 1 (or vice-versa), with probability π. Here, given the genotype for a specific marker and the probability of miscoding (π), the distribution of the observed genotypes can be written as: where m ij and g ij are the observed and true genotypes (m ij , g ij = 0, 1), respectively, for locus j (j = 1, 2) of individual i (i = 1, 2, . . . , n). If a "recombination event" between the loci is observed, this may be due to either a true genetic recombination between them, or to an artifact caused by miscoding. Hereinafter, a "recombination" observed by genotyping the markers will be denoted as the "apparent recombination", to distinguish between observed and "true" recombination events.
The probability of observing an apparent recombination between markers 1 and 2 for individual i can be written as: where s i = |m i1 − m i2 | and r i = |g i1 − g i2 | stand for apparent and real recombination events, respectively; and Pr[r It is easy to realize, therefore, that recombination rates estimated from recombinations observed by genotyping the marker loci, ignoring the possibility of miscoding, would be biased upwards whenever the markers are linked (θ < 0.5) and π > 0. Figure 1 shows the expected apparent recombination rates as function of π, for some selected recombination rate values. It seems that the smaller the genetic recombination rate, the worse the relative bias produced by miscoded genotypes. The variance of the apparent recombination event is equal to: Thus, the variance of apparent recombination events is larger than the variance of the real recombination events whenever the markers are linked (θ < 0.5) and π > 0. Figure 2 shows the variance of the apparent recombination events as a function of π, for some different values of recombination rates. In view of the possibility of miscoding for each marker genotype (i.e. ambiguity about their genotypes), standard methods commonly used for genetic map inferences overestimate the recombination rate between loci (or, in other words, underestimate genetic linkage), and underestimate its precision [15]. For example, the maximum likelihood estimator of the recombination rate between the loci (if the possibility of miscoding is ignored) is: with expectation and variance given by (1) and (2), respectively. In more general situations, we have more than just two marker loci, and the goal is to construct the genetic map, i.e., to order these marker loci and to estimate the genetic distances between them. Again, all inferences are based on recombination events observed (apparent recombinations) between the marker loci. The problem of ignoring miscoding may lead to even worse difficulties, e.g., to the mistaken ordering of the loci, specially with dense maps.

BAYESIAN APPROACH FOR GENETIC MAP CONSTRUCTION
First, we will review a Bayesian approach for map construction when miscoding is not taken into account [9]. Consider the genotype of m markers for the individual i as g i = (g i1 , g i2 , . . . , g im ). In a backcross design, for example, g ij = 0 if the individual i is homozygous for the locus j, and 1 otherwise. The sampling model of g i , assuming the Haldane map function, is given by: where λ is the order of the genetic marker loci and θ j is the recombination rate between the loci j and j+1. Considering a sample of n independent individuals, the likelihood of λ and θ is given by: where G is the (n × m) matrix of marker genotypes, with each row representing one individual, and each column related to one marker locus.
In a Bayesian context, rather than maximizing the likelihood, it is modified by a prior and integrated to produce inference summaries for the unknown components in the model. The prior can be chosen based on earlier studies or information from the literature. Here, we use a prior expressed as: p(λ, θ|τ, α, β) = p(θ|λ, α, β)p(λ|τ), (5) where p(λ|τ) is a probability distribution over the m!/2 different orders for the m markers, τ is a set of prior probabilities of each order, and is the recombination rate between genetic markers j and j + 1. A special case of these prior distributions would be uniform across different gene orders, and The Bayes theorem combines the information from the data and the prior knowledge to produce a posterior distribution over all unknown quantities. In this case, the posterior density of λ and θ is given by: p(λ, θ|G, τ, α, β) ∝ p(G|λ, θ)p(θ|λ, α, β)p(λ|τ). (6) Distribution (6) is intractable analytically but MCMC methods such as the Gibbs sampler and the Metropolis-Hastings algorithm [7,8] can be used to draw samples, from which features of marginal distributions of interest can be inferred.

Fully conditional posterior distributions
The Gibbs sampler draws samples iteratively from conditional posterior distributions deriving from (6). The fully conditional posterior distribution of each recombination rate θ j is: where | is the number of recombination events between the loci j and j + 1. This is the kernel of a Beta distribution with parameters (q j + α j ) and (n − q j + β j ).
The updating for the gene order λ involves moves between a set of models, because for distinct ordering, the recombination rates have different meanings. George et al. [9] discuss a reversible jump algorithm, for which recombination rates are converted into map distances, and reverted to new recombination rates after shifting a randomly selected marker around a pivot marker.
Here, another Metropolis-Hastings [12] scheme is presented for the MCMC updating of λ and θ, simultaneously. A new gene ordering is proposed according to a candidate generator density q(.), and new recombination rates are simulated for this new order, using (7). The Markov chain moves from the current state T = (λ, θ) to T * = (λ * , θ * ) with probability: where p(λ, θ|G, τ, α, β) is the joint conditional posterior distribution of the gene ordering λ and recombination rates θ, given by: Under these circumstances, the choice of q(.) is extremely important for an efficient implementation of the MCMC, especially in situations with a large number of marker loci. A bad choice of q(.) would generate a large number of unlikely orders, or even generate inconsistent orders, in relation to the data set. In order to have a better implementation and mixing of the MCMC, some alternatives for the generation of candidate orders for the Metropolis-Hastings step are described in the Appendix.

Missing data
In practice, some marker genotypes are missing. The missing data can be handled by the MCMC approach, with an additional step for updating each missing genotype based on this fully conditional density. For instance, suppose g ij is missing, the genotype for the j-th marker of the individual i. Its fully conditional distribution is Bernoulli with probability p ij = Pr(g ij = 1|G −ij ) given by: where G −ij refers to all elements in G but g ij , and k = 0, 1. Under the Haldane independence assumption, p(g ij = k|θ, G −ij , τ, α, β) depends just on the recombination rates between the locus j and its flanking neighbors, as well as on the genotypes of these neighbor loci, so it can be written as ).

THE PROBABILITY OF MISCODING GENOTYPES
At present, the methods commonly used for map construction ignore the possibility that molecular (marker) data may be read with error, or the error rate has a fixed and known value, as in Lincoln and Lander [15]. Often, however, there is ambiguity about the genotypes. To address these situations, we introduced a new parameter into the model, the probability π of miscoding a genotype. Now we consider that the matrix G of genotypes is unknown, and that we observe a matrix M of genotypes, possibly with some miscoding. The probability of observing a genotype m ij , i.e. the genotype of locus j for individual i, given that the actual genotype is g ij , may be expressed as: where k 1 and k 2 assume values equal to 0 or 1.
Assuming independence between miscodings in different loci and individuals, and considering that the miscoding rate is constant over the genome, the probability of observing a matrix M of genotypes, given the matrix G of actual genotypes, can be expressed as: where n is the number of individuals, m is the number of marker loci, and t = n i=1 m j=1 |m ij − g ij | is the number of miscoding genotypes in the data set. Note that under these circumstances, M is the observed data, and G is now an auxiliary and non-observed matrix. The joint posterior distribution of all unknowns in the model is written now as the product of (9) by (4), (5) and the prior distribution of π, which gives: λ, θ, π|M, τ, α, β, a, b) Assuming a uniform prior probability distribution for λ; Beta(α j , β j ) as prior for each θ j ; and Beta(a, b) as the prior distribution for π, the expression (10) becomes: |, as already defined, is the number of recombination events between the loci j and j + 1. Note that the dependence of this distribution on λ is rendered implicit by the definition of θ j as the recombination rate between the ordered loci j and j + 1.

Fully conditional posterior distributions
The fully conditional posterior distributions of λ and of each θ j have the same forms as discussed before. In the case of G, its conditional distribution is: Given the independence between the recombination events in different intervals (by the Haldane map function), each element in G can be updated independently. If j = 1, i.e. g ij refers to genotypes at one end of the linkage group, its fully conditional posterior distribution can be written as: where G −i1 represents all the elements in G but g i1 , and similarly for g im .
For genotypes at interior markers in the linkage group, the fully conditional posterior distribution becomes: for j = 2, 3, . . . , m − 1. The conditional distribution of the probability of miscoding π is given by: which is the kernel of a Beta distribution with parameters (a+t) and (b+nm−t).

Example 1
Three data sets were simulated to examine the ability of the model discussed in Section 4 to correctly estimate genetic distances and the probability of miscoding. Each simulation considered 300 individuals with genotypes for 5 loci, denoted as ABCDE. The recombination rates between consecutive loci were assumed to be θ AB = 0.09, θ BC = 0.11, θ CD = 0.05 and θ DE = 0.14. The data sets were generated considering π = 0, 0.02 and 0.04, and 3% of missing data for each.
These data sets were analyzed using models with and without the miscoding parameter (π). An equal probability distribution was adopted as prior for the different loci orders. For each recombination rate, a Uniform (0, 0.5) process was considered as prior distribution. Computations were performed using the IML procedure of SAS [19]. Graphical inspection and the Raftery and Lewis diagnostic [18] for the Gibbs output using CODA [1] were used for assessing convergence to the equilibrium distribution, the joint posterior. A burn-in period of 1 000 iterations was adopted, followed by 60 000 iterations with thinning intervals of 20, based on a lag-correlation study. Hence, 3 000 samples were retained for the post-Gibbs analysis.
For all data sets, the gene order was estimated perfectly by both models, with 100% of the MCMC iterations sampling the order ABCDE. It seems that, up to certain levels, inferences about gene ordering is robust to miscoding genotypes, if these occur at random. As discussed earlier (Sect. 2), the effect of miscoding is larger for smaller genetic distances between loci, such as in fine mapping studies. In these cases, the miscoding may lead to ordering estimated with some positions switched for tightly linked markers, as discussed in the next example.  Table I shows the posterior mean and standard deviation for each recombination rate, for the data set without miscoding. The estimates obtained by each model do not present any relevant difference, so it seems that the introduction of the extra parameter (π) into the model, in situations where there is no miscoding, does not affect the estimated genetic map. In this example, the estimate for π was very close to zero, denoting the ability of the model to recognize situations without miscoding. However, because π = 0 relies on the boundary of the parameter space of π, to test for the absence of miscoding for a particular data set, another approach should be employed, such as comparing both models (with and without miscoding) using some criteria, e.g., the Bayes factor or the likelihood ratio test.
The Bayes factors may be computed by taking ratios between estimates of the marginal densities of the data (after integrating out all parameters). If models are taken as equally probable, a priori, then the Bayes factor gives the ratio between the posterior probabilities of the corresponding models. Here, the marginal densities were estimated by calculating harmonic means of likelihoods evaluated at the posterior draws of the Gibbs output [16], and these are presented in Table I. The Bayes factor (in favor of the model without the miscoding parameter) of 20.5 does not denote important differences between both models for modeling this data set.
The results obtained by both models for the data set with 2% miscoding (π = 0.02) are presented in Table II. As expected, the model that ignores the miscoding problem had estimates biased upwards. When the probability of miscoding was introduced into the model, there was improvement on the estimates. In addition, the probability of miscoding was adequately estimated. For the robust model, all the parameter values were inside a credible set of 0.95 of probability. The Bayes factor of 2.01 × 10 6 , in favor of the model with the miscoding parameter, denotes its greater plausibility, when compared to the model ignoring miscoding genotypes. Table II. True parameter values and posterior means and standard deviations (in parenthesis) of the recombination rates considering the data set with 2% of miscoding genotypes and the two models, with and without the miscoding parameter.  Table III. True parameter values and posterior means and standard deviations (in parenthesis) of the recombination rates considering the data set with 4% of miscoding genotypes and the two models, with and without the miscoding parameter.

Example 2
Thirty data sets were simulated, where half had no miscoding and half had 5% miscoding. Here, our main interest was to examine the performance of the models (with and without miscoding) to correctly estimate the gene order with relatively small data sets, both under situations without miscoding and with high levels of miscoding genotypes. Each data set had 100 individuals with genotypes for five markers; no missing data were considered in this study. The recombination rates between consecutive loci were: θ AB = 0.05, θ BC = 0.18, θ CD = 0.02 and θ DE = 0.07. Prior distributions and computations were similar to those described for the previous simulation study.
For the 15 data sets without miscoding, both models (with and without miscoding) yielded the highest posterior probability for the correct order ABCDE.
The Bayes factor was favorable to the model without miscoding for 8 of the data sets (and favorable to the model with the miscoding parameter for the remaining data sets), but with no expressive or relevant values (ranging from 1.2 to 57.9). For the datasets with 5% miscoding, the Bayes factor was always favorable to the model considering miscoding genotypes, with values that ranged from 2.42 × 10 2 to 8.31 × 10 10 .
The model ignoring the miscoding gave the highest posterior probability for a gene ordering other than ABCDE in 8 of the 15 data sets. In the case of the model with the miscoding parameter, just four data sets had the highest posterior probability for a wrong order. If credibility sets with minimum probability of 0.90 are considered, three of the data sets had the correct gene ordering outside of the set. For the model with the miscoding parameter, just one data set presented a probability set that did not contain the correct order.
The results suggest that the model ignoring the miscoding overstates the precision in relation to the gene ordering, sometimes concentrating posterior probability on the wrong (set of) order(s).

ANALYSIS OF EXPERIMENTAL DATA
The data set refers to the RFLP study with Brassica napus using F1-derived double haploid lines. Materials and methods related to the DNA extraction and a preliminary linkage map construction (using maximum likelihood approach) are presented by Ferreira et al. [6]. These data, combined with phenotypic information (flowering time under one of the three vernalization treatments considered in that study), were also analyzed by Ferreira et al. [5] and by Satagopan et al. [20] for the study on quantitative trait loci.
Here, we focus on the estimation of the probability of miscoding, and also on robust construction of the linkage map for a set of marker loci. To illustrate the methods, we consider the data for 105 progeny and 10 maker loci related to the linkage group 9, for which 9% of the genotypes were missing. For simplicity, the marker loci are denoted here by letters (from A through J), according to the order that was estimated by Ferreira et al. [5].
Prior distributions and computations were similar to those describe in Section 5, for the simulation studies. In this case, a longer burn-in period of 2 000 iterations was adopted, followed by 100 000 iterations with thinning intervals of 20. Hence, 5 000 samples were used for the post-Gibbs analysis.
The miscoding rate for this data set was estimated at a level of approximately 1.5%, with a 95% probability set [0.0055; 0.0266]. The model with the miscoding parameter was much more plausible than the one without it, with the Bayes Factor of 4.91 × 10 6 . Posterior probabilities for different gene ordering, estimated by both models (with and without the miscoding parameter) are presented in Table IV. The most probable order for both models was the  (1) 0.0090 0.0388 (1) Loci orders with posterior probabilities smaller than 0.0050. sequence ABCDEFGHIJ, with approximate posterior probabilities of 0.68 and 0.75, respectively for models with and without the miscoding parameter. Some uncertainty on the order arose with the position of the two first markers (A and B), but very high posterior probabilities were found for the sequence CDEFGHIJ of the other eight loci (respectively 0.85 and 0.93 for the models with and without the miscoding parameter). The recombination rates, as expected, presented higher estimates by the model ignoring the possibility of miscoding (Tab. V). Figure 3 shows the estimated map from these two models, using the inverse of the Haldane map function to the recombination rates drawn at each Gibbs iteration.

CONCLUDING REMARKS
The model discussed in this paper provides an appealing robust alternative for genetic map construction in the presence of non-systematic miscoding genotypes. The MCMC implementation of the Bayesian analysis is straightforward, with just some caution to be addressed in relation to the Metropolis-Hastings step for updating the gene ordering. This approach provides more reliable estimates for subsequent studies that use information on genetic maps, such as quantitative trait loci (QTL) search and marker assisted selection.
High values of miscoding probability estimates, however, should raise concern about the molecular data, and a revaluation of the marker genotypes may be a good approach. In situations with relatively large rates of miscoding, the high frequency of apparent recombinations may not be recognized as the reflect of miscoding genotypes, but due to bigger values of real genetic recombinations. For these cases, a multilocus feasible map function, which assumes interdependence between different marker intervals [17], could be a better alternative to the Haldane map function.
This paper can be extended in various ways to analyze genetic data originated from different designs (e.g. F2 progenies, granddaughter design, etc.). Furthermore, the idea of considering the possibility of miscoding genotypes may be used for QTL analysis as well. The methodology for robust estimation under miscoding genotypes can be adapted to handle multiallelic loci situations. For example, consider that g ij can assume one of t genotypes, denoted as 1, 2, . . . , t. In these cases, the probability of observing a genotype m ij equal to r (r = 1, 2, . . . , t) would be: Pr(m ij = r) = Pr(g ij = r) Pr(m ij = r|g ij = r) + Pr(g ij = r) Pr(m ij = r|g ij = r) = Pr(g ij = r) Pr(m ij = r|g ij = r) + t s=1;s =r π rs Pr(g ij = s), where π rs = Pr(m ij = r|g ij = s) represents the miscoding probability of observing a genotype as r, when the actual genotype is s. The miscoding probabilities π rs could be considered, for example, proportional to the distance of each allele in the gel, which reflects the size of each allele (in base pairs). In this case, the miscoding probabilities would be written as: π rs = Pr(m ij = r|g ij = s) = φ (distance between alleles r and s).
The results found in this work suggest that, unless there is strong reason to believe in the absence of ambiguity about genotypes, it may be safer to use the robust model, which would provide a more reliable estimate of the genetic map.
Towards the completion of this research, we became aware of related studies by Keller [13] under the supervision of G.A. Churchill. Their work confirms the utility of our approach to assess miscoding and improve the estimation of map distances. Some of their methods appear as part of the R/QTL software module of Broman [2].