GATK4.0和全基因組數據分析實踐（下）

05-02

前言

在上一篇文章中我已經用例子仔細跟大家分享了WGS從原始數據到變異數據（Fastq->VCF）的具體執行過程。那麼，在這一篇文章里，我們就來好好談談後續非常重要的一個環節——也是本次實踐分析的最後一個部分——變異的質控。

什麼是質控？我們不妨先給它下個定義：質控的含義和目的是指通過一定的標準，最大可能地剔除假陽性的結果，並儘可能地保留最多的正確數據。有了這麼一個定義之後，我們也就能夠更加清晰地知道接下來該做些什麼了。

在上次的文章里我已經說到在GATK HaplotypeCaller之後，首選的質控方案是GATK VQSR，它通過機器學習的方法利用多個不同的數據特徵訓練一個模型（高斯混合模型）對變異數據進行質控，然而不幸的是使用VQSR需要具備以下兩個條件：

第一，需要一個精心準備的已知變異集，它將作為訓練質控模型的真集。比如，對於我們人來說，就有Hapmap、OMNI，1000G和dbsnp等這些國際性項目的數據，這些可以作為高質量的已知變異集。GATK的bundle主要就是對這四個數據集做了精心的處理和選擇，然後把它們作為VQSR時的真集位點。這裡我強調一個地方：是真集的『位點』而不是真集的『數據』！還請大家多多注意。因為，VQSR並不是用這些變異集里的數據來訓練的，而是用我們自己的變異數據。這個對於剛接WGS的同學來說特別容易搞混，不要因為VQSR中用了那四份變異數據，就以為是用它們的數據來訓練模型。

實際上，這些已知變異集的意義是告訴我們群體中哪些位點存在著變異，如果在其他人的數據里能觀察到落入這個集合中的變異位點，那麼這些被已知集包括的變異就有很大的可能是正確的。也就是說，我們可以從數據中篩選出那些和真集『位點』相同的變異，把它們當作是真實的變異結果。接著，進行VQSR的時候，程序就可以用這個篩選出來的數據作為真集數據來訓練，並構造模型了。

關於VQSR的內在原理，前不久在我的知識星球中我做過簡單的回答（下圖），這裡就不展開了，感興趣的同學看下圖的內容基本上也是足夠的，雖然不詳細，但應該可以幫你建立一個關於VQSR的基本認識。對於希望深入理解演算法細節的同學來說，我的建議是直接閱讀GATK這一部分的代碼。但要注意，類似這樣的過濾演算法實際上還可以用很多不同的機器學習演算法來解決，比如SVM，或者用深度學習來構造這個質控模型也都是OK的。

第二，要求新檢測的結果中有足夠多的變異，不然VQSR在進行模型訓練的時候會因為可用的變異位點數目不足而無法進行。

由於條件1的限制，會導致很多非人的物種在完成變異檢測之後沒法使用GATK VQSR的方法進行質控。而由於條件2，也常常導致一些小panel甚至外顯子測序，由於最後的變異位點不夠，也無法使用VQSR。這個時候，我們就不得不選擇硬過濾的方式來質控了。

那什麼叫做硬過濾呢？所謂硬過濾其實就是通過人為設定一個或者若干個指標閾值（也可以叫數據特徵值），然後把所有不滿足閾值的變異位點採用一刀切掉的方法。

那麼如何執行硬過濾？首先，需要我們確定該用哪些指標來評價變異的好壞。這個非常重要，選擇對了事半功倍，選得不合理，過濾的結果有時還不如不過濾的。如果把這個問題放在從前，我們需要做比較多的嘗試才能確定一些合適的指標，但現在就方便很多了，可以直接使用GATK VQSR所用的指標——畢竟這些指標都是經過精挑細選的。我想這應該不難理解，既然VQSR就是用這些指標來訓練質控模型的，那麼它們就可以在一定程度上描述每個變異的質量，我們用這些指標設置對應的閾值來進行硬過濾也將是合理的。VQSR使用的數據指標有6個（這些指標都在VCF文件的INFO域中，如果不是GATK得到的變異，可能會有所不同，但知道它們的含義之後也是可以自己計算的），分別是：

