We have given bounds on the probability distribution of the number of uncovered gaps when arcs of fixed length are placed randomly on a circle or interval. As discussed in [6], one motivation for such approximations is the issue of computational overflow arising when the exact solution is applied. Typically they involve simple, well-known probability distributions and this aids both computation and further mathematical analysis. Our own motivation in beginning this work was seeing certain quite poor approximations used in practice, both locally and in the literature. For a cautionary example we take s = 10-2 and n = 750, values which arise when 0.3% of the clones in a metagenomic library consisting of 2.5 × 105 40-kilobase fosmid inserts are from the genome of interest, and the genome of interest has length 4 megabases. Here ad-hoc approximations using elementary probability theory can indicate a 95% probability that the library completely covers that genome, while the true probability is 66% (neglecting biologically related bias; calculations are to the nearest integer and are not given).
We stress the similarity between our approximations and those already given in the literature. In [2, 3, 7] simplifying assumptions are made which give a binomial distribution for the number of gaps. Our approximation is a Poisson distribution, and it is well known in probability theory that the binomial distribution converges to the Poisson distribution in certain limits; one such convergence is proved in [8]. Indeed, these different approximations are generally numerically close. The contribution of the present paper is therefore to provide an approximation with error bounds.
We refer the reader to [2] for a simple modification to coverage problems when a certain minimum overlap is required between arcs, and give a worked example in Additional file 2. Note that in [2] the expected number of gaps is calculated approximately; as shown in Additional file 1, the exact value for the expected number of gaps is m = n(1 - s)n-1. This is the parameter of the Poisson distribution in Proposition 1. Our approximating Poisson distribution therefore has the same expectation as the exact distribution, although its variance differs (the exact variance of the number of gaps is given in Additional file 1).
Our results may be applied to the stopping problem. Suppose we desire p% probability that no more than w gaps remain at the end of the undirected part of a shotgun sequencing project. By inverting Proposition 1 (see Additional file 2 for details) we obtain an approximate value for the number of clones, and hence coverage depth, required. Here again, our contribution is not the estimate but rather the lower bound given by Proposition 2 for the probability that no more than w gaps remain (neglecting biologically related bias). Other solutions to the stopping problem have been proposed: see for example [6] and [3] for alternative strategies and further discussion of the stopping problem.
The practical relevance of our approximations clearly depends on the size of the error bound. Figure 1 plots the error bound against coverage depth for arc lengths s = 10-1, 10-2, 10-3 and 10-7 (the curves for s = 10-4, 10-5 and 10-6 are almost indistinguishable from the s = 10-7 curve at this scale). We emphasise that s is the relative arc length – so for genomic applications, s is the actual arc length divided by the length of the genome. These arc lengths are intended to represent the full range of typical genomic projects: for example the smallest, s = 10-7, would correspond to covering the largest known eukaryote genome, the amoeba Chaos chaos [9], with 400 kilobase bacterial artificial chromosome (BAC) inserts. A lookup table of relative arc lengths for recent shotgun sequencing projects is given in [4]. It can be seen that for these relative arc lengths, error bounds on the order of 0.1% are achieved at 10× coverage. Further discussion of Figure 1 is given in Additional file 1; it should also be remarked that the error bound at 5× coverage is considerably larger. For the particular experimental parameters relevant to the user, the spreadsheet in Additional file 3 may be used to obtain values for the approximation and error bounds.
The validity of coverage problems in general for genomic applications depends of course on the extent to which they capture the actual problem. For example, a given genome might contain one gene which is toxic to the E. coli host in a BAC library. Since library fragments will not contain this gene, the corresponding coverage problem is on the interval rather than the circle. With a priori knowledge of an unclonable region we may therefore apply Proposition 3 rather than Proposition 1; without such information we may choose either to neglect this effect or to model it (for example using Propositions 1 and 3 and conditional probability). Pathological cases certainly exist, for example the highly repetitive maize genome for which as many as 80% of arc positions may not be cloned [7]. The approach should therefore be chosen using the best available information. Certain other biases are not so dependent on the particular target genome, for example the inevitable small variations in clone length and position bias. These have been discussed in [3, 4, 10] using empirical data and simulations and the consensus is that they may be neglected. The interested mathematical reader may check that by first conditioning on the fragment lengths, our method gives a lower bound for P(W ≤ w) in Proposition 2 when the arc lengths are random, although we do not pursue this.
Another modelling issue arises in metagenomics, which is the culture-independent study of a mixed community of genomes. In a metagenomic library the number of clones n from a genome in the community is random, having a binomial distribution (in the absence of bias). If the composition of the community is known then, from the central limit theorem of probability, n is well approximated by its average and this value may be used in Proposition 1 (see Additional file 2). Further, since the distribution of n is concentrated at a few values around its mean, it is typically computationally inexpensive to obtain satisfactory bounds for the distribution of the number of gaps by an application of conditional probability.
