轉(zhuǎn)自：
GATK4.0和全基因組數(shù)據(jù)分析實(shí)踐（下）

QualByDepth（QD）

QD是變異質(zhì)量值（Quality）除以覆蓋深度（Depth）得到的比值。這里的變異質(zhì)量值就是VCF中QUAL的值——用來衡量變異的可靠程度，這里的覆蓋深度是這個(gè)位點(diǎn)上所有 含有變異堿基的樣本的覆蓋深度之和，通俗一點(diǎn)說，就是這個(gè)值可以通過累加每個(gè)含變異的樣本（GT為非0/0的樣本）的覆蓋深度（VCF中每個(gè)樣本里面的DP）而得到。舉個(gè)例子：

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

這個(gè)位點(diǎn)是1:1429249，VCF格式，但我把FILTER和INFO的信息省略了，它的變異質(zhì)量值QUAL=1044.77。我們可以從中看到一共有三個(gè)樣本，其中一個(gè)是雜合變異（GT=0/1），一個(gè)純合的非變異（GT=0/0），最后一個(gè)是純合的變異（GT=1/1）。每個(gè)樣本的覆蓋深度都在其各自的DP域上，分別是：63，47和76。按照定義，這個(gè)位點(diǎn)的QD值就應(yīng)該等于質(zhì)量值除以另外兩個(gè)含有變異的樣本的深度之和（排除中間GT=0/0這個(gè)不含變異的樣本），也就是:

QD = 1044.77 / (63+76) = 7.516

QD這個(gè)值描述的實(shí)際上就是單位深度的變異質(zhì)量值，也可以理解為是對(duì)變異質(zhì)量值的一個(gè)歸一化，QD越高一般來說變異的可信度也越高。在質(zhì)控的時(shí)候，相比于QUAL或者DP（深度）來說，QD是一個(gè)更加合理的值。因?yàn)槲覀冎?，原始的變異質(zhì)量值實(shí)際上與覆蓋的read數(shù)目是密切相關(guān)的，深度越高的位點(diǎn)QUAL一般都是越高的，而任何一個(gè)測序數(shù)據(jù)，都不可避免地會(huì)存在局部深度不均的情況，如果直接使用QUAL或者DP都會(huì)很容易因?yàn)楦采w深度的差異而帶來有偏的質(zhì)控結(jié)果。

在上面『執(zhí)行硬過濾』的例子里面，我們看到認(rèn)為好的SNP變異，QD的值不能低于2，但 問題是為什么是2，而不是3或者其它數(shù)值呢？

要回答這個(gè)問題，我們可以通過利用NA12878 VQSR質(zhì)控之后的變異數(shù)據(jù)和原始的變異數(shù)據(jù)來進(jìn)行比較，并把它說明白。

QD密度分布

從這個(gè)圖里，我們可以看到QD的范圍主要集中在0~40之間。同時(shí)，我們可以明顯地看到有兩個(gè)峰值（QD=12和QD=32）。這兩個(gè)峰所反映的恰恰是雜合變異和純合變異的QD值所集中的地方。這里大家可以思考一下，哪一個(gè)是代表雜合變異的峰，哪一個(gè)是代表純合變異的峰呢？

回答是，第一個(gè)峰（QD=12）代表雜合，而第二峰（QD=32）代表純合，為什么呢？因?yàn)閷?duì)于純合變異來說，貢獻(xiàn)于質(zhì)量值的read是雜合變異的兩倍，同樣深度的情況下，QD會(huì)更大。對(duì)于大多數(shù)的高深度測序數(shù)據(jù)來說，QD的分布和上面的情況差不多，因此這個(gè)分布具有一定的代表性。

接著，我們同時(shí)畫出VQSR之后所有 可信變異（FILTER=Pass）和 不可信變異的QD分布圖，如下，淺綠色代表可信變異的QD分布圖，淺紅色代表不可信變異的QD分布圖。

VQSR QD密度分布

你可以看到，大多數(shù)Fail的變異，都集中在左邊的低QD區(qū)域，而且紅波峰恰好是QD=2的地方，這就是為什么硬過濾時(shí)設(shè)置QD>2的原因了。

可是在上面的圖里，我想你也看到了，有很多Fail的變異它的QD還是很高的，有些甚至高于30，通過這樣的硬過濾參數(shù)所得到的結(jié)果中就會(huì)包含這部分本該要過濾掉的壞變異；而同樣的，在低QD（<2）區(qū)域其實(shí)也有一些是好的變異，但是卻被過濾掉了。這其實(shí)也是硬過濾的一大弊端，它不會(huì)像VQSR那樣，通過多個(gè)不同維度的數(shù)據(jù)構(gòu)造合適的高維分類模型，從而能夠更準(zhǔn)確地區(qū)分好和壞的變異，而僅僅是一刀切。

