BACKGROUND
the population mutation rate  remains one of the most fundamental parameters in genetics, ecology, and evolutionary biology  <cit> .
this interest in θ derives from the fact that this parameter measures the effective size  and whole-locus mutation rate  of a population, which are of great importance in understanding its demography and history.
specifically, θ is a compound parameter that is calculated as the product of 2pneμ .
correspondingly, a number of alternative methods are available to estimate θ from a population sample of allelic sequences  <cit> .
these alternative methods range from relatively simple summary statistics  to full coalescent/mutation models.
indeed, the estimation of θ remains central to even the most complex coalescent/mutation models that are otherwise concerned with the determination of other population genetic parameters .
a population sample of sequences is obtained from interbreeding or potentially interbreeding individuals and is therefore usually associated with a small number of mutations  <cit> .
thus, when estimating θ from a population sample, sequence errors can pose a real problem, since their numbers can begin to approach or even surpass those for the mutations  <cit> .
this problem becomes particularly acute when working with error prone data such as expressed sequence tags , low coverage draft sequences, and other such unfinished products  <cit> .
for example, an error rate of one mistake per every  <dig> nucleotides  will make a significant contribution to the observed variation among sequences that differ because of mutations by < <dig> to 2%.
if uncorrected, such errors can lead to an inflated estimate of θ and even erroneous conclusions about the biology of their population  <cit> .
many sequence errors arise as random accidents that occur during the nucleic acid isolation, cloning/amplification, sequencing, and recording phases of a dna sequencing study  <cit> .
as chance events that are rare , each of these random mistakes will most likely be limited to a single sequence, rather than repeated among two or more different ones within the population sample  <cit> .
thus, these random errors will most likely inflate the number of singletons within the dataset  .
in contrast, these rare chance mistakes will make a much smaller contribution to the shared polymorphisms .
this study relies on the simple premise that the random mistakes of error prone sequences can be avoided by ignoring the singletons within their population sample.
this strategy is first used to obtain independently the same new watterson  <cit>  and tajima  <cit>  estimators of θ, which were recently reported by achaz  <cit>  for error prone sequences.
this approach is then implemented in the recent maximum likelihood  model of knudsen and miyamoto  <cit> , which incorporates various factors of experimental error and design with those for coalescence and mutation.
by relying on only shared polymorphisms, these three new approaches allow for more accurate estimates of θ. however, this greater accuracy comes with a cost as singletons due to actual mutations are ignored along with those due to random errors.
to assess this tradeoff, the three new methods are tested against each other and their original predecessors that count singletons with evolutionary simulations and/or analyses of a real population dataset for the california seahare .
these tests document that these new approaches offer reliable and fast alternatives for the determination of θ from error prone sequences.
RESULTS
three new methods for estimating θ from error prone sequences
infinite sites model
as in their original versions, the new watterson, tajima, and knudsen/miyamoto methods rely on the infinite sites model to accommodate a neutral mutation process  <cit> .
the infinite sites model assumes that only a single mutation can occur at any homologous position of the population sample.
thus, each variable site will be represented by only two bases that subdivide the sampled sequences into two non-overlapping subsets consisting of those with the first nucleotide versus the remainder with the second base.
correspondingly, each mutation will map to a specific branch within the sample genealogy, which thereby partitions the sequences into their two non-overlapping subsets.
for example, the shared polymorphism at position  <dig> in figure  <dig> is attributable to a unique mutation along the internal branch that partitions sequences i and ii from iii, iv, and v.
as random errors are treated as rare chance events, they are also modeled in this study along with the mutations by the infinite sites process.
thus, only a single random error or mutation is allowed at any site of the sampled sequences.
in turn, each random error is limited to a single sequence  in contrast to a mutation that can also result in a shared polymorphism .
the reason is that random errors arise during the experimental determination and recording of individual sequences, whereas mutations occur at specific points within the sample genealogy.
thus, a mutation along an internal branch of the genealogy will result in a new base that will be shared by two or more of its descendant sequences.
new watterson estimator
define ti as the length of time  during which there are exactly i ancestors for n sampled sequences.
standard coalescent theory tells us that:
and
<cit> .
the expected total branch length of the genealogy for n sequences  can now be calculated as:
let ns be the observed number of segregating  sites in the dataset.
under the infinite sites model, ns also counts the number of mutations, since each observed variable site is attributable to a single mutation.
the expected number of mutations per locus per unit of branch length is θ/ <dig>  thus, an estimate of θ can be obtained as:
<cit> .
a new watterson estimator that avoids singletons  can now be derived from equation  by adjusting both its numerator and denominator.
the numerator is easily adjusted by counting only the shared polymorphic sites in the original dataset .
although more complicated, the denominator can also be readily adjusted by including in its calculation only the lengths of the internal branches where shared polymorphisms can arise .
let us again consider a point in the genealogy where there are exactly i ancestors for the n sampled sequences.
looking forward in time, the probability that a particular branch of the genealogy is not chosen for the next split  is .
thus, the probability that this branch remains unbroken to the present is:
by combining equations  and , the total length of the external branches where singletons can occur can now be calculated as:
our equation  is equivalent to equation  of fu and li  <cit> , apart from our use of different symbols and terms and of time as scaled by ne generations .
thus, as previously noted by them, the total length of the external branches where a singleton can occur is independent of the original number of sampled sequences.
fu and li  <cit>  also obtained the variance for the total length of the external branches as their equation .
in an asymmetrical genealogy with a basal split of one versus  sampled sequences, a mutation in the internal branch leading to the common ancestor of the  group will result in a singleton within the dataset  .
thus, the length of this  basal branch  must also be accounted for in the adjustment of the denominator for θ'w.
the weighted length of this  basal branch can be calculated as follows.
first, the chance of this branch is determined as the probability of an  asymmetrical topology or equivalently as the probability that either one of the two basal branches for the genealogy remains unbroken to the present.
according to equation , the latter probability is  <dig> .
second, the unweighted length of this branch is then calculated as the length of the time interval t <dig>  which is  <dig> .
thus, the weighted length of the  basal branch is:
by combining equations , , and , θ'w can now be obtained as:
in words, θ'w ignores the singletons  in the original dataset, while restricting the total branch length calculation to only those internal internodes where a shared polymorphism can arise.
apart from our use of different symbols and terms for calculating the denominator, our equation  is equivalent to equation  of achaz  <cit> .
achaz  <cit>  also derived his equation b <dig> for the calculation of the associated variance for n's .
new tajima estimator
let πij denote the number of observed pairwise differences between sampled sequences i and j .
standard coalescent theory tells us that the expected waiting time for these two sequences to coalesce is ~ ne generations.
thus, the expected value of each πij is θ, whereas that for their sum is:
rearranging equation  leads to:
<cit> .
the tajima and watterson estimators have the same expected value of θ, even though the former is based on the total differences between each sampled sequence pair whereas the latter is obtained from the total number of segregating sites within the sample.
their different summaries of the observed variation are the basis of tajima's d that tests for departures from the standard neutral model  <cit>  .
to avoid random errors, equation  can be adjusted in a manner similar to that used to modify equation  of θw.
first, let π'ij count the number of observed pairwise differences between sequences i and j after the removal of all singletons from their dataset.
then, the expected number of singletons can be calculated from equations  and  as the total length of the external branches and  basal branch  multiplied by θ/ <dig> {i.e.,  θ}.
the removal of each singleton from the original dataset reduces the expected sum in equation  by .
combining these results, we obtain:
rearranging equation  then leads to:
our equation  for θ't is equivalent to equation  of achaz  <cit> , except for our use of different symbols and total π'ij  for the population sample.
achaz  <cit>  also derived his equation b <dig> for the calculation of the associated variance for average π'ij ).
new knudsen/miyamoto model
knudsen and miyamoto  developed a full coalescent/mutation model that accounts for three specific factors of experimental error and design:  for random sequence errors,  for unobserved polymorphisms due to missing data, and  for the uncertain assignment of the multiple sequencing reads for a diploid or polyploid individual to its two or more homologues.
their km model uses recursion to calculate an exact probability of the population sample under the standard fisher  <cit>  and wright  <cit>  model for reproduction  and the infinite sites process for both mutations and random sequence errors  <cit> .
their model relies on ml to estimate θ and the expected number of errors per full length sequence .
at the heart of the km model is equation  of knudsen and miyamoto  <cit> , which is reproduced here as:
the parameters and factors of this equation are defined in table  <dig>  as one works backwards in time, the probability that the next observed event is a specific coalescence is given by the first term in the top line before the double sum.
as indicated by this double sum, if indeed a coalescent event occurs, then it will happen between two sampled and/or ancestral alleles with compatible  sequences.
such sequences are referred to as combinable .
m
s
the first term in the bottom line before the single sum is then for the probability that the next observed event is instead a mutation.
as indicated by this sum, if indeed a mutation occurs, then it will happen to a sampled or ancestral sequence with at least one "singleton" .
here, the definition of a "singleton" is expanded to include the derived mutations of the common ancestors for the different groups of related sampled sequences as well as those for the singletons  of the original dataset .
this expanded use of the term acknowledges that a shared polymorphism under the infinite sites model is due to a unique mutation within the common ancestor of those sampled sequences sharing the derived base .
thus, even though they result in shared polymorphisms among the sampled sequences, these derived mutations are ultimately counted as "singletons" as one works backwards in time.
this expanded definition of a "singleton" allows for the economical use of σi alone to track the mutations of both the original shared polymorphisms and singletons.
in the km model, the available region of each sampled sequence is summarized as a closed interval that is scored over the range of .
the |si| factor then quantifies the total amount of sequence available for this sampled allele.
for example, the  score for the leftmost  sampled sequence in figure 4a indicates that it is lacking the initial 20% of the full multiple alignment.
thus, |si| =  <dig>  for this partial, leftmost, sampled sequence.
in turn, as one works backwards in time, the closed intervals for the available regions of the sequences for common ancestors are calculated as the union of the known lengths for their two immediate descendants.
thus, the closed interval of the common ancestor for the two leftmost sampled sequences in figure 4a is  ∪  = , with |si| =  <dig> .
the purpose of Σi |si| in equation  is to track the total available length of all sampled and/or ancestral sequences for the detection of mutations as one works backwards in time.
the corresponding use of |si| in the bottom line of equation  allows for the adjustment of σi due to unobserved polymorphisms resulting from missing data.
the km model uses only equation  when working with error free sequences.
in turn, this model also relies on equation  of knudsen and miyamoto  <cit>  when dealing with error prone sequences.
this additional equation of the km model assumes that the random errors are uniformly distributed along the sampled sequences and that their total number is poisson distributed with an intensity of λ = ε Σi |si|.
the inclusion of Σi |si| in the calculation of λ allows for the adjustment of this error rate to account for incomplete sampled sequences.
in contrast to its predecessor, the km' model uses only equation  when working with either error prone or error free sequences.
as for the new watterson and tajima estimators, the km' model operates by counting only those internal branches of the genealogy where a mutation will result in a shared polymorphism .
the sampled sequences are rescored as unknowns with empty intervals , as is the sequence of the common ancestor for the  basal group of each asymmetrical genealogy .
correspondingly, |si| is then reset to  <dig>  for these sequences.
the km' model ignores the external and basal branches of the genealogy, where a mutation will result in a singleton, by multiplying σi for each sampled and/or  ancestral sequence by |si| =  <dig>  in equation .
in this way, the km' model discounts the singletons in favor of the shared polymorphisms within the dataset without the use of equation  and its associated parameters  of the original km model.
for sampled sequences without missing data, the previous explanation is complete as to how the km' model corrects for the branches of the genealogy where a mutation will result in a singleton.
however, for incomplete sampled sequences, the "above and below" test of the km' model also becomes necessary to correct for the regions of each common ancestral sequence where a mutation in its internal branch will lead to a singleton, rather than to a shared polymorphism, because of this missing information .
in this test, "below" refers to those sampled sequences that belong to the monophyletic group of the common ancestor in question.
"above" then corresponds to the remaining, more distantly related, sampled sequences.
in figure 4b, the two leftmost sampled sequences are direct descendants of the first, leftmost, common ancestor, whereas the other three are not.
thus, these two sampled sequences lie "below," whereas the other three occur "above" the internal branch for this leftmost common ancestor.
the above and below test checks whether comparable information for a region is known for at least two, descendant, sampled sequences below and at least two, more distantly related, sampled sequences above an internal branch in the genealogy.
if not, then a mutation within this region of the ancestral sequence that corresponds to this internal branch will result in a false singleton within the original dataset.
for example, a mutation within the first 20% of the leftmost common ancestor in figure  <dig> will result in a false singleton of the second sampled sequence rather than in a true shared polymorphism of the first and second.
the problem is that the leftmost sampled sequence is missing comparable information for the detection of this mutation as a shared change in the leftmost common ancestor.
the above and below test corrects for such regions of the common ancestral sequences during the rescoring of their closed intervals and |si|.
as the first 20% of the leftmost common ancestor fails the below half of this test, the km' model disregards this region during the recalculation of its closed interval as  and |si| =  <dig>  .
evolutionary simulations and a. californica dataset
evolutionary simulations
to test the three new procedures, evolutionary simulations were conducted according to standard methods  <cit> .
two hundred datasets apiece were simulated for eight or  <dig> sequences of length  <dig> from a single population under the baseline conditions of the standard fw and infinite sites models with θ =  <dig>   <dig>   <dig>  or  <dig>  in addition to these baseline conditions, sequence errors were introduced as four or eight randomly placed changes among the eight and  <dig> sequences, respectively, for an expected ε of  <dig> .
estimates of θ were then obtained for the  <dig> datasets of each tested combination with the fw model and the original and new watterson estimators, tajima estimators, and km and km' models.
as expected, the fw model and original watterson and tajima estimators overestimate θ when the datasets contain random errors .
in these cases, the random errors are counted as mutations, thereby inflating their estimates of θ. furthermore, these overestimations are greater for the watterson estimator than for the tajima estimator.
the reason is that each singleton makes the same contribution as a shared polymorphism to the number of segregating sites  in equation  of θw, but a smaller one to the sum of the pairwise differences in equation  of θt.
specifically, each singleton is limited in equation  to the  pairwise comparisons of the unique sequence to the  other sequences.
thus, the tajima estimator is less vulnerable than the watterson estimator to random errors even though its vulnerability is still significant.
these results are for the standard fw model  and the original and new watterson estimators , tajima estimators , and km and km' models .
they are summarized as the means ± twice their standard deviations for  <dig> simulated datasets apiece for each of the  <dig> different combinations of θ, n , and ε .
mean estimates that are significantly greater than or less than the true θ are highlighted in italics and boldface, respectively.
no results are provided for the fw and km models with θ =  <dig> or  <dig> , since these calculations are too computationally intensive .
in contrast, the km model underestimates θ when the sequences are errorless .
furthermore, this tendency to underestimate θ is also evident  when the sequences contain errors .
in these situations, mutations are sometimes counted as random errors, thereby deflating their estimates of θ. in cases where few to no random errors are expected , one can first perform a likelihood ratio test to determine if ε =  <dig>  <cit> .
if this null hypothesis cannot be rejected, then the km analysis should be restricted to only equation .
however, if ε >  <dig> according to this likelihood ratio test, then equation  of knudsen and miyamoto  <cit>  is also needed and the user of the km model must remain aware that her/his estimate of θ most likely includes some slight downward bias.
unlike their original versions, the new watterson estimator, tajima estimator, and km' model all consistently recover the true θ in the simulations both with and without errors .
furthermore, these three methods also avoid the tendency of the km model to underestimate θ due to unnecessary or "greedy" parameters for random errors.
in particular, the reliance of the km' model on one less equation  and fewer parameters  makes it much simpler and less prone to over-parameterizations than its predecessor.
the mean estimates of θ for the new watterson estimator, tajima estimator, and km' model are also generally associated with greater standard deviations than those for their original versions .
the proximal reason for these increased standard deviations is that the elimination of singletons results in the loss of some actual mutations along with the random errors.
however, this cost of the three new methods appears to be relatively small, given that their standard deviations are never twice as great as those for their original versions .
indeed, the standard deviations for the three new methods are less than or equal to those of their counterparts in four cases .
perhaps more surprising is that the standard deviations for the new watterson and tajima estimators are comparable to those for the km' model, particularly when θ =  <dig> or  <dig>  these similarities in their standard deviations are surprising given that the km' model is a full ml model.
as summary statistics, the original and new watterson and tajima estimators are all computationally fast, requiring much less than  <dig> cpu second on a  <dig>  ghz pentium  <dig> cpu to analyze each of the simulated datasets.
more importantly, the km' model is much faster than the fw and km models as documented by its completed analyses of the more complex datasets .
in contrast, the fw and km models fail to complete their analyses of these more complex datasets due to time and memory constraints.
the km' model is much faster than the fw and km models, because it relies on only shared polymorphic sites.
less variable positions results in faster coalescences and fewer choices as one works back through the coalescent/mutation recursion of equation .
aplysia dataset
the a. californica dataset was the original motivating force behind the development of the km model  <cit> .
thus, this real dataset was also analyzed with the km' model to test further its performance against that of its predecessor.
this dataset consists of  <dig> sequencing reads for six diploid individuals from a laboratory population of the california seahare at the laboratory for marine bioscience, university of florida .
three cloned inserts were sequenced from each individual as a pair of single sequencing passes starting from both ends of an internal segment of  <dig> base pairs for the protein-coding region of the nuclear fmrf gene.
these pairs of passes overlap in the middle for nine sequences, but at most by only  <dig> bases.
thus, these  <dig> essentially single-pass sequences contain many random errors and some missing data and their assignments to the two homologues of each diploid remain uncertain .
correspondingly, the km analysis of this dataset with  <dig> singletons and  <dig> shared polymorphisms required the full use of its factors for experimental error and design .
a full ml analysis of this complex dataset by the km model proved too time and memory consuming and a two step procedure was therefore adopted instead whereby the number of errors was first ml estimated followed by the determination of θ for this ml value  <cit> .
this heuristic approach resulted in an estimate of θkm =  <dig>  for this sample, which still took more than  <dig> hours to complete on the same cpu as for the simulations .
in contrast, a full ml analysis of this dataset by the km' model using the same factors for experimental design is much faster as it took less than  <dig> hour on this cpu to obtain a similar estimate of θ'km =  <dig> .
correspondingly, nucleotide diversity  for this sample is calculated as  =  <dig>  mutations per site.
to summarize, the km and km' models both support similar estimates of θ for this real dataset, but the latter  is much faster as confirmed by its completed, full, ml analysis.
CONCLUSIONS
expressed sequence tags, low coverage draft sequences, and other such unfinished products are known to contain random errors due to chance accidents during their data collection and recording  <cit> .
however, such sequences often constitute the only available allelic information for a population and methods to deal with their random mistakes are therefore needed to obtain accurate estimates of θ from these error prone data.
this study is based on the simple premise that random sequence errors are distributed as singletons.
thus, one can avoid the random mistakes of error prone sequences by ignoring their singletons in favor of their shared polymorphisms .
this strategy is implemented in the new watterson estimator, tajima estimator, and km' model.
these new methods are all accurate and fast according to their evolutionary simulations and analysis of the real complex dataset for a. californica .
these methods come with the cost of increased standard deviations, but this price appears small or even negligible compared to their advantages of significantly improved accuracy and/or computational speed.
obviously, additional evolutionary simulations and applications to real datasets are now needed to evaluate more fully under what conditions the removal of singletons is warranted in light of this tradeoff.
nevertheless, the current successes with the new watterson estimator, tajima estimator, and km' model support our recommendation that these three methods be given serious consideration when estimating θ from error prone sequences.
unlike the watterson and tajima estimators that represent summary statistics, the km and km' models both constitute full ml models that offer the framework for the further incorporation of other experimental and population genetic factors  <cit> .
in particular, such further developments are encouraged for the km' model in light of its current successes in the evolutionary simulations and a. californica analysis .
for example, biased errors within a sample due to the systematic misreading of specific bases during dna sequencing, the postmortem biochemical degradation of ancient dna, and/or other such sources of error can be accommodated by a finite sites process that allows for repeated mistakes as well as mutations at the same sequence positions  <cit> .
likewise, additional factors to account for other population genetic processes such as recombination  can be accommodated by the use of ancestral recombination graphs  <cit> .
inevitably, these more complex versions of the km' model will require the use of sampling based procedures for their implementation , since the current use of direct ml evaluation will remain practical for only the smaller datasets and simpler models  <cit> .
authors' contributions
bk derived the equations, wrote the computer program for the estimation of θ, and conducted the evolutionary simulations.
mmm provided biological interpretations about the results.
both authors contributed to the design and main ideas of this study, wrote and edited the manuscript, and read and approved its final draft.

 the population mutation rate remains one of the most fundamental parameters in genetics, ecology, and evolutionary biology. this interest in θ derives from the fact that this parameter measures the effective size and whole-locus mutation rate of a population, which are of great importance in understanding its demography and history. specifically, θ is a compound parameter that is calculated as the product of 2pneμ. correspondingly, a number of alternative methods are available to estimate θ from a population sample of allelic sequences. these alternative methods range from relatively simple summary statistics to full coalescent/mutation models. indeed, the estimation of θ remains central to even the most complex coalescent/mutation models that are otherwise concerned with the determination of other population genetic parameters. a population sample of sequences is obtained from interbreeding or potentially interbreeding individuals and is therefore usually associated with a small number of mutations. thus, when estimating θ from a population sample, sequence errors can pose a real problem, since their numbers can begin to approach or even surpass those for the mutations. this problem becomes particularly acute when working with error prone data such as expressed sequence tags, low coverage draft sequences, and other such unfinished products. for example, an error rate of one mistake per every nucleotides will make a significant contribution to the observed variation among sequences that differ because of mutations by < to 2%. if uncorrected, such errors can lead to an inflated estimate of θ and even erroneous conclusions about the biology of their population. many sequence errors arise as random accidents that occur during the nucleic acid isolation, cloning/amplification, sequencing, and recording phases of a dna sequencing study. as chance events that are rare, each of these random mistakes will most likely be limited to a single sequence, rather than repeated among two or more different ones within the population sample. thus, these random errors will most likely inflate the number of singletons within the dataset. in contrast, these rare chance mistakes will make a much smaller contribution to the shared polymorphisms. this study relies on the simple premise that the random mistakes of error prone sequences can be avoided by ignoring the singletons within their population sample. this strategy is first used to obtain independently the same new watterson and tajima estimators of θ, which were recently reported by achaz for error prone sequences. this approach is then implemented in the recent maximum likelihood model of knudsen and miyamoto, which incorporates various factors of experimental error and design with those for coalescence and mutation. by relying on only shared polymorphisms, these three new approaches allow for more accurate estimates of θ. however, this greater accuracy comes with a cost as singletons due to actual mutations are ignored along with those due to random errors. to assess this tradeoff, the three new methods are tested against each other and their original predecessors that count singletons with evolutionary simulations and/or analyses of a real population dataset for the california seahare. these tests document that these new approaches offer reliable and fast alternatives for the determination of θ from error prone sequences.  three new methods for estimating θ from error prone sequences infinite sites model as in their original versions, the new watterson, tajima, and knudsen/miyamoto methods rely on the infinite sites model to accommodate a neutral mutation process. the infinite sites model assumes that only a single mutation can occur at any homologous position of the population sample. thus, each variable site will be represented by only two bases that subdivide the sampled sequences into two non-overlapping subsets consisting of those with the first nucleotide versus the remainder with the second base. correspondingly, each mutation will map to a specific branch within the sample genealogy, which thereby partitions the sequences into their two non-overlapping subsets. for example, the shared polymorphism at position in figure is attributable to a unique mutation along the internal branch that partitions sequences i and ii from iii, iv, and v. as random errors are treated as rare chance events, they are also modeled in this study along with the mutations by the infinite sites process. thus, only a single random error or mutation is allowed at any site of the sampled sequences. in turn, each random error is limited to a single sequence in contrast to a mutation that can also result in a shared polymorphism. the reason is that random errors arise during the experimental determination and recording of individual sequences, whereas mutations occur at specific points within the sample genealogy. thus, a mutation along an internal branch of the genealogy will result in a new base that will be shared by two or more of its descendant sequences. new watterson estimator define ti as the length of time during which there are exactly i ancestors for n sampled sequences. standard coalescent theory tells us that: and. the expected total branch length of the genealogy for n sequences can now be calculated as: let ns be the observed number of segregating sites in the dataset. under the infinite sites model, ns also counts the number of mutations, since each observed variable site is attributable to a single mutation. the expected number of mutations per locus per unit of branch length is θ/ thus, an estimate of θ can be obtained as:. a new watterson estimator that avoids singletons can now be derived from equation by adjusting both its numerator and denominator. the numerator is easily adjusted by counting only the shared polymorphic sites in the original dataset. although more complicated, the denominator can also be readily adjusted by including in its calculation only the lengths of the internal branches where shared polymorphisms can arise. let us again consider a point in the genealogy where there are exactly i ancestors for the n sampled sequences. looking forward in time, the probability that a particular branch of the genealogy is not chosen for the next split is. thus, the probability that this branch remains unbroken to the present is: by combining equations and, the total length of the external branches where singletons can occur can now be calculated as: our equation is equivalent to equation of fu and li, apart from our use of different symbols and terms and of time as scaled by ne generations. thus, as previously noted by them, the total length of the external branches where a singleton can occur is independent of the original number of sampled sequences. fu and li also obtained the variance for the total length of the external branches as their equation. in an asymmetrical genealogy with a basal split of one versus sampled sequences, a mutation in the internal branch leading to the common ancestor of the group will result in a singleton within the dataset. thus, the length of this basal branch must also be accounted for in the adjustment of the denominator for θ'w. the weighted length of this basal branch can be calculated as follows. first, the chance of this branch is determined as the probability of an asymmetrical topology or equivalently as the probability that either one of the two basal branches for the genealogy remains unbroken to the present. according to equation, the latter probability is. second, the unweighted length of this branch is then calculated as the length of the time interval t which is. thus, the weighted length of the basal branch is: by combining equations,, and, θ'w can now be obtained as: in words, θ'w ignores the singletons in the original dataset, while restricting the total branch length calculation to only those internal internodes where a shared polymorphism can arise. apart from our use of different symbols and terms for calculating the denominator, our equation is equivalent to equation of achaz. achaz also derived his equation b for the calculation of the associated variance for n's. new tajima estimator let πij denote the number of observed pairwise differences between sampled sequences i and j. standard coalescent theory tells us that the expected waiting time for these two sequences to coalesce is ~ ne generations. thus, the expected value of each πij is θ, whereas that for their sum is: rearranging equation leads to:. the tajima and watterson estimators have the same expected value of θ, even though the former is based on the total differences between each sampled sequence pair whereas the latter is obtained from the total number of segregating sites within the sample. their different summaries of the observed variation are the basis of tajima's d that tests for departures from the standard neutral model. to avoid random errors, equation can be adjusted in a manner similar to that used to modify equation of θw. first, let π'ij count the number of observed pairwise differences between sequences i and j after the removal of all singletons from their dataset. then, the expected number of singletons can be calculated from equations and as the total length of the external branches and basal branch multiplied by θ/ {i.e., θ}. the removal of each singleton from the original dataset reduces the expected sum in equation by. combining these results, we obtain: rearranging equation then leads to: our equation for θ't is equivalent to equation of achaz, except for our use of different symbols and total π'ij for the population sample. achaz also derived his equation b for the calculation of the associated variance for average π'ij ). new knudsen/miyamoto model knudsen and miyamoto developed a full coalescent/mutation model that accounts for three specific factors of experimental error and design: for random sequence errors, for unobserved polymorphisms due to missing data, and for the uncertain assignment of the multiple sequencing reads for a diploid or polyploid individual to its two or more homologues. their km model uses recursion to calculate an exact probability of the population sample under the standard fisher and wright model for reproduction and the infinite sites process for both mutations and random sequence errors. their model relies on ml to estimate θ and the expected number of errors per full length sequence. at the heart of the km model is equation of knudsen and miyamoto, which is reproduced here as: the parameters and factors of this equation are defined in table as one works backwards in time, the probability that the next observed event is a specific coalescence is given by the first term in the top line before the double sum. as indicated by this double sum, if indeed a coalescent event occurs, then it will happen between two sampled and/or ancestral alleles with compatible sequences. such sequences are referred to as combinable. m s the first term in the bottom line before the single sum is then for the probability that the next observed event is instead a mutation. as indicated by this sum, if indeed a mutation occurs, then it will happen to a sampled or ancestral sequence with at least one "singleton". here, the definition of a "singleton" is expanded to include the derived mutations of the common ancestors for the different groups of related sampled sequences as well as those for the singletons of the original dataset. this expanded use of the term acknowledges that a shared polymorphism under the infinite sites model is due to a unique mutation within the common ancestor of those sampled sequences sharing the derived base. thus, even though they result in shared polymorphisms among the sampled sequences, these derived mutations are ultimately counted as "singletons" as one works backwards in time. this expanded definition of a "singleton" allows for the economical use of σi alone to track the mutations of both the original shared polymorphisms and singletons. in the km model, the available region of each sampled sequence is summarized as a closed interval that is scored over the range of. the |si| factor then quantifies the total amount of sequence available for this sampled allele. for example, the score for the leftmost sampled sequence in figure 4a indicates that it is lacking the initial 20% of the full multiple alignment. thus, |si| = for this partial, leftmost, sampled sequence. in turn, as one works backwards in time, the closed intervals for the available regions of the sequences for common ancestors are calculated as the union of the known lengths for their two immediate descendants. thus, the closed interval of the common ancestor for the two leftmost sampled sequences in figure 4a is ∪ =, with |si| =. the purpose of Σi |si| in equation is to track the total available length of all sampled and/or ancestral sequences for the detection of mutations as one works backwards in time. the corresponding use of |si| in the bottom line of equation allows for the adjustment of σi due to unobserved polymorphisms resulting from missing data. the km model uses only equation when working with error free sequences. in turn, this model also relies on equation of knudsen and miyamoto when dealing with error prone sequences. this additional equation of the km model assumes that the random errors are uniformly distributed along the sampled sequences and that their total number is poisson distributed with an intensity of λ = ε Σi |si|. the inclusion of Σi |si| in the calculation of λ allows for the adjustment of this error rate to account for incomplete sampled sequences. in contrast to its predecessor, the km' model uses only equation when working with either error prone or error free sequences. as for the new watterson and tajima estimators, the km' model operates by counting only those internal branches of the genealogy where a mutation will result in a shared polymorphism. the sampled sequences are rescored as unknowns with empty intervals, as is the sequence of the common ancestor for the basal group of each asymmetrical genealogy. correspondingly, |si| is then reset to for these sequences. the km' model ignores the external and basal branches of the genealogy, where a mutation will result in a singleton, by multiplying σi for each sampled and/or ancestral sequence by |si| = in equation. in this way, the km' model discounts the singletons in favor of the shared polymorphisms within the dataset without the use of equation and its associated parameters of the original km model. for sampled sequences without missing data, the previous explanation is complete as to how the km' model corrects for the branches of the genealogy where a mutation will result in a singleton. however, for incomplete sampled sequences, the "above and below" test of the km' model also becomes necessary to correct for the regions of each common ancestral sequence where a mutation in its internal branch will lead to a singleton, rather than to a shared polymorphism, because of this missing information. in this test, "below" refers to those sampled sequences that belong to the monophyletic group of the common ancestor in question. "above" then corresponds to the remaining, more distantly related, sampled sequences. in figure 4b, the two leftmost sampled sequences are direct descendants of the first, leftmost, common ancestor, whereas the other three are not. thus, these two sampled sequences lie "below," whereas the other three occur "above" the internal branch for this leftmost common ancestor. the above and below test checks whether comparable information for a region is known for at least two, descendant, sampled sequences below and at least two, more distantly related, sampled sequences above an internal branch in the genealogy. if not, then a mutation within this region of the ancestral sequence that corresponds to this internal branch will result in a false singleton within the original dataset. for example, a mutation within the first 20% of the leftmost common ancestor in figure will result in a false singleton of the second sampled sequence rather than in a true shared polymorphism of the first and second. the problem is that the leftmost sampled sequence is missing comparable information for the detection of this mutation as a shared change in the leftmost common ancestor. the above and below test corrects for such regions of the common ancestral sequences during the rescoring of their closed intervals and |si|. as the first 20% of the leftmost common ancestor fails the below half of this test, the km' model disregards this region during the recalculation of its closed interval as and |si| =. evolutionary simulations and a. californica dataset evolutionary simulations to test the three new procedures, evolutionary simulations were conducted according to standard methods. two hundred datasets apiece were simulated for eight or sequences of length from a single population under the baseline conditions of the standard fw and infinite sites models with θ = or in addition to these baseline conditions, sequence errors were introduced as four or eight randomly placed changes among the eight and sequences, respectively, for an expected ε of. estimates of θ were then obtained for the datasets of each tested combination with the fw model and the original and new watterson estimators, tajima estimators, and km and km' models. as expected, the fw model and original watterson and tajima estimators overestimate θ when the datasets contain random errors. in these cases, the random errors are counted as mutations, thereby inflating their estimates of θ. furthermore, these overestimations are greater for the watterson estimator than for the tajima estimator. the reason is that each singleton makes the same contribution as a shared polymorphism to the number of segregating sites in equation of θw, but a smaller one to the sum of the pairwise differences in equation of θt. specifically, each singleton is limited in equation to the pairwise comparisons of the unique sequence to the other sequences. thus, the tajima estimator is less vulnerable than the watterson estimator to random errors even though its vulnerability is still significant. these results are for the standard fw model and the original and new watterson estimators, tajima estimators, and km and km' models. they are summarized as the means ± twice their standard deviations for simulated datasets apiece for each of the different combinations of θ, n, and ε. mean estimates that are significantly greater than or less than the true θ are highlighted in italics and boldface, respectively. no results are provided for the fw and km models with θ = or, since these calculations are too computationally intensive. in contrast, the km model underestimates θ when the sequences are errorless. furthermore, this tendency to underestimate θ is also evident when the sequences contain errors. in these situations, mutations are sometimes counted as random errors, thereby deflating their estimates of θ. in cases where few to no random errors are expected, one can first perform a likelihood ratio test to determine if ε =. if this null hypothesis cannot be rejected, then the km analysis should be restricted to only equation. however, if ε > according to this likelihood ratio test, then equation of knudsen and miyamoto is also needed and the user of the km model must remain aware that her/his estimate of θ most likely includes some slight downward bias. unlike their original versions, the new watterson estimator, tajima estimator, and km' model all consistently recover the true θ in the simulations both with and without errors. furthermore, these three methods also avoid the tendency of the km model to underestimate θ due to unnecessary or "greedy" parameters for random errors. in particular, the reliance of the km' model on one less equation and fewer parameters makes it much simpler and less prone to over-parameterizations than its predecessor. the mean estimates of θ for the new watterson estimator, tajima estimator, and km' model are also generally associated with greater standard deviations than those for their original versions. the proximal reason for these increased standard deviations is that the elimination of singletons results in the loss of some actual mutations along with the random errors. however, this cost of the three new methods appears to be relatively small, given that their standard deviations are never twice as great as those for their original versions. indeed, the standard deviations for the three new methods are less than or equal to those of their counterparts in four cases. perhaps more surprising is that the standard deviations for the new watterson and tajima estimators are comparable to those for the km' model, particularly when θ = or these similarities in their standard deviations are surprising given that the km' model is a full ml model. as summary statistics, the original and new watterson and tajima estimators are all computationally fast, requiring much less than cpu second on a ghz pentium cpu to analyze each of the simulated datasets. more importantly, the km' model is much faster than the fw and km models as documented by its completed analyses of the more complex datasets. in contrast, the fw and km models fail to complete their analyses of these more complex datasets due to time and memory constraints. the km' model is much faster than the fw and km models, because it relies on only shared polymorphic sites. less variable positions results in faster coalescences and fewer choices as one works back through the coalescent/mutation recursion of equation. aplysia dataset the a. californica dataset was the original motivating force behind the development of the km model. thus, this real dataset was also analyzed with the km' model to test further its performance against that of its predecessor. this dataset consists of sequencing reads for six diploid individuals from a laboratory population of the california seahare at the laboratory for marine bioscience, university of florida. three cloned inserts were sequenced from each individual as a pair of single sequencing passes starting from both ends of an internal segment of base pairs for the protein-coding region of the nuclear fmrf gene. these pairs of passes overlap in the middle for nine sequences, but at most by only bases. thus, these essentially single-pass sequences contain many random errors and some missing data and their assignments to the two homologues of each diploid remain uncertain. correspondingly, the km analysis of this dataset with singletons and shared polymorphisms required the full use of its factors for experimental error and design. a full ml analysis of this complex dataset by the km model proved too time and memory consuming and a two step procedure was therefore adopted instead whereby the number of errors was first ml estimated followed by the determination of θ for this ml value. this heuristic approach resulted in an estimate of θkm = for this sample, which still took more than hours to complete on the same cpu as for the simulations. in contrast, a full ml analysis of this dataset by the km' model using the same factors for experimental design is much faster as it took less than hour on this cpu to obtain a similar estimate of θ'km =. correspondingly, nucleotide diversity for this sample is calculated as = mutations per site. to summarize, the km and km' models both support similar estimates of θ for this real dataset, but the latter is much faster as confirmed by its completed, full, ml analysis.  expressed sequence tags, low coverage draft sequences, and other such unfinished products are known to contain random errors due to chance accidents during their data collection and recording. however, such sequences often constitute the only available allelic information for a population and methods to deal with their random mistakes are therefore needed to obtain accurate estimates of θ from these error prone data. this study is based on the simple premise that random sequence errors are distributed as singletons. thus, one can avoid the random mistakes of error prone sequences by ignoring their singletons in favor of their shared polymorphisms. this strategy is implemented in the new watterson estimator, tajima estimator, and km' model. these new methods are all accurate and fast according to their evolutionary simulations and analysis of the real complex dataset for a. californica. these methods come with the cost of increased standard deviations, but this price appears small or even negligible compared to their advantages of significantly improved accuracy and/or computational speed. obviously, additional evolutionary simulations and applications to real datasets are now needed to evaluate more fully under what conditions the removal of singletons is warranted in light of this tradeoff. nevertheless, the current successes with the new watterson estimator, tajima estimator, and km' model support our recommendation that these three methods be given serious consideration when estimating θ from error prone sequences. unlike the watterson and tajima estimators that represent summary statistics, the km and km' models both constitute full ml models that offer the framework for the further incorporation of other experimental and population genetic factors. in particular, such further developments are encouraged for the km' model in light of its current successes in the evolutionary simulations and a. californica analysis. for example, biased errors within a sample due to the systematic misreading of specific bases during dna sequencing, the postmortem biochemical degradation of ancient dna, and/or other such sources of error can be accommodated by a finite sites process that allows for repeated mistakes as well as mutations at the same sequence positions. likewise, additional factors to account for other population genetic processes such as recombination can be accommodated by the use of ancestral recombination graphs. inevitably, these more complex versions of the km' model will require the use of sampling based procedures for their implementation, since the current use of direct ml evaluation will remain practical for only the smaller datasets and simpler models. authors' contributions bk derived the equations, wrote the computer program for the estimation of θ, and conducted the evolutionary simulations. mmm provided biological interpretations about the results. both authors contributed to the design and main ideas of this study, wrote and edited the manuscript, and read and approved its final draft. 
