Breaking Lander-Waterman’s Coverage Bound

Lander-Waterman’s coverage bound establishes the total number of reads required to cover the whole genome of size G bases. In fact, their bound is a direct consequence of the well-known solution to the coupon collector’s problem which proves that for such genome, the total number of bases to be sequenced should be O(G ln G). Although the result leads to a tight bound, it is based on a tacit assumption that the set of reads are first collected through a sequencing process and then are processed through a computation process, i.e., there are two different machines: one for sequencing and one for processing. In this paper, we present a significant improvement compared to Lander-Waterman’s result and prove that by combining the sequencing and computing processes, one can re-sequence the whole genome with as low as O(G) sequenced bases in total. Our approach also dramatically reduces the required computational power for the combined process. Simulation results are performed on real genomes with different sequencing error rates. The results support our theory predicting the log G improvement on coverage bound and corresponding reduction in the total number of bases required to be sequenced.


Introduction
Data generated from DNA sequencing machines are growing at an unprecedented rate. Extracting knowledge from these data is extremely tedious and usually requires very powerful computing machines. The main reason is that the volume of data generated for an experiment usually contains redundant data and one needs to pay the price of extracting useful information and removing redundant information at the processing step. As an example, in the whole genome sequencing of Human genome with 100x coverage, each base averagely is present in 100 sequencing reads which means 99 percent of the data is redundant. The first question that comes in mind is whether the volume of data generated by the sequencing machines can be reduced without affecting the overall performance. In this paper, we focus on the whole genome sequencing problem and seek fundamental results on the redundancy level required to obtain the desired result.
The first fundamental result in this area has been due to Lander and Waterman in [1] where they present a lower bound on the number of reads, N, required to assemble the whole genome.
We refer to this bound as coverage bound. The coverage bound states that for a genome of size G and reads of length L, at least N cov ¼ G L ln G L À Á reads are needed such that the whole genome is covered with a probability of no less than 1 − [2]. Therefore, we should have N ! N cov .
For the aforementioned scenario, the total number of bases sequenced by the sequencing machine is NL, which requires to be of the order of G log G, from Lander-Waterman's result. Consequently, the non-reducible redundancy level in such setup will be of the order of log G. However, such result is based on the underlying assumption that sequencing and computing steps are performed independently, i.e., a machine takes samples from the genome and sequences as many reads required to cover the whole genome and then another machine processes the reads to assemble the genome. The overall architecture of such approach for whole genome sequencing is shown in Fig 1(a) and consists of the cascade of two blocks one for sampling and sequencing and one for assembly.
Although such separation between sequencing and processing has been traditionally assumed in the literature, one can question whether such separation is fundamentally optimal with respect to amount of sequenced data which is generated and subsequently processed or not. In other words, can we improve the overall performance of such system by merging the two components? In fact, in order to verify the level of performance improvement achieved by such integration, one needs to answer the following two key questions. First, is there any improvement on lowering the redundancy of generated data by merging the two functions? Second, is it physically possible to build such a machine to perform both functions Whole genome shotgun sequencing (a) and our method platforms (b). In the classic method (a), the processing starts after termination of random sampling and sequencing. In our proposed sequencing platform (b), sampling and processing machines work cooperatively. The results of alignment are fed back to preempt reading redundant data in the sequencing machine.
simultaneously? In the rest of this paper, we will try to answer the first question by proving that one can break the coverage bound of Lander and Waterman and reduce the number of sequenced bases to as low as O(G). It is worth mentioning that even if we may not have access to a machine that can efficiently combine the sampling and computing units, the approach proposed in this paper can still reduce the computational power required to assemble the genome. This is due to the fact that the processing unit processes the data on the fly and if it detects that the information from the remaining part of data is redundant it will stop further processing of that piece of data. As we will show, such early termination can have significant effect on reducing computational complexity of the whole process. Note that the proposed method can control depth of sequencing (for obtaining the desired accuracy) by sequencing only informative bases and reducing over-read bases. Using this capability, also for noisy reads, we can obtain O(G) read bases. In addition, at the worst case, our method obtains an accuracy which is the same as the Lander-Waterman's method.
Answering the second question is beyond the scope of this paper. However, our approach clearly shows that if a sequencing machine can be built that can preempt sequencing at the instant the computation part sees best fit, a much more efficient sequencing machine may actually be obtained.
Our strategy to merge the two functions is shown in Fig 1(b) where the processing machine controls the sequencing machine by blocking sequencing of redundant bases. Hence, we assume that the sequencing machine sequences the DNA fragments base by base and it will stop sequencing a fragment once a blocking command from the processing machine is initiated for that fragment.
In this paper, we only focus on the re-sequencing problem in the processing machine. We first present theoretical results for i.i.d. genome and real genomes with given repeat structures, with noiseless and noisy reads. Our simulations are performed on chr19 of Human genome hg19 for both noiseless and noisy reads. We have shown significant improvement on coverage bound for real genome is achieved by using this method. For processing machine, Meta-aligner [3] scheme is used to align reads on the reference genome.
The longer the read lengths, more reduction on number of read bases will be achieved by our method. A number of Next Generation Sequencing (NGS) methods [4], such as PacBio and Nanopore [5][6][7] already provide reads of several thousand bases long and are suitable candidates for such analysis. Fig 2 shows read length distribution of PacBio technology for the first two read archive in NCBI GenBank SRX533609 (http://www.ncbi.nlm.nih.gov/sra?term= SRX533609.). It can be envisioned that other sequencing methods might also provide longer reads as sequencing technology further advances in that direction in coming years.

Materials and Methods
Conventionally, in the re-sequencing problem, a reference genome and a set of reads from a target genome are available at the processing step. In this framework, the sequencing machine produces reads of length L from N DNA fragments. Due to Lander-Waterman's coverage bound, NL, the total number of bases read by the machine, is required to be O(G log G).
Our method changes this bound by assuming that sequencing can be controlled by a processing machine such that it can terminate sequencing at any base. In other words, the sequencer starts reading bases of one side of a DNA fragment one by one and it will stop as soon as a command is initiated by the processing machine. In order to explain the basic ideas and to prove that sequencing can be performed efficiently, we first analyse our proposed methods on i.i.d. genomes. We extend the results to real genomes where repeats play an important role in their structure.

i.i.d. genomes
In this part, we assume that the reference and target genomes are i.i.d. random sequences of {A, C, G, T} with uniform base probabilities. We also assume that reads are sampled uniformly and independently from the target genome. Consider that reads are noiseless. The key strategy that we use to terminate the sequencing process in a controlled manner is as follows. We divide L by some integer number K 2 {0, . . ., L} and without loss of generality assume that ℓ = L/K is also an integer. We allow the reading machine to read the first ℓ bases of all the DNA fragments. Let R 1 ¼ fR 1 ð'Þ; Á Á Á ; R N ð'Þg denote the set of starting ℓ bases of all the fragments. Here, R i (ℓ) is the first ℓ bases of the i th fragment.
After generating R 1 , all reads are mapped to the reference genome. Some of the the reads can be mapped uniquely to a location on the genome. We call such reads anchored. More precisely, a read R is assumed to be anchored if there is only one location on the genome with Hamming distance no more than α|R| where |R| is the read length and α is some fixed constant.
After mapping, we partition the set R 1 into three disjoint sets: R C 1 the set of reads that are anchored to some location on the genome and in addition, extending them does not increase the coverage, R A 1 the set of reads anchored to some location on the genome and in addition, extending them will increase the coverage, and R F 1 the set of reads that are not anchored in the first step. For a read R i (ℓ) in R C 1 a termination command is initiated to stop further reading of the i th fragment. The union of R A 1 and R F 1 is denoted by R 2 which is the set of fragments where reading process will be continued on them.
Subsequently, the next base of all fragments in R 2 are read and we use the same procedure for mapping and termination. Therefore, at the end of this step, we end up with the set R C 2 of anchored and terminated fragments with length ℓ + 1 and the set R 3 that is used for extension in the next step. In this way, one can proceed up to step L − ℓ + 1 where all the fragments are extended to the maximum length L.
If we denote the set of reads that are uniquely mapped in the algorithm by O, then Our proposed algorithm is then detailed in Algorithm 1.

Algorithm 1
Input: N fragments with size L i of a target genome plus a reference genome of length G. Output: A set of reads O, mapped to the reference genome.

Initiate:
Set R 1 to be the set of all fragments. 1: for k = 1 to L − ℓ + 1 do 2: if k = 1 then 3: Sequence the first ℓ bases of all reads in R k . 4: else 5: Sequence the (ℓ+k − 1)-th base of all reads in R k . 6: end if 7: Map all reads in R k to the reference genome with their last ℓ fragments and Hamming distance d = bαℓc. 8: Add uniquely mapped reads in R k to the set R A . Put the rest of reads in the set R kþ1 . 9: Add reads in R A to O, if by further extensions they will not cover a new base on the reference genome. 10: Add reads in R A to R kþ1 , if by further extensions they will cover new bases on the reference genome.

11: end for
The two parameters of the algorithm, i.e., ℓ and α, should be specified based on the structure of reference genome as well as, G, N and L. One can choose ℓ to be as low as 1. However, the chance of finding uniquely mapped reads is very small for short reads resulting in much higher processing time. On the other hand, if ℓ is chosen to be large, then a lot of reads will overlap after the first step of the algorithm and we will sample a lot of redundant bases. Therefore, an optimal choice of ℓ is desired. The optimal value of ℓ depends on the size and repeat structure of the genome as well as the statistics of the variations between target and reference genomes. In the following, we describe selection of these parameters for noiseless reads.
It can be shown that for a given random DNA sequence of length G, the probability of observing two exact copies of a substring with length ℓ is lower than G 2 4 −ℓ [2]. Hence, a noiseless read of length ℓ > log G from target genome almost certainly can be uniquely mapped to the reference genome. Conversely, a noiseless read of length ℓ < log G from target genome will be mapped to at least two locations on the reference genome. Therefore, for noiseless reads with no variation between target and reference genomes, we can choose ℓ = log G.
The choice of α depends on the mapper quality and allowable Hamming distance between reads and the reference genome. If we assume a perfect mapper, the Hamming distance between a read and its true location is ν|R| where ν is the variation rate between target and reference genomes. Therefore, it suffices to set α = ν. In this scenario, we consider that variation between reference and target genome is negligible.
In order to evaluate the performance of the algorithm, we need to prove that the genome can be completely covered by the reads in O and the number of bases read more than once is small. After the first step of the algorithm, i.e. sequencing ℓ bases of all reads, all reads are anchored to their correct location on the genome with maximum Hamming distance of d = 0. This is due to the preceding discussion in the case of random genomes where duplicate segments are scarce. We distinguish between two cases: 1. Two subsequent reads have common bases, such as i th read and (i + 1) th read in Fig 3. 2. Two subsequent reads do not have a common base, such as j th read and (j + 1) th read in Fig 3. In the first case, we will encounter reads whose extension does not increase the coverage and therefore, their sequencing should be terminated. These kind of reads belong to R. We call the bases common between i th and (i + 1) th read in such case, over-read bases. In the second case, we will encounter reads whose extension will increase the coverage. We put reads of this kind in R A , to be further extended in the subsequent steps of the algorithm. Consequently, sequencing of an aligned read in R A continues until it becomes a member of O.
For further analysis, we assume that starting point of reads are a Poisson point process with rate l ¼ N G . Hence, the inter-arrival locations have independent exponential distributions. Let B i be a random variable representing the number of extra bases read by the sequencing machine for the i th read. If B denotes the total number of over-read bases, then B ¼ To compute EfB j g, note that the next read of the j th read starts x bases after the j th read. If x ℓ, Fig 4. Since x has an exponential distribution, we obtain In the case of i.i.d. genomes, we set ℓ = log G. Hence, the average number of over-read bases of all reads is as follow: For a constant value κ, the order of EfBg becomes O(G) and number of reads becomes N ¼ k G log G . To be able to cover the genome with N reads, we should be able to close gaps of at least G log G/N bases (from the coverage bound). Therefore, length of reads We observe that by tuning κ, one can control maximum read lengths as well as total number of bases read by the machine. For instance, by choosing κ = 1, we obtain EfBg ¼ 0:36G and the maximum read length becomes L % 1000 bps.
Next, we consider noisy reads with the error rate . Moreover, we assume that the target and reference genomes vary in their sequences with the variation rate ν. Before presenting our algorithm for noisy reads, we set a parameter denoted by C for the coverage depth. The coverage depth helps in removing sequencing errors by averaging over several reads. The coverage depth is chosen such that if C reads cover a base then it is possible to correctly recover that base with a given probability P .
To compute C , we choose a base as the target base if it achieves the highest number of votes over reads covering that location. Let I i denote the error event of incorrect calling of the i th base. If the random variable C i denotes the coverage of the i th base, then error occurs if the corresponding base in more reads are incorrect. Thus, the probability of error for the i th base is, Hence, For example, if P = 10 −4 and = 0.07, C becomes 6. For having coverage depth of C , it suffices to change the k th step of Algorithm 1 as follows: we add a read in R A to R kþ1 when its extension will cover a base c times on the reference genome, where c C . Details of the proposed algorithm for noisy reads is presented in Algorithm 2.

Algorithm 2
Input: N fragments with size L i of a target genome with G bases plus a reference genome with the same length. Output: A set of reads O, mapped to the reference genome.

Initiate:
Sequence the first ℓ bases of all the reads in R k . 4: else 5: Sequence the (ℓ+k − 1)-th base of all the reads in R k . 6: end if 7: Map all the reads in R k to the reference genome with their last ℓ bases and Hamming distance d max . 8: Add uniquely mapped reads in R k to the set R A . Put the rest of reads in the set R kþ1 . 9: Add reads in R A to O, if by further extensions they will cover a base more than C times on the reference genome. 10: Add reads in R A to R kþ1 , if by further extensions they will cover a base c times on the reference genome, where c C . 11: end for Because of sequencing errors, mapping fragments of length ℓ with maximum Hamming distance d max to the reference genome is error prone and we first analyze the performance of the alignment procedure. Let T i and F i denote the true and false alignment events for the i th read, respectively. The i th read of length ℓ with a maximum Hamming distance of d max is mapped to its true location with probability and is mapped to a false location on the genome with maximum Hamming distance d max with probability (using the union bound), where P w (ℓ) represents the probability of incorrect alignment of a fragment of length ℓ to another position on the reference genome. Denote, F ¼ [ N i¼1 F i as the false alignment event. Hence, , PfF g tends to zero and all reads are mapped to their true location with PfT i g's. In the worst case, if any read is extracted from each base of the reference genome, i.e. N = G, then P w (ℓ) scales as O 1 . It can be easily verified that for = 0.05, ℓ = 2log G and d max = 8, PfF g tends to zero and PfT i g % 1. Therefore for 0.05, all noisy fragments of length ℓ are uniquely mapped to their correct locations on the reference genome.
In order to compute the number of extra read bases, we use a similar argument as the one used in the noiseless case. To obtain EfB j g for the j th read in this case, assume that the C th read after this read starts x bases further. Again B j ¼ maxf0; ' À xg. Since, x has an Erlang distribution, we obtain where γ(s, x) is the lower incomplete gamma function and for s 2 N is equal to Thus, the average number of over-read bases for all reads in this step is, where k n ¼ N' G . Therefore, for any constant κ n , EfBg becomes O(G). The coverage bound when each base of the genome is covered by at least C reads can be determined as follows, where E c and E c;i are error events that at least one base and the i th base of the genome are not covered by at least C reads, respectively. Therefore, if then each base of the genome is covered by at least C reads almost surely. Subsequently, for a fix ℓ and constant value of κ n , number of reads and read length from Eq (9) become N ¼ k n G ' and L ¼ ' k n ð log G þ ðC À 1Þ log ð log GÞÞ, respectively. Again, there exists a trade-off between the length of reads and number of read bases that can be controlled by κ n .

Real Genomes
In this section, we consider DNA sequencing of real genomes where many repeats are dispersed across the genome. First, we assume that reads are noiseless. Note that, if all the ℓ-mers of the genome are repetitive elements then Algorithm 1 fails in anchoring reads correctly to the reference genome and therefore reading O(G log G) is unavoidable. However, as we will show, the repeat patterns in real genomes allow successful coverage of bases with only O(G) reading bases.
A mosaic model for capturing the repeat structure of the genome is presented in [3]. In this model, the reference genome consists of two types of intervals: repeat and random intervals. These types are defined based on two parameters ℓ and d: representing the fragment length and mismatch factor, respectively. Repeat (random) intervals are consecutive bases where any fragment of length ℓ starting from a base within these intervals can be aligned to some other location(s) (one location) of the genome with maximum Hamming distance d. For the sake of simplicity, we consider only d = 0.
Let us denote the set of all exact repeat intervals of the reference genome as R. Also, assume that a repeat R 2 R has length ℓ R and repeat lengths have the distribution f ℓ . We need to treat reads starting from repeat and random intervals, differently. For this purpose, we consider three starting regions for starting point of a given read, as S 1 , S 2 , and S 3 . We can determine the average number of total over-read bases (i.e. EfBg) as follows, where S i is a random variable that denotes the starting region of the i th read. Hence, In the following, we determine each term of Eq (11). Region S 1 consists of random intervals such that reads starting from random intervals can be anchored to their true locations based on their first ℓ bases. Therefore, we can readily compute the average number of over-read bases in random intervals using Eq (1). More precisely, where the total average number of these reads is The normalized average number of over-read bases with P = 10 −4 . The normalized total number of read bases as a function of read length for different sequencing error rates and P = 10 −4 . Note that, total number of read bases at each error rate tends to its corresponding C . doi:10.1371/journal.pone.0164888.g006

Breaking Lander-Waterman's Coverage Bound
On the other hand, the read starting from a repeat interval can not be anchored unless it contains an ℓ-mer which resides in random interval. Using this fact, c.f. Fig 7, each repeat interval of length ℓ R can be partitioned into two disjoint intervals: 1) Mappable zone: the last min{L − ℓ, ℓ R } bases, 2) Un-mappable zone: the first max{0, ℓ R − L + ℓ} bases. Regions S 2 and S 3 are mappable and un-mappable zones, respectively.
Clearly, reads from un-mappable zones cannot be anchored and therefore need to be read up to length L. Using Eq (1), we compute the average number of over-read bases for each read in un-mappable zones as: where the total average number of these reads is The next step is to compute the number of over-read bases for reads from mappable zones. Consider all mappable regions s m 2 S 2 of length l m , for m ¼ f1; Á Á Á ; jRjg. If the m th repeat interval has length ℓ m , then l m = min{L − ℓ, ℓ m }. At least l m + ℓ bases of each read within s m are being sequenced. Consider the i th read has distance l i from the end of its mappable zone. The average number of over-read bases for the i th read in mappable zones can be determined as: Thus, the average number of total over-read bases for all reads in mappable zones becomes, Define, Therefore, the average number of over-read bases by considering repeat structure in Eq (10) becomes, Given the repeat length distribution (i.e. f ℓ R ) of any real genome, we can determine the EfBg for that genome. The distribution of log f ℓ r for Human genome hg19 is illustrated in Fig  8. Also, Fig 9 shows EfBg=G for whole genome of hg19 and i.i.d. genome. In this simulation, we used ℓ = log G % 30. Results confirm that for real data set, we can read only O(G) bases to assemble the genome.
When reads are contaminated with sequencing errors of rate , we use a proper value of ℓ such that a fragment length of ℓ is aligned to its correct location with a probability close to one. Also, consider the coverage depth, i.e. C , is given by Eq (3). We model the reference genome with Meta-aligner (ℓ, d = 0)-model. Thus, using the same argument as noiseless reads, the average number of over-read bases can be determined similar to Eq (13) by incorporating V C (.) in Eq (6) instead of V 1 (.). Therefore, the average number of over-read bases for real genomes in the presence of noise becomes,

