ARTICLE

TOP ABSTRACT ARTICLE REFERENCES

The complete DNA sequence of a genome is a powerful tool for studying an organism. Biological research in the 21st century will surely require obtaining the sequence of large numbers of important organisms, including many higher animals and plants with large genomes.

Various approaches have been proposed for sequencing large genomes, which broadly fall into two categories: (1) whole-genome shotgun sequencing and (2) clone-based shotgun sequencing.

Shotgun sequencing involves breaking a target into random fragments, sequencing these fragments, and reconstructing the full sequence from these pieces. Shotgun sequencing was invented by Sanger, who applied it to the genome of bacteriophage  (Sanger et al. 1980, 1982). It was subsequently extended to genomes of large recombinant plasmids, large viruses, mitochondria, chloroplasts, and bacteria (Ohyama et al. 1986; Goebel et al. 1990; Oda et al. 1992; Fleischmann et al. 1995). More recently, Weber and Myers (1997) proposed that the approach could be extended to very large genomes, such as the human genome, by making greater use of long-range linking information, for example, sequences at the opposite ends of large-insert clones. There are many potential pitfalls, such as the possibility that the huge number of repeats (e.g., one million copies of the Alu repeat in the human genome) may result in many large-scale sequence misassemblies that are hard to detect and correct (Green 1997).

Clone-based shotgun-sequencing involves obtaining a collection of large-insert clones covering a genome and performing shotgun sequencing on each clone. This approach was used for sequencing the genomes of the yeast Saccharomyces cerevisiae (Oliver et al. 1992; Dujon 1996) and Caenorhabditis elegans (The C. elegans Sequencing Consortium 1998) and is being used for sequencing of mammalian genomes such as the human and mouse. Clone-based sequencing has the advantage that it eliminates the possibility of large-scale misassemblies, because each clone is assembled individually, but it requires that one generate the collection of clones covering the genome. It also has the drawback that there is some redundant sequencing as a result of the overlap of adjacent clones. (There is no workable technique for shotgun sequencing, only the nonoverlapping portion of a clone. Each clone must be subjected to full shotgun sequencing, with the overlapping portion therefore being shotgun-sequenced twice. Such duplication can be avoided in the directed finishing phase, in which gaps and ambiguities in the shotgun sequence are resolved. Such overlaps thus increase the cost for the overlapping region by somewhat less than twofold.)

The current cloning system of choice for clone-based sequencing is the bacterial artifical chromosome (BAC) with inserts of ~200 kb. The discussion below will refer to BACs (rather than generic "large-insert clones"), but the results, of course, are completely general.

Two approaches have been proposed for generating the BAC clones to be sequenced: (1) physical mapping (or "map first, sequence second") and (2) walking (or "map as you go").

Physical mapping involves constructing a complete physical map covering the genome before beginning sequencing (map first, sequence second). One "fingerprints" a BAC library, using a technique such as restriction digestion to characterize each clone, and then attempts to use this information to infer the order of the clones. A "path" of clones is then selected for sequencing. This approach was successfully used to generate physical maps of cosmids used to sequence the S. cerevisiae and C. elegans genomes (Coulson et al. 1986; Olson et al. 1986). It is being applied to the human genome, although the problem is more challenging because of the larger clones sizes (yielding more complex fingerprints), the larger genome size, and the more complex repeat structure.

The walking approach proceeds directly to sequencing without a prior physical map: One starts by sequencing an initial collection of random clones and then "walks" the genome by iteratively selecting minimally overlapping clones. Venter et al. (1996) proposed an efficient method for selecting minimally overlapping clones, which is commonly referred to as "BAC-end sequencing" or, more formally, as "sequence-tagged connectors" (STCs). A BAC library is characterized by sequencing the inserts of each clone at its two ends (using primers located in the two vector arms), and a database of the resulting BAC-end sequences is created. Given a fully sequenced BAC clone C, one can walk by searching the database to identify all overlapping BACs. The location and orientation of each overlapping BAC is immediately apparent from the position of its end sequence within C (Fig. 1). One can also tell whether an overlapping BAC closes a gap in the genomic sequence by whether its other end sequence lies in another fully sequenced BAC clone C'.

