Efficient and Accurate Construction of Genetic Linkage Maps from the Minimum Spanning Tree of a Graph (original) (raw)

Abstract

Genetic linkage maps are cornerstones of a wide spectrum of biotechnology applications, including map-assisted breeding, association genetics, and map-assisted gene cloning. During the past several years, the adoption of high-throughput genotyping technologies has been paralleled by a substantial increase in the density and diversity of genetic markers. New genetic mapping algorithms are needed in order to efficiently process these large datasets and accurately construct high-density genetic maps. In this paper, we introduce a novel algorithm to order markers on a genetic linkage map. Our method is based on a simple yet fundamental mathematical property that we prove under rather general assumptions. The validity of this property allows one to determine efficiently the correct order of markers by computing the minimum spanning tree of an associated graph. Our empirical studies obtained on genotyping data for three mapping populations of barley (Hordeum vulgare), as well as extensive simulations on synthetic data, show that our algorithm consistently outperforms the best available methods in the literature, particularly when the input data are noisy or incomplete. The software implementing our algorithm is available in the public domain as a web tool under the name MSTmap.

Author Summary

Genetic linkage maps are cornerstones of a wide spectrum of biotechnology applications. In recent years, new high-throughput genotyping technologies have substantially increased the density and diversity of genetic markers, creating new algorithmic challenges for computational biologists. In this paper, we present a novel algorithmic method to construct genetic maps based on a new theoretical insight. Our approach outperforms the best methods available in the scientific literature, particularly when the input data are noisy or incomplete.

Introduction

Genetic linkage mapping dates back to the early 20th century when scientists began to understand the recombinational nature and cellular behavior of chromosomes. In 1913 Sturtevant studied the first genetic linkage map of chromosome X of Drosophila melanogaster [1]. Genetic linkage maps began with just a few to several tens of phenotypic markers obtained one by one by observing morphological and biochemical variations of an organism, mainly following mutation. The introduction of DNA-based markers such as restriction fragment length polymorphism (RFLP), randomly amplified polymorphic DNA (RAPD), simple sequence repeats (SSR) and amplified fragment length polymorpshim (AFLP) caused genetic maps to become much more densely populated, generally into the range of several hundred to more than a thousand markers per genome. More recently, the number of markers has surged well above 1,000 in a number of organisms with the adoption of DArT, SFP and especially SNP markers, the latter providing avenues to 100,000 s to millions of markers per genome. In plants, one of the most densely populated maps is that of Brassica napus [2], which was developed from an initial set of 13,551 markers. High density genetic maps facilitate many biological studies including map-based cloning, association genetics and marker assisted breeding. Because they do not require whole genome sequencing and require relatively small expenditures for data acquisition, high density genetic linkage maps are currently of great interest.

A genetic map usually is built using input data composed of the states of loci in a set of meiotically derived individuals obtained from controlled crosses. When an order of the markers is computed from the data, the recombinational distance is also estimated. To characterize the quality of an order, various objective functions have been proposed, e.g., minimum Sum of Square Errors (SSE) [3], minimum number of recombination events (COUNT) [4], Maximum Likelihood (ML) [5], Modified Maximum Likelihood (MML) [6] which tries to incorporate the presence of possible genotype errors into the ML model, maximum Sum of adjacent LOD scores (SALOD) [7], minimum Sum of Adjacent Recombination Fractions (SARF) [8], minimum Product of Adjacent Recombination Fractions (PARF) [9]. Searching for an optimal order with respect to any of these objective functions is computationally difficult. Enumerating all the possible orders quickly becomes infeasible because the total number of distinct orders is proportional to n!, which can be very large even for a small number n of markers.

The connection between the traveling salesman problem and a variety of genomic mapping problem is well known, e.g., for the physical mapping problem [10],[11], the genetic mapping problem [12],[13] and the radiation hybrid ordering problem [14]. Various searching heuristics that were originally developed for the traveling salesman problem, such as simulated annealing [15], genetic algorithms [16], tabu search [17],[18], ant colony optimization, and iterative heuristics such as K-opt and Lin-Kernighan heuristic [19] have been applied to the genetic mapping problem in various computational packages. For example, JoinMap [5] and Tmap [6] implement simulated annealing, Carthagene [12],[20] uses a combination of simulated annealing, tabu search and genetic algorithms, AntMap [21] exploits the ant colony optimization heuristic, [22] is based on genetic algorithms, and [23] takes advantage of evolutionary algorithms. Finally, Record [4] implements a combination of greedy and Lin-Kernighan heuristics.

Most of the algorithms proposed in the literature for genetic linkage mapping find reasonably good solutions. Nonetheless, they fail to identify and exploit the combinatorial structures hidden in the data. Some of them simply start to explore the space of the solutions from a purely random order (see, e.g., [12],[23],[5],[21]), while others start from a simple greedy solution (see, e.g., [4],[3]). Here we show both theoretically and empirically that when the data quality is high, the optimal order can be identified very quickly by computing a minimal spanning tree of the graph associated with the genotyping data. We also show that when the genotyping data is noisy or incomplete, our algorithm consistently constructs better genetic maps than the best available tools in the literature. The software implementing our algorithm is currently available as a web tool under the name MSTmap.

Materials and Methods

We are concerned with genetic markers in the form of single nucleotide polymorphism (SNP), more specifically biallelic SNPs. By convention, the two alternative allelic states are denoted as A and B respectively. The organisms considered here are diploids with two copies of each chromosome, one from the mother and the other from the father. A SNP locus may exist in the homozygous state if the two allele copies are identical, and in the heterozygous state otherwise.

Various population types have been studied in association with genetic mapping, which includes Back Cross (BC1), Doubled Haploid (DH), Haploid (Hap), Recombinant Inbred Line (RIL), advanced RIL, etc. Our algorithm can handle all of the aforementioned population types. For the sake of clarity, in what follows we will concentrate on the DH population (see the section on barley genotyping data for details on DH populations). The application of our method to Hap, advanced RIL and BC1 populations is straightforward. In Supplementary Text S1, we will discuss the extension of our method to the RIL population (see, e.g., [24] for an introduction to RIL populations).

