本文是閱讀一個(gè)做算法的文章，有些收獲和想法庭呜，作為一種記錄
我們?cè)谑褂胓atk call snp的時(shí)候呻袭，常常使用HaplotypeCaller這個(gè)參數(shù)，這個(gè)參數(shù)的目的是檢測(cè)出用來進(jìn)行snp及小indel變異
HaplotypeCaller這個(gè)參數(shù)使用的是預(yù)組裝的方法蜻韭，能提高變異檢測(cè)的準(zhǔn)確度，但是在某種程度上增加了資源的消耗和分析時(shí)長(zhǎng)柿扣。在分析時(shí)肖方，該模塊并不會(huì)在基因組范圍進(jìn)行全局的變異檢測(cè)，而是劃定高變區(qū)間檢測(cè)未状，以下是具體的步驟及原理圖俯画。
一般在重測(cè)序變異檢測(cè)中，需要測(cè)一定深度的reads司草，那么gatk HaplotypeCaller首先根據(jù)參考基因組過一遍艰垂，找到那些高變異的區(qū)間，然后對(duì)區(qū)間內(nèi)的數(shù)據(jù)及基因組進(jìn)行組裝并且預(yù)估單倍型埋虹，再就是根據(jù)該單倍型計(jì)算似然值猜憎，最后判定
1.定義區(qū)間
為了加速程序的處理速度，gatk并不會(huì)對(duì)全基因組范圍內(nèi)的所有位點(diǎn)進(jìn)行變異檢測(cè)搔课，只會(huì)對(duì)高變區(qū)進(jìn)行檢測(cè)胰柑，也就是所謂的Active Regions。此處定義的區(qū)間選擇，來源于比對(duì)結(jié)果柬讨，包括其中的錯(cuò)配崩瓤、插入、缺失以及soft clipping踩官。
計(jì)算時(shí)却桶，會(huì)計(jì)算每個(gè)位點(diǎn)出現(xiàn)突變的概率，這一步會(huì)同時(shí)考慮分型最大似然值以及雜合率蔗牡，之后通過高斯核做卷積來擴(kuò)大概率值肾扰。當(dāng)概率值超過0.002時(shí)，該區(qū)間就會(huì)被定義為高變區(qū)蛋逾。在計(jì)算過程中，程序會(huì)對(duì)毗鄰的高變區(qū)(100bp)進(jìn)行合并處理窗悯，設(shè)定的區(qū)域大小最小50bp区匣，最長(zhǎng)300bp，如果合并區(qū)超過300bp蒋院，會(huì)被截?cái)嘈纬?個(gè)或多個(gè)高變區(qū)亏钩，并且存在overlap的reads在相鄰的兩個(gè)高變區(qū)中都參與運(yùn)算。
2.數(shù)據(jù)組裝
高變區(qū)的reads以及對(duì)應(yīng)的基因組區(qū)域會(huì)被切割成存在overlap的小片段欺旧，若參考基因組對(duì)應(yīng)的片段集合中存在重復(fù)姑丑，會(huì)將短片段的長(zhǎng)度遞增直到?jīng)]有重復(fù)或者達(dá)到最大長(zhǎng)度限制(65nt)，若存在65nt的重復(fù)短序列辞友，組裝停止栅哀。序列切割完成后使用de-Bruijn-like graphs算法進(jìn)行組裝。組裝獲得的邊緣根據(jù)比對(duì)到的reads數(shù)目分配權(quán)重称龙，比對(duì)數(shù)量較少的邊會(huì)從圖形中刪除留拾，以此來對(duì)圖形進(jìn)行精簡(jiǎn)并形成單倍型數(shù)據(jù)。
組裝完成后會(huì)給每個(gè)單倍型生成CIGRA信息鲫尊，以此為依據(jù)進(jìn)行位點(diǎn)的分型痴柔。
這一過程類似于序列比對(duì)：
也就是說在第一步檢測(cè)高變異區(qū)間以后，gatk會(huì)根據(jù)測(cè)序reads的情況重組單倍型疫向，當(dāng)然根據(jù)位點(diǎn)的不同變異會(huì)組裝出許多單倍型咳蔚，那么重組好的單倍型相互之間進(jìn)行比對(duì)，從而篩選出它們之間的變異snp
最后搔驼，將這一些單倍型和參考基因組比對(duì)谈火，這樣就篩選出單倍型之間以及單倍型和參考基因組的snp變異信息了（圖中黃色表格）
3.隱馬爾科夫模型預(yù)估最大似然值
上一步我們說到根據(jù)不同reads的變異信息組裝出許多的單倍型，那么我們就要將這些reads重新比對(duì)會(huì)這些單倍型匙奴，并且構(gòu)成單倍型矩陣（某個(gè)位點(diǎn)的單倍型矩陣）
組裝完成后堆巧，將原有的reads回帖到獲得的單倍型序列，之后根據(jù)比對(duì)結(jié)果預(yù)估在不同單倍型中的最大似然值(likelihood)。具體原理是使用隱馬爾科夫模型(pair-HMM谍肤，paired Hidden Markov Model)對(duì)reads和單倍型之間進(jìn)行最大似然值的預(yù)估
利用馬爾可夫模型進(jìn)行序列比對(duì)
在序列比對(duì)中啦租，我們將Match，Insert和Delete看作為三個(gè)狀態(tài)荒揣，分別用M篷角，I 和D來表示（可參考http://www.reibang.com/writer#/notebooks/42286457/notes/65851344/preview）
比對(duì)所得到的似然值可以構(gòu)成了單倍型矩陣，具體展開為：
那么接下來系任，我們根據(jù)單倍型以及reads比對(duì)到每個(gè)單倍型的情況恳蹲，將該矩陣進(jìn)行轉(zhuǎn)換，轉(zhuǎn)換為突變基因的似然值：
其中俩滥，上面的矩陣為某位點(diǎn)的單倍型矩陣嘉蕾，橫坐標(biāo)為不同組裝的單倍型（1，2霜旧，3错忱，4，可視為不同的單倍型）挂据，縱坐標(biāo)為測(cè)序的reads（1以清，2，3代表reads 1崎逃，2掷倔，3），矩陣對(duì)應(yīng)的值就是每條reads比對(duì)到該位點(diǎn)的似然值个绍，我們選似然值較高的單倍型作為突變后的基因型
如果考慮突變成兩個(gè)堿基（即考慮等位基因突變）
我們通常將每一條reads比對(duì)到每種單倍型上勒葱，得到相應(yīng)的似然值，對(duì)于每一條reads來說選擇似然值最大的兩個(gè)單倍型
即將每一條reads在四個(gè)單倍型中最高的兩個(gè)似然值選出來構(gòu)成兩列的矩陣障贸，根據(jù)該矩陣中每條reads在該位點(diǎn)堿基出現(xiàn)的頻率來判斷突變的類型错森，比方說第一列中2號(hào)單倍型出現(xiàn)頻率最多（綠色一共出現(xiàn)2次），即定為此單倍型為第一備選單倍型篮洁；而第二列中3號(hào)單倍型出現(xiàn)頻率最多（暗橙色一共出現(xiàn)2次）涩维，選該單倍型為第二備選單倍型，而第一備選單倍型在相應(yīng)高突變位點(diǎn)為C堿基袁波，第二備選單倍型在相應(yīng)高突變位點(diǎn)為A堿基瓦阐，這樣就篩選出了兩個(gè)備選堿基組合
4.分型結(jié)果判定
對(duì)高突變區(qū)的重組單倍型完成后，接下來需要進(jìn)行分型檢驗(yàn)篷牌，即檢測(cè)所測(cè)物種突變成什么類型的堿基睡蟋，以及突變?yōu)橐粋€(gè)堿基還是兩個(gè)堿基
確定該物種突變所得的類型以后，與參考基因組對(duì)比枷颊，看看突變情況
上圖表示：
- A/C戳杀，0/0表示突變類型為A/A该面；0/1表示突變類型為A/C；1/1表示突變類型為C/C
- A/G信卡，0/0表示突變類型為A/A隔缀；0/1表示突變類型為A/G；1/1表示突變類型為G/G
其中0/0和1/1表示突變?yōu)?個(gè)堿基傍菇，0/1表示為突變成2個(gè)堿基
根據(jù)上述計(jì)算猾瘸，我們可以得到突變基因矩陣，并且定義如下關(guān)系式子：
其中G為先驗(yàn)的基因型信息,為已知信息,故P(G) = 1
P(G | R)為在測(cè)得的reads的信息條件下,基因型為G的概率,為后驗(yàn)信息
接下來我們根據(jù)等位基因矩陣計(jì)算每一種突變類型的概率值,其中G(C/C)為突變基因矩陣第一列的每個(gè)元素和自己相加再除以2,最后結(jié)果相乘;而G(C/A)為突變基因矩陣每一行的前后兩個(gè)元素相加,再除以2,最后結(jié)果相乘,以此類推,最后將這些基因型的似然值相加,得到0.002116
最后用貝葉斯公式計(jì)算每個(gè)堿基突變類型的概率值,選取最大的作為最合適的突變類型丢习，該例子顯然是C/C為最佳突變類型牵触，即突變?yōu)镃堿基
之后可以將這個(gè)基因型于參考基因組對(duì)比,看看是否發(fā)生變異
參考：http://www.reibang.com/p/d6964cdcf4c6
https://mp.weixin.qq.com/s/wt1DuVPN1qOpIolfwrgfsg
https://gatk.broadinstitute.org/hc/en-us/articles/360035890531?id=4442