View larger version (28K): [in this window] [in a new window]   Figure 1   Overlapping BACs in library with depth d = 12-fold. A typical BAC (shown in black) will tend to have 12 BACs overlapping its right end and 12 BACs overlapping its left end. Each overlapping clone can be identified because one of its end sequences (shown by arrowheads) is contained within the sequence of the target BAC. The precise position of each overlapping clone is inferred by the location and orientation of the end sequence within the sequence of the target BAC.

The purpose of this paper is to analyze the strategy of sequencing a genome by clone-based walking. Although the basic notion of walking is straightforward, key strategic questions have not been addressed: How many independent walks should be performed in parallel? With too few, one will need too many walking steps to cover the genome in a reasonable time. With too many, one may perform too much redundant sequencing as walks "bump" into one another with randomand therefore suboptimaloverlap. How can one maximize the speed of covering the genome while minimizing the amount of redundant sequencing?

The purpose of this paper is to present a mathematical analysis that expresses the amount of redundant sequencing in terms of the initial density of walks and the depth of the BAC library. We also describe and analyze a variant strategy in which one uses a second BAC library with smaller insert clones to close gaps efficiently; this strategy dramatically reduces the amount of redundant sequencing. The results provide direct guidance for planning a genome sequencing project.

Basic Model

We begin by describing the basic model to be analyzed.

BAC Library

Sequencing

i

View larger version (5K): [in this window] [in a new window]   Figure 2   Serial walking of the genome from a single initial clone C0.

1

1

View larger version (3K): [in this window] [in a new window]   Figure 3   Serial walking of the genome from a collection of seed clones (shown in black). Bidirectional walking steps are shown, with steps that close an ocean shown in grey. The oceans are closed with relatively little overlap in the first two cases but with large overlap in the third case [marked with an asterisk (*)].

Mathematical Analysis

Our goal is to derive simple formulas providing a good approximation for the proportion of excess sequencing. Toward this end, we will make certain simplifying assumptions. The reasonableness of the approximations will then be demonstrated by their close agreement with simulations, below.

We start with a collection of nonoverlapping seed clones. Following the terminology of Lander and Waterman (1988), each clone is an "island" followed by an "ocean" to be walked. Walking proceeds outward from each clone in a bidirectional fashion, with overlapping clones recognized on the basis of their end sequence. At each round, one identifies instances in which two islands can be joined by a single clone (readily apparent from the two end sequences) and ensures that a walking step is taken from only one of the two islands (to avoid unnecessary excess sequencing).

Although walking proceeds bidirectionally in practice, it simplifies the discussion to suppose that walking proceeds unidirectionally to the right. Each bidirectional walking step is clearly equivalent to two unidirectional walking steps. Figure 4 shows a typical seed clone C₀ followed by consecutive walking steps C₁, C₂, ..., C_j, where C_j is the first clone that overlaps the seed clone to the right of C₀ (a fact that is readily apparent from its end sequence).

View larger version (6K): [in this window] [in a new window]   Figure 4   Unidirectional walking from a seed clone C0. Each successive clone Ci extends the island, with the redundant sequence shown in grey. The last clone Cj closes the ocean, with redundant sequencing due to both overlap with Cj-1 and with the next seed clone to the right. Each clone has size 1, and the ocean has size X. The number of clones sequenced is j + 1, with each having length 1. The amount of unique sequence obtained (the initial island together with the ocean) is X + 1.

Suppose that the seeds cover proportion  of the genome in clones of length 1. The oceans therefore cover proportion 1  , and the mean ocean size must be  = (1  ) / .

Our first simplifying assumption concerns ocean lengths:
Assumption of exponential oceans (AEO). Nonoverlapping seed clones will be assumed to be chosen such that the resulting ocean lengths follow an exponential distribution (with mean ).

What is the proportion of redundant sequencing entailed in sequencing the island C₀ together with the ocean on its right? Let the random variable J denote the sum of the lengths of C₁, C₂, ..., C_j. (Because the clones have unit length, J = j.) The amount of total sequencing is J + 1, whereas the amount of unique sequence obtained is X + 1 (see Fig. 4). The expected proportion of redundant sequencing is thus
where E( ) denotes the expected value. As noted above, the mean ocean size E(X) is  = (1  ) / . Therefore, we simply need to calculate E(J).

