为什么不能通过 GATK 的 PL 直接计算基因型剂量（Genotype dosage）
原创 黄树嘉 碱基矿工
image
GATK 的 PL 比较特殊，它是不能直接用于基因型剂量（Genotype dosage）的计算的。这次我们就来谈一谈这个问题。
有时候我们需要在项目中用基因型剂量来代替基因型（Genotype），特别是进行低深度（<10x）数据的全基因组关联分析（GWAS）时，就经常会做这个转换。这是因为低深度数据由于样本覆盖深度不足，在个体基因型上往往会存在更多的不确定性甚至Missing的情况，这会降低GWAS的功效（Power）——有时候还十分明显。这时如果能够将数据中的这种不确定性体现出来，是能够有效改善结果的Power的。基因型剂量恰好能描述这类不确定性，它描述的是某一个样本在一个位点上预期的突变碱基（即非参考序列碱基）个数，计算公式很简单，如下：
dosage = Pr(het|data) + 2 * Pr(hom|data)
其中，Pr(het|data) 指的是某一样本在该位点的基因型为杂合突变的后验概率，Pr(hom|data)指的是某一样本在该位点的基因型为纯合突变的后验概率，经过计算你会发现这个 dosage 将确定的基因型转换为了通过概率来描述的突变等位基因（Allelle）个数——并且是一个连续的数值，值域为0到2。
通常来说，你使用Freebayes、GATK、甚至samtools这些工具得到样本的变异结果之后，其对应的结果文件（VCF格式）中通常都有会有一个用来描述每一个样本在三种基因型上的后验概率值，找出来就可以计算了。只不过为了表示上的方便，这些后验概率值通常都会被转化为Phred-scale，一般用 PL 标签记录（如附图）。
image
通过后验概率计算 PL 的公式是：
PL = -10 log (Pr)
需要注意的是，这里的 log 是底数为 10 的对数，Pr 是基因型的后验概率，对于二倍体基因组（比如人类）将样本里每个位点的三种基因型代入进去就可以得到该样本在每一个位点上三种基因型的 PL 分别是多少了，举个例子：
* 参考序列碱基(Reference)：A
* 突变碱基：C
该位点上 AA、AC 和 CC 这三个不同基因型各自的后验概率假设如下：
Pr(AA|data) = 0.001
Pr(AC|data) = 0.010
Pr(CC|data) = 0.989
那么它们各自的 PL 分别是：
-10*log(0.001) = 30
-10*log(0.010) = 20
-10*log(0.989) = 0.05
从这个计算我们可以看出，PL 最小的基因型是最可信的基因型。反过来，也可以根据上面的公式很容易地将 PL 还原为基因型的后验概率。这样一来通过 PL 计算基因型剂量这本身应该是一个很简单的事情，事实上，bcftools 都有直接的计算命令可以使用。那我为什么还要大费周章专门写一篇文章来讨论呢？这个原因就出在GATK上。当你仔细去看 GATK 得到的 PL 时，你会发现事情不对了！你可以看到 GATK HaplotypeCaller 或者经过 GenotypeGVCFs 之后，后验概率最大的那个基因型它的 PL 竟然都是0，这时直接通过 PL 转换计算之后，所有样本的 Genotype dosage 不但是错的，而且还会出现大于 2 的情况——而这是不可能的。那这是怎么回事呢？原来，这是因为 GATK HaplotypeCaller 和GATK GenotypeGVCFs 所得到的 PL 并非是直接由基因型后验概率转换而来，而是经过了一次预处理之后才给出。那为何要做这个预处理呢？还是看上面我给出的例子，Pr(CC|data) 的 PL 等于 0.05，这个数和其他的两个整数放在一起多少显得不够“漂亮”，不够简洁！GATK 为了结果的简洁和清晰，就将三个基因型的后验概率全部除以那个最大的后验概率值，这个预处理在GATK中称之为“归一化(Normalized)”——其实就是将数据按照最大值进行了缩放，然后再计算 PL，并且全部取整。这样我们就明白了，做了这个除法运算之后原来后验概率最大的基因型，其概率值就都变成了1（如：Pr(CC|data)/Pr(CC|data) = 1），而其他的两个值，就相当于是最好基因型的几分之几，或者你也可以理解为最好基因型比其他两个分别好上多少倍。虽然这个计算改变了原来的值，但是却可以提升数据的解析度和可读性。因此，如果直接用现有的计算工具（bcftools +dosage），是一定得不到正确的结果的，这个时候，我们就得自己写程序来解决了。可是该怎么办呢？经过上面的描述之后，你可能也大致清楚了，这里的难点就在于要将原来最好基因型的后验概率值重新计算出来，怎么计算呢? 我还是以上面的例子为基础给大家列一下计算过程，一切就都清楚了：我们假设由 GATK 归一化的 PL 值转换（Phred-scale计算公式的逆运算）得到的基因型“相对”后验概率值（加上相对是为了和后面真正的后验概率值作区分）为: n_Pr(AA|data)，n_Pr(AC|data) 和 n_Pr(CC|data)。我们现在的目标是依据这三个值重新计算出它们真正的基因型后验概率值: Pr(AA|data)、Pr(AC|data) 和Pr(CC|data)。
将 PL 转换为概率值是一定要先做的，之后才能完成后续计算。
假设这三个基因型中后验概率值最大的是 Pr(CC|data) —— 如果你要假设为其它的两个也可以，那么依据上面的讨论我们知道：
n_Pr(AA|data) = Pr(AA|data) / Pr(CC|data)
n_Pr(AC|data) = Pr(AC|data) / Pr(CC|data)
n_Pr(CC|data) = Pr(CC|data) / Pr(CC|data)
将上面三个数学表达式相加，得到：
n_Pr(AA|data) + n_Pr(AC|data) + n_Pr(CC|data) = [Pr(AA|data)+Pr(AC|data)+Pr(CC|data)]/Pr(CC|data)
而我们知道三个基因型的原始后验概率之和一定是等于 1 的，所以，我们就可以得到，Pr(CC|data) 这个最大的后验概率值是：
Pr(CC|data) = 1/(n_Pr(AA|data) + n_Pr(AC|data) + n_Pr(CC|data))
刚好就是 1 除以三个相对后验概率值之和！得到 Pr(CC|data) 之后，剩下的 Pr(AA|data) 和Pr(AC|data) 也同样可以得到。
那么，通过 GATK 的 PL 计算基因型剂量的问题也就解决了：
dosage = Pr(AC|data) + 2 * Pr(CC|data)
最后，我将这个计算转换的过程写成了Python代码，可以直接使用（附图为部分代码示意），完整的代码只分享在知识星球——“达尔文生信星球”****之中了（链接：https://t.zsxq.com/Aay37aE）。不过，我在截取图片的时候，已经将计算dosage的核心代码包含在内了，如果此刻你觉得还不需要加入我的知识星球，那么也可以参考这一段代码去实现你的程序。
What I cannot create, I do not understand.
- Richard P.Feynman(理查德.菲利普斯.费曼)
image
网友评论