QualByDepth（QD）
FisherStrand (FS)
StrandOddsRatio (SOR)
RMSMappingQuality (MQ)
MappingQualityRankSumTest (MQRankSum)
ReadPosRankSumTest (ReadPosRankSum)

指標有了，那麼閾值應該設置為多少？下面我想先給出一個硬過濾的例子，然後再逐個來對其進行分析，以便大家能夠更好地理解變異質控的思路。值得注意的是不同的數據，有不同的情況，它的閾值有時是不同的。不過不用擔心，當你掌握了如何做的思路之後完全有能力根據具體的情況舉一反三。

執行硬過濾

首先是硬過濾的例子，這個過程我都用最新的GATK來完成。GATK 4.0中有一個專門的VariantFiltration模塊（繼承自GATK 3.x），它可以很方便地幫我們完成這個事情。不過，過濾的時候，需要分SNP和Indel這兩個不同的變異類型來進行，它們有些閾值是不同的，需要區別對待。在下面的例子里，我們還是用上一篇文章中最後得到的變異數據（E_coli_K12.vcf.gz）為例子，這是具體的執行命令：

# 使用SelectVariants，選出SNPtime /Tools/common/bin/gatk/4.0.1.2/gatk SelectVariants -select-type SNP -V ../output/E.coli/E_coli_K12.vcf.gz -O ../output/E.coli/E_coli_K12.snp.vcf.gz# 為SNP作硬過濾time /Tools/common/bin/gatk/4.0.1.2/gatk VariantFiltration -V ../output/E.coli/E_coli_K12.snp.vcf.gz --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "Filter" -O ../output/E.coli/E_coli_K12.snp.filter.vcf.gz# 使用SelectVariants，選出Indeltime /Tools/common/bin/gatk/4.0.1.2/gatk SelectVariants -select-type INDEL -V ../output/E.coli/E_coli_K12.vcf.gz -O ../output/E.coli/E_coli_K12.indel.vcf.gz# 為Indel作過濾time /Tools/common/bin/gatk/4.0.1.2/gatk VariantFiltration -V ../output/E.coli/E_coli_K12.indel.vcf.gz --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "Filter" -O ../output/E.coli/E_coli_K12.indel.filter.vcf.gz# 重新合併過濾後的SNP和Indeltime /Tools/common/bin/gatk/4.0.1.2/gatk MergeVcfs -I ../output/E.coli/E_coli_K12.snp.filter.vcf.gz -I ../output/E.coli/E_coli_K12.indel.filter.vcf.gz -O ../output/E.coli/E_coli_K12.filter.vcf.gz# 刪除無用中間文件rm -f ../output/E.coli/E_coli_K12.snp.vcf.gz* ../output/E.coli/E_coli_K12.snp.filter.vcf.gz* ../output/E.coli/E_coli_K12.indel.vcf.gz* ../output/E.coli/E_coli_K12.indel.filter.vcf.gz*

與上一篇文章的目錄邏輯一樣，我們把這些shell命令都寫入到bin目錄下一個名為「variant_filtration.sh」的文件中，然後運行它。最後，只要符合了上面任意一個閾值的變異都會被設置為「Filter」，剩下的會被認為是正常的變異，並標記為「PASS」。流程的最後，我們需要把分開質控的SNP和Indel結果重新合併在一起，然後再把那些不必要的中間文件刪除掉。

在具體的項目中，你如果需要使用硬過濾的策略，這個例子中的參數可以作為參考，特別是對於高深度數據而言。接下來我結合GATK所提供的資料與大家分享如何理解這些指標以及得出這些閾值的思路。

如何理解硬過濾的指標和閾值的計算

考慮到SNP和Indel在判斷指標和閾值方面的思路是一致的，因此沒必要重複說。所以，下面我只以SNP為例子，告訴大家設定閾值的思路。強調一下，為了更具有通用價值，這些閾值是借用NA12878（來自GIAB）的高深度數據進行計算獲得的，所以如果你的數據（或者物種）相對比較特殊（不是哺乳動物），那麼就不建議直接套用了，但可以依照類似的思路去尋找新閾值。

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都會很容易因為覆蓋深度的差異而帶來有偏的質控結果。

在上面『執行硬過濾』的例子裡面，我們看到認為好的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區域，而且紅波峰恰好是QD=2的地方，這就是為什麼硬過濾時設置QD>2的原因了。

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

