轉(zhuǎn)載
GATK4.0和全基因組數(shù)據(jù)分析實(shí)踐（下）
原創(chuàng)：?礦工?堿基礦工?2018-03-23
前言
在上一篇文章中我已經(jīng)用例子仔細(xì)跟大家分享了WGS從原始數(shù)據(jù)到變異數(shù)據(jù)（Fastq->VCF）的具體執(zhí)行過程。那么，在這一篇文章里，我們就來好好談?wù)労罄m(xù)非常重要的一個(gè)環(huán)節(jié)——也是本次實(shí)踐分析的最后一個(gè)部分——變異的質(zhì)控。
什么是質(zhì)控？我們不妨先給它下個(gè)定義：質(zhì)控的含義和目的是指通過一定的標(biāo)準(zhǔn)，最大可能地剔除假陽性的結(jié)果，并盡可能地保留最多的正確數(shù)據(jù)。有了這么一個(gè)定義之后，我們也就能夠更加清晰地知道接下來該做些什么了。
在上次的文章里我已經(jīng)說到在GATK HaplotypeCaller之后，首選的質(zhì)控方案是GATK VQSR，它通過機(jī)器學(xué)習(xí)的方法利用多個(gè)不同的數(shù)據(jù)特征訓(xùn)練一個(gè)模型（高斯混合模型）對(duì)變異數(shù)據(jù)進(jìn)行質(zhì)控，然而不幸的是使用VQSR需要具備以下兩個(gè)條件：
第一，需要一個(gè)精心準(zhǔn)備的已知變異集，它將作為訓(xùn)練質(zhì)控模型的真集。比如，對(duì)于我們?nèi)藖碚f，就有Hapmap、OMNI，1000G和dbsnp等這些國(guó)際性項(xiàng)目的數(shù)據(jù)，這些可以作為高質(zhì)量的已知變異集。GATK的bundle主要就是對(duì)這四個(gè)數(shù)據(jù)集做了精心的處理和選擇，然后把它們作為VQSR時(shí)的真集位點(diǎn)。這里我強(qiáng)調(diào)一個(gè)地方：是真集的『位點(diǎn)』而不是真集的『數(shù)據(jù)』！還請(qǐng)大家多多注意。因?yàn)椋?b>VQSR并不是用這些變異集里的數(shù)據(jù)來訓(xùn)練的，而是用我們自己的變異數(shù)據(jù)。這個(gè)對(duì)于剛接WGS的同學(xué)來說特別容易搞混，不要因?yàn)閂QSR中用了那四份變異數(shù)據(jù)，就以為是用它們的數(shù)據(jù)來訓(xùn)練模型。
實(shí)際上，這些已知變異集的意義是告訴我們?nèi)后w中哪些位點(diǎn)存在著變異，如果在其他人的數(shù)據(jù)里能觀察到落入這個(gè)集合中的變異位點(diǎn)，那么這些被已知集包括的變異就有很大的可能是正確的。也就是說，我們可以從數(shù)據(jù)中篩選出那些和真集『位點(diǎn)』相同的變異，把它們當(dāng)作是真實(shí)的變異結(jié)果。接著，進(jìn)行VQSR的時(shí)候，程序就可以用這個(gè)篩選出來的數(shù)據(jù)作為真集數(shù)據(jù)來訓(xùn)練，并構(gòu)造模型了。
關(guān)于VQSR的內(nèi)在原理，前不久在我的知識(shí)星球中我做過簡(jiǎn)單的回答（下圖），這里就不展開了，感興趣的同學(xué)看下圖的內(nèi)容基本上也是足夠的，雖然不詳細(xì)，但應(yīng)該可以幫你建立一個(gè)關(guān)于VQSR的基本認(rèn)識(shí)。對(duì)于希望深入理解算法細(xì)節(jié)的同學(xué)來說，我的建議是直接閱讀GATK這一部分的代碼。但要注意，類似這樣的過濾算法實(shí)際上還可以用很多不同的機(jī)器學(xué)習(xí)算法來解決，比如SVM，或者用深度學(xué)習(xí)來構(gòu)造這個(gè)質(zhì)控模型也都是OK的。
第二，要求新檢測(cè)的結(jié)果中有足夠多的變異，不然VQSR在進(jìn)行模型訓(xùn)練的時(shí)候會(huì)因?yàn)榭捎玫淖儺愇稽c(diǎn)數(shù)目不足而無法進(jìn)行。
由于條件1的限制，會(huì)導(dǎo)致很多非人的物種在完成變異檢測(cè)之后沒法使用GATK VQSR的方法進(jìn)行質(zhì)控。而由于條件2，也常常導(dǎo)致一些小panel甚至外顯子測(cè)序，由于最后的變異位點(diǎn)不夠，也無法使用VQSR。這個(gè)時(shí)候，我們就不得不選擇硬過濾的方式來質(zhì)控了。
那什么叫做硬過濾呢？所謂硬過濾其實(shí)就是通過人為設(shè)定一個(gè)或者若干個(gè)指標(biāo)閾值（也可以叫數(shù)據(jù)特征值），然后把所有不滿足閾值的變異位點(diǎn)采用一刀切掉的方法。
那么如何執(zhí)行硬過濾？首先，需要我們確定該用哪些指標(biāo)來評(píng)價(jià)變異的好壞。這個(gè)非常重要，選擇對(duì)了事半功倍，選得不合理，過濾的結(jié)果有時(shí)還不如不過濾的。如果把這個(gè)問題放在從前，我們需要做比較多的嘗試才能確定一些合適的指標(biāo)，但現(xiàn)在就方便很多了，可以直接使用GATK VQSR所用的指標(biāo)——畢竟這些指標(biāo)都是經(jīng)過精挑細(xì)選的。我想這應(yīng)該不難理解，既然VQSR就是用這些指標(biāo)來訓(xùn)練質(zhì)控模型的，那么它們就可以在一定程度上描述每個(gè)變異的質(zhì)量，我們用這些指標(biāo)設(shè)置對(duì)應(yīng)的閾值來進(jìn)行硬過濾也將是合理的。VQSR使用的數(shù)據(jù)指標(biāo)有6個(gè)（這些指標(biāo)都在VCF文件的INFO域中，如果不是GATK得到的變異，可能會(huì)有所不同，但知道它們的含義之后也是可以自己計(jì)算的），分別是：
QualByDepth（QD）
FisherStrand (FS)
StrandOddsRatio (SOR)
RMSMappingQuality (MQ)
MappingQualityRankSumTest (MQRankSum)
ReadPosRankSumTest (ReadPosRankSum)
指標(biāo)有了，那么閾值應(yīng)該設(shè)置為多少？下面我想先給出一個(gè)硬過濾的例子，然后再逐個(gè)來對(duì)其進(jìn)行分析，以便大家能夠更好地理解變異質(zhì)控的思路。值得注意的是不同的數(shù)據(jù)，有不同的情況，它的閾值有時(shí)是不同的。不過不用擔(dān)心，當(dāng)你掌握了如何做的思路之后完全有能力根據(jù)具體的情況舉一反三。
執(zhí)行硬過濾
首先是硬過濾的例子，這個(gè)過程我都用最新的GATK來完成。GATK 4.0中有一個(gè)專門的VariantFiltration模塊（繼承自GATK 3.x），它可以很方便地幫我們完成這個(gè)事情。不過，過濾的時(shí)候，需要分SNP和Indel這兩個(gè)不同的變異類型來進(jìn)行，它們有些閾值是不同的，需要區(qū)別對(duì)待。在下面的例子里，我們還是用上一篇文章中最后得到的變異數(shù)據(jù)（E_coli_K12.vcf.gz）為例子，這是具體的執(zhí)行命令：
# 使用SelectVariants，選出SNP
time?/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?"PASS"?\
-O ../output/E.coli/E_coli_K12.snp.filter.vcf.gz
# 使用SelectVariants，選出Indel
time /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?"PASS"?\
-O ../output/E.coli/E_coli_K12.indel.filter.vcf.gz
# 重新合并過濾后的SNP和Indel
time /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目錄下一個(gè)名為“variant_filtration.sh”的文件中，然后運(yùn)行它。最后，只要符合了上面任意一個(gè)閾值的變異都會(huì)被留下來。流程的最后，我們需要把分開質(zhì)控的SNP和Indel結(jié)果重新合并在一起，然后再把那些不必要的中間文件刪除掉。
在具體的項(xiàng)目中，你如果需要使用硬過濾的策略，這個(gè)例子中的參數(shù)可以作為參考，特別是對(duì)于高深度數(shù)據(jù)而言。接下來我結(jié)合GATK所提供的資料與大家分享如何理解這些指標(biāo)以及得出這些閾值的思路。
如何理解硬過濾的指標(biāo)和閾值的計(jì)算
考慮到SNP和Indel在判斷指標(biāo)和閾值方面的思路是一致的，因此沒必要重復(fù)說。所以，下面我只以SNP為例子，告訴大家設(shè)定閾值的思路。強(qiáng)調(diào)一下，為了更具有通用價(jià)值，這些閾值是借用NA12878（來自GIAB）的高深度數(shù)據(jù)進(jìn)行計(jì)算獲得的，所以如果你的數(shù)據(jù)（或者物種）相對(duì)比較特殊（不是哺乳動(dòng)物），那么就不建議直接套用了，但可以依照類似的思路去尋找新閾值。
QualByDepth（QD）
QD是變異質(zhì)量值（Quality）除以覆蓋深度（Depth）得到的比值。這里的變異質(zhì)量值就是VCF中QUAL的值——用來衡量變異的可靠程度，這里的覆蓋深度是這個(gè)位點(diǎn)上所有含有變異堿基的樣本的覆蓋深度之和，通俗一點(diǎn)說，就是這個(gè)值可以通過累加每個(gè)含變異的樣本（GT為非0/0的樣本）的覆蓋深度（VCF中每個(gè)樣本里面的DP）而得到。舉個(gè)例子：
11429249?.CT1044.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è)測(cè)序數(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)行比較，并把它說明白。
首先，我們可以先把所有變異位點(diǎn)的QD值都提取出來，然后畫一個(gè)密度分布圖（Y軸代表的是對(duì)應(yīng)QD值所占總數(shù)的比例，而不是個(gè)數(shù)），看看QD值的總體分布情況（如下圖，來自NA12878的數(shù)據(jù)）。
從這個(gè)圖里，我們可以看到QD的范圍主要集中在0~40之間。同時(shí)，我們可以明顯地看到有兩個(gè)峰值（QD=12和QD=32）。這兩個(gè)峰所反映的恰恰是雜合變異和純合變異的QD值所集中的地方。這里大家可以思考一下，哪一個(gè)是代表雜合變異的峰，哪一個(gè)是代表純合變異的峰呢？
回答是，第一個(gè)峰（QD=12）代表雜合，而第二峰（QD=32）代表純合，為什么呢？因?yàn)?b>對(duì)于純合變異來說，貢獻(xiàn)于質(zhì)量值的read是雜合變異的兩倍，同樣深度的情況下，QD會(huì)更大。對(duì)于大多數(shù)的高深度測(cè)序數(shù)據(jù)來說，QD的分布和上面的情況差不多，因此這個(gè)分布具有一定的代表性。
接著，我們同時(shí)畫出VQSR之后所有可信變異（FILTER=Pass）和不可信變異的QD分布圖，如下，淺綠色代表可信變異的QD分布圖，淺紅色代表不可信變異的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)換而來的值，它要描述的是測(cè)序或者比對(duì)時(shí)對(duì)于只含有變異的read以及只含有參考序列堿基的read是否存在著明顯的正負(fù)鏈特異性（Strand bias，或者說是差異性）。這個(gè)差異反應(yīng)了測(cè)序過程不夠隨機(jī)，或者是比對(duì)算法在基因組的某些區(qū)域存在一定的選擇偏向。如果測(cè)序過程是隨機(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ù)。
與QD一樣，我們先來看一下質(zhì)控前所有變異的FS總體密度分布圖（如下）。很明顯與QD相比，F(xiàn)S的范圍更加的大，從0到好幾百的都有。不過從圖中也可以看出，絕大部分的變異還是在100以下的。
下面這一個(gè)圖則是經(jīng)過VQSR之后，畫出來的FS分布圖。跟上面的QD一樣，淺綠色代表好變異，淺紅色代表壞變異。我們可以看到，大部分好變異的FS都集中在0~10之間，而且壞變異的峰值在60左右的位置上，因此過濾的時(shí)候，我們把FS設(shè)置為大于60。其實(shí)設(shè)置這么高的一個(gè)閾值是比較激進(jìn)的（留下很多假變異），但是從圖中你也可以看到，不過設(shè)置得多低，我們總會(huì)保留下很多假的變異，既然如此我們就干脆選擇盡可能保留更多好的變異，然后祈禱可以通過『執(zhí)行硬過濾』里其他的閾值來過濾掉那些無法通過FS過濾的假變異。
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)表（如下），公式也十分簡(jiǎn)單，我把公式進(jìn)行了簡(jiǎn)單的展開，從中可以清楚地看出，它考慮的其實(shí)就是ALT和REF這兩個(gè)堿基的read覆蓋方向的比例是否有偏，如果完全無偏，那么應(yīng)該等于1。
sor?= (ref_fwd/ref_rev) / (alt_fwd/alt_rev) = (ref_fwd * alt_rev) / (ref_rev * alt_fwd)
OK，那么同樣的，我們先看一下這個(gè)值總體的密度分布情況（如下）?？偟膩碚f，這個(gè)分布明顯集中在0~3之間，這也和我們的預(yù)期比較一致。不過也有比較明顯的長(zhǎng)尾現(xiàn)象，這個(gè)時(shí)候我們也沒必要定下太過明確的閾值預(yù)期，先看VQSR的分布結(jié)果。
下面這個(gè)圖就是在VQSR之后，區(qū)分了好和壞變異之后，SOR的密度分布。很明顯，好的變異基本就在1附近。結(jié)合這個(gè)分布圖，我們?cè)谏厦娴睦永锇阉拈撝刀?基本上也不會(huì)過損失好的變異了，雖然單靠這個(gè)閾值還是會(huì)保留下不少假的變異，但是至少不合理的長(zhǎng)尾部分可以被砍掉。
RMSMappingQuality（MQ）
MQ這個(gè)值是所有比對(duì)至該位點(diǎn)上的read的比對(duì)質(zhì)量值的均方根（先平方、再平均、然后開方，如下公式）。
它和平均值相比更能夠準(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的就可以過濾掉了。
但是，理性告訴我們還是要看一下VQSR的對(duì)比結(jié)果（下圖）。
你會(huì)發(fā)現(xiàn)似乎所有好的變異都緊緊集中在60旁邊了，其它地方就都是假的變異了，所以MQ的閾值設(shè)置為50也是合理的。但是同樣要注意到的地方是，60這個(gè)范圍實(shí)際上依然有假的變異位點(diǎn)在那里，我們把這個(gè)區(qū)域放大來看，如下圖，這里你就會(huì)發(fā)現(xiàn)其實(shí)假變異的密度分布圖也覆蓋到60這個(gè)范圍了。
考慮到篇幅的問題，接下來MappingQualityRankSumTest（MQRankSum）和ReadPOSRankSumTest（ReadPOSRankSum）的閾值設(shè)定原理，我不打算再細(xì)說下去，思路和上面的4個(gè)是完全一樣的。都是通過比較VQSR之后的密度分布圖，最后確定了硬過濾的閾值。
但請(qǐng)不要以為這只是適用于GATK得到的變異，實(shí)際上，只要我們弄懂了這些指標(biāo)選擇的原因和過濾的思路，那么通過任何其他的變異檢測(cè)工具也是依舊可以適用的，區(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。（如下圖）
另外，在哺乳動(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)行變異過濾的。
小結(jié)
雖然本文一直在談?wù)摰氖侨绾巫龊糜策^濾，但不管我們的指標(biāo)和閾值設(shè)置的多么完美，硬過濾的硬傷都是一刀切，它并不會(huì)根據(jù)多個(gè)維度的數(shù)據(jù)自動(dòng)設(shè)計(jì)出更加合理的過濾模式。硬過濾作為一個(gè)簡(jiǎn)單粗暴的方法，不到不得已的時(shí)候不推薦使用，即使使用了，也一定要明白每一個(gè)過濾指標(biāo)和閾值都意味著什么。
最后，我也希望通過這一篇文章能夠完整地為你呈現(xiàn)一個(gè)變異質(zhì)控的思路。
※ ?※ ?※
推薦閱讀
GATK4.0和全基因組數(shù)據(jù)分析實(shí)踐（上）
從零開始完整學(xué)習(xí)全基因組測(cè)序數(shù)據(jù)分析：第4節(jié) 構(gòu)建WGS主流程
從零開始完整學(xué)習(xí)全基因組測(cè)序數(shù)據(jù)分析：第3節(jié) 數(shù)據(jù)質(zhì)控
從零開始完整學(xué)習(xí)全基因組測(cè)序數(shù)據(jù)分析：第1節(jié) 測(cè)序技術(shù)
這是我的知識(shí)星球：解螺旋技術(shù)交流圈，是一個(gè)與讀者朋友們的私人朋友圈，歡迎你的加入。這是知識(shí)星球上第一個(gè)真正與基因組學(xué)和生物信息學(xué)強(qiáng)相關(guān)的圈子，我們的目標(biāo)是共同營(yíng)造一個(gè)高質(zhì)量的組學(xué)知識(shí)圈和人脈圈，讓圈子里的人都能夠受益。圈子從正式開放三個(gè)月以來，主題內(nèi)容已經(jīng)超過200個(gè)，其中精華內(nèi)容38個(gè)，同時(shí)已經(jīng)過兩次漲價(jià)，目前人數(shù)即將達(dá)到140人，屆時(shí)將會(huì)按照規(guī)則漲到117元/年。
微信掃一掃
關(guān)注該公眾號(hào)
