BACKGROUND
next generation sequencing technologies have unprecedented potential for assessing genetic diversity in viral populations  <cit> .
as compared to traditional sanger bulk sequencing, the cost per base of sequencing is low, due in part to the massively parallel nature of these technologies  <cit> .
roche- <dig> is currently the most popular technology to investigate diversity in viral populations, thanks to its superior read length compared to other platforms such as illumina gaiix  <cit> , but both technologies have advantages and disadvantages for analysing viral diversity, suggesting that experimental design should be planned carefully  <cit> .
currently, one run of roche- <dig> sequencing can produce one million reads, with a modal read length of 700nt - enough to sequence five 10kb viral genomes to a coverage of  <dig> x each.
this high coverage facilitates deep sequencing experiments, where reads are treated as a population sample containing information about the underlying genetic structure of the population.
for example, deep sequencing has been used to study bottlenecks in early hcv infection  <cit>  and to identify mutations associated with drug resistance in hiv  <cit> .
a major challenge in deep sequencing experiments involves separating true variants from sequencing errors.
although average error rates are constantly improving, error rates remain highly heterogeneous as they are affected by sequence context  and alignment errors  <cit> .
this means that locally, error rates can still be quite high, leading to false positive snv calls  <cit> .
a number of tools are available for calling snvs within the context of sequencing errors.
however, most approaches are tailored to calling snvs in human resequencing projects and are not suitable for viral deep sequencing data.
when resequencing a human genome, snvs are either heterologous and can be expected to be present in 50% of all reads, or homologous and should be present in 100% of reads.
commonly, a probabilistic approach is employed to decide whether a potential snv is either an error, a homologous snv, or a heterologous snv  <cit> .
in deep sequencing experiments, however, very low frequency snvs, including snvs present in less than 1% of reads, are often of interest  <cit> .
for instance, the presence of low frequency variants in the viral population of a single host may determine the potential for drug resistance and treatment failure  <cit> .
to date, deep sequencing data has been analysed by applying various read quality filters before and after alignment and then further manually filtering the results, according to the researcher’s understanding of the sequencing technology and the individual data set  <cit> .
filtering may also be performed with a software program such as varscan  <cit> , an snv caller which relies on read coverage, base qualities, and variant frequencies to separate true snvs from errors.
an emerging statistical method for detecting snvs in deep sequencing data involves performing a statistical test for each potential snv, taking into account a prior positional error model  <cit> .
this method detected snvs with a population frequency of  <dig> % in control illumina data.
however, the requirement for a prior positional error model relies on the availability of multiple replicates, which is not always feasible.
a more recent approach to snv detection in viral populations, implemented in the program v-phaser, utilises the covariance between observed snvs present on the same read to increase sensitivity  <cit> .
the power of covariance was also harnessed in a recent study of hcv early evolution, where snv calling was achieved by parsing haplotypes reconstructed using probabilistic clustering performed by the program shorah  <cit> .
shorah uses a bayesian approach to group reads into haplotypes.
both sensitivity and specificity are improved by assuming that the consensus sequence for each cluster is the true haplotype sequence; all other variants present in the reads within a cluster are considered to be errors  <cit> .
shorah has been shown to accurately reconstruct variants with a population frequency of  <dig> %, and to efficiently remove random pcr errors  <cit> .
error correction via probabilistic clustering will fail, however, when errors themselves appear clustered.
this can occur in regions of low diversity characterised by an absence of true variants, or when sequencing errors are systematic .
the potential for systematic errors to confound clustering may explain why an empirical k-mer based algorithmic approach was found to be more efficient than shorah at removing false haplotypes, although both approaches detected all true haplotypes  <cit> .
there is growing evidence that systematic errors may exhibit a strand bias, where a particular error is more likely to occur in reads traversing the genome in one direction as opposed to the other  <cit> .
for instance, in a recent analysis addressing errors from control illumina data using prior positional error models, error rates between forward and reverse reads were uncorrelated  <cit> .
because of this, several snv detection methods consider forward reads and reverse reads separately, using some combination of the two obtained p-values to inform the final snv calling decision  <cit> .
researchers have also used strand bias directly to filter out false positive snvs.
for example, in a deep sequencing analysis of hiv- <dig> quasispecies, any potential snvs with more than  <dig> fold difference in the forward and reverse read count were considered to be systematic errors  <cit> , while in the hcv study mentioned above, raw results from shorah were further cleaned by removing any potential snvs that showed a difference of more than 10% between the number of forward and reverse reads  <cit> .
the latest version of samtools formalises this approach, reporting a p-value based on an exact test of strand bias  <cit> .
a two tailed fisher’s exact test is also included in the recently released snv caller lofreq  <cit> .
here, we show that in addition to its primary function of haplotype reconstruction, shorah may also be used to reliably call snvs.
this is demonstrated through analysis of control roche- <dig> data from both hcv and hiv samples, and by development of a statistical argument in support of snv calling via probabilistic clustering.
as strand bias has been observed in systematic errors retained during clustering  <cit> , we also investigate the benefits of combining probabilistic clustering with a statistical test of strand bias.
RESULTS
our snv calling method is based on two assumptions:  it is more powerful to compare several variants that co-occur on the same read simultaneously than individually; and  the frequency of sequencing errors, but not of true snvs, is often higher on reads sequenced in one direction than the other.
while there is empirical evidence for the second assumption, the first one can be motivated and made precise by a simplified statistical model of snv detection .
in this model, two haplotypes of frequencies f and  <dig> − f that differ at d sites co-occurring within an observed read, can be distinguished in the presence of a per site error probability ϵ only if  f > ϵd.
for d =  <dig>  this result confirms the intuitive notion that a single snv can be detected only if the error rate is smaller than the snv frequency.
however, if the haplotypes differ at two sites , then the probability that these two variants are both errors is ϵ <dig>  and the limit of haplotype detection decreases to f > ϵ <dig>  in general, considering multiple variant sites at once by comparing haplotypes increases snv calling accuracy .
this finding motivates analysing entire reads or large segments of reads in order to include as many snvs as possible, as is done in read clustering methods such as shorah.
to investigate real life performance of snv calling via probabilistic clustering, we analysed roche- <dig> data from two distinct control experiments.
for the first experiment, the same library of cloned moderate-diversity hcv fragments mixed at equal quantities was subjected to two independent roche- <dig> sequencing runs   <cit> .
for the second experiment, one round of roche- <dig> sequencing was performed on a high diversity pcr amplified set of hiv fragments mixed according to a geometric series  <cit> .
snvs were called using either probabilistic clustering , quality score based filtering , or a per site quality score based statistical test   <cit> .
a statistical test of strand bias was applied to the results from shorah and varscan, utilising either a beta-binomial  model of forward read distribution for true snvs, or fisher’s exact test .
our experimental design allowed us to assess the performance of probabilistic clustering, and variation in the true snv strand bias distribution, both between samples with different diversity levels and between individual sequencing runs.
alignment and clone frequencies
for hcv run  <dig> and hcv run  <dig>   <dig> % and  <dig> % of reads aligned uniquely to the reference, giving an average read depth of  <dig> reads and  <dig> reads, respectively.
for the hiv data,  <dig> % of reads aligned, giving an average read depth of  <dig> reads.
true haplotype frequency estimates are given in table  <dig>  and table  <dig> .
for all samples, the frequencies of at least some haplotypes differed substantially from their intended frequencies.
as the estimated frequencies for hcv run  <dig> and hcv run  <dig> are very similar, it is likely that deviation from intended frequencies is a result of mixing errors, rather than bias introduced during library preparation and sequencing.
lists of true snvs for both hcv and hiv samples are provided in the additional files  <dig> and  <dig>
raw snv calls
for hcv run  <dig> and hcv run  <dig>  both shorah and varscan detected all true snvs, i.e.
the false negative rate was zero , while lofreq called all but one true snv in each run.
varscan, however, also called up to  <dig>  and  <dig>  false positives, more than seven times as many as shorah.
more errors were incorrectly identified as snvs for hcv run  <dig> compared to hcv run  <dig> by both varscan and shorah, indicating a higher error rate for hcv run  <dig>  varscan’s recall was greater than shorah’s for the hiv data, with 93% of true snvs identified, compared to 81% for shorah.
however, shorah’s recall was greater than lofreq’s, which only identified 64% of true snvs.
for the hiv sample, all of shorah’s false negatives had true population frequencies under  <dig> %, while one of varscan’s false negatives had a true population frequency of  <dig> % .
upon inspection of the alignment and the varscan file, this appeared to be the result of confusion in reporting positions with more than one variant allele.
most of lofreq’s false negatives were very low frequency , with three having true population frequencies of under 6%.
varscan’s superior recall was accompanied by a high number of false positive calls – over the reference 766bp considered,  <dig> false positive snvs were also called.
by comparison, shorah only made two false positive calls .
overall, lofreq exhibited an exceptionally low false positive rate.
for the hcv runs, one unique false positive was identified in each run, while no false positives were called by lofreq for the hiv data.
snv calling statistics for shorah, varscan, and lofreq.
for shorah and varscan, snv calls without the strand bias test  and using different values of the beta-binomial dispersion parameter σ in the strand bias test are given, with σ =  <dig> corresponding to a binomial forward read distribution.
for lofreq, the results of applying our strand bias test are absent as this software does not report forward and reverse strand counts.
for all snv calling methods, the results of applying fisher’s exact test to the raw output are also given .
reported statistics include true positives , i.e., genomic sites with a variant matching a known true variant, false positives , i.e., genomic sites with a variant that is not a known true variant, and false negatives , i.e., known variants which are not identified by the relevant snv calling method.
individual genomic sites may contribute to both true positives and false positives.
recall ) and precision ) are also reported.
strand bias test
in an attempt to improve shorah’s performance under conditions of sub-optimal diversity or clustered systematic errors, we applied a statistical test of strand bias to the raw snv calls.
for each potential snv detected, we test the hypothesis that the distribution of forward reads carrying this snv follows a beta-binomial distribution, with mean value equal to the overall forward read distribution at that position.
the overdispersion of the forward read distribution is controlled by the parameter σ, which we estimate separately for each sequencing run using maximum likelihood  methods based on the observed forward read distribution for all true snvs within the sample .
ml estimates of σ for hcv run  <dig>  hcv run  <dig>  and hiv samples were  <dig> ,  <dig> , and  <dig>  respectively.
to avoid overfitting, we applied tests using all estimated values of σ to each dataset, also performing a test with σ= <dig> .
for comparison with the strand bias test implemented in other software, we also performed a fisher’s exact test of strand bias.
snv calling statistics for strand bias tests applied to the raw results of both shorah and varscan are given in table  <dig>  the results of lofreq’s own implementation of fisher’s exact test are also given.
applying a strand bias test to the raw snv calls for both shorah and varscan improved precision in all cases, except for hiv snv calls made using shorah with σ =  <dig> , where precision remained unaltered compared to the raw shorah output.
in general, choice of σ did not have a substantial impact on precision.
however, smaller values of σ were more likely to result in reduced recall, as was fisher’s exact test.
recall was also compromised by lofreq’s implementation of fisher’s exact test.
this was particularly evident for the hiv data, and is consistent with an overdispersed forward read distribution for true snvs, justifying our choice of a beta-binomial forward read distribution.
for the beta-binomial test, the reduction in recall was minimised in all cases by using σ =  <dig> .
to further investigate the appropriateness of a beta-binomial based strand bias test , compared to either a binomial based test or fisher’s exact test, we generated precision recall curves for each test applied to snv calls from shorah and varscan .
in all cases, across the full range of recall values obtained the beta-binomial test gave equivalent or higher precision than either of the other tests, confirming its suitability.
for the hcv samples, the strand bias test was particularly good at reducing errors in results from shorah, with more than 70% of errors eliminated.
after the strand bias test, less than six errors for every  <dig> reference sites remained for hcv run  <dig>  while less than four errors for every  <dig> reference sites remained for hcv run  <dig>  the strand bias test also eliminated 48% of errors in varscan’s output, however as varscan’s raw output contained more errors than shorah, this still left a minimum of  <dig> errors for every  <dig> reference sites.
lofreq’s strand bias test successfully eliminated the single false positives originally identified in hcv run  <dig> and hcv run  <dig>
figure  <dig> presents a detailed analysis of snv calling by frequency for varscan and shorah, for snvs with frequencies less than 15%.
lofreq results are excluded from this analysis, due to its near perfect precision for all data sets.
figure 2a shows that for the high diversity hiv data, perfect recall could be achieved for snvs with population frequency greater than around  <dig> %.
any reduction in recall due to the strand bias test was limited to snvs with frequencies below this detection limit.
thus, although consideration of both recall and precision indicates the strand bias test does not improve on the raw hiv shorah results, it is not significantly detrimental.
the lack of improvement is due in part to the extremely high precision of  <dig>  obtained by probabilistic clustering alone .
for the hiv varscan results, although the strand bias test did increase precision, more than  <dig> errors for every  <dig> reference sites remained.
figure 2b indicates that this is the result of low frequency false positives.
restricting snv calls to those covered by a minimum of  <dig> reads therefore rectified this problem to some extent , however over  <dig> errors for every  <dig> reference sites were retained.
recall was also adversely affected, due to the presence of very low frequency true snvs in the hiv sample.
using a minimum read coverage of  <dig> had no effect on hiv shorah results.
figures 2b to 2d describe the number of false positives observed for each data set, binned in frequency increments of 1%.
in general, shorah calls less false positives than varscan, particularly for very low frequency snvs .
in most cases, the strand bias test was able to reduce or eliminate false positives across the range of frequencies observed .
the effectiveness of the strand bias test was most notable for shorah snv calls from the moderate diversity hcv samples.
discussion
we have shown that snv calling via probabilistic clustering is a powerful technique allowing detection of snvs at frequencies lower than the corresponding sequencing error rate, as long as sample diversity is high  and errors are random.
when these conditions are not met, a statistical test of strand bias improves snv calling precision.
the distribution of forward reads covering a true snv at any one position could be expected to follow a binomial distribution, with mean value equal to the fraction of forward reads at that position multiplied by the total number of reads covering the variant.
however, our results show that more accurate snv calling is achieved if the expected number of forward reads is modelled using a beta-binomial read distribution, compared to either a binomial distribution or applying fisher’s exact test.
variation in the number of false positives between the hcv runs indicated error rates differed markedly, even though all samples were sequenced using roche- <dig> titanium technology.
notably, the ability to utilise empirical run specific error rates was one of the motivating factors behind lofreq’s development.
a similar result has been reported for illumina sequencing as well, with individual runs sequenced on the same machine having dissimilar error profiles  <cit> .
snv calling precision after strand bias testing was relatively robust to choice of dispersion factor σ. using a very low value of σ, however, compromised snv recall, thus we would recommend choosing a value of σ closer to our upper estimate of  <dig>  when applying our test to real data.
for this reason, σ is set to  <dig>  by default in our updated version of shorah .
for the hcv samples, which featured a moderate number of relatively high frequency true snvs, the beta-binomial strand bias test provided a substantial reduction in error rates with no cost to true positive recall.
this was true when applied to both the raw shorah output and the raw varscan output.
however, both the relative reduction in error and absolute precision attained was greatest when the test was applied to the raw shorah output.
unfortunately, it was not possible to apply our beta-binomial strand bias test to lofreq’s raw snv calls, as variant counts by forward and reverse read are not given in lofreq’s output.
lofreq’s own fisher’s exact test of strand bias successfully removed all false positives, however it also removed a number of true positives.
given the results reported here, we believe that a beta-binomial strand bias test may also remove these false positives, with substantially less cost to true positive recall.
for the hiv sample, the raw shorah output gave the best results.
this can be explained by the extremely strong performance of error correction via probabilistic clustering in this case.
the ability to correctly identify true haplotypes via probabilistic clustering  is a function of the frequency of the haplotypes, the number of variant positions within the haplotypes, and the error rate of the sequencing technology .
whenever there is more than one true variant within a haplotype under consideration, it becomes possible to detect snvs with population frequencies below the error rate of the sequencing technology.
consideration of hiv false positives called using varscan indicates that locally, error rates may exceed 5%.
for the hiv data, we achieved a detection limit well below this local error rate using shorah alone, validating our argument.
figures 2a and 2b show that perfect accuracy  can be achieved by applying a detection limit of 3%, while a detection limit of only  <dig> % gives perfect recall with only two false positives.
another implication of our model is that accuracy should improve as the number of snvs within a window of sites considered simultaneously increases, as demonstrated by the fact that shorah performs better on the intrinsically high diversity hiv sample than on the moderate diversity hcv samples.
our argument is further supported by varscan’s performance.
although for the hiv sample varscan detected more true positives than shorah, the corresponding extremely large number of low frequency false positives within the varscan results suggests that this is the result of only considering one site at a time and indiscriminately calling low frequency snvs below the associated detection limit .
lofreq’s approach of discriminating snvs individually by inferring an error probability based on the quality score assigned to individual nucleotides clearly overcomes this problem, although at the expense of recall.
while our strand bias test could not improve shorah’s raw snv calls for the hiv sample, we consider the slight reduction in recall with σ =  <dig>   to be an acceptable cost, compared to the potential benefits of the test.
it is usually not possible to know, a priori, how diverse a deep sequenced population may be.
even within highly diverse populations, genomic regions under purifying selection may have sub-optimal levels of diversity, allowing errors to appear clustered as true snvs.
our strand bias test can help remove such errors, while leaving snv calls in highly diverse regions relatively untouched.
our strand bias test could also be used to provide a level of confidence in results.
if snv calls are unaffected by applying a strand bias test, as in the case of our hiv sample, then it is likely that diversity was sufficient for probabilistic clustering to perform at its best.
applying our strand bias test to raw varscan snv calls also improved precision in all cases.
however, false positive snv calls could not be reduced to a biologically useful level.
the relative advantage of applying a strand bias test to shorah’s output compared to varscan’s output can be understood by considering how each of these programs identifies snvs.
varscan distinguishes errors from true snvs by examining the individual base quality scores of reads covering the potential snv.
varscan is therefore reliant on accurate assignment of quality scores during the sequencing platform’s base calling pipeline.
if an error nucleotide is falsely assigned a high quality score, then varscan will retain it.
although lofreq is also dependent on quality scores assigned by the sequencing platform, it has some resilience to this effect as quality information from all aligned nucleotides at a position is considered when making a call; errors at a given position would have to systematically be assigned high quality scores for lofreq to incorrectly call an snv.
in contrast to varscan, shorah is independent of quality score information.
shorah works by removing any potential snvs that are randomly distributed among reads, rather than clustering with nearby potential snvs - regardless of whether they occurred during library preparation, pcr amplification, or read synthesis.
errors retained during clustering are likely to be systematic sequencing errors, which tend to exhibit a strand bias.
methods
statistical analysis of snv calling based on probabilistic clustering
snv calling from shorah output involves parsing a set of ‘haplotypes’ for variants with respect to a reference sequence.
these haplotypes are generated by deriving consensus sequences from groups of probabilistically clustered aligned reads.
because only the consensus sequence within each cluster is considered, snvs will be correctly called, provided reads are correctly clustered so that the haplotypes can be identified from the centroids of the clusters.
the detection limit for true snvs therefore depends on the detection limit for haplotypes.
the detection limit for haplotypes is the haplotype frequency at which observed variants are equally likely to be errors  or biological variants .
for example, if two haplotypes differ at only one site , then the detection limit is intuitively equal to the error rate, ϵ. if, however, the haplotypes differ at two sites , then it takes two specific sequencing errors to mutate a read originating from one haplotype into a read that is equal to the other haplotype.
this occurs with a probability of the order ϵ <dig>  now, ϵ2< ϵ, allowing haplotypes  to be detected at frequencies lower than the error rate provided two sites are considered simultaneously .
to investigate this property more formally, we study the probability of correctly assigning a read to either of two haplotypes: a minor haplotype with frequency f containing the same variants as the read under consideration, and a major haplotype  with frequency  <dig> − f .
to simplify our analysis, we consider only those sites where these two haplotypes differ, as other sites are uninformative when choosing which haplotype to assign the read to.
according to this model, the major and minor haplotypes may be represented by the binary sequences 0 =  <dig> .
.
.
<dig> and 1 =  <dig> .
.
.
<dig> of length d, where d is the hamming distance between these two haplotypes.
now, let the random variable h =  denote the true haplotype, such that p =  <dig> − f and p = f .
observation errors are assumed to occur independently with uniform per site probability p = ϵ, where each random variable ei indicates the occurrence of an error at site i. the observed read, y, is the true haplotype h subject to observation error: yi =  <dig> if  =  or , and yi =  <dig> otherwise.
the full model is defined as follows:
h∼bernoullif,h∈0¯,1¯ei∼bernoulliϵ,ei∈ <dig> ,i= <dig> ...,dyi∼hi+eimod <dig> i= <dig> ...,d.
in particular, p  = ϵd and p  = d.
in order to assign the read correctly, we need to decide whether the observed read, y =  <dig>  is more likely to originate from the minor haplotype h =  <dig> without error or from the major haplotype h =  <dig> through  errors.
this decision represents a crucial step in the iterative probabilistic assignment of reads to haplotypes within shorah – if the constituent reads of a low frequency true haplotype are just as likely to originate from a higher frequency haplotype through sequencing errors, then they may be incorrectly assigned to the higher frequency haplotype.
in this case, the low frequency haplotype falls below the limit of haplotype detection, and may not be represented in the final set of clustered reads.
based on this model, the assignment will be correct if:
ph=0¯|y=1¯<ph=1¯|y=1¯.
using bayes’ theorem and eqs.
<dig>  this inequality becomes ϵd < df, or
ϵ<1+1−ff1/d− <dig>
for noise levels ϵ below this upper bound, we would assign the reads correctly.
for d =  <dig>  the condition is ϵ < f, and for ϵ and f small, it is approximately ϵd < f. thus for a given per nucleotide sequencing error rate, the frequency limit for haplotype detection decreases as the number of true snvs observed within a read increases.
in reality, the process of haplotype inference is of course more complicated than presented here – our simplified argument is designed to merely demonstrate in a quantitative way that by observing the co-occurrence of variant sites, more power can be achieved when calling low frequency haplotypes  than by considering variant sites independently.
strand bias test
we assume that at any given genomic position, the number of forward reads carrying a variant follows a beta-binomial distribution.
more formally, let fi and ri, respectively, denote the total number of forward and reverse reads covering position i. the fraction of forward reads is μi = fi/.
we define nib as the total number of reads  covering the variant base b at sequence position i. the number of forward reads carrying variant b at position i is modelled by the beta-binomial random variable xib ~ betabin such that  <dig> ≤ xib ≤ nib and
pxib=x=Γnib+1Γx+1Γnib−x+1×Γ1/σΓx+μi/σΓnib+1−μi/σ−xΓnib+1/σΓμi/σΓ1−μi/σ
the expected value of xib is e = nibμi, and σ is a dispersion parameter, defined such that
varxib=nibμi1−μi1+σnib−11+σ.
the dispersion parameter σ determines the spread of the distribution.
we assume that, unlike nib and μi, σ does not change between each variant within a sample.
both nib and μi can be easily estimated directly from the data for each individual potential snv under question.
we employed ml estimation to obtain σ estimates for each sample, by considering the k true snvs xi1b <dig> …,xikbk.
the likelihood function is:
lσ=∏j=1kpxijbj=xijbj.
ml estimation was performed in r with the package mle <dig>  using the default nelder and mead optimisation method and a start value of one.
.
our null hypothesis, h <dig>  is that for any given true snv, the number of forward reads follows a beta-binomial distribution with parameters μ, n, and σ. the null hypothesis h <dig> is rejected if, under h <dig>  the probability to observe data as extreme or more extreme than that observed is smaller than a threshold p. that is, given an observation xib, we reject h <dig> if min {Φ,  <dig> − Φ } × 2 < p, where Φ = p is the beta-binomial cumulative distribution function , with parameters estimated as described above.
if, for a potential snv, the null hypothesis is rejected, then the potential snv is considered to be an error.
as we do not care which direction the strand bias occurs in, a two-sided test was performed by multiplying the smallest tail of Φ by two.
we then applied the benjamini-hochberg correction for multiple testing, giving the final p-value p. if p <  <dig> , we reject the potential variant as an error.
to investigate the effect of σ on the strand bias test , strand bias tests with each estimate of σ were applied to snv calling results obtained from shorah and varscan for all data sets.
a test with σ =  <dig> was also performed, which models a binomial forward read distribution.
for comparison, two-sided fisher’s exact tests were also applied, based on counts of the forward and reverse reads, again with a benjamini-hochberg correction for multiple testing.
it was not possible to apply our strand bias test to snv calls from lofreq, as lofreq’s output does not record the number of forward and reverse reads covering a variant.
lofreq does however report a p-value based on a fisher’s exact test of strand bias.
thus we were able to report the results of a fisher’s exact test of strand bias applied to raw lofreq snv calls, again applying a benjamini-hochberg correction for multiple testing within r and rejecting any variants for which p <  <dig> .
to further compare the performance of fisher’s exact test, a binomial strand bias test, and our beta-binomial test , precision recall curves were constructed within r for each test applied to shorah or varscan snv calls from hcv run  <dig> and hcv run  <dig>  this was done by calculating the precision and recall for p-value cut offs varying between zero and one  and plotting the results.
sample sequencing, alignment, and frequencies
for hcv e1/e <dig> samples, a mixture of four different 2324bp fragments , each covering the 3′ end of the hcv core to the 5′ end of the p <dig> region was cloned into the pgem-t easy vector and sequenced using roche- <dig> flx titanium as described previously   <cit> .
to investigate variation in σ between sequencing runs, a second round of roche- <dig> flx titanium sequencing was performed on this same library of hcv clone mixtures .
for the hiv gag/pol-gene fragments,  <dig> pcr products from subtype-b clinical isolates were propagated within the pcrii-topo vector, excised, and subjected to a single round pcr before roche- <dig> flx titanium sequencing, as previously described  <cit> .
prior to alignment, reads were trimmed using a sliding window of 8nt.
the trim was made at the start of the first 8nt window from either end with an average quality greater than  <dig>  for all samples, alignment was performed using bwasw  <cit> , with default parameters.
for the hcv samples, alignment was performed against the entire pgem genome with an hcv 1a strain insert, corresponding to the inserted fragments  <cit> .
for the hiv sample, alignment was performed against the appropriate region of the hxb <dig> hiv reference genome, from base  <dig> to  <dig>
to estimate the true frequency of each fragment within the mixture, the known sequences of each fragment were used to locate snv positions uniquely identifying each fragment.
fragment frequencies were then calculated by averaging the frequency of all snvs unique to a fragment haplotype, then normalising the results so that the sum of all haplotype frequencies was one.
snv calling
shorah  <dig>  was run on all three alignments described above.
for each run, α was set to  <dig> , window size was set to  <dig> with a step size of  <dig>  and the number of iterations to  <dig> times the window coverage for each window.
to call an snv from the raw shorah output, locally reconstructed windows were parsed for variant positions.
for an snv to be called, the same variant was required to be found in at least two of three overlapping windows, with a posterior probability greater than  <dig> .
this approach was compared to single site snv calling, as implemented in lofreq and varscan.
lofreq is a single site variant caller that uses per nucleotide quality score information to inform a statistical test for each potential snv.
lofreq is tailored to calling low frequency variants in deep sequencing data, and is thus ideally suited to identifying variants in viral populations.
snv calling was performed using the script lofreq_snpcaller.py  with default options.
varscan was also employed due to its ability to report the required read statistics for the strand bias test, and its single site mode of snv calling.
however, we acknowledge that varscan is intended for use on diploid resequencing data rather than deep sequencing of viral populations, and do not imply that results reported here are indicative of performance under normal usage conditions.
varscan’s pileup2snp command  was used to call snvs.
default values for quality filtering were used.
minimum variant frequency was set to zero, to allow comparison with shorah.
for all variant callers employed, we limited called variants to those covered by more than one read.
to allow a fair comparison, snv calling using lofreq and varscan was only attempted in areas of the genome covered by three shorah windows.
this corresponds to positions  <dig> to  <dig> of the hcv fragments , and positions  <dig> to  <dig> of the hiv fragments .
in total,  <dig> true snvs were analysed for the hcv samples, and  <dig> true snvs were analysed for the hiv sample.
an attempt to call snvs using v-phaser was made, as this program combines joint snv calling with quality score based methods.
however, the program crashed when used on any of our samples, possibly due to the high coverage of samples.
correspondence with the authors did not resolve the issue, thus the comparison was not able to be made.
snv calling was also attempted using an empirical k-mer based algorithm   <cit> .
after running for more than  <dig> cpu hours, this software did not produce a result, at which point this comparison was also abandoned.
snv calling with deepsnv  <cit>  and using the approach described by flaherty et al.
<cit>  was not attempted due to the lack of control or multiple samples.
for shorah and varscan, raw results, as well as results subjected to our strand bias test, are reported.
for the purposes of this paper, all strand bias tests were performed in r as described above.
this was done to facilitate an exact comparison between varscan and shorah results.
however, in conjunction with this manuscript we have released a new version  <dig>  of shorah with integrated snv calling and strand bias testing, which does not require r.
CONCLUSIONS
we have demonstrated that by combining probabilistic clustering with a statistical test of strand bias, highly accurate snv calling can be achieved.
this is because the strengths and weaknesses of these two methods are complementary: probabilistic clustering relies on errors occurring randomly and independently, while strand bias tests detect systematic errors with atypical read direction distributions.
using a statistical model and through analysis of control data we investigated the conditions under which probabilistic clustering can be used to call snvs.
we show that the limit of snv detection is ϵd < f, where ϵ is the rate of randomly occurring sequencing errors, d is the number of true variants observed within a read length, and f is the frequency of the haplotype to which the variants belong.
in general, the snv detection limit decreases as true diversity increases.
this is consistent with our empirical analysis of snv calling from control data with various diversity levels.
importantly, for snv calling via probabilistic clustering to have an advantage over approaches that consider individual variants one at a time, more than one true variant must be contained within a read length.
to insure against situations where diversity is inadequate or errors are non-random, we suggest a strand bias test be applied to snv calls generated via probabilistic clustering.
the test we propose uses a beta-binomial model of forward read distribution for true snvs.
by comparing our test to the performance of either fisher’s exact test or a test based on a binomial model of forward read distribution, we show that the forward read ratio of true snvs is overdispersed.
our analysis indicates that as long as this overdispersion is taken into account, a statistical test of strand bias will not affect recall.
this means such a test can be safely applied and will not adversely affect accuracy, even in situations where no improvement in precision is obtained.
while our analysis is presented in the context of snv calling using shorah, the insights into snv calling via probabilistic clustering should be applicable to any method that considers multiple variant sites simultaneously.
furthermore, our strand bias test has the potential to improve snv calling precision for any variant detection pipeline.
availability
the snv caller and strand bias test have been integrated into shorah  <dig> , available at http://www.cbgethz.ch/software/shorah.
haplotypes, aligned reads, and reference sequences used in this study may be obtained from https://wiki-bsse.ethz.ch/display/shorah/data.
contact: shorah@bsse.ethz.ch.
competing interests
the authors declare that they have no competing interest.
authors’ contributions
km implemented the strand bias test, analysed control data, and in conjunction with oz developed the software shorah v <dig> .
nb developed the snv calling model.
km wrote the manuscript, with editing by fl, oz, and nb.
nb and fl conceived of the overall study.
rb and oz provided control data.
all authors read and approved the final manuscript.
supplementary material
additional file 1
tab delimited list of true snvs for hcv haplotypes.
format: position \t reference nucleotide \t variant nucleotide \t frequency run  <dig> \t frequency run  <dig>  to facilitate analysis with supplied reference sequences and alignments, position is given with respect to the start of the vector; subtract  <dig> to get position with respect to the hcv fragment.
frequencies are calculated based on the estimated true frequency of haplotypes .
click here for file
additional file 2
tab delimited list of true snvs for hiv haplotypes.
format: position \t reference nucleotide \t variant nucleotide \t frequency.
frequencies are calculated based on the estimated true frequency of haplotypes .
click here for file
acknowledgements
km was supported by an australian postgraduate award.
fl and rb are both nhmrc training fellows.
we acknowledge support from ach <dig> and nhmrc.
part of this work was supported by the swiss national science foundation under grant number cr32i <dig>  <dig> .

 next generation sequencing technologies have unprecedented potential for assessing genetic diversity in viral populations. as compared to traditional sanger bulk sequencing, the cost per base of sequencing is low, due in part to the massively parallel nature of these technologies. roche- is currently the most popular technology to investigate diversity in viral populations, thanks to its superior read length compared to other platforms such as illumina gaiix, but both technologies have advantages and disadvantages for analysing viral diversity, suggesting that experimental design should be planned carefully. currently, one run of roche- sequencing can produce one million reads, with a modal read length of 700nt - enough to sequence five 10kb viral genomes to a coverage of x each. this high coverage facilitates deep sequencing experiments, where reads are treated as a population sample containing information about the underlying genetic structure of the population. for example, deep sequencing has been used to study bottlenecks in early hcv infection and to identify mutations associated with drug resistance in hiv. a major challenge in deep sequencing experiments involves separating true variants from sequencing errors. although average error rates are constantly improving, error rates remain highly heterogeneous as they are affected by sequence context and alignment errors. this means that locally, error rates can still be quite high, leading to false positive snv calls. a number of tools are available for calling snvs within the context of sequencing errors. however, most approaches are tailored to calling snvs in human resequencing projects and are not suitable for viral deep sequencing data. when resequencing a human genome, snvs are either heterologous and can be expected to be present in 50% of all reads, or homologous and should be present in 100% of reads. commonly, a probabilistic approach is employed to decide whether a potential snv is either an error, a homologous snv, or a heterologous snv. in deep sequencing experiments, however, very low frequency snvs, including snvs present in less than 1% of reads, are often of interest. for instance, the presence of low frequency variants in the viral population of a single host may determine the potential for drug resistance and treatment failure. to date, deep sequencing data has been analysed by applying various read quality filters before and after alignment and then further manually filtering the results, according to the researcher’s understanding of the sequencing technology and the individual data set. filtering may also be performed with a software program such as varscan, an snv caller which relies on read coverage, base qualities, and variant frequencies to separate true snvs from errors. an emerging statistical method for detecting snvs in deep sequencing data involves performing a statistical test for each potential snv, taking into account a prior positional error model. this method detected snvs with a population frequency of % in control illumina data. however, the requirement for a prior positional error model relies on the availability of multiple replicates, which is not always feasible. a more recent approach to snv detection in viral populations, implemented in the program v-phaser, utilises the covariance between observed snvs present on the same read to increase sensitivity. the power of covariance was also harnessed in a recent study of hcv early evolution, where snv calling was achieved by parsing haplotypes reconstructed using probabilistic clustering performed by the program shorah. shorah uses a bayesian approach to group reads into haplotypes. both sensitivity and specificity are improved by assuming that the consensus sequence for each cluster is the true haplotype sequence; all other variants present in the reads within a cluster are considered to be errors. shorah has been shown to accurately reconstruct variants with a population frequency of %, and to efficiently remove random pcr errors. error correction via probabilistic clustering will fail, however, when errors themselves appear clustered. this can occur in regions of low diversity characterised by an absence of true variants, or when sequencing errors are systematic. the potential for systematic errors to confound clustering may explain why an empirical k-mer based algorithmic approach was found to be more efficient than shorah at removing false haplotypes, although both approaches detected all true haplotypes. there is growing evidence that systematic errors may exhibit a strand bias, where a particular error is more likely to occur in reads traversing the genome in one direction as opposed to the other. for instance, in a recent analysis addressing errors from control illumina data using prior positional error models, error rates between forward and reverse reads were uncorrelated. because of this, several snv detection methods consider forward reads and reverse reads separately, using some combination of the two obtained p-values to inform the final snv calling decision. researchers have also used strand bias directly to filter out false positive snvs. for example, in a deep sequencing analysis of hiv- quasispecies, any potential snvs with more than fold difference in the forward and reverse read count were considered to be systematic errors, while in the hcv study mentioned above, raw results from shorah were further cleaned by removing any potential snvs that showed a difference of more than 10% between the number of forward and reverse reads. the latest version of samtools formalises this approach, reporting a p-value based on an exact test of strand bias. a two tailed fisher’s exact test is also included in the recently released snv caller lofreq. here, we show that in addition to its primary function of haplotype reconstruction, shorah may also be used to reliably call snvs. this is demonstrated through analysis of control roche- data from both hcv and hiv samples, and by development of a statistical argument in support of snv calling via probabilistic clustering. as strand bias has been observed in systematic errors retained during clustering, we also investigate the benefits of combining probabilistic clustering with a statistical test of strand bias.  our snv calling method is based on two assumptions: it is more powerful to compare several variants that co-occur on the same read simultaneously than individually; and the frequency of sequencing errors, but not of true snvs, is often higher on reads sequenced in one direction than the other. while there is empirical evidence for the second assumption, the first one can be motivated and made precise by a simplified statistical model of snv detection. in this model, two haplotypes of frequencies f and − f that differ at d sites co-occurring within an observed read, can be distinguished in the presence of a per site error probability ϵ only if f > ϵd. for d =  this result confirms the intuitive notion that a single snv can be detected only if the error rate is smaller than the snv frequency. however, if the haplotypes differ at two sites, then the probability that these two variants are both errors is ϵ and the limit of haplotype detection decreases to f > ϵ in general, considering multiple variant sites at once by comparing haplotypes increases snv calling accuracy. this finding motivates analysing entire reads or large segments of reads in order to include as many snvs as possible, as is done in read clustering methods such as shorah. to investigate real life performance of snv calling via probabilistic clustering, we analysed roche- data from two distinct control experiments. for the first experiment, the same library of cloned moderate-diversity hcv fragments mixed at equal quantities was subjected to two independent roche- sequencing runs. for the second experiment, one round of roche- sequencing was performed on a high diversity pcr amplified set of hiv fragments mixed according to a geometric series. snvs were called using either probabilistic clustering, quality score based filtering, or a per site quality score based statistical test. a statistical test of strand bias was applied to the results from shorah and varscan, utilising either a beta-binomial model of forward read distribution for true snvs, or fisher’s exact test. our experimental design allowed us to assess the performance of probabilistic clustering, and variation in the true snv strand bias distribution, both between samples with different diversity levels and between individual sequencing runs. alignment and clone frequencies for hcv run and hcv run % and % of reads aligned uniquely to the reference, giving an average read depth of reads and reads, respectively. for the hiv data, % of reads aligned, giving an average read depth of reads. true haplotype frequency estimates are given in table  and table . for all samples, the frequencies of at least some haplotypes differed substantially from their intended frequencies. as the estimated frequencies for hcv run and hcv run are very similar, it is likely that deviation from intended frequencies is a result of mixing errors, rather than bias introduced during library preparation and sequencing. lists of true snvs for both hcv and hiv samples are provided in the additional files and raw snv calls for hcv run and hcv run both shorah and varscan detected all true snvs, i.e. the false negative rate was zero, while lofreq called all but one true snv in each run. varscan, however, also called up to and false positives, more than seven times as many as shorah. more errors were incorrectly identified as snvs for hcv run compared to hcv run by both varscan and shorah, indicating a higher error rate for hcv run varscan’s recall was greater than shorah’s for the hiv data, with 93% of true snvs identified, compared to 81% for shorah. however, shorah’s recall was greater than lofreq’s, which only identified 64% of true snvs. for the hiv sample, all of shorah’s false negatives had true population frequencies under %, while one of varscan’s false negatives had a true population frequency of %. upon inspection of the alignment and the varscan file, this appeared to be the result of confusion in reporting positions with more than one variant allele. most of lofreq’s false negatives were very low frequency, with three having true population frequencies of under 6%. varscan’s superior recall was accompanied by a high number of false positive calls – over the reference 766bp considered, false positive snvs were also called. by comparison, shorah only made two false positive calls. overall, lofreq exhibited an exceptionally low false positive rate. for the hcv runs, one unique false positive was identified in each run, while no false positives were called by lofreq for the hiv data. snv calling statistics for shorah, varscan, and lofreq. for shorah and varscan, snv calls without the strand bias test and using different values of the beta-binomial dispersion parameter σ in the strand bias test are given, with σ = corresponding to a binomial forward read distribution. for lofreq, the results of applying our strand bias test are absent as this software does not report forward and reverse strand counts. for all snv calling methods, the results of applying fisher’s exact test to the raw output are also given. reported statistics include true positives, i.e., genomic sites with a variant matching a known true variant, false positives, i.e., genomic sites with a variant that is not a known true variant, and false negatives, i.e., known variants which are not identified by the relevant snv calling method. individual genomic sites may contribute to both true positives and false positives. recall ) and precision ) are also reported. strand bias test in an attempt to improve shorah’s performance under conditions of sub-optimal diversity or clustered systematic errors, we applied a statistical test of strand bias to the raw snv calls. for each potential snv detected, we test the hypothesis that the distribution of forward reads carrying this snv follows a beta-binomial distribution, with mean value equal to the overall forward read distribution at that position. the overdispersion of the forward read distribution is controlled by the parameter σ, which we estimate separately for each sequencing run using maximum likelihood methods based on the observed forward read distribution for all true snvs within the sample. ml estimates of σ for hcv run hcv run and hiv samples were,, and respectively. to avoid overfitting, we applied tests using all estimated values of σ to each dataset, also performing a test with σ=. for comparison with the strand bias test implemented in other software, we also performed a fisher’s exact test of strand bias. snv calling statistics for strand bias tests applied to the raw results of both shorah and varscan are given in table  the results of lofreq’s own implementation of fisher’s exact test are also given. applying a strand bias test to the raw snv calls for both shorah and varscan improved precision in all cases, except for hiv snv calls made using shorah with σ =, where precision remained unaltered compared to the raw shorah output. in general, choice of σ did not have a substantial impact on precision. however, smaller values of σ were more likely to result in reduced recall, as was fisher’s exact test. recall was also compromised by lofreq’s implementation of fisher’s exact test. this was particularly evident for the hiv data, and is consistent with an overdispersed forward read distribution for true snvs, justifying our choice of a beta-binomial forward read distribution. for the beta-binomial test, the reduction in recall was minimised in all cases by using σ =. to further investigate the appropriateness of a beta-binomial based strand bias test, compared to either a binomial based test or fisher’s exact test, we generated precision recall curves for each test applied to snv calls from shorah and varscan. in all cases, across the full range of recall values obtained the beta-binomial test gave equivalent or higher precision than either of the other tests, confirming its suitability. for the hcv samples, the strand bias test was particularly good at reducing errors in results from shorah, with more than 70% of errors eliminated. after the strand bias test, less than six errors for every reference sites remained for hcv run while less than four errors for every reference sites remained for hcv run the strand bias test also eliminated 48% of errors in varscan’s output, however as varscan’s raw output contained more errors than shorah, this still left a minimum of errors for every reference sites. lofreq’s strand bias test successfully eliminated the single false positives originally identified in hcv run and hcv run figure  presents a detailed analysis of snv calling by frequency for varscan and shorah, for snvs with frequencies less than 15%. lofreq results are excluded from this analysis, due to its near perfect precision for all data sets. figure 2a shows that for the high diversity hiv data, perfect recall could be achieved for snvs with population frequency greater than around %. any reduction in recall due to the strand bias test was limited to snvs with frequencies below this detection limit. thus, although consideration of both recall and precision indicates the strand bias test does not improve on the raw hiv shorah results, it is not significantly detrimental. the lack of improvement is due in part to the extremely high precision of obtained by probabilistic clustering alone. for the hiv varscan results, although the strand bias test did increase precision, more than errors for every reference sites remained. figure 2b indicates that this is the result of low frequency false positives. restricting snv calls to those covered by a minimum of reads therefore rectified this problem to some extent, however over errors for every reference sites were retained. recall was also adversely affected, due to the presence of very low frequency true snvs in the hiv sample. using a minimum read coverage of had no effect on hiv shorah results. figures 2b to 2d describe the number of false positives observed for each data set, binned in frequency increments of 1%. in general, shorah calls less false positives than varscan, particularly for very low frequency snvs. in most cases, the strand bias test was able to reduce or eliminate false positives across the range of frequencies observed. the effectiveness of the strand bias test was most notable for shorah snv calls from the moderate diversity hcv samples. discussion we have shown that snv calling via probabilistic clustering is a powerful technique allowing detection of snvs at frequencies lower than the corresponding sequencing error rate, as long as sample diversity is high and errors are random. when these conditions are not met, a statistical test of strand bias improves snv calling precision. the distribution of forward reads covering a true snv at any one position could be expected to follow a binomial distribution, with mean value equal to the fraction of forward reads at that position multiplied by the total number of reads covering the variant. however, our results show that more accurate snv calling is achieved if the expected number of forward reads is modelled using a beta-binomial read distribution, compared to either a binomial distribution or applying fisher’s exact test. variation in the number of false positives between the hcv runs indicated error rates differed markedly, even though all samples were sequenced using roche- titanium technology. notably, the ability to utilise empirical run specific error rates was one of the motivating factors behind lofreq’s development. a similar result has been reported for illumina sequencing as well, with individual runs sequenced on the same machine having dissimilar error profiles. snv calling precision after strand bias testing was relatively robust to choice of dispersion factor σ. using a very low value of σ, however, compromised snv recall, thus we would recommend choosing a value of σ closer to our upper estimate of when applying our test to real data. for this reason, σ is set to by default in our updated version of shorah. for the hcv samples, which featured a moderate number of relatively high frequency true snvs, the beta-binomial strand bias test provided a substantial reduction in error rates with no cost to true positive recall. this was true when applied to both the raw shorah output and the raw varscan output. however, both the relative reduction in error and absolute precision attained was greatest when the test was applied to the raw shorah output. unfortunately, it was not possible to apply our beta-binomial strand bias test to lofreq’s raw snv calls, as variant counts by forward and reverse read are not given in lofreq’s output. lofreq’s own fisher’s exact test of strand bias successfully removed all false positives, however it also removed a number of true positives. given the results reported here, we believe that a beta-binomial strand bias test may also remove these false positives, with substantially less cost to true positive recall. for the hiv sample, the raw shorah output gave the best results. this can be explained by the extremely strong performance of error correction via probabilistic clustering in this case. the ability to correctly identify true haplotypes via probabilistic clustering is a function of the frequency of the haplotypes, the number of variant positions within the haplotypes, and the error rate of the sequencing technology. whenever there is more than one true variant within a haplotype under consideration, it becomes possible to detect snvs with population frequencies below the error rate of the sequencing technology. consideration of hiv false positives called using varscan indicates that locally, error rates may exceed 5%. for the hiv data, we achieved a detection limit well below this local error rate using shorah alone, validating our argument. figures 2a and 2b show that perfect accuracy can be achieved by applying a detection limit of 3%, while a detection limit of only % gives perfect recall with only two false positives. another implication of our model is that accuracy should improve as the number of snvs within a window of sites considered simultaneously increases, as demonstrated by the fact that shorah performs better on the intrinsically high diversity hiv sample than on the moderate diversity hcv samples. our argument is further supported by varscan’s performance. although for the hiv sample varscan detected more true positives than shorah, the corresponding extremely large number of low frequency false positives within the varscan results suggests that this is the result of only considering one site at a time and indiscriminately calling low frequency snvs below the associated detection limit. lofreq’s approach of discriminating snvs individually by inferring an error probability based on the quality score assigned to individual nucleotides clearly overcomes this problem, although at the expense of recall. while our strand bias test could not improve shorah’s raw snv calls for the hiv sample, we consider the slight reduction in recall with σ = to be an acceptable cost, compared to the potential benefits of the test. it is usually not possible to know, a priori, how diverse a deep sequenced population may be. even within highly diverse populations, genomic regions under purifying selection may have sub-optimal levels of diversity, allowing errors to appear clustered as true snvs. our strand bias test can help remove such errors, while leaving snv calls in highly diverse regions relatively untouched. our strand bias test could also be used to provide a level of confidence in results. if snv calls are unaffected by applying a strand bias test, as in the case of our hiv sample, then it is likely that diversity was sufficient for probabilistic clustering to perform at its best. applying our strand bias test to raw varscan snv calls also improved precision in all cases. however, false positive snv calls could not be reduced to a biologically useful level. the relative advantage of applying a strand bias test to shorah’s output compared to varscan’s output can be understood by considering how each of these programs identifies snvs. varscan distinguishes errors from true snvs by examining the individual base quality scores of reads covering the potential snv. varscan is therefore reliant on accurate assignment of quality scores during the sequencing platform’s base calling pipeline. if an error nucleotide is falsely assigned a high quality score, then varscan will retain it. although lofreq is also dependent on quality scores assigned by the sequencing platform, it has some resilience to this effect as quality information from all aligned nucleotides at a position is considered when making a call; errors at a given position would have to systematically be assigned high quality scores for lofreq to incorrectly call an snv. in contrast to varscan, shorah is independent of quality score information. shorah works by removing any potential snvs that are randomly distributed among reads, rather than clustering with nearby potential snvs - regardless of whether they occurred during library preparation, pcr amplification, or read synthesis. errors retained during clustering are likely to be systematic sequencing errors, which tend to exhibit a strand bias. methods statistical analysis of snv calling based on probabilistic clustering snv calling from shorah output involves parsing a set of ‘haplotypes’ for variants with respect to a reference sequence. these haplotypes are generated by deriving consensus sequences from groups of probabilistically clustered aligned reads. because only the consensus sequence within each cluster is considered, snvs will be correctly called, provided reads are correctly clustered so that the haplotypes can be identified from the centroids of the clusters. the detection limit for true snvs therefore depends on the detection limit for haplotypes. the detection limit for haplotypes is the haplotype frequency at which observed variants are equally likely to be errors or biological variants. for example, if two haplotypes differ at only one site, then the detection limit is intuitively equal to the error rate, ϵ. if, however, the haplotypes differ at two sites, then it takes two specific sequencing errors to mutate a read originating from one haplotype into a read that is equal to the other haplotype. this occurs with a probability of the order ϵ now, ϵ2< ϵ, allowing haplotypes to be detected at frequencies lower than the error rate provided two sites are considered simultaneously. to investigate this property more formally, we study the probability of correctly assigning a read to either of two haplotypes: a minor haplotype with frequency f containing the same variants as the read under consideration, and a major haplotype with frequency − f. to simplify our analysis, we consider only those sites where these two haplotypes differ, as other sites are uninformative when choosing which haplotype to assign the read to. according to this model, the major and minor haplotypes may be represented by the binary sequences 0 = ... and 1 = ... of length d, where d is the hamming distance between these two haplotypes. now, let the random variable h =  denote the true haplotype, such that p =  − f and p = f. observation errors are assumed to occur independently with uniform per site probability p = ϵ, where each random variable ei indicates the occurrence of an error at site i. the observed read, y, is the true haplotype h subject to observation error: yi = if  =  or, and yi = otherwise. the full model is defined as follows: h∼bernoullif,h∈0¯,1¯ei∼bernoulliϵ,ei∈,i=...,dyi∼hi+eimod i=...,d. in particular, p = ϵd and p  = d. in order to assign the read correctly, we need to decide whether the observed read, y =  is more likely to originate from the minor haplotype h =  without error or from the major haplotype h =  through errors. this decision represents a crucial step in the iterative probabilistic assignment of reads to haplotypes within shorah – if the constituent reads of a low frequency true haplotype are just as likely to originate from a higher frequency haplotype through sequencing errors, then they may be incorrectly assigned to the higher frequency haplotype. in this case, the low frequency haplotype falls below the limit of haplotype detection, and may not be represented in the final set of clustered reads. based on this model, the assignment will be correct if: ph=0¯|y=1¯<ph=1¯|y=1¯. using bayes’ theorem and eqs. this inequality becomes ϵd < df, or ϵ<1+1−ff1/d− for noise levels ϵ below this upper bound, we would assign the reads correctly. for d =  the condition is ϵ < f, and for ϵ and f small, it is approximately ϵd < f. thus for a given per nucleotide sequencing error rate, the frequency limit for haplotype detection decreases as the number of true snvs observed within a read increases. in reality, the process of haplotype inference is of course more complicated than presented here – our simplified argument is designed to merely demonstrate in a quantitative way that by observing the co-occurrence of variant sites, more power can be achieved when calling low frequency haplotypes than by considering variant sites independently. strand bias test we assume that at any given genomic position, the number of forward reads carrying a variant follows a beta-binomial distribution. more formally, let fi and ri, respectively, denote the total number of forward and reverse reads covering position i. the fraction of forward reads is μi = fi/. we define nib as the total number of reads covering the variant base b at sequence position i. the number of forward reads carrying variant b at position i is modelled by the beta-binomial random variable xib ~ betabin such that ≤ xib ≤ nib and pxib=x=Γnib+1Γx+1Γnib−x+1×Γ1/σΓx+μi/σΓnib+1−μi/σ−xΓnib+1/σΓμi/σΓ1−μi/σ the expected value of xib is e = nibμi, and σ is a dispersion parameter, defined such that varxib=nibμi1−μi1+σnib−11+σ. the dispersion parameter σ determines the spread of the distribution. we assume that, unlike nib and μi, σ does not change between each variant within a sample. both nib and μi can be easily estimated directly from the data for each individual potential snv under question. we employed ml estimation to obtain σ estimates for each sample, by considering the k true snvs xi1b …,xikbk. the likelihood function is: lσ=∏j=1kpxijbj=xijbj. ml estimation was performed in r with the package mle using the default nelder and mead optimisation method and a start value of one.. our null hypothesis, h is that for any given true snv, the number of forward reads follows a beta-binomial distribution with parameters μ, n, and σ. the null hypothesis h is rejected if, under h the probability to observe data as extreme or more extreme than that observed is smaller than a threshold p. that is, given an observation xib, we reject h if min {Φ, − Φ } × 2 < p, where Φ = p is the beta-binomial cumulative distribution function, with parameters estimated as described above. if, for a potential snv, the null hypothesis is rejected, then the potential snv is considered to be an error. as we do not care which direction the strand bias occurs in, a two-sided test was performed by multiplying the smallest tail of Φ by two. we then applied the benjamini-hochberg correction for multiple testing, giving the final p-value p. if p < , we reject the potential variant as an error. to investigate the effect of σ on the strand bias test, strand bias tests with each estimate of σ were applied to snv calling results obtained from shorah and varscan for all data sets. a test with σ = was also performed, which models a binomial forward read distribution. for comparison, two-sided fisher’s exact tests were also applied, based on counts of the forward and reverse reads, again with a benjamini-hochberg correction for multiple testing. it was not possible to apply our strand bias test to snv calls from lofreq, as lofreq’s output does not record the number of forward and reverse reads covering a variant. lofreq does however report a p-value based on a fisher’s exact test of strand bias. thus we were able to report the results of a fisher’s exact test of strand bias applied to raw lofreq snv calls, again applying a benjamini-hochberg correction for multiple testing within r and rejecting any variants for which p < . to further compare the performance of fisher’s exact test, a binomial strand bias test, and our beta-binomial test, precision recall curves were constructed within r for each test applied to shorah or varscan snv calls from hcv run and hcv run this was done by calculating the precision and recall for p-value cut offs varying between zero and one and plotting the results. sample sequencing, alignment, and frequencies for hcv e1/e samples, a mixture of four different 2324bp fragments, each covering the 3′ end of the hcv core to the 5′ end of the p region was cloned into the pgem-t easy vector and sequenced using roche- flx titanium as described previously. to investigate variation in σ between sequencing runs, a second round of roche- flx titanium sequencing was performed on this same library of hcv clone mixtures. for the hiv gag/pol-gene fragments, pcr products from subtype-b clinical isolates were propagated within the pcrii-topo vector, excised, and subjected to a single round pcr before roche- flx titanium sequencing, as previously described. prior to alignment, reads were trimmed using a sliding window of 8nt. the trim was made at the start of the first 8nt window from either end with an average quality greater than for all samples, alignment was performed using bwasw, with default parameters. for the hcv samples, alignment was performed against the entire pgem genome with an hcv 1a strain insert, corresponding to the inserted fragments. for the hiv sample, alignment was performed against the appropriate region of the hxb hiv reference genome, from base to to estimate the true frequency of each fragment within the mixture, the known sequences of each fragment were used to locate snv positions uniquely identifying each fragment. fragment frequencies were then calculated by averaging the frequency of all snvs unique to a fragment haplotype, then normalising the results so that the sum of all haplotype frequencies was one. snv calling shorah was run on all three alignments described above. for each run, α was set to, window size was set to with a step size of and the number of iterations to times the window coverage for each window. to call an snv from the raw shorah output, locally reconstructed windows were parsed for variant positions. for an snv to be called, the same variant was required to be found in at least two of three overlapping windows, with a posterior probability greater than. this approach was compared to single site snv calling, as implemented in lofreq and varscan. lofreq is a single site variant caller that uses per nucleotide quality score information to inform a statistical test for each potential snv. lofreq is tailored to calling low frequency variants in deep sequencing data, and is thus ideally suited to identifying variants in viral populations. snv calling was performed using the script lofreq_snpcaller.py with default options. varscan was also employed due to its ability to report the required read statistics for the strand bias test, and its single site mode of snv calling. however, we acknowledge that varscan is intended for use on diploid resequencing data rather than deep sequencing of viral populations, and do not imply that results reported here are indicative of performance under normal usage conditions. varscan’s pileup2snp command was used to call snvs. default values for quality filtering were used. minimum variant frequency was set to zero, to allow comparison with shorah. for all variant callers employed, we limited called variants to those covered by more than one read. to allow a fair comparison, snv calling using lofreq and varscan was only attempted in areas of the genome covered by three shorah windows. this corresponds to positions to of the hcv fragments, and positions to of the hiv fragments. in total, true snvs were analysed for the hcv samples, and true snvs were analysed for the hiv sample. an attempt to call snvs using v-phaser was made, as this program combines joint snv calling with quality score based methods. however, the program crashed when used on any of our samples, possibly due to the high coverage of samples. correspondence with the authors did not resolve the issue, thus the comparison was not able to be made. snv calling was also attempted using an empirical k-mer based algorithm. after running for more than cpu hours, this software did not produce a result, at which point this comparison was also abandoned. snv calling with deepsnv and using the approach described by flaherty et al. was not attempted due to the lack of control or multiple samples. for shorah and varscan, raw results, as well as results subjected to our strand bias test, are reported. for the purposes of this paper, all strand bias tests were performed in r as described above. this was done to facilitate an exact comparison between varscan and shorah results. however, in conjunction with this manuscript we have released a new version of shorah with integrated snv calling and strand bias testing, which does not require r.  we have demonstrated that by combining probabilistic clustering with a statistical test of strand bias, highly accurate snv calling can be achieved. this is because the strengths and weaknesses of these two methods are complementary: probabilistic clustering relies on errors occurring randomly and independently, while strand bias tests detect systematic errors with atypical read direction distributions. using a statistical model and through analysis of control data we investigated the conditions under which probabilistic clustering can be used to call snvs. we show that the limit of snv detection is ϵd < f, where ϵ is the rate of randomly occurring sequencing errors, d is the number of true variants observed within a read length, and f is the frequency of the haplotype to which the variants belong. in general, the snv detection limit decreases as true diversity increases. this is consistent with our empirical analysis of snv calling from control data with various diversity levels. importantly, for snv calling via probabilistic clustering to have an advantage over approaches that consider individual variants one at a time, more than one true variant must be contained within a read length. to insure against situations where diversity is inadequate or errors are non-random, we suggest a strand bias test be applied to snv calls generated via probabilistic clustering. the test we propose uses a beta-binomial model of forward read distribution for true snvs. by comparing our test to the performance of either fisher’s exact test or a test based on a binomial model of forward read distribution, we show that the forward read ratio of true snvs is overdispersed. our analysis indicates that as long as this overdispersion is taken into account, a statistical test of strand bias will not affect recall. this means such a test can be safely applied and will not adversely affect accuracy, even in situations where no improvement in precision is obtained. while our analysis is presented in the context of snv calling using shorah, the insights into snv calling via probabilistic clustering should be applicable to any method that considers multiple variant sites simultaneously. furthermore, our strand bias test has the potential to improve snv calling precision for any variant detection pipeline. availability the snv caller and strand bias test have been integrated into shorah, available at http://www.cbgethz.ch/software/shorah. haplotypes, aligned reads, and reference sequences used in this study may be obtained from https://wiki-bsse.ethz.ch/display/shorah/data. contact: shorah@bsse.ethz.ch. competing interests the authors declare that they have no competing interest. authors’ contributions km implemented the strand bias test, analysed control data, and in conjunction with oz developed the software shorah v. nb developed the snv calling model. km wrote the manuscript, with editing by fl, oz, and nb. nb and fl conceived of the overall study. rb and oz provided control data. all authors read and approved the final manuscript. supplementary material additional file 1 tab delimited list of true snvs for hcv haplotypes. format: position \t reference nucleotide \t variant nucleotide \t frequency run \t frequency run to facilitate analysis with supplied reference sequences and alignments, position is given with respect to the start of the vector; subtract to get position with respect to the hcv fragment. frequencies are calculated based on the estimated true frequency of haplotypes. click here for file additional file 2 tab delimited list of true snvs for hiv haplotypes. format: position \t reference nucleotide \t variant nucleotide \t frequency. frequencies are calculated based on the estimated true frequency of haplotypes. click here for file acknowledgements km was supported by an australian postgraduate award. fl and rb are both nhmrc training fellows. we acknowledge support from ach and nhmrc. part of this work was supported by the swiss national science foundation under grant number cr32i.