Algorithm Coverage Analysis
In this section, we compute the number of gaped bases in the reference genome when the proposed algorithms are used. First, consider i.i.d. genomes. We show that the whole genome is in fact covered by the noiseless reads when using Algorithm 1. Suppose that all reads of length ℓ are aligned to the reference genome. Let E denote the error event which is the event that a base is not covered in our algorithm. Let us denote the event of not covering the i th base by E i . Thus, From union bound, we have for arbitrary i 2 {1, Á Á Á, G}. Define the set S i , as the set consisting of starting points of reads that are aligned to the reference genome with less than L bases before the i th base location.
Since the nearest read in S i to the i th base does not overlap with other reads from its right hand side, this read will be extended in subsequent steps of the algorithm and at the end, the i th base will be covered by this read. Hence, S i must be an empty set and E i occurs when no read's starting point is located less than L bases before the i th base. This condition is the same as the coverage bound condition. Thus, Consequently, if the number of reads N and the reads' length L satisfy the coverage condition in [1] (i.e. NL ! G log G), the sequence is completely covered by the reads in our method for noiseless reads and i.i.d. genomes.
For noisy reads, we show the perfect coverage of reference genome when noisy reads are used in Algorithm 2. For this purpose, the same argument as noiseless reads is considered. Let E n denote the error event which is the event that a base is not covered by at least C reads in our algorithm. Also, denote error event for the i th base by E n;i . Therefore, E n ¼ [ G i¼1 E n;i and E n;i occurs when less than C read's starting points are located within L bases before the i th base. This condition is the same as the coverage bound condition with a given C in Eq (8). Thus, if number of reads N and reads' length L satisfy the coverage condition for noisy reads in Eq (9), the sequence is completely covered by at least C noisy reads in our method as well. Now consider a real genome. We must determine how many bases are covered with Algorithm 1 or Algorithm 2 (based on noiseless or noisy reads). For this purpose, let E k denote the error event that the k th base is not covered by reads with the proposed algorithms. We only consider coverage in this section, therefore, we use C = 1 for noisy reads. Based on the base location and its neighboring repeat intervals within the genome, this base is classified into two different classes as shown in Fig 11. Using these two classes, different sub-classes for locating random and repeat intervals can be modelled. Note that similar to i.i.d. genome, locating one read within distance of L bases before a given base is sufficient for covering that base. In the following, we determine probability of the E k for each class. In the proposed analysis, we assume that each fragment of length ℓ is mapped uniquely to the reference genome with probability p t . Also, we denote the number of reads within a random interval of length l as N ðlÞ.
1. Class A: Assume d 1 and d 2 bases within distance of L bases before and after the k th base are in random interval, respectively. If a read has a fragment within a random interval, it can be mapped to the reference genome uniquely. If d 1 + d 2 ! L, we divide the interval of length L before the k th base to three parts: 1) all bases with distance [L, d 1 ] from base k, 2) all bases with distance [d 1 , L − d 2 ] from base k, and 3) the remaining bases of the random interval with distance [L − d 2 , 0] from base k. Divide the first part to LÀ d 1 ' disjoint sub-intervals of length ℓ. If any read starts within the j th sub-interval, the number of fragments of that read within the random interval is d 1 ' þ j À 1. If any read starts within the second part, the number of fragments of that read within the random interval is L ' . In addition, divide the third part to LÀ d 2 ' sub-intervals of length ℓ such that if any read starts within the j th sub-interval, the number of fragments of that read within the random interval is L ' À j. Thus, The same result is obtained when d 1 + d 2 < L, i.e., 2. Class B: Assume d 1 and d 2 bases within distance of L bases before and after the k th base are in repeat interval, respectively. If d 1 + d 2 ! L, we divide the interval of length L before the k th base to three parts similar to class A. Thus, The same result is obtained when d 1 + d 2 < L, i.e., Breaking Lander-Waterman's Coverage Bound We are interested in error probabilities of some special cases. These cases of interest are illustrated as sub-classes in Fig 11. The error probability of each sub-class is determined in the following.
• Sub-class (I): This sub-class can be modelled with the first class with d 1 = d 2 = L. Thus, : ð22Þ • Sub-class (II): Consider that d bases within distance of L bases after the k th base is in random interval. This sub-class can be modelled with the first class with d 1 = L and d = d 2 < L. Thus, • Sub-class (III): Consider that d bases within distance of L bases before the k th base is in random interval. This sub-class can be modelled with the first class with d 2 = L and d = d 1 < L. Thus, the error probability of this class is the same as the sub-class (II) except that L − d bases within distance of L before the k th base exist in random interval.
• Sub-class (IV): Consider that d bases within distance of L bases before the k th base is in repeat interval. This sub-class can be modelled with the second class with d 2 = L and d = d 1 < L. Thus, • Sub-class (V): Consider that d bases within distance of L bases after the k th base is in repeat interval. This sub-class can be modelled with the second class with d 1 = L and d = d 2 < L. Thus, the error probability of this class is the same as the sub-class (IV) except that L − d bases exist within distance of L before the k th base in repeat interval.
• Sub-class (VI): This sub-class can be modelled with the second class with d 1 = d 2 = L. Thus, the error probability of this sub-class is 1.
Thus, by considering repeat structure of the genome, we can determine the probability of coverage for the genome. We consider Meta-aligner (30, 0)-model for the reference genome. We classify bases of Human genome hg19 and determine probability of gap for each class using Eqs (18)-(21). The average probabilities of gap, i.e. E k fPfE k gg, for different values of l ¼ lL= log G and read lengths (L) are shown in Figs 12 and 13. These probabilities of gap are shown for p t = 1 and p t = 0.7 (dotted line). Note that, the coverage bound shows that using reads of length L ! log G andl ! 1, all bases of an i.i.d. genome are covered with probability almost one.

