美文网首页进击的GATK科研信息学BioStat
HaplotypeCaller是怎么计算PL的?(VCF解密)

HaplotypeCaller是怎么计算PL的?(VCF解密)

作者: 徐广惠_6f76 | 来源:发表于2019-05-17 14:23 被阅读36次
贝聿铭先生设计的卢浮宫

写在前面

PL是HaplotypeCaller等GATK变异检测软件在sample层面给出的注释,一般记录在VCF文件的FORMAT/sample一栏中。它代表了根据某位点变异情况给出的基因型判定不正确的可能性。

计算方法

计算PL的基本公式如下:
PL=−10∗logP(Data|Genotype)
公式右面的P(Data | Genotype)指的是:根据观察到的数据 D 计算出的基因型为 G 的条件概率(想知道P(D|G)是怎么计算的请戳这儿)。将P(D|G)取log值并乘以-10后将其转换为Phred-scale格式即为PL,然后将所有基因型的PL进行归一化,使得最有可能的基因型PL为0。

举个例子

假如某位点的参考等位基因是A,现在观察到一个read在该位点为T,并且我们现在有HaplotypeCaller根据这个read计算出来的各基因型的条件概率P(D|G)(当然如果有多个read,它们的贡献将会相加):

Alleles
Reference: A
Read: T

Conditional probabilities calculated by HC
P(AA | Data) = 0.000001
P(AT | Data) = 0.000100
P(TT | Data) = 0.010000

计算初始PL值:

我们要分别计算基因型为0/0, 0/1 和 1/1时的初始PL值,根据前述的公式,结果如下:

Genotype A/A A/T T/T
Raw PL -10 * log(0.000001) = 60 -10 * log(0.000100) = 40 -10 * log(0.010000) = 20

可以发现P(D|G)最低的基因型在转换为PL后变为最大值。这是在我们的意料之中的,因为PL指的就是这个基因型是不正确的概率。一个基因型的Raw PL值越小,代表越有可能是真实的基因型。

标准化:

在将PL值写入VCF之前,还要对Raw PL进行一个小小的计算:对所有的PL值进行标准化,使得最小的PL值为零。很简单,取所有基因型的PL值和最小PL值之间的差值即可。

Genotype A/A A/T T/T
Normalized PL 60 - 20 = 40 40 - 20 = 20 20 - 20 = 0

我们从中可以发现,最终的PL值和原始的P(G|D)之间的关系:我们取的原始P(G|D)之间依次相差100倍,最终的PL之间最终相差20,也就是100取Phred-scale后的值。这样我们扫一眼VCF里面PL值的大小,就可以方便的比较各个基因型之间可能性的差异了。

PL和GQ之间的联系和区别

GQ的值其实就是PL次小值和最小值之间的差异,由于PL的最小值总是0,所以GQ就是PL的次小值。在上面的例子中,GQ = 20 - 0 = 20。需要注意的是,为了节省计算空间,在GATK中,GQ值最大为99,就算实际计算出的GQ大于99,VCF中也只会记录为99哦。

参考资料

https://software.broadinstitute.org/gatk/documentation/article?id=5913

相关文章

网友评论

    本文标题:HaplotypeCaller是怎么计算PL的?(VCF解密)

    本文链接:https://www.haomeiwen.com/subject/mfbroqtx.html