當(dāng)你理解了上面有關(guān)QD的計(jì)算和閾值選擇的過程之后，要弄懂后面的指標(biāo)和閾值也就容易了，因?yàn)橛玫囊捕际峭瑯拥乃悸贰?/p>

FisherStrand（FS）

FS是一個(gè)通過Fisher檢驗(yàn)的p-value轉(zhuǎn)換而來的值，它要描述的是測序或者比對(duì)時(shí)對(duì)于只含有變異的read以及只含有參考序列堿基的read是否存在著明顯的正負(fù)鏈特異性（Strand bias，或者說是差異性）。這個(gè)差異反應(yīng)了測序過程不夠隨機(jī)，或者是比對(duì)算法在基因組的某些區(qū)域存在一定的選擇偏向。如果測序過程是隨機(jī)的，比對(duì)是沒問題的，那么不管read是否含有變異，以及是否來自基因組的正鏈或者負(fù)鏈，只要是真實(shí)的它們就都應(yīng)該是比較均勻的，也就是說，不會(huì)出現(xiàn)鏈特異的比對(duì)結(jié)果，F(xiàn)S應(yīng)該接近于零。

這里多說一句，在VCF的INFO中有時(shí)除了FS之外，有時(shí)你還會(huì)看到SB或者SOR。它們實(shí)際上是從不同的層面對(duì)鏈特異的現(xiàn)象進(jìn)行描述。只不過SB給出的是原始各鏈比對(duì)數(shù)目，而FS則是對(duì)這些數(shù)據(jù)做了精確Fisher檢驗(yàn)；SOR原理和FS類似，但是用了不同的統(tǒng)計(jì)檢驗(yàn)方法計(jì)算差異程度，它更適合于高覆蓋度的數(shù)據(jù)。

FS 密度分布

下面這一個(gè)圖則是經(jīng)過VQSR之后，畫出來的FS分布圖。跟上面的QD一樣，淺綠色代表好變異，淺紅色代表壞變異。我們可以看到，大部分好變異的FS都集中在0~10之間，而且壞變異的峰值在60左右的位置上，因此過濾的時(shí)候，我們把FS設(shè)置為大于60。其實(shí)設(shè)置這么高的一個(gè)閾值是比較激進(jìn)的（留下很多假變異），但是從圖中你也可以看到，不過設(shè)置得多低，我們總會(huì)保留下很多假的變異，既然如此我們就干脆選擇盡可能保留更多好的變異，然后祈禱可以通過『執(zhí)行硬過濾』里其他的閾值來過濾掉那些無法通過FS過濾的假變異。

VQSR 密度分布

StrandOddsRatio（SOR）

關(guān)于SOR在上面講到FS的時(shí)候，我就在注釋里提及過了。它同樣是對(duì)鏈特異（Strand bias）的一種描述，但是從上面我們也可以看到FS在硬過濾的時(shí)候并不是非常給力，而且由于很多時(shí)候read在外顯子區(qū)域末端的覆蓋存在著一定的鏈特異（這個(gè)區(qū)域的現(xiàn)象其實(shí)是正常的），往往只有一個(gè)方向的read，這個(gè)時(shí)候該區(qū)域中如果有變異位點(diǎn)的話，那么FS通常會(huì)給出很差的分值，這時(shí)SOR就能夠起到比較好的校正作用了。計(jì)算SOR所用的統(tǒng)計(jì)檢驗(yàn)方法也與FS不同，它用的是symmetric odds ratio test，數(shù)據(jù)是一個(gè)2×2的列聯(lián)表（如下），公式也十分簡單，我把公式進(jìn)行了簡單的展開，從中可以清楚地看出，它考慮的其實(shí)就是ALT和REF這兩個(gè)堿基的read覆蓋方向的比例是否有偏，如果完全無偏，那么應(yīng)該等于1。

SOR列聯(lián)表

sor = (ref_fwd/ref_rev) / (alt_fwd/alt_rev) = (ref_fwd * alt_rev) / (ref_rev * alt_fwd)

OK，那么同樣的，我們先看一下這個(gè)值總體的密度分布情況（如下）?？偟膩碚f，這個(gè)分布明顯集中在0~3之間，這也和我們的預(yù)期比較一致。不過也有比較明顯的長尾現(xiàn)象，這個(gè)時(shí)候我們也沒必要定下太過明確的閾值預(yù)期，先看VQSR的分布結(jié)果。

SOR 密度分布

下面這個(gè)圖就是在VQSR之后，區(qū)分了好和壞變異之后，SOR的密度分布。很明顯，好的變異基本就在1附近。結(jié)合這個(gè)分布圖，我們在上面的例子里把它的閾值定為3基本上也不會(huì)過損失好的變異了，雖然單靠這個(gè)閾值還是會(huì)保留下不少假的變異，但是至少不合理的長尾部分可以被砍掉。