Results and Discussion
In this section, we propose simulation results in two cases: simulated reads from the chr19 of Human genome hg19 and real reads from Human genome hg19 obtained from 454 technology.

Benchmark
In the simulated reads case, the chr19 of Human genome hg19 is used as the reference and noisy reads are also extracted uniformly from the reference genome and are mapped to this genome. We consider sequencing error rates of = {0, 5, 10}% with 90% mismatches and 10% indels. Errors are added in an i.i.d. manner. For two cases, we use Meta-aligner [3] to align reads with their two fragments of length ℓ and not using all their bases. We consider only the first stage Meta-aligner. Since with a mismatch percentage of α and read length of ℓ, there are α × ℓ bases altered in each read on the average, we allow Meta-aligner to align reads to the reference genome with a distance of dα × ℓe. For both cases, with reference size of G bases, we use ℓ = log G % 30. Therefore, N = (30 + 4(C − 1))G/L reads are randomly generated from each reference genome for any read length L and C . Also, we consider C 0 = 1, C 5 = 6, C 10 = 8.

Simulated Reads
In each simulation of the simulated reads case, we present number of aligned reads by Metaaligner. Reports of Meta-aligner show its robustness to sequencing errors such that it aligns many reads almost correctly at its first stage. It should be noted that we need to read 2ℓ bases at the first step of Meta-aligner which increases the number of over-read bases. However, most of the mapped reads are located on the genome correctly.
We first determine number of mapped reads at end of the first stage of Meta-aligner. Fig 14  shows fraction of mapped reads for = {0, 0.05, 0.1}. Results show that most of reads are mapped uniquely to the reference genome.
In Fig 15, the total number of bases not covered by mapped reads (also known as genome gaps) for various sequencing error rates after the first stage of Meta-aligner is presented. By increasing read length, a higher fraction of repeats are bridged by reads and the gap fraction decreases. Also, Fig 16 shows   that using the proposed method, the reference genome is gradually covered by the reads. Consequently, since the remaining bases of the genome are located within long repeat regions, they will be covered by the reads.