it is then used to modify the recent, full, maximum-likelihood model of knudsen and miyamoto, which incorporates various factors for experimental error and design with those for coalescence and mutation. in light of these results, we recommend the use of these three new methods for the determination of θ from error prone sequences. however, its accurate estimation can be seriously compromised when working with error prone data such as expressed sequence tags, low coverage draft sequences, and other such unfinished products. these new methods are all accurate and fast according to evolutionary simulations and analyses of a real complex population dataset for the california seahare. this strategy is implemented under an infinite sites model that focuses on only the internal branches of the sample genealogy where a shared polymorphism can arise. this approach is first used to derive independently the same new watterson and tajima estimators of θ, as recently reported by achaz for error prone sequences. the population mutation rate remains one of the most fundamental parameters in genetics, ecology, and evolutionary biology. thus, one can avoid these random errors by ignoring the singletons within a dataset. this study is premised on the simple idea that a random sequence error due to a chance accident during data collection or recording will be distributed within a population dataset as a singleton. in particular, we advocate the new maximum likelihood model as a starting point for the further development of more complex coalescent/mutation models that also account for experimental error and design.
 the population mutation rate remains one of the most fundamental parameters in genetics, ecology, and evolutionary biology. however, its accurate estimation can be seriously compromised when working with error prone data such as expressed sequence tags, low coverage draft sequences, and other such unfinished products. this study is premised on the simple idea that a random sequence error due to a chance accident during data collection or recording will be distributed within a population dataset as a singleton. thus, one can avoid these random errors by ignoring the singletons within a dataset.  this strategy is implemented under an infinite sites model that focuses on only the internal branches of the sample genealogy where a shared polymorphism can arise. this approach is first used to derive independently the same new watterson and tajima estimators of θ, as recently reported by achaz for error prone sequences. it is then used to modify the recent, full, maximum-likelihood model of knudsen and miyamoto, which incorporates various factors for experimental error and design with those for coalescence and mutation. these new methods are all accurate and fast according to evolutionary simulations and analyses of a real complex population dataset for the california seahare.  in light of these results, we recommend the use of these three new methods for the determination of θ from error prone sequences. in particular, we advocate the new maximum likelihood model as a starting point for the further development of more complex coalescent/mutation models that also account for experimental error and design. 