VQSR SOR密度分布

RMSMappingQuality（MQ）

MQ這個(gè)值是所有比對(duì)至該位點(diǎn)上的read的比對(duì)質(zhì)量值的均方根（先平方、再平均、然后開方，如下公式）。

RMS 計(jì)算公式

它和平均值相比更能夠準(zhǔn)確地描述比對(duì)質(zhì)量值的離散程度。而且有意思的是，如果我們的比對(duì)工具用的是bwa mem，那么按照它的算法，對(duì)于一個(gè)好的變異位點(diǎn)，我們可以預(yù)期，它的MQ值將等于60。

下面是所有未過濾的變異位點(diǎn)上MQ的密度分布圖?；旧暇椭辉?0的地方存在一個(gè)很瘦很高的峰?？梢哉f這是目前為止這幾個(gè)指標(biāo)中圖形最為規(guī)則的了，在這個(gè)圖上，我們甚至就可以直接定出MQ的閾值了，比如所有小于50的就可以過濾掉了。

MQ 密度分布

但是，理性告訴我們還是要看一下VQSR的對(duì)比結(jié)果（下圖）。

VQSR 密度分布

你會(huì)發(fā)現(xiàn)似乎所有好的變異都緊緊集中在60旁邊了，其它地方就都是假的變異了，所以MQ的閾值設(shè)置為50也是合理的。但是同樣要注意到的地方是，60這個(gè)范圍實(shí)際上依然有假的變異位點(diǎn)在那里，我們把這個(gè)區(qū)域放大來看，如下圖，這里你就會(huì)發(fā)現(xiàn)其實(shí)假變異的密度分布圖也覆蓋到60這個(gè)范圍了。

VQSR 密度分布

考慮到篇幅的問題，接下來MappingQualityRankSumTest（MQRankSum）和ReadPOSRankSumTest（ReadPOSRankSum）的閾值設(shè)定原理，我不打算再細(xì)說下去了 ，思路和上面的4個(gè)是完全一樣的。都是通過比較VQSR之后的密度分布圖，最后確定了硬過濾的閾值。

但請(qǐng)不要以為這只是適用于GATK得到的變異，實(shí)際上，只要我們弄懂了這些指標(biāo)選擇的原因和過濾的思路，那么通過任何其他的變異檢測工具也是依舊可以適用的，區(qū)別就在于GATK幫我們把這些要用的指標(biāo)算好了。

同樣地，這些指標(biāo)也不是一成不變的，可以根據(jù)實(shí)際的情況換成其他，或者我們自己重新計(jì)算。

Ti/Tv處于合理的范圍

Ti/Tv的值是物種在與自然相互作用和演化過程中在基因組上留下來的一個(gè)統(tǒng)計(jì)標(biāo)記，在物種中這個(gè)值具有一定的穩(wěn)定性。因此，一般來說，在完成了以上的質(zhì)控之后，還會(huì)看一下這些變異位點(diǎn)Ti/Tv的值是多少，以此來進(jìn)一步確定結(jié)果的可靠程度。

Ti（Transition）指的是嘌呤轉(zhuǎn)嘌呤，或者嘧啶轉(zhuǎn)嘧啶的變異位點(diǎn)數(shù)目，即A<->G或C<->T；
Tv（Transversion）指的則是嘌呤和嘧啶互轉(zhuǎn)的變異位點(diǎn)數(shù)目，即A<->C，A<->T，G<->C和G<->T。（如下圖）

Ti和Tv轉(zhuǎn)換關(guān)系

另外，在哺乳動(dòng)物基因組上C->T的轉(zhuǎn)換比較多，這是因?yàn)榛蚪M上的胞嘧啶C在甲基化的修飾下容易發(fā)生C->T的轉(zhuǎn)變。

說了這么多，Ti/Tv的比值應(yīng)該是多少才是正常的呢？如果沒有 選擇壓力的存在，Ti/Tv將等于0.5，因?yàn)閺母怕噬现vTv將是Ti的兩倍。但現(xiàn)實(shí)當(dāng)然不是這樣的，比如對(duì)于人來說，全基因組正常的Ti/Tv在2.1左右，而外顯子區(qū)域是3.0左右，新發(fā)的變異（Novel variants）則在1.5左右。

最后多說一句，Ti/Tv是一個(gè)被動(dòng)指標(biāo)，它是對(duì)最后質(zhì)控結(jié)果的一個(gè)反應(yīng)，我們是不能夠在一開始的時(shí)候使用這個(gè)值來進(jìn)行變異過濾的。

作者：黃樹嘉
鏈接：http://m.itdecent.cn/p/ff8204ae7ebf
來源：簡書
簡書著作權(quán)歸作者所有，任何形式的轉(zhuǎn)載都請(qǐng)聯(lián)繫作者獲得授權(quán)並註明出處。