by combining probabilistic clustering with a statistical test of strand bias, snvs may be called from deeply sequenced viral populations with high accuracy. finally, we show that in situations where probabilistic clustering retains false positive snvs, applying a strand bias test based on a beta-binomial model of forward read distribution can improve precision, with negligible cost to true positive recall. such experiments harness the high coverage afforded by next generation sequencing protocols by treating sequencing reads as a population sample. we demonstrate these findings on two control data sets, showing that snv detection is more reliable on a high diversity human immunodeficiency virus sample as compared to a moderate diversity sample of hepatitis c virus. current protocols are characterised by high false positive rates, with results requiring time consuming manual checking. by statistical modelling, we show that if multiple variant sites are considered at once, snvs can be called reliably from high coverage viral deep sequencing data at frequencies lower than the error rate of the sequencing technology, and that snv calling accuracy increases as true sequence diversity within a read length increases.
 deep sequencing is a powerful tool for assessing viral genetic diversity. such experiments harness the high coverage afforded by next generation sequencing protocols by treating sequencing reads as a population sample. distinguishing true single nucleotide variants from sequencing errors remains challenging, however. current protocols are characterised by high false positive rates, with results requiring time consuming manual checking.  by statistical modelling, we show that if multiple variant sites are considered at once, snvs can be called reliably from high coverage viral deep sequencing data at frequencies lower than the error rate of the sequencing technology, and that snv calling accuracy increases as true sequence diversity within a read length increases. we demonstrate these findings on two control data sets, showing that snv detection is more reliable on a high diversity human immunodeficiency virus sample as compared to a moderate diversity sample of hepatitis c virus. finally, we show that in situations where probabilistic clustering retains false positive snvs, applying a strand bias test based on a beta-binomial model of forward read distribution can improve precision, with negligible cost to true positive recall.  by combining probabilistic clustering with a statistical test of strand bias, snvs may be called from deeply sequenced viral populations with high accuracy. 