BACKGROUND
rare variants are responsible for a large portion of the heritability of some common complex human diseases  <cit> .
genome-wide association studies have focused on the contribution of variants of low minor allele frequency , or of rare variants   <cit> .
the functional and evolutionary impacts of rare variants have been reported; therefore, large-scale screening for disease-associated rare variants is becoming increasingly important  <cit> .
one major application of rare variants genotyping is in screens for rare genetic disorders such as tay–sachs disease and thalassemia  <cit> .
because of the extremely low frequency of rare variants, sample sizes must be large enough to guarantee efficient observations.
the cost of dna sequencing has dropped dramatically with the introduction of the massively parallel sequencing technologies.
however, identifying rare variant carriers by sequencing individual samples separately remains prohibitively expensive because of the huge challenge of constructing sequencing libraries for thousands of samples  <cit> .
barcoding has been developed as a powerful approach to cost-effectively determine the genotype of each individual  <cit> .
to further reduce the cost of large-scale screens for rare variant carriers, several techniques based on the group testing theory  <cit>  and compressed sensing  <cit>  to construct overlapping pool sequencing strategies have been used.
these strategies have helped decrease the sequencing times for rare variant carrier identification and further lower the cost  <cit> .
because a large number of samples are pooled together and then sequenced, overlapping pool sequencing uses fewer pools to identify rare variant carriers among large numbers of individuals.
thus, overlapping pool sequencing can vastly reduce the time and cost for dna library preparation.
some overlapping pool sequencing programs return true/false values after testing each pool; this scheme was adopted by erlich et al.
<cit> , prabhu and pe’er  <cit> , and cao et al.
<cit> , who used the number of normal and variant reads in each pool to determine whether a pool contained carriers.
however, the quantitative information in the sequencing data, which can be used to estimate the percentage of variant chromosomes in a pool, is missed in these methods.
quantitative group testing is an alternative scheme that takes the quantitative information into account, thus rare variant carriers can be identified efficiently  <cit> .
we propose an efficient random overlapping pool sequencing strategy with quantitative group testing for the identification of rare variant carriers using massively parallel sequencing data.
because of the excellent performance of random designs in classic group testing  <cit> , we employed a random k-set pool design  <cit>  to mix samples.
the parameters of the random k-set pool design can be selected appropriately according to an indicative probability value.
based on a depth model for pooled sequencing, we calculated the optimal cut-off of the number of reads containing variants to distinguish pools containing variant carriers from those that do not.
using the quantitative information contained in the sequencing data, we designed a heuristic bayesian decoding algorithm to identify variant carriers accurately.
compared with the dna sudoku algorithm  <cit>  and compressed sequencing  <cit> , our method required less data throughput.
finally, we applied our method to identify variant carriers among  <dig> simulated escherichia coli strains using simulated pools and publicly available illumina sequencing data.
the results showed that our method could successfully identify carriers for  <dig> – <dig> % of the variants with frequencies ranging from  <dig>  to  <dig> %.
methods
random k-set pool design
random k-set pool designs were first proposed by bruno et al.
<cit>  for efficient dna clone library screening.
in such a design, each clone is pooled in exact k pools that are chosen with equal chance.
random k-set pool designs are easy to specify for any number of pools and are known to be efficient in classic group testing, but they have not been used in real situations, partly because of the presence of different sample sets with identical test sets that are indistinguishable when the test results are qualitative  <cit> .
however, this problem can be overcome by quantitative tests such as sequencing experiments.
for n samples containing d positive samples, the basic objective of a random k-set pool design is to identify all the positive samples with a small number of pooled tests.
in such a design, each sample occurs in exact k pools, and a pool is defined as positive only when it contains at least one positive sample; otherwise, it is defined as negative.
for a random k-set pool design with t pools, hwang  <cit>  calculated the probability that a given set of i pools is a negative one ) and the expected number of negative pools .
negi=∑h=it−1h−it−ih−it−hkd/tkd×i∈n_min,n_max
e=∑i=n_minn_maxitinegi
where n_min and n_max are the minimum and maximum number of negative pools, respectively, and h is a temporary variable.
n_max = t - k, and n_min =  <dig>  or t - dk .
to evaluate the performance of random designs, barillot et al.
<cit>  proposed that the number of unresolved negative samples  can be taken as a criterion.
an unresolved negative sample is defined as a negative sample that occurred only in positive pools, as a result, its status can only be confirmed by testing it separately.
negative samples that are contained in at least one negative pool can confidently be determined as negative; therefore, hwang  <cit>  calculated the expectation ) and probability distributions ) for the number of unresolved negative samples.
en¯=n−d∑i=n_minn_maxtinegit−ik−d/tk−d
pn¯=j=∑i=n_minn_maxtinegit−ik−dj×) to evaluate the performance of random k-set designs in quantitative group testing.
pi is the probability that positive pools are more than the sum of unresolved negative samples and real positive samples; therefore, designs with high pis will perform better than designs with low pis.
pi=∑i=p_minp_maxtt−inegt−i∑j=0i−dpn¯=j
where p_min and p_max are the minimum and maximum number of positive pools, respectively, p_min = t - n_max, and p_max = t - n_min.
the derivation of eq.
is given in appendix  <dig>
optimal cut-off value for pooled sequencing
for overlapping pool sequencing, the dna samples are mixed and then sequenced.
samples from variant carriers are treated as positive and samples from normal individuals are treated as negative.
to distinguish positive pools containing variant carriers from negative pools consisting of normal individuals, the cut-off for the number of reads containing variants must be selected carefully to declare whether a pool contains carriers or not.
ideally, the cut-off value must guarantee that the minimum error rates are obtained, including false-positive and false-negative rates.
the variation of sequencing depth distribution is significantly greater than the mean  <cit> ; therefore, in recent studies, negative binomial distribution rather than poisson distribution has been used to model sequencing depth.
following the study reported by miller et al.
<cit> , we used a negative binomial model to estimate the sequencing depth distribution.
given the average sequencing depth d, the depth that represents the number of times a nucleotide is read follows a negative binomial distribution nbdr− <dig> r where r is the variance/mean ratio; r is related to sequencing platforms and genomes and can be estimated from sequencing results ).
pdepth=i=nbi;dr− <dig> r
for a pool with n samples, given sequencing depth d and average sequencing error rate perror, the probabilities pnv and ppv that nv reads containing variants are observed in negative pools and positive pools, respectively, are given by eqs.
and .
for a locus sequenced i times, the number of sequencing errors follows a binomial distribution bin.
pnvnv=∑i=nv∞nbi;dr− <dig> rbinnv;i,perror
ppvnv=∑x=0nv∑i=x∞nbbinx;i,perror×∑j=nv−x∞nbbinj−nv+x;j,perror
similarly, the probabilities pnn and ppn that nn reads without variants are observed in negative pools and positive pools, respectively, are given by eqs.
and .
pnnnn=∑i=nn∞nbi;dr− <dig> rbini−nn;i,perror
ppnnn=∑x=0nn∑i=x∞nbbini−x;i,perror×∑j=nn−x∞nbbinnn−x;j,perror
where p is the percentage of variant chromosomes in the pool.
in a positive pool with n diploid dna samples and c heterozygous variant carriers, ignoring the bias in mixing samples, the percentage of the variant chromosomes is p = c/, while for haploid samples p = c/n.
the derivations of eqs.
– are given in appendix  <dig>
given pnv and ppv, the formula for the false-positive rate fp and false-negative rate fn in classifying pools with a cut-off value t can be constructed  and ).
fp=∑i=t∞pnvi
fn=∑i=0t−1ppvi
the optimal cut-off t can be defined as the value that minimizes the maximum values of fn and fp.
t=argminmaxfn,fp|t∈ <dig> d
decoding algorithm
our decoding procedure involves two steps.
in the first step, we identify negative pools according to the sequencing results and cut-off values for each pool.
samples that participate in any negative pools are classified as from normal individuals.
then, separate the real positive samples from unresolved negative samples according to the quantitative information in the sequencing results.
the probability of observing the sequencing results under the exact set of variant carriers should be greater than taking other set of samples as variant carriers.
assuming a is the set of variant carriers, the probability that the sequencing result o is observed is given by eq.
.
po|a=∏i=1tpoiv,oin|a
where oiv and oin are the number of reads with and without variants in the ith pool.
given that a is the set of variant carriers, if the ith pool is negative, then p = pnvpnn; otherwise, p = ppvppn.
for the second step, we designed a bayesian decoding algorithm based on eq.
.
to identify variant carriers among haploid samples, after excluding resolved negative samples, the rest of the samples form a set a0 = {s <dig> …,sc}.
first, assuming that all the samples in a <dig> are variant carriers, we calculate the probability of observing the sequencing results and denote it as pmax_ <dig>  next, replace one positive sample in set a <dig> with a negative sample in turn and repeat the probability calculation.
denote the derived set that results in the maximum probability as a <dig> and the corresponding probability as pmax_ <dig>  replace a <dig> with a <dig> and repeat this step until ac and pmax_c are obtained.
finally, the set ai that results in the maximum corresponding probability pmax_i is defined as the set of variant carriers.
these steps are written as algorithm  <dig>
two kinds of positive samples need to be considered while identifying variant carriers among diploid samples: heterozygous carriers and homozygous carriers.
first, suppose that there are only heterozygous variant carriers; this is analogous to finding variant carriers among haploid samples.
then we present algorithm  <dig> which is very similar to algorithm  <dig> to recognize heterozygous and homozygous variant carriers.
our decoding procedure to identify variant carriers among diploid samples is summarized in algorithm  <dig>  first, we suppose that only heterozygous variant carriers exist and run algorithm  <dig> to find the set of variant carriers c <dig>  then, run algorithm  <dig> to recognize heterozygous and homozygous variant carriers..
RESULTS
optimal cut-off value
to approximate the sequencing depth distribution for data obtained by using illumina sequencing platforms, miller et al.
<cit>  found that the negative binomial distribution with the variance/mean ratio r of  <dig> performed much better than the poisson distribution.
therefore, r was set to  <dig> in our simulation unless otherwise stated.for a pool consisting of  <dig> diploid samples, we calculated the false-positive and false-negative probabilities with different cut-off values when only one heterozygous variant carrier was allowed in the positive pool.
the average sequencing error rate and whole depth were set to  <dig>  and 600×, respectively .
the results verified the importance of selecting an appropriate cut-off value.
with smaller cut-off values, the probability of misclassifying negative pools as positive is high.
with higher cut-off values, some positive pools will be misclassified.
both will lower the speed and accuracy of decoding.
from the results we can infer that the optimal cut-off value is 14; here both the false-negative and false-positive probabilities are very low .in most studies, because of the rarity of positive samples, the number of positive and negative pools is unequal.
therefore, selecting a cut-off value based on the expected number of errors in classifying pools is a more practical approach.
for instance, when there are  <dig> negative and  <dig> positive pools mentioned above, the optimal cut-off value is  <dig> .
in the following simulation experiments, we adopted this kind of scheme unless otherwise stated.
as mentioned, the variance/mean ratio r is related to the sequencing platforms and genomes.
because the observed variation is significantly greater than the mean of the sequencing depth, r is larger than  <dig>  different values for r yield distinct distributions.
the pooling design mentioned above consisting of  <dig> negative and  <dig> positive pools was used to estimate the effect of r on our methods.
we calculated the least depth required for each pool to make the expected number of errors in classifying pools smaller than  <dig> by increasing the depth gradually .
from the results, we can see that our method performed even better for smaller r, which required less data throughput.
performance of the pipeline
to evaluate the performance of our method, we conducted substantial simulations to identify four heterozygous variant carriers among  <dig> haploid samples.
<dig> replicates were done for each pair of design parameters: pool number t and column weight k. the pooling matrix was designed by collecting random binary vectors with length t and weight k, meaning that each sample was mixed in k of t pools.
identical vectors were deleted and the steps were repeated until  <dig> vectors were obtained to form the matrix.
we used the random function in perl to simulate the number of reads with and without variants in pooled sequencing.
sequencing error and mixing bias were added to the simulation procedure to bring it closer to a real situation.
sequencing error follows a binomial distribution in sequencing results, and in the simulation the average sequencing error rate was set as  <dig> .
mixing bias is caused by the practical difficulty of mixing exactly equal amounts of dna samples.
based on the study conducted by shental et al.
<cit> , a random variable following the gaussian distribution was added to each non-zero element of the pooling matrix to simulate mixing bias.
the standard deviation of the gaussian distribution was  <dig> , reflecting up to 5% average noise in the mixed quantities of each sample.
after simulating the pooled sequencing results containing the sequencing errors and mixing bias, the genotypes of the  <dig> samples were reconstructed using our decoding algorithm.
the correct decoding rates, namely the percentages of simulations that identified all the variant carriers correctly, were determined .
the results showed that either a too large or too small k negatively influenced the correct decoding rates.
moreover, a large k meant more pooling procedures and increased experimental costs.
therefore, a proper column weight k is critical for conducting experiments successfully and efficiently.
cost-effective overlapping pool sequencing
the column weight k denotes the mixing times for each sample in a random k-set pool design.
for a given number of pools, a k that is too large or too small will lower the decoding accuracy.
we designed an indicative probability pi, which reflects the performance of random k-set designs that could be used to choose the optimal column weight k.
we calculated the correct decoding rates for different k under the condition that  <dig> pools were allowed to identify four heterozygous variant carriers among  <dig> diploid samples by conducting  <dig> replicates for each k. next, pi was computed based on eq.
and the results are shown in figure  <dig>  a strong correlation was observed between pi and the correct decoding rate , especially before the correct decoding rate reached the saturation point .
the pi values and correct decoding rates for identifying variant carriers were also obtained under different scenarios .
all the scenarios showed strong correlations between pi and the correct decoding rate before the correct decoding rate reached the saturation point.
for a given pool number t, we defined the optimal k as the minimum that obtains the maximum pi value, which could maximize the correct decoding rate.
designs with optimal k require fewer pools or lower sequencing depth.
in practice, the optimal k is selected by calculating the pi value without the need to conduct simulations, thereby greatly reducing the computational time required.
next, we conducted a series of simulated overlapping pool sequencing experiments with 20– <dig> pools and  <dig> – <dig> × overall sequencing data throughput .
one thousand replicates were conducted for each scenario, and the column weight was set as the optimal value .
the correct decoding rates were low when few pools or data throughput were used.
however, adequate pools and data throughput achieved higher accuracy but increased the cost, which conflicted with our motivation in this study.
there is a trade-off between the number of pools and data throughput.
hence, numerous simulations need to be performed to verify whether a pool number and data throughput pair can succeed in achieving high accuracy .
clearly, the optimal design parameters should be selected based on the whole cost of the sequencing experiment.
for a given population with  <dig> diploid individuals containing one heterozygous variant carrier, we generated several candidate designs in which over 95% of the simulations correctly identified the variant carrier .
the sequencing region was set to 30 mb to fit the human exome sequencing project  <cit> .
the cost of a sequencing experiment includes library construction and data production.
using the cost model from our previous work  <cit> , we inferred that the lowest cost design was design ii in table  <dig>  compared with sequencing separately, which requires sequencing depths of  <dig> × for each sample to obtain correct decoding rates over 95%, our method can save at least 50% of the cost.
with the same procedure, we generated the most cost-effective designs for variants with different frequencies and different sequencing region sizes .
for smaller sequencing regions and variants with lower frequencies, there are greater cost reductions with our method compared with those for larger regions and variants with higher frequencies.
a
aaverage value from five simulations.
gb is short for gigabases.
data throughput is the sequencing depth multiplied by the length of the sequencing region.
bsequencing separately is the strategy when each sample is sequenced independently.
all the candidate designs can identify the variant carrier correctly in 95% of the simulations.
the costs were estimated at $ <dig> for one library preparation and $ <dig> for  <dig> gb of data.
the sequencing region for haploid samples was set as 5 mb to fit the average length of the bacterial genome.
the sequencing region for diploid samples was set as 30 mb to fit the human exome sequencing.
comparisons with current methods
in  <dig>  benefiting from the chinese remainder theorem, erlich et al.
<cit>  put forward the dna sudoku design for overlapping pool sequencing.
a pattern consistency decoding algorithm was also developed by erlich et al.
<cit>  to identify variant carriers with the dna sudoku design.
in  <dig>  shental et al.
<cit>  developed a method called compressed sequencing to identify rare variants and their carriers by borrowing techniques from compressed sensing.
two designs were proposed in compressed sequencing: one used pools with a random half of the samples and the other used pools with sizes equal to the square root of the number of samples.
we compared the performance of our method in identifying rare variant carriers with the performances of these two methods.
to identify variant carriers in  <dig> diploid samples, the dna sudoku design with parameter d0 =  <dig> was employed that required  <dig> pools.
to maintain consistency, only  <dig> pools were allowed for the random k-set pool design and compressed sequencing.
since the expected number of positive and negative pools was not clear for the dna sudoku design, the cut-off value for the number of reads containing variants to declare a pool to be positive was set based on the false-negative and false-positive rates, and not on the expected number of errors in the classifying pools.
with  <dig> pools, we computed the least sequencing data throughput required for all the methods by increasing the depth gradually, until 95% of the simulations identified all the carriers correctly for various percentages of heterozygous variant carriers .
our method performed better than both the designs in compressed sequencing.
the advantages of our method were significant with large numbers of variant carriers.
the performance of the dna sudoku design was similar to our method when the number of variant carriers was no more than two, but it did not perform well for variants with higher frequencies because of the limited efficiency of the pattern consistency decoding algorithm.
for these cases, more pools are required for the dna sudoku design than for both our method and compressed sequencing.
the dna sudoku design is hard to specify for any number of pools.
therefore, we compared only the performance of compressed sequencing with that of our method to identify four heterozygous variant carriers among  <dig> diploid samples by using the same amounts of pools and sequencing throughput .
our method performed better for most scenarios, especially when the sequencing throughput was limited.
simulation experiment
we applied our method to identify variant carriers among  <dig> simulated e. coli strains.
illumina sequencing reads of two e. coli strains were downloaded from genbank’s short read archive  and bgi’s ftp site .
we treated the o157:h <dig> strain as the variant carrier and the o104:h <dig> strain as the normal sample.
bowtie <dig> .
<dig>  <cit>  was used to map the o157:h <dig> reads to the o104:h <dig> genome, and samtools  <dig> .
<dig>  <cit>  was used to call single base mutations.
because the mean depth was 134× for o157:h <dig>  mutations with depths lower than  <dig> or higher than  <dig> were removed to control the quality; the remaining  <dig> mutations were used in the analysis.
we conducted three simulation experiments to validate the ability of our method to identify carriers of variants with frequencies ranging from  <dig>  to  <dig> %.
based on the results in table  <dig>  we designed the pooling matrix and sequencing depth so that 95% of the simulations correctly identified the variant carriers.
next, pooled sequencing was conducted by selecting reads randomly from the data set and mixing them in silico.
considering up to 5% average noise in the dna quantities of each sample in the pooling procedure, the number of reads for each sample was revised with a random coefficient following a gaussian distribution to simulate reality.
bowtie was used to map pooled reads, and perl scripts were used to count the reads with and without variants that were mapped at the loci of variants.
after the decoding procedure, variant carriers could be identified correctly for  <dig> – <dig> % variants.
this result was consistent with the design capability .
CONCLUSIONS
here, an efficient method that harnesses random k-set pool designs and massively parallel sequencing technologies to identify rare variant carriers is presented.
the parameters of the random k-set pool design can be selected appropriately depending on an indicative probability.
according to the depth model for pooled sequencing, the optimal cut-off value to separate negative pools from positive pools was designed.
taking advantage of the quantitative information in the sequencing results, a heuristic bayesian decoding algorithm to identify the variant carriers was developed.
compared with the dna sudoku design and compressed sequencing, our method showed potential advantages, especially in decreasing the required data throughput.
finally, we applied our method to identify variant carriers among  <dig> simulated e. coli strains using simulated pools and illumina sequencing data.
our method successfully identified variant carriers at reduced experimental costs.
for the accurate identification of variant carriers, the sequencing depth and pool number must be adequate to overcome sequencing errors and mixing bias.
considering the trade-off between the pool number and data throughput, substantial simulations need to be performed to verify whether a design is capable of identifying all the variant carriers correctly.
because the overall cost of overlapping pool sequencing stems from the pooling procedure, library construction, and data production, the optimal design depends on the whole cost.
our decoding algorithm identifies the variant carriers by maximizing the posterior probability, and does not depend too much on the rarity of variants.
therefore, our approach can succeed even for low frequency variants.
furthermore, the sequencing qualities that indicate the sequencing error probabilities could be integrated into the calculation of the posterior probability in the decoding procedure to improve the accuracy.
compared with compressed sequencing, our decoding procedure was very time-consuming because of the substantial calculation of the posterior probability.
this will be improved in future work.
further improvement could be made with a reasonable depth model.
although in many studies negative binomial distribution rather than poisson distribution has been used to fit the sequencing depth, numerous different models exist.
we could not determine which model fit the depth distribution best because, in previous studies, these models have not been compared.
additionally, different sequencing procedures and platforms, such as exome sequencing and whole genome sequencing, produce distinct depth distributions.
we aim to employ a better depth model to improve the performance of our method.
our method has the advantage over compressed sequencing because required data throughput is reduced.
however, because each sample is sequenced multiple times, the required data throughput is still substantial.
third-generation sequencing technologies  <cit> , which significantly reduce the cost for data production, may help to overcome this drawback.
we expect that our method could be applied not only in sequencing experiments but also in other fields as long as the pooled experimental results contain quantitative information about the number of positive samples.
appendix 1: derivations of eq.
and eqs.
–
eq.
: the indicative probability pi is the probability that positive pools are more than the sum of unresolved negative samples and real positive samples.
if, np is the number of positive pools, n¯ is the number of unresolved negative samples, and d is the number of positive samples, then pi can be written as a.
pi=∑i=p_minp_maxpnp=ipn¯+d≤i
where p_min and p_max are the minimum and maximum number of positive pools, respectively.
because np=i indicates that there are t - i negative pools, p can be formulated as a. because pn¯+d≤i = pn¯≤i−d, pn¯+d≤i can be formulated as a. after integrating a–a, pi can be formulated as a.
pnp=i=tt−inegt−i
pn¯+d≤i=∑j=0i−dpn¯=j
pi=∑i=p_minp_maxtt−inegt−i∑j=0i−dpn¯=j
eq.
and eq.
: these equations define the probabilities that nv reads containing variants are observed in a negative pool ), and nn reads without variants are observed in a negative pool ), respectively.
briefly, pnv can be written as a.
pnvnv=∑i=nv∞pipenv|i
where p is the probability that i reads are obtained, and pe is the probability that nv errors occur among these i reads.
because the depth follows a negative binomial distribution and sequencing errors follow a binomial distribution, these two probabilities can be formulated as a and a. in a, d and r are the mean depth of coverage for pooled sequencing and the variance/mean ratio, respectively.
in a, perror is the mean sequencing error rate.
pi=nbi;dr− <dig> r
penv|i=binnv;i,perror
after integrating a–a, pnv can be formulated as a.
pnvnv=∑i=nv∞nbi;dr− <dig> rbinnv;i,perror
the derivation of the formula for pnn ) is similar to the derivation for pnv.
pnnnn=∑i=nn∞nbi;dr− <dig> rbini−nn;i,perror
eq.
and eq.
: these equations define the probability that nv reads containing variants are observed in a positive pool ) and nn reads without variants are observed in a positive pool ), respectively.
the observations of a variant in a positive pool consist of two parts: real variants from variant chromosomes, and false variants resulting from sequencing errors.
briefly, ppv can be written as a where pn stands for the probability that x reads containing variants stemming from the sequencing results of normal chromosomes, and pp denotes the probability that o - x reads contain variants from variant chromosomes.
ppvnv=∑x=0nvpnxppnv−x
by applying a similar procedure to the one used to obtain a and a, pn and pp can be formulated as a and a. the only difference is the mean sequencing depth of coverage.
because the percentages of variant chromosomes and normal chromosomes are p and  <dig> - p, respectively, the mean depths of coverage for sequencing variant chromosomes and normal chromosomes are pd and d, respectively.
pnx=∑i=x∞nbbinx;i,perror
ppo−x=∑j=o−x∞nbbinj−o+x;j,perror
in the same way, ppv can be obtained by integrating a–a, which is shown as a.
ppvnv=∑x=0nv∑i=x∞nbbinx;i,perror×∑j=nv−x∞nbbinj−nv+x;j,perror
similarly, ppn can be obtained as shown in a.
ppnnn=∑x=0nn∑i=x∞nbbini−x;i,perror×∑j=nn−x∞nbbinnn−x;j,perror
competing interests
the authors declare that they have no competing interests.
authors’ contributions
ccc, cl, and xs developed the method.
ccc performed the experiments and wrote the manuscript.
xs revised the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1: figure s1
least depth to make the expected number of errors in classifying pools smaller than  <dig>  figure s <dig>  correlation between the pi value and the correct decoding rate for different scenarios.
figure s <dig>  comparison of the correct decoding rate between our method and compressed sequencing using the second kind of design.
table s <dig>  optimal column weight for various numbers of pools to identify four variant carriers among  <dig> samples.
table s <dig>  least data throughput required to achieve a 95% correct decoding rate in the identification of heterozygous variant carriers among  <dig> diploid samples under the condition that only  <dig> pools are allowed.
click here for file
acknowledgements
this work was supported by the national basic research program of china  and the national natural science foundation of china .

 rare variants are responsible for a large portion of the heritability of some common complex human diseases. genome-wide association studies have focused on the contribution of variants of low minor allele frequency, or of rare variants. the functional and evolutionary impacts of rare variants have been reported; therefore, large-scale screening for disease-associated rare variants is becoming increasingly important. one major application of rare variants genotyping is in screens for rare genetic disorders such as tay–sachs disease and thalassemia. because of the extremely low frequency of rare variants, sample sizes must be large enough to guarantee efficient observations. the cost of dna sequencing has dropped dramatically with the introduction of the massively parallel sequencing technologies. however, identifying rare variant carriers by sequencing individual samples separately remains prohibitively expensive because of the huge challenge of constructing sequencing libraries for thousands of samples. barcoding has been developed as a powerful approach to cost-effectively determine the genotype of each individual. to further reduce the cost of large-scale screens for rare variant carriers, several techniques based on the group testing theory and compressed sensing to construct overlapping pool sequencing strategies have been used. these strategies have helped decrease the sequencing times for rare variant carrier identification and further lower the cost. because a large number of samples are pooled together and then sequenced, overlapping pool sequencing uses fewer pools to identify rare variant carriers among large numbers of individuals. thus, overlapping pool sequencing can vastly reduce the time and cost for dna library preparation. some overlapping pool sequencing programs return true/false values after testing each pool; this scheme was adopted by erlich et al., prabhu and pe’er, and cao et al., who used the number of normal and variant reads in each pool to determine whether a pool contained carriers. however, the quantitative information in the sequencing data, which can be used to estimate the percentage of variant chromosomes in a pool, is missed in these methods. quantitative group testing is an alternative scheme that takes the quantitative information into account, thus rare variant carriers can be identified efficiently. we propose an efficient random overlapping pool sequencing strategy with quantitative group testing for the identification of rare variant carriers using massively parallel sequencing data. because of the excellent performance of random designs in classic group testing, we employed a random k-set pool design to mix samples. the parameters of the random k-set pool design can be selected appropriately according to an indicative probability value. based on a depth model for pooled sequencing, we calculated the optimal cut-off of the number of reads containing variants to distinguish pools containing variant carriers from those that do not. using the quantitative information contained in the sequencing data, we designed a heuristic bayesian decoding algorithm to identify variant carriers accurately. compared with the dna sudoku algorithm and compressed sequencing, our method required less data throughput. finally, we applied our method to identify variant carriers among simulated escherichia coli strains using simulated pools and publicly available illumina sequencing data. the results showed that our method could successfully identify carriers for – % of the variants with frequencies ranging from to %. methods random k-set pool design random k-set pool designs were first proposed by bruno et al. for efficient dna clone library screening. in such a design, each clone is pooled in exact k pools that are chosen with equal chance. random k-set pool designs are easy to specify for any number of pools and are known to be efficient in classic group testing, but they have not been used in real situations, partly because of the presence of different sample sets with identical test sets that are indistinguishable when the test results are qualitative. however, this problem can be overcome by quantitative tests such as sequencing experiments. for n samples containing d positive samples, the basic objective of a random k-set pool design is to identify all the positive samples with a small number of pooled tests. in such a design, each sample occurs in exact k pools, and a pool is defined as positive only when it contains at least one positive sample; otherwise, it is defined as negative. for a random k-set pool design with t pools, hwang calculated the probability that a given set of i pools is a negative one ) and the expected number of negative pools. negi=∑h=it−1h−it−ih−it−hkd/tkd×i∈n_min,n_max e=∑i=n_minn_maxitinegi where n_min and n_max are the minimum and maximum number of negative pools, respectively, and h is a temporary variable. n_max = t - k, and n_min =  or t - dk. to evaluate the performance of random designs, barillot et al. proposed that the number of unresolved negative samples can be taken as a criterion. an unresolved negative sample is defined as a negative sample that occurred only in positive pools, as a result, its status can only be confirmed by testing it separately. negative samples that are contained in at least one negative pool can confidently be determined as negative; therefore, hwang calculated the expectation ) and probability distributions ) for the number of unresolved negative samples. en¯=n−d∑i=n_minn_maxtinegit−ik−d/tk−d pn¯=j=∑i=n_minn_maxtinegit−ik−dj×) to evaluate the performance of random k-set designs in quantitative group testing. pi is the probability that positive pools are more than the sum of unresolved negative samples and real positive samples; therefore, designs with high pis will perform better than designs with low pis. pi=∑i=p_minp_maxtt−inegt−i∑j=0i−dpn¯=j where p_min and p_max are the minimum and maximum number of positive pools, respectively, p_min = t - n_max, and p_max = t - n_min. the derivation of eq. is given in appendix optimal cut-off value for pooled sequencing for overlapping pool sequencing, the dna samples are mixed and then sequenced. samples from variant carriers are treated as positive and samples from normal individuals are treated as negative. to distinguish positive pools containing variant carriers from negative pools consisting of normal individuals, the cut-off for the number of reads containing variants must be selected carefully to declare whether a pool contains carriers or not. ideally, the cut-off value must guarantee that the minimum error rates are obtained, including false-positive and false-negative rates. the variation of sequencing depth distribution is significantly greater than the mean; therefore, in recent studies, negative binomial distribution rather than poisson distribution has been used to model sequencing depth. following the study reported by miller et al., we used a negative binomial model to estimate the sequencing depth distribution. given the average sequencing depth d, the depth that represents the number of times a nucleotide is read follows a negative binomial distribution nbdr− r where r is the variance/mean ratio; r is related to sequencing platforms and genomes and can be estimated from sequencing results ). pdepth=i=nbi;dr− r for a pool with n samples, given sequencing depth d and average sequencing error rate perror, the probabilities pnv and ppv that nv reads containing variants are observed in negative pools and positive pools, respectively, are given by eqs. and. for a locus sequenced i times, the number of sequencing errors follows a binomial distribution bin. pnvnv=∑i=nv∞nbi;dr− rbinnv;i,perror ppvnv=∑x=0nv∑i=x∞nbbinx;i,perror×∑j=nv−x∞nbbinj−nv+x;j,perror similarly, the probabilities pnn and ppn that nn reads without variants are observed in negative pools and positive pools, respectively, are given by eqs. and. pnnnn=∑i=nn∞nbi;dr− rbini−nn;i,perror ppnnn=∑x=0nn∑i=x∞nbbini−x;i,perror×∑j=nn−x∞nbbinnn−x;j,perror where p is the percentage of variant chromosomes in the pool. in a positive pool with n diploid dna samples and c heterozygous variant carriers, ignoring the bias in mixing samples, the percentage of the variant chromosomes is p = c/, while for haploid samples p = c/n. the derivations of eqs. – are given in appendix given pnv and ppv, the formula for the false-positive rate fp and false-negative rate fn in classifying pools with a cut-off value t can be constructed and ). fp=∑i=t∞pnvi fn=∑i=0t−1ppvi the optimal cut-off t can be defined as the value that minimizes the maximum values of fn and fp. t=argminmaxfn,fp|t∈ d decoding algorithm our decoding procedure involves two steps. in the first step, we identify negative pools according to the sequencing results and cut-off values for each pool. samples that participate in any negative pools are classified as from normal individuals. then, separate the real positive samples from unresolved negative samples according to the quantitative information in the sequencing results. the probability of observing the sequencing results under the exact set of variant carriers should be greater than taking other set of samples as variant carriers. assuming a is the set of variant carriers, the probability that the sequencing result o is observed is given by eq.. po|a=∏i=1tpoiv,oin|a where oiv and oin are the number of reads with and without variants in the ith pool. given that a is the set of variant carriers, if the ith pool is negative, then p = pnvpnn; otherwise, p = ppvppn. for the second step, we designed a bayesian decoding algorithm based on eq.. to identify variant carriers among haploid samples, after excluding resolved negative samples, the rest of the samples form a set a0 = {s …,sc}. first, assuming that all the samples in a are variant carriers, we calculate the probability of observing the sequencing results and denote it as pmax_ next, replace one positive sample in set a with a negative sample in turn and repeat the probability calculation. denote the derived set that results in the maximum probability as a and the corresponding probability as pmax_ replace a with a and repeat this step until ac and pmax_c are obtained. finally, the set ai that results in the maximum corresponding probability pmax_i is defined as the set of variant carriers. these steps are written as algorithm two kinds of positive samples need to be considered while identifying variant carriers among diploid samples: heterozygous carriers and homozygous carriers. first, suppose that there are only heterozygous variant carriers; this is analogous to finding variant carriers among haploid samples. then we present algorithm which is very similar to algorithm to recognize heterozygous and homozygous variant carriers. our decoding procedure to identify variant carriers among diploid samples is summarized in algorithm first, we suppose that only heterozygous variant carriers exist and run algorithm to find the set of variant carriers c then, run algorithm to recognize heterozygous and homozygous variant carriers..  optimal cut-off value to approximate the sequencing depth distribution for data obtained by using illumina sequencing platforms, miller et al. found that the negative binomial distribution with the variance/mean ratio r of performed much better than the poisson distribution. therefore, r was set to in our simulation unless otherwise stated.for a pool consisting of diploid samples, we calculated the false-positive and false-negative probabilities with different cut-off values when only one heterozygous variant carrier was allowed in the positive pool. the average sequencing error rate and whole depth were set to and 600×, respectively. the results verified the importance of selecting an appropriate cut-off value. with smaller cut-off values, the probability of misclassifying negative pools as positive is high. with higher cut-off values, some positive pools will be misclassified. both will lower the speed and accuracy of decoding. from the results we can infer that the optimal cut-off value is 14; here both the false-negative and false-positive probabilities are very low.in most studies, because of the rarity of positive samples, the number of positive and negative pools is unequal. therefore, selecting a cut-off value based on the expected number of errors in classifying pools is a more practical approach. for instance, when there are negative and positive pools mentioned above, the optimal cut-off value is. in the following simulation experiments, we adopted this kind of scheme unless otherwise stated. as mentioned, the variance/mean ratio r is related to the sequencing platforms and genomes. because the observed variation is significantly greater than the mean of the sequencing depth, r is larger than different values for r yield distinct distributions. the pooling design mentioned above consisting of negative and positive pools was used to estimate the effect of r on our methods. we calculated the least depth required for each pool to make the expected number of errors in classifying pools smaller than by increasing the depth gradually. from the results, we can see that our method performed even better for smaller r, which required less data throughput. performance of the pipeline to evaluate the performance of our method, we conducted substantial simulations to identify four heterozygous variant carriers among haploid samples. replicates were done for each pair of design parameters: pool number t and column weight k. the pooling matrix was designed by collecting random binary vectors with length t and weight k, meaning that each sample was mixed in k of t pools. identical vectors were deleted and the steps were repeated until vectors were obtained to form the matrix. we used the random function in perl to simulate the number of reads with and without variants in pooled sequencing. sequencing error and mixing bias were added to the simulation procedure to bring it closer to a real situation. sequencing error follows a binomial distribution in sequencing results, and in the simulation the average sequencing error rate was set as. mixing bias is caused by the practical difficulty of mixing exactly equal amounts of dna samples. based on the study conducted by shental et al., a random variable following the gaussian distribution was added to each non-zero element of the pooling matrix to simulate mixing bias. the standard deviation of the gaussian distribution was, reflecting up to 5% average noise in the mixed quantities of each sample. after simulating the pooled sequencing results containing the sequencing errors and mixing bias, the genotypes of the samples were reconstructed using our decoding algorithm. the correct decoding rates, namely the percentages of simulations that identified all the variant carriers correctly, were determined. the results showed that either a too large or too small k negatively influenced the correct decoding rates. moreover, a large k meant more pooling procedures and increased experimental costs. therefore, a proper column weight k is critical for conducting experiments successfully and efficiently. cost-effective overlapping pool sequencing the column weight k denotes the mixing times for each sample in a random k-set pool design. for a given number of pools, a k that is too large or too small will lower the decoding accuracy. we designed an indicative probability pi, which reflects the performance of random k-set designs that could be used to choose the optimal column weight k. we calculated the correct decoding rates for different k under the condition that pools were allowed to identify four heterozygous variant carriers among diploid samples by conducting replicates for each k. next, pi was computed based on eq. and the results are shown in figure  a strong correlation was observed between pi and the correct decoding rate, especially before the correct decoding rate reached the saturation point. the pi values and correct decoding rates for identifying variant carriers were also obtained under different scenarios. all the scenarios showed strong correlations between pi and the correct decoding rate before the correct decoding rate reached the saturation point. for a given pool number t, we defined the optimal k as the minimum that obtains the maximum pi value, which could maximize the correct decoding rate. designs with optimal k require fewer pools or lower sequencing depth. in practice, the optimal k is selected by calculating the pi value without the need to conduct simulations, thereby greatly reducing the computational time required. next, we conducted a series of simulated overlapping pool sequencing experiments with 20– pools and – × overall sequencing data throughput. one thousand replicates were conducted for each scenario, and the column weight was set as the optimal value. the correct decoding rates were low when few pools or data throughput were used. however, adequate pools and data throughput achieved higher accuracy but increased the cost, which conflicted with our motivation in this study. there is a trade-off between the number of pools and data throughput. hence, numerous simulations need to be performed to verify whether a pool number and data throughput pair can succeed in achieving high accuracy. clearly, the optimal design parameters should be selected based on the whole cost of the sequencing experiment. for a given population with diploid individuals containing one heterozygous variant carrier, we generated several candidate designs in which over 95% of the simulations correctly identified the variant carrier. the sequencing region was set to 30 mb to fit the human exome sequencing project. the cost of a sequencing experiment includes library construction and data production. using the cost model from our previous work, we inferred that the lowest cost design was design ii in table  compared with sequencing separately, which requires sequencing depths of × for each sample to obtain correct decoding rates over 95%, our method can save at least 50% of the cost. with the same procedure, we generated the most cost-effective designs for variants with different frequencies and different sequencing region sizes. for smaller sequencing regions and variants with lower frequencies, there are greater cost reductions with our method compared with those for larger regions and variants with higher frequencies. a aaverage value from five simulations. gb is short for gigabases. data throughput is the sequencing depth multiplied by the length of the sequencing region. bsequencing separately is the strategy when each sample is sequenced independently. all the candidate designs can identify the variant carrier correctly in 95% of the simulations. the costs were estimated at $ for one library preparation and $ for gb of data. the sequencing region for haploid samples was set as 5 mb to fit the average length of the bacterial genome. the sequencing region for diploid samples was set as 30 mb to fit the human exome sequencing. comparisons with current methods in benefiting from the chinese remainder theorem, erlich et al. put forward the dna sudoku design for overlapping pool sequencing. a pattern consistency decoding algorithm was also developed by erlich et al. to identify variant carriers with the dna sudoku design. in shental et al. developed a method called compressed sequencing to identify rare variants and their carriers by borrowing techniques from compressed sensing. two designs were proposed in compressed sequencing: one used pools with a random half of the samples and the other used pools with sizes equal to the square root of the number of samples. we compared the performance of our method in identifying rare variant carriers with the performances of these two methods. to identify variant carriers in diploid samples, the dna sudoku design with parameter d0 =  was employed that required pools. to maintain consistency, only pools were allowed for the random k-set pool design and compressed sequencing. since the expected number of positive and negative pools was not clear for the dna sudoku design, the cut-off value for the number of reads containing variants to declare a pool to be positive was set based on the false-negative and false-positive rates, and not on the expected number of errors in the classifying pools. with pools, we computed the least sequencing data throughput required for all the methods by increasing the depth gradually, until 95% of the simulations identified all the carriers correctly for various percentages of heterozygous variant carriers. our method performed better than both the designs in compressed sequencing. the advantages of our method were significant with large numbers of variant carriers. the performance of the dna sudoku design was similar to our method when the number of variant carriers was no more than two, but it did not perform well for variants with higher frequencies because of the limited efficiency of the pattern consistency decoding algorithm. for these cases, more pools are required for the dna sudoku design than for both our method and compressed sequencing. the dna sudoku design is hard to specify for any number of pools. therefore, we compared only the performance of compressed sequencing with that of our method to identify four heterozygous variant carriers among diploid samples by using the same amounts of pools and sequencing throughput. our method performed better for most scenarios, especially when the sequencing throughput was limited. simulation experiment we applied our method to identify variant carriers among simulated e. coli strains. illumina sequencing reads of two e. coli strains were downloaded from genbank’s short read archive and bgi’s ftp site. we treated the o157:h strain as the variant carrier and the o104:h strain as the normal sample. bowtie. was used to map the o157:h reads to the o104:h genome, and samtools. was used to call single base mutations. because the mean depth was 134× for o157:h mutations with depths lower than or higher than were removed to control the quality; the remaining mutations were used in the analysis. we conducted three simulation experiments to validate the ability of our method to identify carriers of variants with frequencies ranging from to %. based on the results in table  we designed the pooling matrix and sequencing depth so that 95% of the simulations correctly identified the variant carriers. next, pooled sequencing was conducted by selecting reads randomly from the data set and mixing them in silico. considering up to 5% average noise in the dna quantities of each sample in the pooling procedure, the number of reads for each sample was revised with a random coefficient following a gaussian distribution to simulate reality. bowtie was used to map pooled reads, and perl scripts were used to count the reads with and without variants that were mapped at the loci of variants. after the decoding procedure, variant carriers could be identified correctly for – % variants. this result was consistent with the design capability.  here, an efficient method that harnesses random k-set pool designs and massively parallel sequencing technologies to identify rare variant carriers is presented. the parameters of the random k-set pool design can be selected appropriately depending on an indicative probability. according to the depth model for pooled sequencing, the optimal cut-off value to separate negative pools from positive pools was designed. taking advantage of the quantitative information in the sequencing results, a heuristic bayesian decoding algorithm to identify the variant carriers was developed. compared with the dna sudoku design and compressed sequencing, our method showed potential advantages, especially in decreasing the required data throughput. finally, we applied our method to identify variant carriers among simulated e. coli strains using simulated pools and illumina sequencing data. our method successfully identified variant carriers at reduced experimental costs. for the accurate identification of variant carriers, the sequencing depth and pool number must be adequate to overcome sequencing errors and mixing bias. considering the trade-off between the pool number and data throughput, substantial simulations need to be performed to verify whether a design is capable of identifying all the variant carriers correctly. because the overall cost of overlapping pool sequencing stems from the pooling procedure, library construction, and data production, the optimal design depends on the whole cost. our decoding algorithm identifies the variant carriers by maximizing the posterior probability, and does not depend too much on the rarity of variants. therefore, our approach can succeed even for low frequency variants. furthermore, the sequencing qualities that indicate the sequencing error probabilities could be integrated into the calculation of the posterior probability in the decoding procedure to improve the accuracy. compared with compressed sequencing, our decoding procedure was very time-consuming because of the substantial calculation of the posterior probability. this will be improved in future work. further improvement could be made with a reasonable depth model. although in many studies negative binomial distribution rather than poisson distribution has been used to fit the sequencing depth, numerous different models exist. we could not determine which model fit the depth distribution best because, in previous studies, these models have not been compared. additionally, different sequencing procedures and platforms, such as exome sequencing and whole genome sequencing, produce distinct depth distributions. we aim to employ a better depth model to improve the performance of our method. our method has the advantage over compressed sequencing because required data throughput is reduced. however, because each sample is sequenced multiple times, the required data throughput is still substantial. third-generation sequencing technologies, which significantly reduce the cost for data production, may help to overcome this drawback. we expect that our method could be applied not only in sequencing experiments but also in other fields as long as the pooled experimental results contain quantitative information about the number of positive samples. appendix 1: derivations of eq. and eqs. – eq.: the indicative probability pi is the probability that positive pools are more than the sum of unresolved negative samples and real positive samples. if, np is the number of positive pools, n¯ is the number of unresolved negative samples, and d is the number of positive samples, then pi can be written as a. pi=∑i=p_minp_maxpnp=ipn¯+d≤i where p_min and p_max are the minimum and maximum number of positive pools, respectively. because np=i indicates that there are t - i negative pools, p can be formulated as a. because pn¯+d≤i = pn¯≤i−d, pn¯+d≤i can be formulated as a. after integrating a–a, pi can be formulated as a. pnp=i=tt−inegt−i pn¯+d≤i=∑j=0i−dpn¯=j pi=∑i=p_minp_maxtt−inegt−i∑j=0i−dpn¯=j eq. and eq.: these equations define the probabilities that nv reads containing variants are observed in a negative pool ), and nn reads without variants are observed in a negative pool ), respectively. briefly, pnv can be written as a. pnvnv=∑i=nv∞pipenv|i where p is the probability that i reads are obtained, and pe is the probability that nv errors occur among these i reads. because the depth follows a negative binomial distribution and sequencing errors follow a binomial distribution, these two probabilities can be formulated as a and a. in a, d and r are the mean depth of coverage for pooled sequencing and the variance/mean ratio, respectively. in a, perror is the mean sequencing error rate. pi=nbi;dr− r penv|i=binnv;i,perror after integrating a–a, pnv can be formulated as a. pnvnv=∑i=nv∞nbi;dr− rbinnv;i,perror the derivation of the formula for pnn ) is similar to the derivation for pnv. pnnnn=∑i=nn∞nbi;dr− rbini−nn;i,perror eq. and eq.: these equations define the probability that nv reads containing variants are observed in a positive pool ) and nn reads without variants are observed in a positive pool ), respectively. the observations of a variant in a positive pool consist of two parts: real variants from variant chromosomes, and false variants resulting from sequencing errors. briefly, ppv can be written as a where pn stands for the probability that x reads containing variants stemming from the sequencing results of normal chromosomes, and pp denotes the probability that o - x reads contain variants from variant chromosomes. ppvnv=∑x=0nvpnxppnv−x by applying a similar procedure to the one used to obtain a and a, pn and pp can be formulated as a and a. the only difference is the mean sequencing depth of coverage. because the percentages of variant chromosomes and normal chromosomes are p and - p, respectively, the mean depths of coverage for sequencing variant chromosomes and normal chromosomes are pd and d, respectively. pnx=∑i=x∞nbbinx;i,perror ppo−x=∑j=o−x∞nbbinj−o+x;j,perror in the same way, ppv can be obtained by integrating a–a, which is shown as a. ppvnv=∑x=0nv∑i=x∞nbbinx;i,perror×∑j=nv−x∞nbbinj−nv+x;j,perror similarly, ppn can be obtained as shown in a. ppnnn=∑x=0nn∑i=x∞nbbini−x;i,perror×∑j=nn−x∞nbbinnn−x;j,perror competing interests the authors declare that they have no competing interests. authors’ contributions ccc, cl, and xs developed the method. ccc performed the experiments and wrote the manuscript. xs revised the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1: figure s1 least depth to make the expected number of errors in classifying pools smaller than figure s correlation between the pi value and the correct decoding rate for different scenarios. figure s comparison of the correct decoding rate between our method and compressed sequencing using the second kind of design. table s optimal column weight for various numbers of pools to identify four variant carriers among samples. table s least data throughput required to achieve a 95% correct decoding rate in the identification of heterozygous variant carriers among diploid samples under the condition that only pools are allowed. click here for file acknowledgements this work was supported by the national basic research program of china and the national natural science foundation of china.
