QualByDepth（QD）
QD是變異質量值（Quality）除以覆蓋深度（Depth）得到的比值。這里的變異質量值就是VCF中QUAL的值——用來衡量變異的可靠程度，這里的覆蓋深度是這個位點上所有 含有變異堿基的樣本的覆蓋深度之和，通俗一點說，就是這個值可以通過累加每個含變異的樣本（GT為非0/0的樣本）的覆蓋深度（VCF中每個樣本里面的DP）而得到。舉個例子：
1 1429249 . C T 1044.77 . . GT:AD:DP:GQ:PL 0/1:48,15:63:99:311,0,1644 0/0:47,0:47:99:392,0,0 1/1:0,76:76:99:3010,228,0
這個位點是1:1429249，VCF格式，但我把FILTER和INFO的信息省略了，它的變異質量值QUAL=1044.77。我們可以從中看到一共有三個樣本，其中一個是雜合變異（GT=0/1），一個純合的非變異（GT=0/0），最后一個是純合的變異（GT=1/1）。每個樣本的覆蓋深度都在其各自的DP域上，分別是：63，47和76。按照定義，這個位點的QD值就應該等于質量值除以另外兩個含有變異的樣本的深度之和（排除中間GT=0/0這個不含變異的樣本），也就是:
QD = 1044.77 / (63+76) = 7.516
QD這個值描述的實際上就是單位深度的變異質量值，也可以理解為是對變異質量值的一個歸一化，QD越高一般來說變異的可信度也越高。在質控的時候，相比于QUAL或者DP（深度）來說，QD是一個更加合理的值。因為我們知道，原始的變異質量值實際上與覆蓋的read數目是密切相關的，深度越高的位點QUAL一般都是越高的，而任何一個測序數據，都不可避免地會存在局部深度不均的情況，如果直接使用QUAL或者DP都會很容易因為覆蓋深度的差異而帶來有偏的質控結果。
在上面『執(zhí)行硬過濾』的例子里面，我們看到認為好的SNP變異，QD的值不能低于2，但 問題是為什么是2，而不是3或者其它數值呢？
要回答這個問題，我們可以通過利用NA12878 VQSR質控之后的變異數據和原始的變異數據來進行比較，并把它說明白。
首先，我們可以先把所有變異位點的QD值都提取出來，然后畫一個密度分布圖（Y軸代表的是對應QD值所占總數的比例，而不是個數），看看QD值的總體分布情況（如下圖，來自NA12878的數據）。
從這個圖里，我們可以看到QD的范圍主要集中在0~40之間。同時，我們可以明顯地看到有兩個峰值（QD=12和QD=32）。這兩個峰所反映的恰恰是雜合變異和純合變異的QD值所集中的地方。這里大家可以思考一下，哪一個是代表雜合變異的峰，哪一個是代表純合變異的峰呢？
回答是，第一個峰（QD=12）代表雜合，而第二峰（QD=32）代表純合，為什么呢？因為對于純合變異來說，貢獻于質量值的read是雜合變異的兩倍，同樣深度的情況下，QD會更大。對于大多數的高深度測序數據來說，QD的分布和上面的情況差不多，因此這個分布具有一定的代表性。
接著，我們同時畫出VQSR之后所有 可信變異（FILTER=Pass）和 不可信變異的QD分布圖，如下，淺綠色代表可信變異的QD分布圖，淺紅色代表不可信變異的QD分布圖。
你可以看到，大多數Fail的變異，都集中在左邊的低QD區(qū)域，而且紅波峰恰好是QD=2的地方，這就是為什么硬過濾時設置QD>2的原因了。
可是在上面的圖里，我想你也看到了，有很多Fail的變異它的QD還是很高的，有些甚至高于30，通過這樣的硬過濾參數所得到的結果中就會包含這部分本該要過濾掉的壞變異；而同樣的，在低QD（<2）區(qū)域其實也有一些是好的變異，但是卻被過濾掉了。這其實也是硬過濾的一大弊端，它不會像VQSR那樣，通過多個不同維度的數據構造合適的高維分類模型，從而能夠更準確地區(qū)分好和壞的變異，而僅僅是一刀切。
當你理解了上面有關QD的計算和閾值選擇的過程之后，要弄懂后面的指標和閾值也就容易了，因為用的也都是同樣的思路。
FisherStrand（FS）
FS是一個通過Fisher檢驗的p-value轉換而來的值，它要描述的是測序或者比對時對于只含有變異的read以及只含有參考序列堿基的read是否存在著明顯的正負鏈特異性（Strand bias，或者說是差異性）。這個差異反應了測序過程不夠隨機，或者是比對算法在基因組的某些區(qū)域存在一定的選擇偏向。如果測序過程是隨機的，比對是沒問題的，那么不管read是否含有變異，以及是否來自基因組的正鏈或者負鏈，只要是真實的它們就都應該是比較均勻的，也就是說，不會出現(xiàn)鏈特異的比對結果，F(xiàn)S應該接近于零。
與QD一樣，我們先來看一下質控前所有變異的FS總體密度分布圖（如下）。很明顯與QD相比，F(xiàn)S的范圍更加的大，從0到好幾百的都有。不過從圖中也可以看出，絕大部分的變異還是在100以下的。這里多說一句，在VCF的INFO中有時除了FS之外，有時你還會看到SB或者SOR。它們實際上是從不同的層面對鏈特異的現(xiàn)象進行描述。只不過SB給出的是原始各鏈比對數目，而FS則是對這些數據做了精確Fisher檢驗；SOR原理和FS類似，但是用了不同的統(tǒng)計檢驗方法計算差異程度，它更適合于高覆蓋度的數據。
下面這一個圖則是經過VQSR之后，畫出來的FS分布圖。跟上面的QD一樣，淺綠色代表好變異，淺紅色代表壞變異。我們可以看到，大部分好變異的FS都集中在0~10之間，而且壞變異的峰值在60左右的位置上，因此過濾的時候，我們把FS設置為大于60。其實設置這么高的一個閾值是比較激進的（留下很多假變異），但是從圖中你也可以看到，不過設置得多低，我們總會保留下很多假的變異，既然如此我們就干脆選擇盡可能保留更多好的變異，然后祈禱可以通過『執(zhí)行硬過濾』里其他的閾值來過濾掉那些無法通過FS過濾的假變異。
StrandOddsRatio（SOR）
關于SOR在上面講到FS的時候，我就在注釋里提及過了。它同樣是對鏈特異（Strand bias）的一種描述，但是從上面我們也可以看到FS在硬過濾的時候并不是非常給力，而且由于很多時候read在外顯子區(qū)域末端的覆蓋存在著一定的鏈特異（這個區(qū)域的現(xiàn)象其實是正常的），往往只有一個方向的read，這個時候該區(qū)域中如果有變異位點的話，那么FS通常會給出很差的分值，這時SOR就能夠起到比較好的校正作用了。計算SOR所用的統(tǒng)計檢驗方法也與FS不同，它用的是symmetric odds ratio test，數據是一個2×2的列聯(lián)表（如下），公式也十分簡單，我把公式進行了簡單的展開，從中可以清楚地看出，它考慮的其實就是ALT和REF這兩個堿基的read覆蓋方向的比例是否有偏，如果完全無偏，那么應該等于1。
sor = (ref_fwd/ref_rev) / (alt_fwd/alt_rev) = (ref_fwd * alt_rev) / (ref_rev * alt_fwd)
OK，那么同樣的，我們先看一下這個值總體的密度分布情況（如下）?？偟膩碚f，這個分布明顯集中在0~3之間，這也和我們的預期比較一致。不過也有比較明顯的長尾現(xiàn)象，這個時候我們也沒必要定下太過明確的閾值預期，先看VQSR的分布結果。
下面這個圖就是在VQSR之后，區(qū)分了好和壞變異之后，SOR的密度分布。很明顯，好的變異基本就在1附近。結合這個分布圖，我們在上面的例子里把它的閾值定為3基本上也不會過損失好的變異了，雖然單靠這個閾值還是會保留下不少假的變異，但是至少不合理的長尾部分可以被砍掉。
RMSMappingQuality（MQ）
MQ這個值是所有比對至該位點上的read的比對質量值的均方根（先平方、再平均、然后開方，如下公式）。
它和平均值相比更能夠準確地描述比對質量值的離散程度。而且有意思的是，如果我們的比對工具用的是bwa mem，那么按照它的算法，對于一個好的變異位點，我們可以預期，它的MQ值將等于60。
下面是所有未過濾的變異位點上MQ的密度分布圖?；旧暇椭辉?0的地方存在一個很瘦很高的峰。可以說這是目前為止這幾個指標中圖形最為規(guī)則的了，在這個圖上，我們甚至就可以直接定出MQ的閾值了，比如所有小于50的就可以過濾掉了。
但是，理性告訴我們還是要看一下VQSR的對比結果（下圖）。
你會發(fā)現(xiàn)似乎所有好的變異都緊緊集中在60旁邊了，其它地方就都是假的變異了，所以MQ的閾值設置為50也是合理的。但是同樣要注意到的地方是，60這個范圍實際上依然有假的變異位點在那里，我們把這個區(qū)域放大來看，如下圖，這里你就會發(fā)現(xiàn)其實假變異的密度分布圖也覆蓋到60這個范圍了。
考慮到篇幅的問題，接下來MappingQualityRankSumTest（MQRankSum）和ReadPOSRankSumTest（ReadPOSRankSum）的閾值設定原理，我不打算再細說下去了 ，思路和上面的4個是完全一樣的。都是通過比較VQSR之后的密度分布圖，最后確定了硬過濾的閾值。
但請不要以為這只是適用于GATK得到的變異，實際上，只要我們弄懂了這些指標選擇的原因和過濾的思路，那么通過任何其他的變異檢測工具也是依舊可以適用的，區(qū)別就在于GATK幫我們把這些要用的指標算好了。
同樣地，這些指標也不是一成不變的，可以根據實際的情況換成其他，或者我們自己重新計算。
Ti/Tv處于合理的范圍
Ti/Tv的值是物種在與自然相互作用和演化過程中在基因組上留下來的一個統(tǒng)計標記，在物種中這個值具有一定的穩(wěn)定性。因此，一般來說，在完成了以上的質控之后，還會看一下這些變異位點Ti/Tv的值是多少，以此來進一步確定結果的可靠程度。
Ti（Transition）指的是嘌呤轉嘌呤，或者嘧啶轉嘧啶的變異位點數目，即A<->G或C<->T；
Tv（Transversion）指的則是嘌呤和嘧啶互轉的變異位點數目，即A<->C，A<->T，G<->C和G<->T。（如下圖）
另外，在哺乳動物基因組上C->T的轉換比較多，這是因為基因組上的胞嘧啶C在甲基化的修飾下容易發(fā)生C->T的轉變。
說了這么多，Ti/Tv的比值應該是多少才是正常的呢？如果沒有 選擇壓力的存在，Ti/Tv將等于0.5，因為從概率上講Tv將是Ti的兩倍。但現(xiàn)實當然不是這樣的，比如對于人來說，全基因組正常的Ti/Tv在2.1左右，而外顯子區(qū)域是3.0左右，新發(fā)的變異（Novel variants）則在1.5左右。
最后多說一句，Ti/Tv是一個被動指標，它是對最后質控結果的一個反應，我們是不能夠在一開始的時候使用這個值來進行變異過濾的。
作者：黃樹嘉
鏈接：http://www.itdecent.cn/p/ff8204ae7ebf
來源：簡書
簡書著作權歸作者所有，任何形式的轉載都請聯(lián)繫作者獲得授權並註明出處。
