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

QualByDepth（QD）

QD是變異質(zhì)量值（Quality）除以覆蓋深度（Depth）得到的比值。這里的變異質(zhì)量值就是VCF中QUAL的值——用來(lái)衡量變異的可靠程度人芽，這里的覆蓋深度是這個(gè)位點(diǎn)上所有 含有變異堿基的樣本的覆蓋深度之和智什，通俗一點(diǎn)說(shuō)婚度，就是這個(gè)值可以通過(guò)累加每個(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越高一般來(lái)說(shuō)變異的可信度也越高鸟蟹。在質(zhì)控的時(shí)候乌妙，相比于QUAL或者DP（深度）來(lái)說(shuō)，QD是一個(gè)更加合理的值建钥。因?yàn)槲覀冎捞僭希嫉淖儺愘|(zhì)量值實(shí)際上與覆蓋的read數(shù)目是密切相關(guān)的，深度越高的位點(diǎn)QUAL一般都是越高的熊经，而任何一個(gè)測(cè)序數(shù)據(jù)泽艘，都不可避免地會(huì)存在局部深度不均的情況，如果直接使用QUAL或者DP都會(huì)很容易因?yàn)楦采w深度的差異而帶來(lái)有偏的質(zhì)控結(jié)果镐依。

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

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

QD密度分布

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

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

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

VQSR QD密度分布

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

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

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

FisherStrand（FS）

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

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

FS 密度分布

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

VQSR 密度分布

StrandOddsRatio（SOR）

關(guān)于SOR在上面講到FS的時(shí)候涮俄，我就在注釋里提及過(guò)了。它同樣是對(duì)鏈特異（Strand bias）的一種描述尸闸，但是從上面我們也可以看到FS在硬過(guò)濾的時(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)表（如下），公式也十分簡(jiǎn)單曼氛，我把公式進(jìn)行了簡(jiǎn)單的展開(kāi)豁辉，從中可以清楚地看出，它考慮的其實(shí)就是ALT和REF這兩個(gè)堿基的read覆蓋方向的比例是否有偏舀患，如果完全無(wú)偏，那么應(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è)值總體的密度分布情況（如下）现使〉统祝總的來(lái)說(shuō)，這個(gè)分布明顯集中在0~3之間碳锈，這也和我們的預(yù)期比較一致顽冶。不過(guò)也有比較明顯的長(zhǎng)尾現(xiàn)象，這個(gè)時(shí)候我們也沒(méi)必要定下太過(guò)明確的閾值預(yù)期售碳，先看VQSR的分布結(jié)果强重。

SOR 密度分布

下面這個(gè)圖就是在VQSR之后，區(qū)分了好和壞變異之后贸人，SOR的密度分布间景。很明顯，好的變異基本就在1附近艺智。結(jié)合這個(gè)分布圖倘要，我們?cè)谏厦娴睦永锇阉拈撝刀?基本上也不會(huì)過(guò)損失好的變異了，雖然單靠這個(gè)閾值還是會(huì)保留下不少假的變異十拣，但是至少不合理的長(zhǎng)尾部分可以被砍掉封拧。

VQSR SOR密度分布

RMSMappingQuality（MQ）

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

RMS 計(jì)算公式

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

下面是所有未過(guò)濾的變異位點(diǎn)上MQ的密度分布圖非竿⊥勺牛基本上就只在60的地方存在一個(gè)很瘦很高的峰『熘可以說(shuō)這是目前為止這幾個(gè)指標(biāo)中圖形最為規(guī)則的了承匣，在這個(gè)圖上，我們甚至就可以直接定出MQ的閾值了锤悄，比如所有小于50的就可以過(guò)濾掉了韧骗。

MQ 密度分布

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

VQSR 密度分布

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

VQSR 密度分布

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

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

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

Ti/Tv處于合理的范圍

Ti/Tv的值是物種在與自然相互作用和演化過(guò)程中在基因組上留下來(lái)的一個(gè)統(tǒng)計(jì)標(biāo)記缸濒，在物種中這個(gè)值具有一定的穩(wěn)定性足丢。因此粱腻，一般來(lái)說(shuō)，在完成了以上的質(zhì)控之后斩跌，還會(huì)看一下這些變異位點(diǎn)Ti/Tv的值是多少绍些，以此來(lái)進(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)變吏口。

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

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

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