當你理解了上面有關QD的計算和閾值選擇的過程之後，要弄懂後面的指標和閾值也就容易了，因為用的也都是同樣的思路。

FisherStrand（FS）

FS是一個通過Fisher檢驗的p-value轉換而來的值，它要描述的是測序或者比對時對於只含有變異的read以及只含有參考序列鹼基的read是否存在著明顯的正負鏈特異性（Strand bias，或者說是差異性）。這個差異反應了測序過程不夠隨機，或者是比對演算法在基因組的某些區域存在一定的選擇偏向。如果測序過程是隨機的，比對是沒問題的，那麼不管read是否含有變異，以及是否來自基因組的正鏈或者負鏈，只要是真實的它們就都應該是比較均勻的，也就是說，不會出現鏈特異的比對結果，FS應該接近於零。

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

與QD一樣，我們先來看一下質控前所有變異的FS總體密度分布圖（如下）。很明顯與QD相比，FS的範圍更加的大，從0到好幾百的都有。不過從圖中也可以看出，絕大部分的變異還是在100以下的。

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

StrandOddsRatio（SOR）

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

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

OK，那麼同樣的，我們先看一下這個值總體的密度分布情況（如下）。總的來說，這個分布明顯集中在0~3之間，這也和我們的預期比較一致。不過也有比較明顯的長尾現象，這個時候我們也沒必要定下太過明確的閾值預期，先看VQSR的分布結果。

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

RMSMappingQuality（MQ）

MQ這個值是所有比對至該位點上的read的比對質量值的均方根（先平方、再平均、然後開方，如下公式）。

它和平均值相比更能夠準確地描述比對質量值的離散程度。而且有意思的是，如果我們的比對工具用的是bwa mem，那麼按照它的演算法，對於一個好的變異位點，我們可以預期，它的MQ值將等於60。

下面是所有未過濾的變異位點上MQ的密度分布圖。基本上就只在60的地方存在一個很瘦很高的峰。可以說這是目前為止這幾個指標中圖形最為規則的了，在這個圖上，我們甚至就可以直接定出MQ的閾值了，比如所有小於50的就可以過濾掉了。

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

你會發現似乎所有好的變異都緊緊集中在60旁邊了，其它地方就都是假的變異了，所以MQ的閾值設置為50也是合理的。但是同樣要注意到的地方是，60這個範圍實際上依然有假的變異位點在那裡，我們把這個區域放大來看，如下圖，這裡你就會發現其實假變異的密度分布圖也覆蓋到60這個範圍了。

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

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

同樣地，這些指標也不是一成不變的，可以根據實際的情況換成其他，或者我們自己重新計算。

Ti/Tv處於合理的範圍

Ti/Tv的值是物種在與自然相互作用和演化過程中在基因組上留下來的一個統計標記，在物種中這個值具有一定的穩定性。因此，一般來說，在完成了以上的質控之後，還會看一下這些變異位點Ti/Tv的值是多少，以此來進一步確定結果的可靠程度。

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

另外，在哺乳動物基因組上C->T的轉換比較多，這是因為基因組上的胞嘧啶C在甲基化的修飾下容易發生C->T的轉變。

說了這麼多Ti/Tv的比值應該是多少才是正常的呢？如果沒有選擇壓力的存在，Ti/Tv將等於0.5，因為從概率上講Tv將是Ti的兩倍。但現實當然不是這樣的，比如對於人來說，全基因組正常的Ti/Tv在2.1左右，而外顯子區域是3.0左右，新發的變異（Novel variants）則在1.5左右。

最後多說一句，Ti/Tv是一個被動指標，它是對最後質控結果的一個反應，我們是不能夠在一開始的時候使用這個值來進行變異過濾的。

小結

雖然本文一直在談論的是如何做好硬過濾，但不管我們的指標和閾值設置的多麼完美，硬過濾的硬傷都是一刀切，它並不會根據多個維度的數據自動設計出更加合理的過濾模式。硬過濾作為一個簡單粗暴的方法，不到不得已的時候不推薦使用，即使使用了，也一定要明白每一個過濾指標和閾值都意味著什麼。

最後，我也希望通過這一篇文章能夠完整地為你呈現一個變異質控的思路。

我的微信公眾號：解螺旋的礦工 歡迎關注更及時了解更多信息。