See also
fastq_filter command
Quality scores
Average Q is a bad idea!
Setting an e.e. threshold for long reads
Edgar & Flyvbjerg 2014 paper on expected errors and
Bayesian assembly
Phred scores
Next-generation reads specify a Phred score for each base, also known as a Quality or Q score. The Q score is an integer, typically in the range 2 to 40. Q indicates the probability that the base call is incorrect (P_e). For example, Q=2 means that the error probability is 63%, so the machine is reporting that the base is more likely to be wrong than right, while Q=20 corresponds to an error probability of 1%.
Expected number of errors
Assume that the Phred scores are good, in other words assume that the machine makes good estimates of the error probabilities. Now we can measure the overall quality of a read by calculating the expected number of errors in the read. The expected value of a quantity under a probability distribution is a technical term in statistics, so let me explain what that means in this case.
Take the error probabilities according to the Q scores in a given read, and imagine a very large number of reads containing errors that occur with those probabilities. For example, if the first base in the read has Q=20, then 1% of the reads in our collection will have errors in that position. Some of the reads will have no errors, some will have one error, and so on. The Q scores can't tell us exactly how many errors there are in a given read. However, the Q scores can tell us what the average number of errors would be in that large collection, and this is the definition of the expected number of errors. Because we're taking an average, the value is real rather than integer, and can be less than one.
Calculating the expected number of errors and most probable number of errors
It turns out that it is easy to calculate the expected number of errors: it is the sum of the error probabilities. The most probable number of errors (E*) is also easy to calculate. First calculate E = expected errors = sum P_e. Then round down to the nearest integer, and this is the most probable number of errors. In symbols: E* = floor(E). For example, if E < 1, then the most probable number of errors is zero. Proofs of these results are given in Edgar & Flybjerg (2014).
Quality filtering by expected errors
Small E means high quality, large E means low quality. We can therefore filter low quality reads by setting a maximum number of expected errors (E_max) and discarding reads with E > E_max. This is done by running the fastq_filter command.
Choosing the maximum expected error threshold
A natural choice is E_max = 1, because the most probable number of errors of a filtered read is zero. Though of course, we expect some of them to have one or more errors (this can't be entirely avoided, however stringent the filter, because Q scores are probabilities, not certainties). If you want to filter more stringently, you might choose something like E_max = 0.5 or E_max = 0.25.
Overlapping paired reads
If you have paired reads that overlap, then these can be exploited to improve the error filtering. The paired read merger (assembler) should report re-calculated Phred scores in the overlapping region to reflect that we have two independent observations of the same bases. If the forward and reverse read agree on a base, then our confidence that the base is correct should increase. The Q score should be higher, corresponding to a lower error probability. Conversely, if the base calls disagree, then the re-calculated Phred score should be lower than both the original base calls. The Q score in the merged read (contig) should be calculated as the posterior probability according to Bayes theorem where the P_e's from the forward and reverse read are the prior probabilities.
Quality filtering by maximum expected errors should be performed as a post-processing step after reads have been merged in order to take advantage of the re-calculated Phred scores, which should be better predictions of the true error rates.
It is important to use the USEARCH fastq_mergepairs command for read merging because according to my tests, other read mergers do not correctly calculate the posterior Q scores, including PANDAseq, COPE, PEAR, fastq-join and FLASH. Some of these mergers also have high rates of false positive and incorrect merges, including especially PANDAseq and the make.contigs() command in mothur.
OK, the concept is nice, but how well does this work in practice?
For Illumina reads, this approach works very well, because the Illumina Q scores are pretty good (much better than 454). The accuracy of the Q scores varies with the sequencing machine, base-calling software and probably other factors such as reagents. Amplicon sequencing tends to have less accurate Q scores than shotgun, so you can't assume that results on a PhiX library will be similar to results for amplicon reads. The Q scores may systematically under- or over-estimate the error probabilities in different runs. Errors do have some tendency to correlate, i.e. to cluster in certain reads, so they're not perfectly independent as assumed by the theory. Also, Q scores can only predict sequencer error per se, they cannot predict errors introduced before sequencing, e.g. during PCR. But despite these caveats, the approach works well in practice, better than anything else I've tested.
Long reads
See this FAQ: Should I increase the expected error threshold for long reads?