article: Measuring the overlap between two populations is, in principle, straightforward. Upon fully sampling both populations, the number of shared objects—species, taxonomical units, or gene variants, depending on the context—can be directly counted. In practice, however, only a fraction of each population’s objects are likely to be sampled due to stochastic data collection or sequencing techniques. Although methods exists for quantifying population overlap under subsampled conditions, their bias is well documented and the uncertainty of their estimates cannot be quantified. Here we derive and validate a method to rigorously estimate the population overlap from incomplete samples when the total number of objects, species, or genes in each population is known, a special case of the more general β-diversity problem that is particularly relevant in the ecology and genomic epidemiology of malaria. By solving a Bayesian inference problem, this method takes into account the rates of subsampling and produces unbiased and Bayes-optimal estimates of overlap. In addition, it provides a natural framework for computing the uncertainty of its estimates, and can be used prospectively in study planning by quantifying the tradeoff between sampling effort and uncertainty.
  Of the diverse multigene families of P. falciparum, the var family is the most heavily studied because of its direct links to both malaria’s duration of infection and its virulence [14–17]. Each parasite genome contains a repertoire of ∼ 60 hypervariable and mutually distinct var genes, but repertoires differ between parasites, evolving rapidly through recombination and reassortment. Recent studies of P. falciparum epidemiology and evolution have generated insights by comparing of the sets of genomic var repertoires between parasites [18–23]. Indeed, since var repertoires are, themselves, under selection, theory suggests that if a human population has been exposed to particular var genes, then repertoires containing those var genes will have a lower fitness than repertoires that are entirely unrecognized by local hosts, shaping the var population structure [21,22,24,25]. Methods by which we estimate the extent to which var repertoires overlap are therefore important, particularly as studies of the population genetics and genetic epidemiology of malaria’s antigens become more sophisticated and data rich. However, as with estimates of β-diversity in ecology, traditional estimates of overlap between var repertoires also suffer bias due to subsampling. Due to the massive diversity and recombinant structure of var genes, researchers are restricted to using degenerate PCR primers targeting a small “tag” sequence within a particular var domain called DBLα [26]. Due to their experimental accessibility, DBLα tags have been widely used to study the structure and function of var genes [17,18,21,26–30]. Still, these PCR techniques generate a random sample of 60 or fewer unique tag sequences from each parasite. This means that experimental measurements of repertoire overlap are performed using stochastic subsamples whose empirical overlap may fluctuate from experiment to experiment (Fig 1), motivating the three questions answered by this paper: First, how can we estimate the true overlap between repertoires when we can only measure the overlap between samples from repertoires? Second, how can we quantify the uncertainty around our repertoire overlap estimates? Third, what are the implications of uncertainty for the design and budgeting of var repertoire studies? In the malaria literature, repertoire overlap is most commonly computed using the Sørenson-Dice coefficient where it is often called pairwise type sharing [18]. When PCR methods have produced na and nb tags from parasites a and b, respectively, and when a sequence-level comparison has found nab tags are shared by both repertoires then repertoire similarity is computed using the coefficient S ˚ in Eq (1). When na and nb are nearly 60, the performance of S ˚ is excellent. For instance, when two parasites are completely different, nab = 0, so S ˚ = 0; when two parasites are identical, and both repertoires have been fully sampled, nab = na = nb, so S ˚ = 1. However, when na or nb is smaller (as is overwhelmingly the case in existing studies [18–23]) S ˚ is conservative and systematically underestimates the true overlap between repertoires [7]. In this manuscript, I introduce a method that estimates repertoire overlap using Bayesian inference. By modeling the stochastic process by which repertoires are sampled, I show that this method produces unbiased a posteriori estimates of true repertoire overlap. I then show how the Bayesian framework can be used to estimate uncertainty and produce error bars which represent credible intervals, a Bayesian analog of confidence intervals. These methods are then used to reevaluate past results which used the Sørenson-Dice coefficient S ˚. Finally, in the case of P. falciparum, since each successful PCR amplification randomly samples just one of 60 available tags, I extend the Bayesian approach to compute the tradeoff between increasing sampling and decreasing the uncertainty of overlap estimates. These calculations allow the cost of sampling to be weighed against scientific confidence, illustrating the use of this statistical framework for planning and budgeting experiments. Open-source code and a web tool are freely available (see Acknowledgements).
  Suppose that there are two P. falciparum parasites, each with a repertoire of 60 var types. Our goal is to estimate the true repertoire overlap s (were we to fully sample each parasite) from the knowledge that na samples from parasite a and nb samples from parasite b share nab types. Due to the fact that the underlying sampling process is stochastic (Fig 1), our secondary goal is to quantify the uncertainty in the method’s estimates. Both goals can be met by writing down the process that creates the data in the first place. Therefore, in what follows, we will at first assume that the true overlap s is fixed, model the process of generating data via stochastic sampling, and use that model to compute a likelihood. We will then use Bayes’ Rule to compute the posterior probability for each value of s, given the evidence in the data and the likelihood computed in the first step. Consider the following sampling process, written in the slightly more rigid and generic language of a probability textbook. Suppose that there are s special objects among a total of N objects. We draw n objects uniformly at random without replacement. The number of special objects chosen during this sampling procedure will be distributed according to a hypergeometric distribution, which we write as H (s, N, n). First, with this definition in mind, consider drawing na var genes from parasite a’s 60 total. Of the 60 total, suppose that exactly s are considered special because they are also shared by parasite b. The number of shared sequences that are captured by sequencing parasite a will be a random variable S a = H (s, 60, n a). Depending on the luck of the draw, this number could be as small as zero, or as high as s or na (whichever is smaller). Now consider drawing nb var genes from parasite b’s 60 total, in which exactly sa are special because they are shared by both parasites and were actually drawn during the sequencing of parasite a. This process is identical in construction to the process for sampling parasite a, but with sa special sequences instead of s, and so the number of shared sequences that are captured after sequencing both parasites will be H (s a, 60, n b). Substituting the random variable Sa for a fixed value sa, which we derived in the paragraph above, yields a hypergeometric inside a hypergeometric, which means that the probability of a particular number of shared sequences in the samples nab is given by these sequential (or nested) hypergeometric distributions, P (n a b ∣ n a, n b, s) ∼ H (H (s, 60, n a), 60, n b). (2) Reassuringly, one can switch the order in which the imagined sampling took place, first sequencing parasite b and then sequencing parasite a, or sequencing them both at once, and show that these are mathematically equivalent. In practice, we want to go the other direction, and estimate s from our empirical measurements of na, nb, and nab. Since the distributions above allow us to compute the likelihood of empirical observations, given s, we use Bayes’ rule to formulate the posterior distribution for s, P (s ∣ n a, n b, n a b) = P (n a b ∣ n a, n b, s) P (s) P (n a b), (3) where P (s) is the prior distribution for overlap. In practice, we generally wish to remain agnostic about the level of overlap s and therefore we consider an uninformative prior P (s) ∼ unif[0,60], i. e. P (s) = 1 61. Using the law of total probability to rewrite the denominator, and canceling the factors of 1 61, we get P (s ∣ n a, n b, n a b) = P (n a b ∣ n a, n b, s) ∑ s ′ = 0 60 P (n a b ∣ n a, n b, s ′). (4) Each term on the right hand side of Eq (4) can now be computed directly from the nested hypergeometric distributions in Eq (2) as follows. To generate a specific empirical overlap nab, two things must have happened in succession and independently of each other: first, sa of the original s shared sequences must have been sampled; and second, nab of the intermediate sa shared sequences must then have been sampled. We therefore multiply these two hypergeometric probabilities. However, because this sequential process may occur for any value of the intermediate variable sa, we sum over all possible values of sa, P (n a b ∣ n a, n b, s) = ∑ s a = 0 60 P (n a b ∣ n b, s a) P (s a ∣ n a, s). (5) Thus, computing the posterior probability that the true overlap was s, given the empirical overlap between samples, is given by substituting Eq (5) into Eq (4), yielding P (s ∣ n a, n b, n a b) = ∑ s a = 0 60 P (n a b ∣ n b, s a) P (s a ∣ n a, s) ∑ s ′ = 0 60 ∑ s a = 0 60 P (n a b ∣ n b, s a) P (s a ∣ n a, s ′). (6) The term P (s ∣ na, nb, nab) is a posterior distribution over s, meaning that it tells us the probability for each value of s, given the evidence provided by the actual data. While this equation appears notation-heavy, its inference requires only calls to the hypergeometric probability distribution. To illustrate this graphically, the posterior distribution is plotted for na = 47, nb = 32, and nab = 20 in Fig 2. The posterior distribution can now be used (i) to estimate the true value of s, and (ii) to quantify the uncertainty of that estimate. First, our estimate for the true value of s, which we call s ^, is the expected value of the posterior, s ^ = ∑ s = 0 60 s P (s ∣ n a, n b, n a b). (7) This value is typically (in 99. 85% of all possible cases) larger than the estimate provided by S ˚ (Fig 2). The framework here is easily extended to repertoire sizes other than 60, generalizing to applications beyond var genes. Suppose that populations a and b have total sizes of Na and Nb, and without loss of generality, assume that Na ≤ Nb. The values Na and Nb need only be substituted into Eqs (6) and (7), with P (s) = N a+1- 1. This is done explicitly in Eqs. (1) and (2) in S1 Text, but not shown here for conciseness (see S1 Text). The posterior distribution provides a convenient way to quantify the uncertainty associated with an estimate s ^. Intuitively, if the posterior is sharply peaked around s ^, then our confidence in s ^ is high; if the posterior is broadly distributed then our confidence in s ^ is low. Making use of the Bayesian construction once more, we compute a credible interval by finding the range of s values that account for 90% of the posterior probability (Fig 2). Due to the fact that the posterior distribution is a discrete distribution over only 61 values, it is possible (indeed, highly probable) that no interval will contain exactly 90% of the probability. Nevertheless, we define a conservative equal-tailed 90% credible interval [smin, smax] as the smallest index smin and the largest index smax for which ∑ s = s max 60 P (s ∣ n a, n b, n a b) ≥ 0. 05 ∑ s = 0 s min P (s ∣ n a, n b, n a b) ≥ 0. 05. (8) We first demonstrate that the s ^ computed in Eq (7) produces accurate estimates by simulating the sampling process with known s and evaluating our ability to accurately recover it. Specifically, for each simulation, we consider two var repertoires a and b, of 60 genes each, and specify a priori that they share exactly s sequences. We then choose the number of samples taken from each, na and nb respectively, and draw from each repertoire uniformly at random, without replacement. These draws are compared to compute the number of empirically shared sequences nab. Eq (7) is used to compute the Bayesian repertoire overlap (BRO) estimate s ^, while Eq (1) is used to compute S ˚ using the same data. These estimates are then compared to the true value of s to evaluate accuracy. Varying the values of s, na, and nb allows us to quantify the performance of BRO and S ˚ in a variety of realistic sampling scenarios. Fig 3 shows the results of this simulation for sampling rates of 30,40, and 50 genes, with two independent simulations at each value of s. Intuitively, both BRO and S ˚ are more accurate when na and nb are larger. However, the two methods’ behaviors are fundamentally different. When na and nb are below 60, BRO provides estimates that are distributed around the true overlap, with variance decreasing as sampling rates increase. In contrast, S ˚ systematically underestimates the true overlap, while also showing decreasing variance as sampling rates increase [7]. For realistic sampling rates, BRO provides estimates centered at the true value, while S ˚ provides estimates centered below the true value. These general patterns hold even when na ≠ nb or when total repertoire sizes are unequal (S1 Fig). Credible intervals, which visually show uncertainty in each estimate, can also be easily computed from the simulations described above. For each simulation, Eq (8) uses the posterior distribution over s to produce error bars around the point estimate s ^, shown for sampling rates of 30,40, and 50 in Fig 4. This illustrates the substantial reduction in uncertainty that comes with increased sampling rates. While all simulations shown here use na = nb, this is by no means required (S1 Fig), and in real data scenarios, is rare. We now show how the methods of this paper can be used in practical contexts by applying them to data from three published studies. In particular, this reanalysis highlights the impact of variation in sampling rates across studies, which creates variable bias in S ˚ calculations and produces misleading results. However, we also show that while using BRO in place of S ˚ sidesteps bias problems, the ability to quantify uncertainty with error bars highlights new problems. In short, the conclusions of previous studies may be worth reevaluating. In 2007, Barry et al. introduced S ˚, which they referred to as pairwise type sharing, in an analysis of var data from Amele, Papua New Guinea [18]. In 2010, Albrecht et al. included Barry’s data in a broader analysis of var data from Ariquemes, Brazil [19] which also included sequences from a study by Bull et al. from Kilifi, Kenya for comparison [27]. Each one of these studies, individually, sequenced parasite isolates to a particular target depth, yet the studies varied in their coverage of repertoires. Since the bias of S ˚ depends on the number of samples (Fig 3; see also [7]), the variation of sampling rates across study populations means that different populations are biased downward by different amounts. Albrecht et al. conveniently provide var type data from all three studies, from which we can rebuild their first figure which shows a S ˚ comparison of five populations (Fig 5; left). Overlaps between pairs of parasites can then be recomputed using BRO (Fig 5; middle). The conclusions drawn from these two figures differ substantially. First, according to S ˚, identical clones from Ariquemes share only around 30 sequences with themselves, illustrating the downward bias produced by subsampling—clones ought to share all of their genes with their genetically identical siblings. Indeed, the reanalysis using BRO finds over 75% of overlap estimates to be greater than 50 (and over 50% over 55), far closer to what is expected. Second, the inter-clone overlap and inter-isolate overlap distributions in Ariquemes appear to be similar and overlapping through the lens of S ˚. However, the recalculation using BRO shifts the clones’ distribution dramatically upward but leaves the isolates’ distribution more or less untouched. This is due to the dramatic difference in var coverage: the average number of sequences per clone is n ¯ = 26. 5 while for isolates it is n ¯ = 45. 8, meaning that relatively different amounts of bias are inherited from S ˚ (illustrated in simulations in Fig 3). Finally, the distributions from Brazil (n ¯ = 17. 3) and Amele (n ¯ = 15. 6) also shift dramatically upward when the bias of S ˚ is removed (Fig 5; left, middle). However, this does not necessarily mean that they should be reinterpreted. For each pairwise comparison, Eq (8) allows us to compute the width of the credible interval, smax − smin + 1, quantifying our uncertainty in each estimate. Due to low average coverage, the uncertainty of estimates in the Amele dataset tends to be extremely large (Fig 5; right), with the majority of estimates showing an uncertainty greater than 30 sequences (50% overlap). For comparison, estimates from Ariquemes clones (n ¯ = 26. 5) are also shown, whose dramatically lower uncertainty enables more confident conclusions to be drawn. There are two main methodological findings that result from using rigorous and unbiased methods. First, the boxplots of Fig 5 clearly illustrate that sampling rates can have a dramatic impact on findings, reinforcing the simulation results of Fig 3. Second, uncertainty is an issue when n ¯ is too small, and datasets with low sampling rates may have such wide error bars that their estimates should not be trusted, as shown in the histograms of Fig 5, reinforcing the simulation results of Fig 4. Additional sequencing efforts come at a cost, however, and so in the next subsection we use the methods of this paper to quantify the tradeoff between increased sequencing efforts and decreased uncertainty. In the previous section, the reanalysis of published results shows clearly that the number of samples per parasite has a dramatic impact on the uncertainty (and therefore the interpretability) of painstakingly collected parasite sequence data. Naturally, increasing the sampling rates, na and nb, decreases the uncertainty in s ^, our estimate of s (Fig 4). However, additional samples cost time, effort, and money. Complicating matters, generating additional var sequences may or may not increase na, since the previously sequenced var tags may be redundantly sequenced. Thus, there is a stochastic tradeoff between increased laboratory effort and decreased uncertainty about repertoire overlap, which we now calculate. To obtain var tags, the DNA is PCR amplified using degenerate primers that are designed to universally capture all var genes with DBLα domains. This product is then cloned into a vector that allows single products to integrate, and these vectors are then transformed into bacteria and plated such that each colony contains one vector and one insert (see e. g. [21] for detailed methods, but see also [25] which uses a different pipeline based on next-generation sequencing). Therefore, among a large number of colonies, there are likely to be multiple colonies with the same var gene while some genes may not be covered by any colony. How many colonies should be separated and sequenced in order to get an accurate estimate of the repertoire overlap between two parasites? Put more formally, if we repeatedly perform an experiment in which we sequence c colonies each from two parasites and estimate their overlap s ^, how much more accurate will s ^ become if we increase c? To answer this question, we split it into two parts. First, if we sequence c colonies, how many unique var genes n are we likely to have sampled? Second, what implications will this have for our repertoire overlap estimates, discussed in the previous section? The first question can be answered by considering a process in which there are k = 60 distinct sequences in total and we draw c of them, one at a time, independently and with replacement. For a fixed c, we can compute the probability mass function for the number of distinct sequences by a straightforward recursion: At any point during the process of drawing sequences, if n distinct sequences have already been drawn, then the probability of drawing an already-discovered sequence is n/k, making the probability of drawing a new sequence 1 − n/k. Each draw is independent of the previous draws, so the incremental accumulation of distinct sequences can be written as a Markov chain with transition matrix π whose non-zero entries are π n → n = n k and π n → n + 1 = 1 - n k. (9) Initially, zero sequences have been drawn (c = 0), making n = 0 with probability 1. For each additional sequence drawn, the probability distribution over the number of distinct sequences evolves according to the transition matrix π, so that after c draws the distribution over distinct sequences is given by the entries of the vector x, x = x 0 T π c, (10) where x0 is initial condition vector of zeros, except for the entry corresponding to the state n = 0, which equals one. This allows us to analytically compute the distribution of the number of unique var genes sampled by a PCR process with c colonies. In other words, we now have a map between laboratory efforts c and the distribution of actual unique var genes sampled, and we write this as P (n ∣ c). A variant of this problem was previously considered with the goal of computing the value of c that would cover at least 60% of each repertoire [31]. Although those calculations can be shown to produce incorrect estimates, Eq (10) can be used to solve that problem variant as well. More widely, this general problem has been charmingly named the coupon collector’s problem by statisticians. The second question focuses on the implications of Eq (10), and specifically requires that we quantify how an increase in sequencing efforts c affects the noisy distribution of estimates s ^. Intuitively, for low c, both na and nb will tend to be small, leading to broad distributions of s ^ around the correct value of s. Similarly, as c grows very large, we expect the distribution of s ^ to concentrate on exactly s. This distribution, P (s ^ ∣ s, c), can be computed by integrating the distribution of estimates, conditioned on particular data, over the probability distribution of having produced those data, conditioned on c and s. Symbolically, the distribution of estimators s ^, given true overlap s and colonies c is given by P (s ^ ∣ s, c) =∑ n a, n b, n a b { P (s ^ ∣ n a, n b, n a b) ×P (n a b ∣ n a, n b, s) P (n b ∣ c) P (n a ∣ c) }. (11) P (s ^ ∣ n a, n b, n a b) is the probability of getting a particular estimate s ^, given information about coverage and overlap. In fact, this is a distribution concentrated at a single point, i. e. , a Dirac δ function, since each triple (na, nb, nab) maps to exactly one point estimate s ^. As a result, this term tells us the locations at which there will be probability mass, while the remaining terms in Eq (11) tell us how much mass there will be at those locations. In other words, this distribution is a discrete probability distribution, and we have written down a fancy form of it above. By aggregating into bins, this distribution can be conveniently visualized as a histogram, which shows how the uncertainty of estimators depends on the true overlap s and the sequencing effort c. Fig 6 shows the effect of increasing sequencing efforts from a half plate (c = 48) to a full 96-well plate (c = 96) and beyond. These calculations succinctly quantify intuition: additional laboratory efforts lead to higher accuracy guarantees. The calculations and distributions in this section show how the Bayesian framework in this manuscript can also be used to plan sequencing studies and estimate study costs. If a desired downstream analysis of repertoire overlap requires results that are accurate to within a particular number of shared sequences, BRO methods can easily specify the sequencing efforts needed.
  This manuscript places the estimation of overlap between fixed-size repertoires or populations from incomplete samples on firm statistical ground. While myriad indices of β-diversity for presence-absence data exist [2], they implicitly treat species counts as complete, leading to bias. Notable exceptions which embrace imperfect sampling exist [7], but require species abundance data to compute. Here, we clearly define a stochastic process for fixed-size repertoires that generates sample presence/absence data, opening the door to rigorous Bayesian inference. In particular, Eq (7) provides point estimates of true repertoire overlap, while Eq (8) provides error bars and uncertainty estimates via credible intervals. If desired, improved estimates of s ^ can be plugged directly into any of the dozens of presence/absence measures of similarity reviewed in Ref. [2]. Figs 3 and 4 show the consistency and accuracy of these calculations across simulated sampling regimes in which the correct answer is known. Bayesian repertoire overlap (BRO) is also useful in real-data scenarios, when the correct answer is unknown. By revisiting previously published studies of the var genes of P. falciparum [18,19,21,27], we showed that switching from the Sørenson-Dice coefficient S ˚ (called pairwise type sharing in the malaria literature) to BRO leads to different conclusions (Fig 5 left, middle) or high uncertainty (Fig 5 right). In particular, these reanalyses point to a clear recommendation for the design of future malaria studies: the number of unique var sequences per isolate should be at least 30. Since each additional PCR product may not contribute an additional unique sequence, we again used the Bayesian framework to translate increased PCR efforts to decreased uncertainty (Fig 6). Accuracy requirements can now be weighed against laboratory costs during the planning of studies. While BRO clearly outperforms S ˚ in practical contexts, it is also more cumbersome to compute. Indeed, S ˚ can be calculated on the back of an envelope while Eq (7) requires a computer, or at least a lot more envelopes. However, as it turns out, there are only around 77,500 possible combinations of na, nb, and nab, which means that a lookup table of every conceivable s ^ value can be computed on a laptop in minutes and attached to an email. Links to open-source code and a convenient web tool can be found in the Acknowledgements. The models introduced in this paper are as correct as their assumptions, which we now revisit. During the construction of the Bayesian repertoire overlap, we assumed that our prior distribution P (s) was uniform, meaning that we treated each possible level of overlap as equally likely. This is easily defensible in practice, as any other choice would introduce unacceptable bias. We also assumed, when computing the tradeoff between sequencing effort and uncertainty, that each sequence in each repertoire was just as likely to have been sampled, which may or may not be true, for two distinct reasons. First, due to the fact that sequences are obtained using degenerate primers, the effects of primer bias may cause some sequences to be amplified more often than others. Second, a single parasite genome might have multiple copies of the exact same var gene, or might have distinct var genes whose DBLα tags are nevertheless identical. This scenario is arguably more likely among South American genomes whose overall var diversity is lower. Experimentally, the probability that a sequence is observed will be scaled upward by its genomic multiplicity, but the scaling may be non-linear since PCR protocols include many rounds of amplification, magnifying the deviations from uniformity. Fully addressing either of these possibilities would require that we modify the probabilities in both the coupon collector’s problem and the repertoire subsampling processes, and then use Monte Carlo methods to numerically compute posterior distributions. New sampling protocols for var genes, based on next-generation sequencing methods [25], may or may or may not meet the assumptions of the estimator presented here. To use the hypergeometric distribution, we require that, if an entire sampling protocol were to be technically replicated many times, that eventually each member of the repertoire would be observed with equal probability. In other words, while fluctuations in any particular set of observations are expected, with technical replication those fluctuations must eventually even out, approaching uniformity. Thus, the issues of primer bias and gene multiplicity violate the assumption of uniformity, but the magnification of random initial fluctuations, e. g. , by the repeated amplification rounds of PCR, do not. Could deviations from the modeling assumption of uniformity could be inferred from the data themselves? If so, this idea could in principle be applied to cloning-based methods and next-generation methods alike. This is an interesting direction for future work, and could draw from advances in abundance-based estimators for β-diversity [7,12], or could incorporate explicit knowledge of the effects of protocols and pipelines, in order to mathematically undo their effects. More practically, the assumption that the var repertoire size is 60 makes the methods of this paper less useful in the context of complex infections with multiple parasite genomes [23,27]. In cases where the multiplicity of infection is known, overlap estimates could be computed using generalizations of the statistics in this paper, computing overlap between infections (instead of between parasites). This would be complicated by possible overlap of parasite repertoires within each infection, but repertoires tend to be quite different in areas of high transmission so the methods herein may be approximately correct. Nevertheless, development of more sophisticated methods would be especially useful in the context of var-based epidemiological studies. This paper focuses on malaria’s var genes, and assumes a total repertoire size of 60, but mathematically relaxing this assumption (S1 Text) broadens the set of possible applications. First, within studies of malaria’s var genes, the total repertoire size fluctuates slightly from parasite to parasite. As larger whole-genome datasets become available, this information can be incorporated directly as a prior over the distribution of repertoire sizes, improving estimates further. This opens the door to the analysis of Plasmodium spp. multigene antigen families such as rif and stevor [32,33], or more general studies of β-diversity in multigene families in which population sizes are fixed or their size distributions have been sampled [34]. Second, outside of malaria, improved estimators may also be useful in comparing, for instance, the genomic archives of antigen-encoding vsg genes used by Trypanosoma brucei for immune evasion [35]. Finally, the mathematics of this paper need not be applied to genetics or even within ecology; large-data applications like Facebook and other online social networks compute the sizes of intersections of sets—how many interests or friends do two individuals have in common? —but use only subsamples of data to decrease computation time in distributed computing settings [13]. Finally, this work presents a Bayesian approach to inference of β-diversity under particular assumptions, which contrasts the vast majority of indices, coefficients, and metrics to date which remain non-probabilistic [2]. Changing the underlying assumptions of the Bayesian repertoire overlap method, or the statistics of the sampling process, would lead to additional estimators for other common cases. Across applications, unbiased estimation combined with the quantification of uncertainty will allow for more reliable results and better prospective study design.
