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...
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...
Comment