GATK 的 PL 比較特殊，它是不能直接用于基因型劑量（Genotype dosage）的計算的障斋。這次我們就來談一談這個問題刚盈。
有時候我們需要在項目中用基因型劑量來代替基因型（Genotype），特別是進行低深度（<10x）數(shù)據(jù)的全基因組關聯(lián)分析（GWAS）時峻呛，就經(jīng)常會做這個轉換栏妖。這是因為低深度數(shù)據(jù)由于樣本覆蓋深度不足乱豆，在個體基因型上往往會存在更多的不確定性甚至Missing的情況，這會降低GWAS的功效（Power）——有時候還十分明顯吊趾。這時如果能夠將數(shù)據(jù)中的這種不確定性體現(xiàn)出來宛裕，是能夠有效改善結果的Power的瑟啃。基因型劑量恰好能描述這類不確定性揩尸，它描述的是某一個樣本在一個位點上預期的突變堿基（即非參考序列堿基）個數(shù)蛹屿，計算公式很簡單，如下：
dosage = Pr(het|data) + 2 * Pr(hom|data)
其中疲酌，Pr(het|data) 指的是某一樣本在該位點的基因型為雜合突變的后驗概率蜡峰，Pr(hom|data)指的是某一樣本在該位點的基因型為純合突變的后驗概率，經(jīng)過計算你會發(fā)現(xiàn)這個 dosage 將確定的基因型轉換為了通過概率來描述的突變等位基因（Allelle）個數(shù)——并且是一個連續(xù)的數(shù)值朗恳，值域為0到2湿颅。
通常來說，你使用Freebayes粥诫、GATK油航、甚至samtools這些工具得到樣本的變異結果之后，其對應的結果文件（VCF格式）中通常都有會有一個用來描述每一個樣本在三種基因型上的后驗概率值怀浆，找出來就可以計算了谊囚。只不過為了表示上的方便，這些后驗概率值通常都會被轉化為Phred-scale执赡，一般用 PL 標簽記錄（如附圖）镰踏。
過后驗概率計算 PL 的公式是：
PL = -10 log (Pr)
需要注意的是，這里的 log 是底數(shù)為 10 的對數(shù)沙合，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 最小的基因型是最可信的基因型。反過來究履，也可以根據(jù)上面的公式很容易地將 PL 還原為基因型的后驗概率滤否。這樣一來通過 PL 計算基因型劑量這本身應該是一個很簡單的事情，事實上最仑，bcftools 都有直接的計算命令可以使用藐俺。那我為什么還要大費周章專門寫一篇文章來討論呢？這個原因就出在GATK上泥彤。當你仔細去看 GATK 得到的 PL 時紊搪，你會發(fā)現(xiàn)事情不對了！你可以看到 GATK HaplotypeCaller 或者經(jīng)過 GenotypeGVCFs 之后全景，后驗概率最大的那個基因型它的 PL 竟然都是0，這時直接通過 PL 轉換計算之后牵囤，所有樣本的 Genotype dosage 不但是錯爸黄，而且還會出現(xiàn)大于 2 的情況——而這是不可能的滞伟。那這是怎么回事呢？原來炕贵，這是因為 GATK HaplotypeCaller 和GATK GenotypeGVCFs 所得到的 PL 并非是直接由基因型后驗概率轉換而來梆奈，而是經(jīng)過了一次預處理之后才給出。那為何要做這個預處理呢称开？還是看上面我給出的例子亩钟，Pr(CC|data) 的 PL 等于 0.05，這個數(shù)和其他的兩個整數(shù)放在一起多少顯得不夠“漂亮”鳖轰，不夠簡潔清酥！GATK 為了結果的簡潔和清晰，就將三個基因型的后驗概率全部除以那個最大的后驗概率值蕴侣，這個預處理在GATK中稱之為“歸一化”——其實就是將數(shù)據(jù)按照最大值進行了縮放焰轻，然后再計算 PL，并且全部取整昆雀。這樣我們就明白了辱志，做了這個除法運算之后原來后驗概率最大的基因型狞膘，其概率值就都變成了1（如：Pr(CC|data)/Pr(CC|data) = 1），而其他的兩個值挽封，就相當于是最好基因型的幾分之幾已球，或者你也可以理解為最好基因型比其他兩個分別好上多少倍场仲。雖然這個計算改變了原來的值，但是卻可以提升數(shù)據(jù)的解析度和可讀性渠缕。因此鸽素，如果直接用現(xiàn)有的計算工具（bcftools +dosage），是一定得不到正確的結果的亦鳞，這個時候馍忽，我們就得自己寫程序來解決了⊙嗖睿可是該怎么辦呢遭笋？經(jīng)過上面的描述之后，你可能也大致清楚了徒探，這里的難點就在于要將原來最好基因型的后驗概率值重新計算出來瓦呼，怎么計算呢? 我還是以上面的例子為基礎給大家列一下計算過程，一切就都清楚了：我們假設由 GATK 歸一化的 PL 值轉換（Phred-scale計算公式的逆運算）得到的基因型“相對”后驗概率值（加上相對是為了和后面真正的后驗概率值作區(qū)分）為: n_Pr(AA|data)测暗，n_Pr(AC|data) 和 n_Pr(CC|data)央串。我們現(xiàn)在的目標是依據(jù)這三個值重新計算出它們真正的基因型后驗概率值: Pr(AA|data)磨澡、Pr(AC|data) 和Pr(CC|data)。
將 PL 轉換為概率值是一定要先做的质和，之后才能完成后續(xù)計算稳摄。
假設這三個基因型中后驗概率值最大的是 Pr(CC|data) —— 如果你要假設為其它的兩個也可以，那么依據(jù)上面的討論我們知道：
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)
將上面三個數(shù)學表達式相加饲宿，得到：
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代碼支竹，可以直接使用（附圖為部分代碼示意）旋廷，完整的代碼只分享在知識星球——“達爾文生信星球”之中了。不過礼搁，我在截取圖片的時候饶碘，已經(jīng)將計算dosage的核心代碼包含在內了，如果此刻你覺得還不需要加入我的知識星球馒吴，那么也可以參考這一段代碼去實現(xiàn)你的程序扎运。
如果喜歡更多的生物信息和組學文章，搜索并關注我的微信公眾號“堿基礦工”（ID: helixminer）
你還可以讀
這是我的知識星球：『達爾文生信星球』（原名：解螺旋技術交流圈）饮戳，是一個我與讀者朋友們的私人朋友圈豪治，如今已有超過400人在星球中一起學習和交流。我有9年前沿而完整的生物信息學扯罐、NGS領域的科研經(jīng)歷负拟，在該領域發(fā)有多篇Nature、Cell級別的科學文章歹河，我希望借助這個知識星球可以與更多的志同道合者溝通和交流掩浙，同時也把自己的一些微薄經(jīng)驗分享給更多對組學感興趣的伙伴們。
這是知識星球上第一個與基因組學和生物信息學強相關的圈子秸歧，也是官方評定的優(yōu)秀星球。如今已經(jīng)累計超過1100個主題键菱，希望能夠借此營造一個高質量的組學知識圈和人脈圈，通過提問、彼此分享芭梯、交流經(jīng)驗、心得等玖喘，促進彼此更好地學習生信知識累奈，共同提升基因組數(shù)據(jù)分析和解讀的能力急但。
在這里你可以結識到全國優(yōu)秀的基因組學和生物信息學專家，同時可以分享你的經(jīng)驗波桩、見解和思考，有問題也可以向我提問和星球里的星友們提問镐躲。