Building a genetic map is a three-step process. First, one has to partition the markers into linkage groups, each of which usually corresponds to a chromosome (sometimes multiple linkage groups can reside on the same chromosome if they are far apart). More specifically, a linkage group is defined as a group of loci known to be physically connected, that is, they tend to act as a single group (except for recombination of alleles by crossing-over) in meiosis instead of undergoing independent assortment. The problem of assigning markers to linkage groups is essentially a clustering problem. Second, given a set of markers in the same linkage group, one needs to determine their correct order. Third, the genetic distances between adjacent markers have to be estimated. Before we describe the algorithmic details, the next section is devoted to a discussion on the input data and our optimization objectives.

Genotyping Data and Optimization Objective Functions

The doubled haploid individuals (a set collectively denoted by N) are genotyped on the set M of markers, i.e., the state of each marker is determined. The genotyping data are collected into a matrix of size m_×_n, where m = |M| and n = |N|. Each entry in corresponds to a marker and individual pair, which is also called an observation. Due to how DH mapping populations are produced (please refer to section on barley genotyping data for details), each observation can exist in two alternative states, namely homozygous A or homozygous E, which are denoted as A and B respectively. The case where there is missing data will be discussed later in this manuscript.

For a pair of markers l 1, l 2 ∈ M and an individual c ∈ N, we say that c is a recombinant with respect to l 1 and l 2 if c has genotype A on l 1 and genotype B on l 2 (or vice versa). If l 1 and l 2 are in the same linkage group, then a recombinant is produced if an odd number of crossovers occurred between the paternal chromosome and the maternal chromosome within the region spanned by l 1 and l 2 during meiosis. We denote with P i,j the probability of a recombinant event with respect to a pair of markers (li,lj). P i,j varies from 0.0 to 0.5 depending on the distance between li and lj At one extreme, if li and lj belong to different LGs, then P i,j = 0.5 because alleles at li and lj are passed down to next generation independently from each other. At the other extreme, when the two markers li and lj are so close to each other that no recombination can occur between them, then P i,j = 0.0. Let (li,lj) and (lp,lq) be two pairs of markers on the same linkage group. We say that the pair (li,lj) is enclosed in the pair (lp,lq) if the region of the chromosome spanned by li and lj is fully contained in the region spanned by lp and lq. A fundamental law in genetics is that if (li,lj) is enclosed in (lp,lq) then P _i,j_≤P p,q.

As mentioned in the Introduction, a wide variety of objective functions have been proposed in the literature to capture the quality of the order (SSE, COUNT, ML, MML, SALOD, SARF, PARF, etc.). With the exception of SSE and MML, the rest of the objective functions listed above can be decomposed into a simple sum of terms involving only pairs of markers. Thus, we introduce a weight function w: M_×_M_→ℜ to be defined on pairs of markers. The function w is said to be semi-linear if w(i, j)≤_w(p, q) for all (li,lj) enclosed in (lp,lq). For example, if we have three markers in order {l 1,l 2,l 3} and an associated weight function w that satisfies semi-linearity, we have w(1,3)≥w(1,2) and w(1,3)≥w(2,3) since (l 1,l 2) and (l 2,l 3) are enclosed in (l 1,l 3), but it is not necessary the case that w(1,3) = w(1,2)+w(2,3). The concept of semi-linearity is essential for the development of our marker ordering algorithm as explained below.

For example, the function w(i, j) = P i,j is semi-linear. Another commonly used weight function is wlp(i, j) = log(P i,j). Since the logarithm function is monotone, then wlp(i, j) is also semi-linear. A more complicated weight function is wml(i, j) = −[P _i,j_log(P i,j)+(1−P i,j)log(1−P i,j)]. It is relatively easy to verify that wml(i, j) is a monotonically increasing function of P i,j when 0≤P _i,j_≤0.5, and therefore wml is also semi-linear. Observe that all these weight functions are functions in P i,j. Although the precise value of P i,j is unknown, we can compute their estimates from the total number of recombinants in the input genotyping data. For DH populations, the total number of recombinants in N with respect to the pair (li,lj) can be easily determined by computing the number di ,j of positions in which row and row do not match, which corresponds to the Hamming distance between and . It is easy to prove that di ,j/n corresponds to the maximum likelihood estimate (MLE) for P i,j. When we replace P i,j by its maximum likelihood estimate di ,j/n, we obtain the following approximate weight functions: _wp_′(i, j) = di ,j/n, _wlp_′ (i, j) = log(di ,j/n), and .

Our optimization objective is to identify a minimum weight traveling salesman path with respect to either of the aforementioned approximated weight functions, which will be discussed in further details below. We should mention that if _wp_′ is used as the weight function, then our optimization objective is equivalent to the SARF or COUNT objective functions (up to a constant). If instead _wlp_′ is used, then our optimization objective is equivalent to the logarithm of the PARF objective function (up to a constant). Lastly, if _wml_′ is employed, our objective function is equivalent to the negative of the logarithm of the ML objective function as being employed in [3],[5],[12],[20] (again, up to a constant). Unless otherwise noted, _wp_′ is the objective function being employed in the rest of this paper. The experimental results will show that the specific choice of objective function does not have a significant impact on the quality of the final map. In particular, both functions _wp_′ and _wml_′ produce very accurate final maps.

Clustering Markers into Linkage Groups

First observe that when two markers li and lj belong to two different linkage groups, then P i,j = 0.5 and consequently di ,j will be large with high probability. More precisely, let li and lj be two markers that belong to two different LGs, and let di ,j be the Hamming distance between and . Then, where δ<0.5. The proof of this bound can be found in Supplementary Text S1.

In order to cluster the markers into linkage groups, we construct a complete graph G(M, E) over the set of all markers. We set the weight of an edge (li, lj) ∈ E to the pairwise distance di ,j between li and lj. As shown in Theorem 1 of Supplementary Text S1, if two markers belong to different LGs, then the distance between them will be large with high probability. Once a small probability is chosen by the user (default is = 0.00001, in general one should choose <0.0001.), we can determine δ by solving the equation −2(n/2−δ)2/n = log_e_. We then remove all the edges from G(M, E) whose weight is larger than or equal to δ. The resulting graph will break up into connected components, each of which is assigned to a linkage group.