Real Reads
In this section, we obtain experimental results for the total number of read bases required to cover the whole human genome in the case of real reads which is published by Roche 454 technology. The data-set can be downloaded from NCBI with accession number SRR003161 (http://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR003161). This data-set consists of N = 31,243,790 reads with average length of 577 bps from the Human genome hg19 (this amounts to approximately 6× coverage). We use the first step of Meta-aligner to map this read-set to the Human genome hg19. Using ℓ = 30 and d = 1,25,999,865 of reads (%83.22%) are handled and the total number of read bases and the fraction of gaps are illustrated in Fig  19. Results show that only 0.6 × G bases are over-read for this read-set. Due to our strategy, the remaining fragments are completely sequenced by the sequencing machine. However, these reads cannot be mapped uniquely to the reference genome which implies they are either from repeat regions of the genome or contaminated by many errors. We use Bowtie2 to map these reads to decrease the gap fraction. Bowtie2 in default mode handles 1,671,647 (%5.35%) reads and as a result, number of gap bases decreases from 0.13 × G to 0.089 × G bases. The main reason for this gap fraction is that reading 6× coverage of the genome is not enough for covering the whole genome. Also, a small fraction of this gap is due to un-mappable reads. As these simulation results show, increasing sequencing error rate leads to an increase in the number of over-read bases. This is due to the fact that each read can not be aligned to the genome at shorter lengths and its length is increased iteratively. In addition, it should be noted that genomes with a larger percentage of repeat patterns naturally lead to a greater level of over-read bases (comparison of chr19 with i.i.d. genome). In such scenarios, less number of reads are uniquely aligned to the genome due to the ambiguity caused by repeating patterns.

Conclusion
Lander and Waterman have presented the coverage bound based on random sampling of i.i.d. DNA sequence. After sampling, read fragments are sent to the processing part. Under such model, the coverage bound shows that minimum number of reads required for covering the whole genome is N % G log G/L. Equivalently, NL % G log G bases are required to cover the whole genome. In our method, sequencing and processing are combined such that first all fragments are sequenced up to ℓ bases, for a properly chosen value of ℓ, and then the processor maps the fragments that are uniquely mapped to the reference genome. Un-mapped reads with non-overlapping reads from the right hand side at the first step are sent back to sequencer for extension to next bases. This procedure is repeated until the process reaches the maximum read length L. As shown in the paper, through use of such approach, the number of bases read in the sequencing part reduces to O(G) bases, a reduction by a log G factor in comparison with Lander-Waterman coverage bound.
We present theoretical results for i.i.d. and real genomes with noiseless and noisy reads. Also, we have simulated our method for chr19 of Human genome hg19 with different sequencing error rates. Simulation results support the validity of the proposed algorithm and demonstrate our improvement on coverage bound for real genomes. For future work, one may expand our algorithm to derive more efficient alignment algorithms in terms of complexity and precision. Also, this method may be extended to Denovo sequencing.

Author Contributions
Conceptualization: DN SAM.