summary: Understanding when two populations are composed of similar species is important for ecologists, epidemiologists, and population geneticists, and in principle it is easy: just sample the two populations, compare the sets of species identified in each, and count how many appear in both populations. In practice, however, this is difficult because sampling methods typically produce only a random subset of the total population, leaving current population overlap estimates biased. Knowing only the number of shared members between two of these partial population samples, this paper shows how we can nevertheless estimate the true overlap between the full populations, when those full populations’ sizes are known. Using Bayesian statistics, we can also compute credible intervals to produce error bars. We show that using this unbiased approach has a dramatic impact on the conclusions one might draw from previously published studies in the malaria literature, which used simple but biased methods. Because the method in this paper quantifies the tradeoff between sampling effort and uncertainty, we also show how to compute the number of samples required to ensure high-confidence results, which may be useful for planning future studies or budgeting lab reagents and time.
keywords: medicine and health sciences
  parasite groups
  ecology and environmental sciences
  plasmodium
  tropical diseases
  parasitic diseases
  parasitic protozoans
  parasitology
  apicomplexa
  probability distribution
  mathematics
  protozoans
  molecular biology techniques
  population biology
  research and analysis methods
  disease ecology
  malarial parasites
  artificial gene amplification and extension
  epidemiology
  population ecology
  molecular biology
  probability theory
  eukaryota
  ecology
  polymerase chain reaction
  biology and life sciences
  physical sciences
  malaria
  organisms
title: Bayes-optimal estimation of overlap between populations of fixed size