A proper choice of appears critical in our clustering algorithm. In practice, however, this is not such a crucial issue because the recombinant probability between nearby markers on the same linkage group is usually very small (usually less than 0.05 in dense genetic maps). According to our experience, our algorithm is capable of determining the correct number of LGs for a fairly large range of values of (see Results and Discussion).

Ordering Markers in Each Linkage Groups

Let us assume now that all markers in M belong to the same linkage group, and that M has been preprocessed so that di ,_j_>0 for all i, j ∈ M. The excluded markers for which di ,j = 0 are called co-segregating markers, and they identify regions of chromosomes that do not recombine. In practice, we coalesce co-segregating markers into bins, where each bin is uniquely identified by any one of its members. Let G(M, E) be an edge-weighted complete undirected graph on the set of vertices M, and let w be one of the weight functions defined above. A traveling salesman path (TSP) Γ in G is a path that visits every marker/vertex once and only once. The weight w(Γ) of a TSP Γ is the sum of the weights of the edges on Γ.

The main theoretical insight behind our algorithm is the following. When w is semi-linear, the minimum weight TSP of G corresponds to the correct order of markers in M. Furthermore when the minimum spanning tree (MST) of G is unique, the minimum weight TSP of G (and thus, the correct order) can be computed by a simple MST algorithm (such as Prim's algorithm). Details of these mathematical facts (with proofs) are given in Supplementary Text S1.

We now turn our attention to the problem of finding a minimum weight TSP in G with respect to one of the approximate weight functions. When the data are clean and n is large, the maximum likelihood estimates di ,j/n will be close to the true probabilities P i,j. Consequently it is reasonable to expect that those approximate weight functions will be also semi-linear, or “almost” semi-linear. Although only in the former case our theory (in particular, Lemma 1 in Supplementary Text S1) guarantees that the minimum weight TSP will correspond to the true order of the markers, in our simulations the order is recovered correctly in most instances. In order to find the minimum weight TSP, we first run Prim's algorithm on G to compute the optimum spanning tree, which takes O(n_log_n). If the MLEs are accurate so that the approximate weight function is semi-linear, our theory (in particular Lemma 2 in Supplementary Text S1) ensures that the MST is a TSP.

In practice, due to noise in the genotyping data or due to an insufficient number of individuals, the spanning tree may not be a path – but hopefully “very close” to a path. This is exactly what we observed when running MST algorithm on both real data and noise-free simulated data – the MST produced is always “almost” a path. In Results and Discussion we compute the fraction ρ of the total number of markers in the linkage group that belong to the longest path of the MST. The closer is ρ to 1.0, the closer is the MST to a path. Table 1 on the barley datasets and Figure 1 on simulated data show that ρ is always very close to 1.0 when the data is noise-free.

Table 1. Summary of the clustering results for the barley data sets.

Data set	# markers (# bins)	# LGs	Sizes of the LGs
OWB	1562(509)	7	168(65), 235(73), 255(91), 211(60), 278(89), 202(64), 213(67)	0.9978
SM	1270(396)	8	148(49), 217(57), 242(63), 130(49), 225(80), 122(40), 183(57), 3(1)	0.9971
MB	1652(443)	8	215(60), 279(72), 246(77), 141(39), 299(74), 219(54), 248(65), 5(2)	1.0000

Figure 1. Average ρ (rho) for thirty runs on simulated data for several choices of the error rates (and no missing data).

The variable n represents the number of individuals, and m represents the number of markers.

When a tree is not a path, we proceed as follows. First, we find the longest path in the MST, hereafter referred to as the backbone. The nodes that do not belong to the path will be first disconnected from it. Then, the disconnected nodes will be re-inserted into the backbone one by one. Each disconnected node is re-inserted at the position which incurs the smallest additional weight to the backbone. The path obtained at the end of this process is our initial solution, which might not be locally optimal.

Once the initial solution is computed, we apply three heuristics that iteratively perform local perturbations in an attempt to improve the current TSP. First, we apply the commonly-used K-opt (K = 2 in this case) heuristic. We cut the current path into three pieces, and try all the possible rearrangements of the three pieces. If any of the resulting paths has less total weight, it will be saved. This heuristic is illustrated in Figure 2-C1 . This procedure is repeated until no further improvement is possible. In the second heuristic, we try to relocate each node in the path to all the other possible positions. If this relocation reduces the weight, the new path will be saved. The second heuristic is illustrated in Figure 2-C2 .

Figure 2. An illustration of the MST-based algorithm.

(A) The MST obtained for a synthetic example; the MST is not a TSP yet; the backbone of the MST is shown with dotted edges. (B) An initial TSP obtained from the backbone (see text for details). The dotted edges represent marker pairs in the wrong order. Several local improvement operations are applied to further improve the TSP, namely 2-OPT (C1), node-relocation (C2) and block-optimize (C3). The final TSP is shown in (D).

In our experiments, we observed that K-opt or node relocation may get stuck in local optima if a block of nodes have to be moved as a whole to a different position in order to further improve the TSP. In order to work around this limitation, we designed a third local optimization heuristic, which is called block-optimize. The heuristic works as follows. We first partition the current TSP into blocks consisting of consecutive nodes. Let l 1, l 2,…, lm be the current TSP. We will place li and li +1 in the same block if (1) w(i, i+1)≤w(i, j) for all i+1<j_≤_m and (2) w(i, i+1)≤w(k, i+1) for all 1≤i. Intuitively, the partitioning of the nodes into blocks reflects the fact that the order between the nodes within a block is stable and should be fixed, while the order among the blocks needs to be further explored. After partitioning the current TSP into blocks, we then carry out the K-opt and node relocation heuristics again by treating a block as a single node. The last heuristic, block-optimize, is illustrated in Figure 2-C3 .

We apply the 2-opt heuristic, the relocation heuristic and the block-optimize heuristic iteratively until none can further reduce the weight of the path. The resulting TSP represents our final solution. A sketch of our ordering algorithm is presented as Algorithm 1 in Supplementary Text S1.

Dealing with Missing Data

In our discussion so far, we assumed no missing genotypes. This assumption is not very realistic in practice. As it turns out, it is common to have missing data about the state of a marker. Our simulations shows that missing observations do not have as much negative impact on the accuracy of the final map as do genotype errors. Thus, it appears beneficial to leave uncertain genotypes as missing observations rather than arbitrarily calling them one way or the other.

We deal with missing observations via an Expectation Maximization (EM) algorithm. Observe that if we knew the order of the markers (or, bins, if we have co-segregating markers), the process of imputing the missing data would be relatively straightforward. For example, suppose we knew that marker l 3 immediately follows marker l 2, and that l 2 immediately follows marker l 1. Let us denote with Pˆ_i,j_ the estimate of the recombinant probabilities between markers li and lj. Let us assume that for an individual c the genotype at locus l 2 is missing, but the genotypes at loci l 1 and l 3 are available. Without loss of generality, let us suppose that they are both A. Then, the posterior probability for the genotype at locus l 2 in individual c is

and P{genotype in c at l 2 is B} = 1−P{genotype in c at l 2 is A}. This posterior probability is the best estimate for the genotype of the missing observation. Similarly, one can compute the posterior probabilities for different combinations of the genotypes at loci l 1 and l 2.

In order to deal with uncertainties in the data and unify the computation with respect to missing and non-missing observations, we replace each entry in the genotype matrix that used to contain symbols A/B with a probability. The probability stored in now represents the confidence that we have about marker li in individual cj of being in state A. For the known observations, the probabilities are fixed to be 1 or 0 depending whether the genotype observed is A or B, respectively. The probabilities for the missing observations will be initially set to 0.5.

Our EM algorithm works as follows. We first compute a reasonably good initial order of the markers by ignoring the missing data. To do so, we compute the normalized pairwise distance di ,j as di ,j = xn/_n_′, where _n_′ is the number of individuals having non-missing genotypes at both loci li and lj, x is the number of individuals having different genotypes at loci li and lj among the _n_′ individuals being considered, and n is the total number of individuals. With the normalized pairwise distances, we rely on the function Order (Supplementary Text S1, Algorithm 1) to compute an initial order.

After an initial order has been computed, we iteratively execute an E-step followed by an M-step. In the E-step, given the current order of the markers, we adjust the estimate for a missing observation at marker on individual cj as follows

	(1)

where is the marker immediately preceding in the most recent ordering, is the marker immediately following , and La ,b,c is the likelihood of the event at the three consecutive loci. The right hand side of Equation (1) is simply the posterior probability of the missing observation being A. The quantity La ,b,c is straightforward to compute. For example, , where , are the MLEs for and respectively, and and are the pairwise distance computed in the previous M step or the initialization step when the initial order is computed. In the case where the missing observation is at the beginning or at the end of the map, the above estimates must be adjusted slightly.

Following the E-step, we execute an M-step. We need to re-compute the pairwise distances according to the new estimates of the missing data. Given that now contains probabilities, the expected pairwise distance between two markers li and lj can be computed as follows

	(2)

With the updated pairwise distances, we use the function Order again to compute a new order of the markers.

An E-step is followed by another M-step, and this iterative process continues until the marker order converges. In our experimental evaluations, the algorithm converges quickly, usually in less than ten iterations. The pseudo-code for the EM algorithm is presented as Algorithm 2 in Supplementary Text S1.

We should mention that our EM algorithm is significantly different from the EM algorithms employed in MapMaker [25] or Carthagene [12]. The EM algorithms used in MapMaker and Carthagene are not used to determine the order, but rather to estimate the recombination probabilities between adjacent markers in the presence of missing data. In MSTmap, our EM method deals with missing data in a way which is very tightly coupled with the problem of finding the best order of the markers.

Detecting and Removing Erroneous Data

As commonly observed in the literature (see, e.g., [26],[27]), with conventional mapping software such as JoinMap, Carthagene or Record, the existence of genotyping errors can have a severe impact on the quality of the final maps. With even a relatively small amount of errors, the order of the markers can be compromised. Therefore, it is necessary to detect erroneous genotype data.

In practice, genotype errors do not distribute evenly across markers. Usually a few “bad markers” tend to be responsible for the majority of errors. Removing those bad markers is relatively easy because they will appear isolated from the other markers in terms of Hamming distance di ,j. We can simply look for markers which are more than a certain distance (a parameter specified by the user, default is 15 cM) away from all other markers. Bad markers are deleted completely from the dataset.

Residual sources of genotyping errors are more difficult to deal with. Given that in practice missing observations have much less negative impact on the quality of the map than errors, our strategy is to identify suspicious data and treat them as missing observations. When doing so, however, we should be careful not to introduce too many missing observations.

In high density genetic mapping, a genotype error usually manifests itself as a singleton (or a double cross-over) under a reasonably accurate ordering of the markers. A singleton is a SNP locus whose state is different from both the SNP marker immediately before and after. An example of a singleton is illustrated in Figure 3. A reasonable strategy to deal with genotyping errors is to iteratively remove singletons by treating them as missing observations and then refine the map by running the ordering algorithm. The problem of this strategy is that at the beginning of this process the number of errors might be high and the marker orders are not very accurate. As a consequence, the identified singletons might be false positives.

Figure 3. An example of a singleton (double crossover).

Each row refers to an individual and each column refers to a marker locus. Given the current order, the entry (c 1, l 7) appears to be a possible error because its state differs from both its immediately preceding and following markers.

We deal with this problem by taking into consideration the neighborhood of a marker instead of just looking at the immediately preceding and following ones. Along the lines of the approach proposed in SMOOTH [26], we define

	(3)

where i is a marker and j is an individual. The quantity estimates the state of a locus given its distance to other markers (and their states). The estimate is a weighted average of the information from all other markers, and the weight is proportional to the inverse of the square of the distance. One can approximate by considering only a small set (say 8) of the closest markers to compute the estimate. When is far from , then the observation (i, j) is regarded as suspicious and is treated as a missing observation. In our implementation, we consider an observation suspicious when .

In our iterative process (1) we detect possible errors using , (2) we refine the map by calling Algorithm 3 (Supplementary Text S1), (3) we estimate the missing data and (4) we re-compute the distances di ,j according to Equation (2). The number of iterations should depend on the quality of the data. If the original data are noisy, more iterations are needed. We propose an adaptive approach to dynamically determine when to stop the iterative process. Let X be the total number of suspicious observations that have been detected so far plus the total number of cross-overs still present in the latest order. Observe that an error usually result in two cross-overs (refer to Figure 3 for an example). By treating an error as a missing observation, the total number of suspicious observations will increase by one, but the total number of crossovers will decrease by two. Overall, the quantity X will decrease by one in the next iteration. On the other hand, if an observation is indeed correct but is mistakenly treated as a missing, X will increase by one in the next iteration. Based on this analysis, we stop the iterative process as soon as the quantity X begins to increase.

In passing, we should mention that Equation (3) can also be used to estimate missing data. According to our experimental studies, it gives comparable performance to the EM algorithm we proposed in the previous section. Our complete algorithm, which incorporates clustering markers into linkage groups, missing data estimation, and error detection, is presented in Supplementary Text S1 as Algorithm 3. We named our algorithm MSTmap since the initial orders of the markers are inferred from the MST of a graph.

Computing the Genetic Distances

Mapping functions are used to convert the recombination probabilities r to a genetic distance D that reflects the actual frequency of crossovers. They correct for undetected double crossovers and cross over interference. The Haldane mapping function [28] assumes that crossovers occur independently and thus do not adjust for interference, while the Kosambi mapping function [29] adjusts for crossover interference assuming that one crossover inhibits another nearby. The Haldane distance function is defined as , whereas the Kosambi distance function is . Both functions are defined for 0<r<0.5. When the crossover interference is not known, Kosambi should be used by default. If the frequency of crossover is low or in the case of high density maps when the distance between adjacent markers is low, either of them can be safely used [30].

Results/Discussion

We implemented our algorithm in C++ and carried out extensive evaluations on both real data and simulated data. The software is available in the public domain at the address http://www.cs.ucr.edu/~yonghui/mstmap.html.

The four tools benchmarked here were run on relatively fast computers by contemporary standards. MSTmap was run on a Linux machine with 32 1.6 GHz Intel Xeon processors and 64 GB memory, Carthagene was executed on a Linux machine with a dual-core 2 GHz Intel Xeon processor and 3 GB memory whereas Record and JoinMap were both run on a Windows XP machine with a dual-core 3 GHz Intel Pentium processor and 3 GB of main memory. We had to use different platforms because some of these tools are platform-specific (i.e., Record and JoinMap only run on Windows, Carthagene only runs on Linux). Note that MSTMap was run on the platform with the slowest CPU. The fact that MSTMap was run on a machine with multiple CPUs and large quantity of main memory did not create an unfair advantage. MSTMap is single threaded and thus it exploits only one CPU. The space complexity of MSTMap is O(n 2), where n is the number of markers per linkage group. Under our simulation studies, n is less than 500, which translates in about of 1 GB memory which is a relatively small amount.

Barley Genotyping Data

The real genotyping data come from an ongoing genetic mapping project for barley (Hordeum vulgare) (see http://barleycap.org/ and http://www.agoueb.org/ for more details on this project). In total we made use of three mapping populations, all of which are DH populations. Doubled haploid (DH) technology refers to the use of microspore or anther culture (ovary culture in some species) to obtain haploid embryos and subsequently double the ploidy level. Briefly, a DH population is prepared as follows. Let M be the set of markers of interest. Pick two highly inbred (fully homozygous) parents p 1 and p 2. We assume that the parents p 1 and p 2 are homozygous for every marker in M (those markers that are heterozygous in either p 1 or p 2 are simply excluded from consideration), and the same marker always has different allelic states in the two parents (those markers having the same allelic state in both parents are also excluded from M). By convention, we use symbol A to denote the allelic states appearing in p 1 and B to denote the allelic states appearing in p 2. Parent p 1 is crossed with parent p 2 to produce the first generation, called F1. The individuals in the F1 generation are heterozygous for every marker in M, with one chromosome being all A and the other chromosome being all B. Gametes produced by meiosis from the F1 generation are fixed in a homozygous state by doubling the haploid chromosomes to produce a doubled haploid individual. The ploidy level of haploid embryo could be doubled by chemical (example colchicine) treatment to obtain doubled haploid plants with 100% homozygosity. This technology is available in some crops to speed up the breeding procedure (see, e.g., [31]).

The first mapping population is the result of crossing Oregon Wolfe Barley Dominant with Oregon Wolfe Barley Recessive (see http://barleyworld.org/oregonwolfe.php). The Oregon Wolfe Barley (OWB) data set consists of 1,562 markers genotyped on 93 individuals. The second mapping population is the result of a cross between Steptoe with Morex (see http://wheat.pw.usda.gov/ggpages/SxM/), which consists of 1,270 markers genotyped on 92 individuals. It will be referred to as the SM dataset from here on. The third mapping population is the result of a cross between Morex with Barke recently developed by Nils Stein and colleagues at the Leibniz Institute of Plant Genetics and Crop Plant Research (IPK), which contains 1,652 markers on 93 individuals. This latter dataset will be referred to as MB in our discussion. The genotypes of SNPs for the above data sets were determined via an Illumina GoldenGate Assay. Very few of the genotypes are missing. The three mapping populations together contain only 51 missing genotype calls out of the total of 417,745. The three barley data sets are expected to contain seven LGs, one for each of the seven barley chromosomes.

Synthetic Genotyping Data

The simulated data set is generated according to the following procedure (which is identical to the one used in [4]). First four parameters are chosen, namely the number m of markers to place on the genetic map, the number n of individuals to genotype, the error rate η and the missing rate γ. Following the choice of m, a “skeleton” map is produced, according to which the genotypes for the individuals will be generated. The markers on the skeleton map are spaced at a distance of 0.5 centimorgan plus a random distance according to a Poisson distribution. On average, the adjacent markers are 1 centimorgan apart from each other. The genotypes for the individuals are then generated as follows. The genotype at the first marker is generated at random with probability 0.5 of being A and probability 0.5 of being B. The genotype at the next marker depends upon the genotype at the previous marker and the distance between them. If the distance between the current marker and the previous marker is x centimorgan, then with probability x/100, the genotype at the current locus is the opposite of that at the previous locus, and with probability 1−x/100 the two genotypes are the same. Finally, according to the specified error rate and missing rate, the current genotype is flipped to purposely introduce an error or is simply deleted to introduce a missing observation. Following this procedure, various datasets for a wide range of choices of the parameters were generated.

Evaluation of the Clustering Algorithm

First, we evaluated the effectiveness of our clustering algorithm on the three datasets for barley. Since the genome of barley consists of seven chromosome pairs, we expected the clustering algorithm to produce seven linkage groups. Using the default value for ε, our algorithm produced seven linkage groups for the OWB data set, eight linkage groups for the MB data set and eight linkage groups for the SM data set. The same results can be obtained in a rather wide range of values of . For example, for any choice of ∈ [0.000001,0.0001] the OWB data set is always clustered into seven LGs. The smallest linkage group in the SM data set contains three markers in a single bin. The smallest linkage group in the MB data set contains five markers in two bins. By comparing the three maps with each other, we determined that these small isolated linkage groups in the SM and MB populations are at the telomere far away from the rest of the markers on the same chromosome. The result of the clustering algorithm is summarized in Table 1. We also compared our clusters with those produced by JoinMap. The clusters were identical.

Evaluation of the Quality of the Minimum Spanning Trees

In this second evaluation step, we verified that on real and simulated data, the MSTs produced by MSTmap are indeed very close to TSPs. This experimental evaluation corroborates the fact that the MSTs produced are very good initial solutions. Here, we computed the fraction ρ of the total number of bins/vertices in the linkage group that belong to the longest path (backbone) of the MST. The closer is ρ to 1, the closer is the MST to a path.

Table 1 shows that on the barley data sets, the average value for ρ for the seven linkage groups (not including the smallest LG in the SM data set) is always very close to 1. Indeed, 19 of the 21 MSTs are paths. The remaining 2 MSTs are all very close to paths, with just one node hanging off the backbone. When the MSTs generated by our algorithm are indeed paths, the resulting maps are guaranteed to be optimal, thus increasing the confidence in the correctness of the orders obtained.

On the simulated dataset with no genotyping errors, ρ is again close to one (see Figure 1) for both n = 100 and n = 200 individuals. When the error rate is 1%, the ratio drops sharply to about 0.6. This is due to the fact that the average distance between nearby markers is only one centimorgan. One percent error introduces an additional distance of two centimorgans which is likely to move a marker around in its neighborhood. The value for ρ for error rates up to 15% are computed and are presented in Figure 1. At 15% error rate, the backbone contains only about 40% of the markers. However, this relatively short backbone is still very useful in obtaining a good map since it can be thought as a sample of the markers in their true order. Also, observe that increasing the number of individuals will slightly increase the length of the backbone. However, the ratio remains the same irrespective of the number of markers we include on a map (data not shown).

Evaluation of the Error Detection Algorithm

Third, we evaluated the accuracy and the effectiveness of the error detection algorithm. Synthetic datasets with a known map and various choices of error rate and missing rate were generated. We ran our tool on each dataset, and by comparing the map produced by MSTMap with the true map we collected a set of relevant statistics.

Given a map produced by MSTmap we define a pair of markers to be erroneous if their relative order is reversed when compared to the true map. The number E of erroneous marker pairs ranges from 0 to m(_m_−1)/2, where m is the number of markers. We have E = 0 when the two maps are identical, and E = m(_m_−1)/2 when one map is the reverse of the other. Since the orientation of the map is not important in this context, whenever _E_>m(_m_−1)/4, one of the maps is flipped and E is recomputed. Notice that E is more sensitive to global reshuffling than to local reshuffling. For example, assume that the true order is the identity permutation. The value of E for the following order is m(_m_−1)/4, whereas E for the order 2,1,4,3,6,5,…,m,_m_−1 is m(_m_−1)/2. For reasonably large m, m(_m_−1)/2 is much smaller than m(_m_−1)/4. The fact that E is more sensitive to global reshuffling is a desirable property since biologists are often more interested in the correctness of the global order of the markers than the local order.

The number of erroneous marker pairs conveys the overall quality of the map produced by MSTmap, however E depends on the number m of markers. The larger is m, the larger E will be. Sometimes it is useful to normalize E by taking the transformation 1−(4_E_/(m(_m_−1))). The resulting statistic is essentially the Kendall's τ statistic. The τ statistic ranges from 0 to 1. The closer is the statistic to 1, the more accurate the map is. We will present the τ statistic along with the E statistic when it is necessary.

The next three statistics we collected are the percentage of true positives, the percentage of false positives, and the percentage of false negatives, which are denoted as %t_ pos, %tf_ pos and %f_ neg respectively. For each dataset, the list of suspicious observations identified by MSTmap is compared with the list of true erroneous observations that were purposely added when the data was first generated. The value of %t_ pos is the number of suspicious observations that are truly erroneous divided by the total number nm of observations. The value %f _ pos is the number suspicious observations that are in fact correct divided by the total number of observations. Likewise, %f_ neg is the number of erroneous observations that are not identified by MSTmap. The three performance metrics are intended to capture the overall accuracy of the error detection scheme. Finally, we collected the running time on each data set.

Table 2 summarizes the statistics for n = 100, m = 100, when the error rate and the missing rate range from 0% to 15%. An inspection of the table reveals that irrespective of the choice of η and γ our error detection scheme is able to detect most of the erroneous observations without introducing too many false positives. When the input data are noisy, the quality of the final maps with error detection is significantly better than those without. However, if the input data are clean (corresponding to rows in the table where η = 0), the quality of the maps with error detection deteriorates slightly. Results for other choices of m and n are presented in Table 4 and Table 5. Similar conclusions can be drawn.

Table 2. Summary of the accuracy and effectiveness of our error detection scheme for m = 100, n = 100 and various choices of η and γ.

Table 4. Comparison between MSTmap and Record for m = 300,500 and n = 100.

Table 5. Comparison between MSTmap and Record for m = 100,300,500 and n = 200.

In Table 2, we also compare the quality of the final maps under different objective functions. The objective functions (SARF) and (ML) give very similar results. Similar results are observed for other choices of n and m (data not shown).

Evaluation of the Accuracy of the Ordering

In the fourth and final evaluation, we use simulated data to compare our tool against several commonly used tools including JoinMap [5], Carthagene [12] and Record [4]. JoinMap is a commercial software that is widely used in the scientific community. It implements two algorithms for genetic map construction, one is based on regression [3] whereas the other based on maximum likelihood [5]. Our experimental results for JoinMap are obtained with the “maximum likelihood based algorithm” since it is orders of magnitude faster than the “regression based algorithm” (see the manual of JoinMap for more details). Due to the fact that JoinMap is GUI-based (non-scriptable), we were able to collect statistics for only a few datasets. Carthagene and Record on the other hand are both scriptable, which allows us to carry out more extensive comparisons. However, due to the slowness of Carthagene (when n = 300, it takes more than several hours to finish), we applied it only to small data sets (n = 100). The most complete comparison was carried out between MSTmap and Record.

As we have done in the previous subsection, we employ the number of erroneous pairs to compare the quality of the maps obtained by different tools. The results for n = 100 and m = 100 are summarized in Table 3. A more thorough comparison of MSTmap and Record is presented in Table 4 and Table 5. Several observations are in order. First, MSTmap constructs significantly better maps than the other tools when the input data are noisy. When the data are clean and contain many missing observations (i.e., η = 0 and γ is large), Carthagene produces maps which are slightly more accurate than those by MSTmap. However, if we knew the data were clean, by turning off the error-detection in MSTmap we would obtain maps of comparable quality to Carthagene in a much shorter running time. Please refer to the “_wp_′ no err” column for the E statistics of MSTmap when the error detection feature is turned off. Second, Carthagene appears to be better than Record when the data are clean (η = 0). When the data are noisy, Record constructs more accurate maps than Carthagene. Third, MSTmap and Record are both very efficient in terms of running time, and they are much faster than Carthagene. A clearer comparison of the running time between MSTmap and Record is presented in Figure 4. The figure shows that MSTmap is even faster than Record when the data set contains no errors. However as the input data set becomes noisier, the running time for MSTmap also increases. This is because our adaptive error detection scheme needs more iterations to identify erroneous observations, and consequently takes more time. However, this lengthened execution does pay off with a significantly more accurate map. Fourth and last, Table 3, 4 and 5 show that the overall quality of the maps produced by MSTmap is usually very high. In most scenarios, the τ statistic is greater than 0.99.

Table 3. Comparison between MSTmap, JoinMap, Carthagene and Record for n = 100 and m = 100.

Figure 4. Running time of MSTmap and Record with respect to error rate or missing rate or error and missing rate.

Every point in the graph is an average of 30 runs. The lines “missing only” correspond to data sets with no error (η = 0, γ is on the _x_-axis). Similarly, lines “error only” correspond to data sets with no missing (γ = 0, η is on the _x_-axis), and lines “error and missing” correspond to data sets with equal missing rate and error rate (η = γ is on the _x_-axis).

An extensive comparison of MSTmap and Record for other choices of m and n is presented in Table 4 and Table 5. Notice that even without error detection, MSTmap is more accurate than Record.

We have also compared MSTmap, Record, JoinMap and Carthagene on real genotyping data for the barley project. We carried out several rounds of cleaning the input data after inspecting the output from MSTmap (in particular, we focused on the list of suspicious markers and genotype calls reported by MSTmap), then the data set was fed into MSTmap, Record, JoinMap and Carthagene. The results show that the genetic linkage maps obtained by MSTmap and JoinMap are identical in terms of marker orders. MSTmap, Carthagene and Record differ only at the places where there are missing observations. At those locations, MSTmap groups markers in the same bin, while Carthagene and Record split them into two or more bins (at a very short distance, usually less than 0.1 cm).

Conclusion

We presented a novel method to cluster and order genetic markers from genotyping data obtained from several population types including doubled haploid, backcross, haploid and recombinant inbred line. The method is based on solid theoretical foundations and as a result is computationally very efficient. It also gracefully handles missing observations and is capable of tolerating some genotyping errors. The proposed method has been implemented into a software tool named MSTMap, which is freely available in the public domain at http://www.cs.ucr.edu/~yonghui/mstmap.html. According to our extensive studies using simulated data, as well as results obtained using a real data set from barley, MSTMap outperforms the best tools currently available, particularly when the input data are noisy or incomplete.

The next computational challenge ahead of us involves the problem of integrating multiple maps. Nowadays, it is increasingly common to have multiple genetic linkage maps for the same organism, usually from a different set of markers obtained with a variety of genotyping technologies. When multiple genetic linkage maps are available for the same organism it is often desirable to integrate them into one single consensus map, which incorporates all the markers and ideally is consistent with each individual map. The problem of constructing a consensus map from multiple individual maps remains a computationally challenging and interesting research topic.

Supporting Information

Text S1

Supplementary Text: Efficient and Accurate Construction of Genetic Linkage Maps.

(0.08 MB PDF)

Acknowledgments

We thank the anonymous reviewers for valuable comments that helped improve the manuscript.

Footnotes

The authors have declared that no competing interests exist.

This project was supported in part by NSF IIS-0447773, NSF DBI-0321756 and USDA CSREES Barley-CAP (visit http://barleycap.org/ for more informations on this project). Funding was used to collect genotyping data and to support graduate student Yonghui Wu and post-doc Prasanna Bhat.

References

• 1.Sturtevant AH. The linear arrangement of six sex-linked factors in drosophila, as shown by their mode of association. J Exp Zool. 1913;14:43–59. [Google Scholar]
• 2.Sun Z, Wang Z, Tu J, Zhang J, Yu F, et al. An ultradense genetic recombination map for Brassica napus, consisting of 13551 srap markers. Theor Appl Genet. 2007;114:1305–1317. doi: 10.1007/s00122-006-0483-z. [DOI] [PubMed] [Google Scholar]
• 3.Stam P. Construction of integrated genetic linkage maps by means of a new computer package: Joinmap. The Plant Journal. 1993;3:739–744. [Google Scholar]
• 4.Os HV, Stam P, Visser RGF, Eck HJV. RECORD: a novel method for ordering loci on a genetic linkage map. Theor Appl Genet. 2005;112:30–40. doi: 10.1007/s00122-005-0097-x. [DOI] [PubMed] [Google Scholar]
• 5.Jansen J, de Jong AG, van Ooijen JW. Constructing dense genetic linkage maps. Theor Appl Genet. 2001;102:1113–1122. [Google Scholar]
• 6.Cartwright DA, Troggio M, Velasco R, Gutin A. Genetic mapping in the presence of genotyping errors. Genetics. 2007;174:2521–2527. doi: 10.1534/genetics.106.063982. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 7.Weeks D, Lange K. Preliminary ranking procedures for multilocus ordering. Genomics. 1987;1:236–42. doi: 10.1016/0888-7543(87)90050-4. [DOI] [PubMed] [Google Scholar]
• 8.Falk CT. Preliminary ordering of multiple linked loci using pairwise linkage data. Genet Epidemiol. 1992;9:367–375. doi: 10.1002/gepi.1370090507. [DOI] [PubMed] [Google Scholar]
• 9.Wilson SR. A major simplification in the preliminary ordering of linked loci. Genet Epidemiol. 1988;5:75–80. doi: 10.1002/gepi.1370050203. [DOI] [PubMed] [Google Scholar]
• 10.Alizadeh F, Karp RM, Newberg LA, Weisser DK. Physical mapping of chromosomes: a combinatorial problem in molecular biology. Proceedings of SODA. 1993. pp. 371–381.
• 11.Alizadeh F, Karp RM, Weisser DK, Zweig G. Physicalmapping of chromosomes using unique probes. Proceedings of SODA. 1994. pp. 489–500.
• 12.Schiex T, Gaspin C. CARTHAGENE: Constructing and joining maximum likelihood genetic maps. Proceedings of ISMB. 1997. pp. 258–267. [PubMed]
• 13.Liu B. The gene ordering problem: an analog of the traveling sales man problem. Plant Genome 1995 [Google Scholar]
• 14.Ben-Dor A, Chor B, Pelleg D. RHO–radiation hybrid ordering. Genome Res. 2000;10:365–378. doi: 10.1101/gr.10.3.365. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 15.Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by simulated annealing. Science. 1983;220:671–680. doi: 10.1126/science.220.4598.671. [DOI] [PubMed] [Google Scholar]
• 16.Goldberg DE. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley Professional; 1989. [Google Scholar]
• 17.Glover F. Tabu search-part I. ORSA Journal on Computing. 1989;1:190–206. [Google Scholar]
• 18.Glover F. Tabu search-part II. ORSA Journal on Computing. 1990;2:4–31. [Google Scholar]
• 19.Lin S, Kernighan B. An effective heuristic algorithm for the traveling sales man problem. Operation research. 1973;21:498–516. [Google Scholar]
• 20.de Givry S, Bouchez M, Chabrier P, Milan D, Schiex T. CARTHAGENE: multipopulation integrated genetic and radiation hybrid mapping. Bioinformatics. 2004;21:1703–1704. doi: 10.1093/bioinformatics/bti222. [DOI] [PubMed] [Google Scholar]
• 21.Iwata H, Ninomiya S. AntMap: constructing genetic linkage maps using an ant colony optimization algorithm. Breeding Science. 2006;56:371–377. [Google Scholar]
• 22.Gaspin C, Schiex T. Genetic algorithms for genetic mapping. Lect Notes Comput Sci. 1998;1363:145–155. [Google Scholar]
• 23.Mester D, Ronin Y, Minkov D, Nevo E, Korol A. Constructing large-scale geneticmaps using an evolutionary strategy algorithm. Lect Notes Comput Sci. 1998;1363:145–155. [Google Scholar]
• 24.Broman KW. The genomes of recombinant inbred lines. Genetics. 2005;169:1133–1146. doi: 10.1534/genetics.104.035212. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 25.Lander ES, Green P. Construction of multilocus genetic linkage maps in humans. PNAS. 1987;84:2363–2367. doi: 10.1073/pnas.84.8.2363. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 26.van Os H, Stam P, Visser RGF, van Eck HJ. Smooth: a statistical method for successful removal of genotyping errors from high-density genetic linkage data. Theor Appl Genets. 2005;112:187–194. doi: 10.1007/s00122-005-0124-y. [DOI] [PubMed] [Google Scholar]
• 27.Lincoln SE, Lander ES. Systematic detection of errors in genetic linkage data. Genomics. 1992;14:604–610. doi: 10.1016/s0888-7543(05)80158-2. [DOI] [PubMed] [Google Scholar]
• 28.Haldane JBS. The combination of linkage values and the calculation of distance between loci of linked factors. J Genet. 1919;8:299–309. [Google Scholar]
• 29.Kosambi DD. The estimation of map distances from recombination values. Ann Eugen. 1944;12:172–175. [Google Scholar]
• 30.Liu B. Statistical Genomics: Linkage mapping and QTL analysis. CRC Press; 1998. [Google Scholar]
• 31.Liu W, Zheng MY, Polle EA, Konzak CF. Highly efficient doubled-haploid production in wheat (Triticum aestivum L.) via induced microspore embryogenesis. Crop Science. 2002;42:686–692. [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Text S1

Supplementary Text: Efficient and Accurate Construction of Genetic Linkage Maps.

(0.08 MB PDF)