Calculating E(J) precisely is complicated, but it is possible to make a good estimate by using a second simplifying assumption:
Assumption of constant overlap (ACO). Each clone C_i will be assumed to overlap the previous clone C_i-1 by exactly the expected amount of overlap, 1 / d, and thus to extend the sequence by exactly the expected amount, 1  (1 / d).

ACO greatly simplifies the mathematical analysis, owing to an elegant property of exponential distributions: If each island is extended by a constant amount, then the lengths of the remaining oceans (i.e., those that are not closed) follow the same exponential distribution as the initial oceans. This stability property allows the walking process to be analyzed by a simple recursion.

Proposition 1

d,

d,

(1d¹)/

d,

k(1  d¹)/

Proof

d,

d,

d,

d,

d,

d,

Proposition 2

Proof

1

1

Results

We now apply the results to study the problem of sequencing a genome by walking. Figure 5A shows the proportion R of redundant sequencing, as a function of the initial mean ocean length  and the BAC library depth d. If the genome is seeded to leave oceans of average length  = 1 (corresponding to initial coverage  = 50%), the proportion of redundant sequencing ranges from 32% to 29% as the BAC library depth ranges from d = 15 to d = . If the genome is seeded more sparsely so that  = 2 and  = 33%, the redundant sequencing R ranges from 23% to 18% over the same range of BAC library depth. If the oceans are still larger ( = 3 and  = 25%), the redundant sequencing R ranges from 19% to 13%.

View larger version (3637K): [in this window] [in a new window]   Figure 5   Proportion of excess sequencing. (A) Curves show the proportion of excess sequencing, R(d,) as a function of mean initial ocean size  for various library depths d = 10, 15, 25, 50, and . (B) Curves show the proportion of excess sequencing attributable to inefficient ocean closure, R*(d,). This is obtained by subtracting the proportion of excess sequencing for an optimal tiling path in a library of depth d, which is 1 / (d  1). That is, R*(d,) = R(d,)  (1 / d  1).

It is useful to compare the results with the situation of serially walking the entire genome with a minimal tiling path (which corresponds to ocean size   ). As noted above, the proportion of redundant sequencing for a minimal tiling path is 1/(d  1). By subtracting this quantity from R, we find the additional redundant sequencing R* caused by inefficient closure of oceans. Figure 5B shows the corresponding graph for R*.

The graph shows that the redundant sequencing caused by inefficient closure of oceans depends primarily on the mean ocean size but much less on the library depth d. This makes intuitive sense, because the inefficiency arises primarily from closing small oceans with clones of unit length. The availability of more clones in a deep library thus makes only a modest improvement. The component R* is decreased by only a few percentage points by going from a library with d = 10 to d = .

Figure 6 shows the number of unidirectional walking steps required to cover 90% of the genome. The number of unidirectional steps is equal to slightly more than twice the mean ocean length. Because actual walking proceeds bidirectionally, it is thus useful to restate the observation in these terms: The number of bidirectional walking steps to cover 90% of the genome is approximately equal to the mean ocean length. (The precise correspondence depends on the library depth d. For example, the number of bidirectional walking steps is ~1.1 times the mean ocean length when d = 15.)

View larger version (50K): [in this window] [in a new window]   Figure 6   Walking steps. The number of unidirectional walking steps required the cover 90% of the genome, as a function of mean initial ocean size  for various library depths d = 10, 15, 25, 50, and . In actual practice, walking occurs in both directions. The number of bidirectional walking steps is thus half as many.

If the number of steps shown in Figure 6 is doubled, the proportion of genome sequenced rises to 98%-99%. [More precisely, the proportion remaining uncovered is 1% × (1  )^-1.] Thus, the number of bidirectional walking steps to cover 98% of the genome is approximately equal to twice the mean ocean length.

From Figures 5 and 6, one can readily assess the increase in redundant sequencing entailed in covering the vast majority of the genome in a given number of steps.

Using Smaller Clones to Close Gaps

