Any suggestion for calculating overall Phred scale quality score for a sequence?

[sulicon]

Hi all,
I have some consensus nucleotides and ptential variant sites with quality scores (phred scale) generated by mpileup and bcftools. We are working on some protein analysis, so I'm trying to get the most reliable sequences from the alignment, and corresponding scores indicating how reliable a sequence should be.

I have met a problem to calculate the quality score for multiple nucleotides.
Suppose we have two nts with Phred scores a and b, then I think the probability that both of them are CORRECT should be:
(1-10**-(a/10)) * (1-10**-(b/10))
The problem is, when the score is high (say, 200), it generates an error in the calculation (the error rate is 0.0, can't be used to calculate the Phred score...)

I think there should be an approximate way to calculate the overall Phred score. For example (c is the overall phred score):
a = 30, b = 30 then c = 26.99
a = 60, b = 60, then c = 56.99
a = 100, b = 100 then c = 96.99
a = 150, b = 150, then c = 146.99

It seems that the overall score is always about a - 3.01 when a == b and a > 30 . But I can't figure out the general formula...

Any suggestion is appreciated.

--------------------
update
0.301 is the value of log10(2). The approximate value could be got by ignoring 10**-(a/10) * 10**-(b/10). A similar approximation could be made when a != b, I think...

[tonybolger]

As a quick hack, if you can live with the error of 'double counting' multi-nt errors, you can just factor out the 'weakest' quality. So assuming both a and b can't happen:

Perror(a) = 10**-(a/10)
Perror(b) = 10**-(b/10)

Perror (a or b) = Perror(a)+Perror(b)
(this double-counts if both happen at once)

Then just use:
m=min(a,b)
Perror (a or b) = Perror(m)*(Perror(a-m)+Perror(b-m))

If a and b are very different, the additional terms will effectively zero anyway, meaning that the combined error is almost exactly equal to the most likely error, which makes sense.

And you don't actually need to calculate Perror(m), you just keep it in Q space, and add it back in later (since multiplication of probabilities = addition of Q values).

You can also apply this approach selectively, say if the minimum quality score is at least say 30, where double errors will be very rare.

[sulicon]

Thanks a lot The key is skipping double errors when the quality score is high enough, as you mentioned in your last sentence.

[srasdk]

If you painstakingly follow the math, the resulting formula would be:

c = 0.5(a+b)
- 10log10( 10**(0.05*(a-b)) + 10**(0.05*(b-a)) - 10**(-0.05*(a+b)))

As you can see, the first term is an average.
The second term reduces the sum, since it will always be positive (a and b are positive making the sum under log > 1).
In your example (a=b and a = 20) the value under log 20 reduces to:
2-10**(-2) ~= 2. The higher 'a' the closer you are to equality
So yes, your observations are correct
No, it will not be the same for a != b or when you start using low quality values.

This formula is safe for all phred values.