with the simulated pools and publicly available illumina sequencing data, our method correctly identified the variant carriers for – % variants with the variant frequency ranging from to %. then, using the quantitative information contained in the sequencing results, we designed a heuristic bayesian probability decoding algorithm to identify variant carriers. based on quantitative group testing, we propose an efficient overlapping pool sequencing strategy that allows the efficient recovery of variant carriers in numerous individuals with much lower costs than conventional methods. quantitative group testingrandom k-set pool designoverlapping pool sequencingrare variants we used random k-set pool designs to mix samples, and optimized the design parameters according to an indicative probability. our method performed better than the published dna sudoku design and compressed sequencing, especially in reducing the required data throughput and cost. finally, we conducted in silico experiments to find variant carriers among simulated escherichia coli strains. genome-wide association studies have revealed that rare variants are responsible for a large portion of the heritability of some complex human diseases. although the massively parallel sequencing technologies have greatly reduced the cost of dna sequencing, the identification of rare variant carriers by large-scale re-sequencing remains prohibitively expensive because of the huge challenge of constructing libraries for thousands of samples. recently, several studies have reported that techniques from group testing theory and compressed sensing could help identify rare variant carriers in large-scale samples with few pooled sequencing experiments and a dramatically reduced cost. using the number of reads, variant carriers could be identified precisely even though samples were randomly selected and pooled. based on a mathematical model of sequencing depth distribution, an optimal threshold was selected to declare a pool positive or negative. this highlights the increasing importance of detecting and screening for rare variants.
 genome-wide association studies have revealed that rare variants are responsible for a large portion of the heritability of some complex human diseases. this highlights the increasing importance of detecting and screening for rare variants. although the massively parallel sequencing technologies have greatly reduced the cost of dna sequencing, the identification of rare variant carriers by large-scale re-sequencing remains prohibitively expensive because of the huge challenge of constructing libraries for thousands of samples. recently, several studies have reported that techniques from group testing theory and compressed sensing could help identify rare variant carriers in large-scale samples with few pooled sequencing experiments and a dramatically reduced cost.  based on quantitative group testing, we propose an efficient overlapping pool sequencing strategy that allows the efficient recovery of variant carriers in numerous individuals with much lower costs than conventional methods. we used random k-set pool designs to mix samples, and optimized the design parameters according to an indicative probability. based on a mathematical model of sequencing depth distribution, an optimal threshold was selected to declare a pool positive or negative. then, using the quantitative information contained in the sequencing results, we designed a heuristic bayesian probability decoding algorithm to identify variant carriers. finally, we conducted in silico experiments to find variant carriers among simulated escherichia coli strains. with the simulated pools and publicly available illumina sequencing data, our method correctly identified the variant carriers for – % variants with the variant frequency ranging from to %.  using the number of reads, variant carriers could be identified precisely even though samples were randomly selected and pooled. our method performed better than the published dna sudoku design and compressed sequencing, especially in reducing the required data throughput and cost. quantitative group testingrandom k-set pool designoverlapping pool sequencingrare variants 