The major inefficiency in walking arises from instances in which a clone of length L = 1 must be sequenced to close a small ocean; most of the sequencing is redundant. This observation suggests a simple improvement: Use a second BAC library with smaller inserts to close small oceans.

Specifically, suppose that we have two BAC libraries in which the clones have been end-sequenced: our original library B with insert size L = 1 and depth d, as well as a second library B' with insert size L' < 1 and depth d'. We will continue to assume that walking from each seed clone proceeds unidirectionally to the right (but see below concerning this point). At each walking step, we would first search our database to see if there is a clone from B' that spans the remaining ocean (based on its end sequences), and if none is found, we would select the minimally overlapping clone from B (which either closes the ocean or extends the walk). In this manner, we would aim to select the smallest clone capable of closing the ocean.

We can adapt the mathematical analysis above to calculate the proportion of redundant sequencing. [To simplify the statement of the result, we will consider only the case in which clones from B typically extend the walk farther than clones from B'; that is, 1  (1 / d) > L'-(1 / d'). This will be true unless L' is close to 1, in which case the second library adds little anyway.]

Proposition 3

d,

0,/L'

Proof

We proceed as in the proof of proposition 1. Part i again follows by a simple application of Wald's equation. For part ii, we calculate E(J) by considering a single walking step and applying recursion. With probability p_d',(/L'), the ocean can be closed with clone from library B', resulting in a total clone length of L'. With probability p_d,  p_d',(/L'), the ocean can be closed with a clone from library B but not B', resulting in a total clone length of 1. With probability 1  p_d,, the walking step cannot close the ocean and remaining ocean having the same exponential distribution. It thus follows by recursion that E(J) = [p_d',(/L')]. L' + [p_d,  p_d',(/L')] 1 + (1  p_d,) [1 + E(J)]. The result follows by solving for E(J).

How much efficiency is gained by using a second library B' with smaller insert clones? Figure 7 shows the proportion of redundant sequencing when an initial library B with depth d = 15 is supplemented with a second library B' with clones whose inserts are half as long (L' = 1/2) at various depths d'.

View larger version (38K): [in this window] [in a new window]   Figure 7   Proportion of excess sequencing, for walking with two BAC libraries. The first library B has clones of insert size L = 1 and fixed depth d = 15, whereas the second library B' has clones of insert size L' = 1/2 and various depths d'. Walking is performed as described in the text. The curves show the proportion of excess sequencing, as a function of mean initial ocean size  and for various library depths (d = 10, 15, 25, 50, and ) for the second library B'.

The savings are substantial, for example, decreasing the proportion of redundant sequencing from ~23%-14% for mean ocean length 2 and from ~18%-12% for mean ocean length 3. The redundant sequencing is considerably closer to the best possible level obtained with a minimal tiling path, ~7.1% for a library with d = 15 (indicated by the broken line in Fig. 7). Using a library with half-size clones roughly halves the redundant sequencing R' caused by inefficient closure of oceans.

Notably, the second library B' does not need to have high depth to have a major impact. A library with depth d = 5 is only slightly less effective than a library with infinite depth. This makes intuitive sense, because even relatively low depth assures the existence of clones suitable for covering very small oceans.

The results above concern a second library B' with insert size L' = 1/2. This choice of insert size L' is nearly optimal for minimizing the amount of redundant sequencing R. Specifically, one can show by straightforward calculus that the value of L' that minimizes R only changes from 0.55 to 0.50 as the mean ocean length  goes from 1 to .

The analysis above assumes that walking from each seed clone proceeds unidirectionally to the right. The fact that walking is actually bidirectional slightly complicates the situation: If the clones for a given round of walking were all sequenced simultaneously, we would sometimes cover oceans inefficiently by walking with large clones at both ends, when using one large and one small clone would suffice (Fig. 8). The additional excess sequencing that would result from such occurences can be shown to be

[Briefly, the first term is the excess contribution to J from sequencing a large clone instead of a short clone; the second term is the proportion of oceans that should properly be closed with one large and one small clone (these are precisely the oceans that would be closed by unidirectional walking with an odd number of large clones followed by a small clone), and the denominator ( + 1) occurs as in the calculation of R in 2(i) above.] For typical values of interest (e.g., d  15, d'  5, L = 1/2, and = 1-5), the additional excess sequencing would be in the range of 3%.

View larger version (4K): [in this window] [in a new window]   Figure 8   Schematic illustration of the situation of an ocean being closed by two large clones (B and B'), when one large and one small clone (B and B'') would suffice.

It is possible to do much better than this by exploiting the fact that the BAC clones will not all be sequenced simultaneously. As shotgun sequence is obtained for each large clone B, one can automatically check the sequence to see whether any subsequent large clone B' can be replaced by an ocean-closing small clone B'' (as in Fig. 8). Provided that B' did not pass through shotgun sequencing simultaneously with B, one can thereby avoid the cost of sequencing a larger clone when a smaller clone will suffice. (One may have incurred the cost of preparing a small-insert library from B', but this is small relative to the cost of sequencing.) The additional excess sequencing is therefore only , where < 1 is the proportion of clones in a given walking step that are processed in parallel through shotgun sequencing. The proportion  depends on the workflow of the sequencing operationspecifically, on the amount of work in process (WIP) in the shotgun sequencing phase. Because the shotgun sequencing phase is relatively rapid (the elapsed time for picking, growing, preparing, and sequencing the small-insert clones is typically on the order of a week), the WIP tends to be relatively small. Even if the proportion of clones simultaneously passing through shotgun sequencing is as high as  = 10%, the additional excess sequencing owing to the occasional failure to use a small clone is only  = 0.3%. In short, the cost is small and does not substantially impact the advantage of using a smaller library. The results derived under the assumption of unidirectional walking are thus not far off.

Optimizing BAC Library Depth

What is the optimal depth of a BAC library to be used for sequencing a genome by walking? One can decrease the cost of redundant sequencing by using deeper BAC libraries, but one incurs the cost of end-sequencing a larger number of BAC clones. Beyond some point, there are diminishing returns to increasing the depth of the BAC library.

The optimal library depth depends on the relative costs of sequencing entire BAC clones versus sequencing BAC ends. With current laboratory procedures and economics, this ratio is in the neighborhood of  = 1000:1. (Roughly 4000 shotgun sequencing reactions are needed to sequence a 200-kb BAC. Two sequencing reactions are required to sequence its ends, although these reactions are considerably more expensive because sequencing directly from a BAC template requires higher concentrations of reagents.) Given the cost ratio , the optimal library depth is the value d that minimizes the sum of the costs of BAC-end sequencing (proportional to d) and redundant BAC sequence [proportional to R(d,)].

Figure 9 shows the optimal library depth when using a single BAC library (curve denoted d*). The optimal depth increases with the initial mean ocean length  or equivalently decreases with the initial proportion  of the genome covered. The optimal depth d ranges from 22.5 to 30.8 as  increases from 1 to 8, approaching an asymptotic limit of 32.8 as   . (One can show that the optimal library density approaches  + 1 as   , by using straightforward calculus to minimize the asymptotic total cost.) The fact that the optimal d is smaller for small  makes intuitive sense, because a dense library offers more limited advantage when most oceans are fairly small.

View larger version (38K): [in this window] [in a new window]   Figure 9   Optimal library depth to minimize total cost, assuming that the ratio of the cost of obtaining the complete sequence of a BAC to the cost of sequencing both ends of a BAC is  = 1000:1. The top curve (marked d*) shows the optimal depth when using a single library. The two lower curves show the optimal depths d and d' when using two library B and B', as described in the text. The optimal depths depend on the mean initial ocean size . The depth d' of the small insert library was constrained to be at least 5, to ensure that the existence of a properly situated, small-insert clone at the end of the vast majority of islands, consistent with the assumption ACO, described in the text.

Figure 9 also shows the optimal library depths d and d' when using two BAC libraries with respective insert sizes of 1 and 1/2. The optimal value of d of the large insert library is only slightly lower than in the previous case, whereas the optimal value of d' is in the neighborhood of 6-7 for relevant values of .

The optima are relatively broad. Overall, it seems reasonable to use libraries B and B' with respective depths of d = 20-25 and d' = 6-7.

Seeding the Genome

We next turn to the issue of generating a collection of nonoverlapping seed clones. The analysis above assumes that seed clones are chosen such that initial ocean sizes are exponentially distributed (AEO). How close can one come to this in practice?

The most straightforward experimental approach is to select seeds sequentially from a list of random clones. Each clone is selected and sequenced, subject only to the condition that it does not overlap (as seen from its end sequences) with a previously selected seed clone. In this fashion, one can progressively seed the genome with random, nonoverlapping clones.

This stochastic process is known in the literature as the "parking process" because it is equivalent to cars of constant size (clones) sequentially choosing random parking spots along a long street (genome). Each car is permitted to park in its chosen parking spot provided that it is not blocked by a previously parked car; otherwise, it must drive away.

The resulting "parking distribution" is relevant to many physical processes and has been well studied (Krapivsky 1992). Suppose that one has processed a list of random clones comprising a sublibrary covering the genome to depth t. The expected proportion of genome covered with seed clones will be (t), and the distribution of oceans sizes will be (x,t), where


The parking distribution of ocean sizes is very close to the exponential distribution for low coverage (t) and is still reasonably close for coverage (t) = 50% (Fig. 10A,B).

View larger version (1918K): [in this window] [in a new window]   Figure 10   Comparison of parking and exponential distributions. (A) Probability density and (B) cumulative probability distribution, for  = 0.3 and  = 0.5. The kinked curves are the parking distribution; the smooth curves are the exponential distribution with the same mean.

The expected proportion of excess sequencing under the parking distribution can be calculated as in proposition 2. The difference  between the excess sequencing expected under the exponential distribution and the parking distributions is shown in Figure 11. The parking distribution entails slightly less excess sequencing than the exponential distribution (i.e.,  > 0), because it has fewer extremely small intervals (as see from Fig. 10A). The difference  is quite small over the range of interest, being ~0.6% for = 1 and decreasing as  increases.

View larger version (17K): [in this window] [in a new window]   Figure 11   Difference between exponential and parking distributions. Graph shows difference  between Rexp, the proportion of excess sequencing under exponential distribution, and Rparking, the proportion of excess sequencing under parking distribution. That is,  = Rexp  Rparking. Two curves are shown, corresponding to d = 15 and d = 50. The difference  is seen to be largely insensitive to d.

It should be noted that the parking process cannot fully cover the genome: As the coverage increases, the remaining spaces become smaller, and the proportion of cars turned away increases. When the remaining spaces are all smaller than a car length, no more cars can park. This happens at the so-called "jamming limit," which occurs at  = 0.747597...sequential selection of strictly nonoverlapping clones thus cannot provide seed coverage beyond this point.

Fortunately, we are not interested in seeding to >50% because the cost of excess sequencing begins to skyrocket for  > 50%. Seeding the genome to 50% requires starting with a random list of clones covering the genome to depth 1.15 [because (1.15) ~0.5].

We briefly mention some other approaches to seeding the genome but do not analyze them in detail.

Simultaneous Seeding

Contig-Based Seeding

Simulations

The simple formulas above are premised on two simplifying assumptions (AEO and ACO).We tested whether the formulas provided reasonable approximations by performing extensive simulations.

The simulations were performed as follows: Libraries were generated by randomly selecting starting points for BAC clones (L = 2 × 10⁵ bp, L' = 1 ×10⁵ bp) from a mammalian genome (3 × 10⁹ bp) using a uniform distribution to achieve the desired depths d and d'. An initial selection of seed clones was generated as in the parking problem, by considering the clones in a random order and accepting successive nonoverlapping clones until the desired proportion  of the genome was reached. [We confirmed that the simulation-yielded ocean sizes closely fit the parking distribution (not shown).] Walking was then performed by selecting minimally overlapping BACs as described above. For each choice of parameters, a total of 40 simulations were performed, and the results were averaged. Simulations were performed for parameters throughout the range of interest, at the parameters d  {10, 15, 25, 50, } and  = {1.0, 1.5, 2.0, 2.5, ..., 7.0}. We calculated R_sim(d,), the average proportion of redundant sequencing over the simulations.

We also performed simulations mimicking the ACO assumption, in which the seed clones were selected as above, but walking steps then occurred under the assumption that the minimally overlapping clone overlapped by exactly 1 / d. We calculated R_sim,ACO(d,), the average proportion of redundant sequencing over these simulations.

We compared the simulation results to R_parking(d,), the excess sequence under the parking distribution predicted by the formula in proposition 2. The difference (d,) = R_parking(d,)  R_sim(d,) is shown in Figure 12A. The difference is negligible (on the order of the S.E.M. from the simulations), except for small d and small . The discrepancy for small d and  is readily seen to be due to ACO, by examining the difference '(d,w) = R_parking(d,)  R_sim,ACO(d,) shown in Figure 12B. This difference is negligible over the entire range.

View larger version (2018K): [in this window] [in a new window]   Figure 12   Difference between formulas and simulations. Plots show the differences  between the proportion of excess sequencing predicted under the parking distribution (Rparking) and the proportion of excess sequencing observed in (A) simulations in which the seed clones and walking clones are chosen from a randomly distributed library (Rsim) and (B) simulations in which seed clones are chosen from a randomly distributed library, but walking steps are then made under the ACO assumption (Rsim,ACO). Results are shown for various values of d and . () d = 10; () d = 15; () d = 25; () d = 50; (gray ) d = .

We also studied the proportion of genome covered after k steps (k = 1, 2, ..., 10), comparing the prediction under the exponential distribution [proposition 1(iii)] to the results seen in the simulations above. In the range of interest, the predicted proportion of genome covered at each stage was within 0.5% of the results from the simulations (not shown). Finally, we simulated the situation of two libraries, in which d and  were as above and the smaller library has depth d' = 5. Again, the difference  was extremely small (not shown).

In summary, the simple formulas derived under the assumption of AEO and ACO are sufficiently close for practical purposes. The differences are quite small and have no impact on the key conclusions.

Conclusion

Sequencing a genome by walking is an inherently serial process. One seeds the genome to an initial coverage  and then engages in sequential rounds of walking.

Completing the task in a reasonable amount of time requires a high degree of parallelism, that is, a high density of seed clones from which walks proceed outward in both directions. The number of such walking steps to cover 90% of the genome is roughly equal to the initial mean ocean size [ = (1  ) / ], whereas the number of steps to cover 98% is roughly twice as large.

There is a clear trade-off between the number of walking steps and the cost of redundant sequencing. As the mean ocean size  grows from 1 to 2 to 3, the proportion of redundant sequencing decreases from 32% to 23% to 19%. (These values correspond to a library with depth d = 15; they are about two percentage points lower for libraries with d = 25.)

One solution is to decrease the cycle time required for each walking step, making it feasible to take more serial walking steps within the desired time frame. Because the time required to prepare and validate high-quality shotgun libraries represents a significant component of the cycle time, one might develop efficient and inexpensive methods to prepare shotgun libraries from all clones in the BAC library in advance of sequencing. (The storage requirements are modest, because each shotgun library is stored as a ligation mixture in a single test tube until required.) One could then rapidly initiate each consecutive walking step. Implementing such an approach would require significant streamlining or automation of existing procedures for preparation of shotgun libraries but may be feasible.

A complementary solution is to use a second BAC library with smaller clones, with roughly half the insert size. This approach dramatically improves efficiency. For the cases involving a 15-fold library above, the proportion of redundant sequencing is 18%, 14%, and 12%. This is much closer to the best possible result obtainable with a 15-fold library of ~7.1% (= 1/14). Importantly, the vast majority of the efficiency is obtained using a second BAC library having relatively modest depth, for example, d' = 5.

The strategy could, of course, be generalized. One could use multiple BAC libraries each with distinct insert sizes or a single BAC library with extremely variable inserts sizes. Mathematical analyses of such strategies would be of interest, although it is notable that the simple two-tiered strategy above already eliminates the majority of the excess redundant sequencing due to inefficent closure of oceans.

We should note that our analysis ignores certain biological and experimental issues:

Variable Insert Size

Cloning Bias

Missing End Sequences

Repeat Sequences at BAC Ends