2
$\begingroup$

i am trying to calculate mean Phred scores for my sequencing data, but i feel not very comfortable about it. There are actually two ways of calculating. (I just use an existing sample) giving: 3 reads with a mean phred score of 20, 10 and 3 (meaning Error rates of 1%, 10% and 50%, respectively)

Way 1 to calculate: Just take the arithmetic mean of the q scores:

(20+10+3)/3 = 11

Way 2: Take the arithmetic mean of the error rates:

(0,5+0,1+0,01)/3 ~= 0.2

And then calculate the q score by

-10*log10(0.2) = 6.99

(https://gigabaseorgigabyte.wordpress.com/2017/06/26/averaging-basecall-quality-scores-the-right-way/)

The Author of this Blog claims, that Way 2 is the correct version, but i still do not get why.

Is somebody able to answer this question?

Thanks a lot

$\endgroup$
0
5
$\begingroup$

Not sure if my explanation is any good but let's try...

First, let's convince ourselves that converting phred to probability is the right thing to do (i.e. your Way 2).

Imagine a read of length 10000 bp with half of the bases with Q=10 (or P=0.1) and the other half with Q=20 (or P=0.01). The number of wrong calls you expect is therefore 550:

n = 10000
(n * 0.01 * 0.5) + (n * 0.1 * 0.5) # -> 550 wrong calls

So the overall error rate is 550/10000 = 0.055 which corresponds to Q=12.6. If you follow your Way 2 you do get this answer:

((n * 0.5 * 0.01) + (n * 0.5 * 0.1))/n = 0.055
-10*log10(0.055) = 12.6 

If you follow your Way 1 you get Q=15, which means an overoptimistic error rate of ~0.03 or ~316 wrong calls. Why is that?

Because phred is on a non-linear, logarithmic, scale: a unit change towards the high values is less important than a unit change towards the lower end. For example, across 10000 bases a change from Q=30 to Q=31 means a change from 10 to ~8 wrong calls - a difference of 2 only. But a change from Q=10 to Q=11 means 1000 to ~800 wrong calls: the difference is a substantial ~200 calls.

Having said this, the straight average of phred scores is not "wrong". It measures the overall phred score but it does not reflect the expected number of wrong calls, which is probably what you want in this case.

$\endgroup$
1
  • $\begingroup$ Thanks alot for this clear answer! This just means, that for calculation of the expected wrong calls it is necessary to have a look first on the sequencing technique. If it is expected to have low/medium output reads (like a not that good MinION sequencing) it is necessary to go with way 2. For higher quality reads (like most of Illumina reads), both methods can be considered, but just to be safe, way 2 is the right one to choose. Thanks alot for the explanation $\endgroup$ – Don Castanio May 20 at 13:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.