GATK--數據預處理，質控，檢測變異

• 版權聲明：本文源自 解螺旋的礦工， 由 XP 整理髮表，共 13781 字。
• 轉載請註明：從零開始完整學習全基因組測序（WGS）數據分析：第4節 構建WGS主流程 | Public Library of Bioinformatics
• 轉載地址：https://www.plob.org/article/11698.html

WGS數據分析的目的是準確檢測出每一個樣本（這裏特指人）基因組中的變異集合，也就是人與人之間存在差別的那些DNA序列。我把整個分析過程按照它們實際要完成的功能，將其分紅了三個大的模塊：

• 原始數據質控
• 數據預處理
• 變異檢測

這或許和不少人看到的WGS分析流程，在結構梳理上有些差別（好比GATK的最佳實踐），但過程當中的各個步驟和所要完成的事情是如出一轍的。

0.準備階段

在開始以前，咱們須要作一些準備工做，主要是部署好相關的軟件和工具。咱們在這個WGS數據分析過程當中用到的全部軟件都是開源的，它們的代碼所有都可以在github上找到，具體以下：

• BWA（Burrow-Wheeler Aligner）: 這是最權威，使用最廣的NGS數據比對軟件，目前已經更新到0.7.16版本；
• Samtools: 是一個專門用於處理比對數據的工具，由BWA的做者（lh3）所編寫；
• Picard: 它是目前最著名的組學研究中心-Broad研究所開發的一款強大的NGS數據處理工具，功能方面和Samtools有些重疊，但更多的是互補，它是由java編寫的，咱們直接下載最新的.jar包就好了。
• GATK: 一樣是Broad研究所開發的，是目前業內最權威、使用最廣的基因數據變異檢測工具。值得注意的是，目前GATK有3.x和4.x兩個不一樣的版本，代碼在github上也是分開的。4.x是今年新推出的，在覈心算法層面並沒太多的修改，但使用了新的設計模式，作了不少功能的整合，是更適合於大規模集羣和雲運算的版本，後續GATK團隊也將主要維護4.x的版本，並且它的代碼是100%開源的，這和3.x只有部分開源的狀況不一樣。看得出GATK今年的此次升級是爲了應對接下來愈來愈多的大規模人羣測序數據而作出的改變，但現階段4.x版本還不穩定，真正在使用的人和機構其實也還很少。短時間來看，3.x版本還將在業內繼續使用一段時間；其次，3.x對於絕大分部的分析需求來講是徹底足夠的。咱們在這裏也以GATK3.8（最新版本）做爲流程的重要工具進行分析流程的構建。

事實上，對於構造WGS分析流程來講，以上這個四個工具就徹底足夠了。它們的安裝都很是簡單，除了BWA和Samtools由C編寫的，安裝時須要進行編譯以外，另外兩個只要保證系統中的java是1.8.x版本及以上的，那麼直接下載jar包就可使用了。操做系統方面推薦linux（集羣）或者Mac OS。

1.原始數據質控

數據的質控，因爲我已經在上一節的文章中講的比較詳細了，所以在本篇中就再也不進行詳細的討論了。並且質控的處理方法都是比較一致的，基本不須要爲特定的分析作定製化的改動，所以，咱們能夠把它做爲WGS主流程以外的一環。但仍是再強調一下，數據質控的地位一樣重要，否則我也沒必要專門爲其單獨寫一篇完整的文章。

2.數據預處理

序列比對

先問一個問題：爲何須要比對？

咱們已經知道NGS測序下來的短序列（read）存儲於FASTQ文件裏面。雖然它們本來都來自於有序的基因組，但在通過DNA建庫和測序以後，文件中不一樣read之間的先後順序關係就已經所有丟失了。所以，FASTQ文件中緊挨着的兩條read之間沒有任何位置關係，它們都是隨機來自於本來基因組中某個位置的短序列而已。

所以，咱們須要先把這一大堆的短序列捋順，一個個去跟該物種的 參考基因組【注】比較，找到每一條read在參考基因組上的位置，而後按順序排列好，這個過程就稱爲測序數據的比對。這 也是核心流程真正意義上的第一步，只有完成了這個序列比對咱們纔有下一步的數據分析。

【注】參考基因組：指該物種的基因組序列，是已經組裝成的完整基因組序列，常做爲該物種的標準參照物，好比人類基因組參考序列，fasta格式。

序列比對本質上是一個尋找最大公共子字符串的過程。你們若是有學過生物信息學的話，應該或多或少知道BLAST，它使用的是動態規劃的算法來尋找這樣的子串，但在面對巨量的短序列數據時，相似BLAST這樣的軟件實在太慢了！所以，須要更加有效的數據結構和相應的算法來完成這個搜索定位的任務。

咱們這裏將用於流程構建的BWA就是其中最優秀的一個，它將BW(Burrows-Wheeler)壓縮算法和後綴樹相結合，可以讓咱們以較小的時間和空間代價，得到準確的序列比對結果。

如下咱們就開始流程的搭建。

首先，咱們須要爲參考基因組的構建索引——這實際上是在爲參考序列進行Burrows Wheeler變換（wiki: 塊排序壓縮），以便可以在序列比對的時候進行快速的搜索和定位。

 

 

1	$ bwa index human.fasta

以咱們人類的參考基因組（3Gb長度）爲例，這個構造過程須要消耗幾個小時的時間（通常3個小時左右）。完成以後，你會看到相似以下幾個以human.fasta爲前綴的文件：

 

 

1 2 3 4 5 6

. ├── human.fasta.amb ├── human.fasta.ann ├── human.fasta.bwt ├── human.fasta.pac └── human.fasta.sa

這些就是在比對時真正須要被用到的文件。這一步完成以後，咱們就能夠將read比對至參考基因組了：

 

 

1	$ bwa mem -t 4 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' /path/to/human.fasta read_1.fq.gz read_2.fq.gz > sample_name.sam

大夥若是之前沒使用過這個比對工具的話，那麼可能不明白上面參數的含義。咱們這裏調用的是bwa的mem比對模塊，在解釋這樣作以前，咱們不妨先看一下bwa mem的官方用法說明，它就一句話：

 

 

1	Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]

其中，[options]是一系列可選的參數，暫時很少說。這裏的 < idxbase>要輸入的是參考基因組的BW索引文件，咱們上面經過bwa index構建好的那幾個以human.fasta爲前綴的文件即是；< in1.fq>和 [in2.fq]輸入的是質控後的fastq文件。但這裏輸入的時候爲何會須要兩個fq（in1.fq和in2.fq）呢？咱們上面的例子也是有兩個：read_1.fq.gz和read_2.fq.gz。這是由於這是雙末端測序（也稱Pair-End）的狀況，那什麼是「雙末端測序」呢？這兩個fq之間的關係又是什麼？這個我須要簡單解析一下。

咱們已經知道NGS是短讀長的測序技術，一次測出來的read的長度都不會太長，那爲了儘量把一個DNA序列片斷儘量多地測出來，既然測一邊不夠，那就測兩邊，因而就有了一種從被測DNA序列兩端各測序一次的模式，這就被稱爲雙末端測序（Pair-End Sequencing，簡稱PE測序）。以下圖是Pair-End測序的示意圖，中間灰色的是被測序的DNA序列片斷，左邊黃色帶箭頭和右邊藍色帶箭頭的分別是測序出來的read1和read2序列，這裏假定它們的長度都是100bp。雖然不少時候Pair-End測序仍是沒法將整個被測的DNA片斷徹底測通，可是它依然提供了極其有用的信息，好比，咱們知道每一對的read1和read2都來自於同一個DNA片斷，read1和read2之間的距離是這個DNA片斷的長度，並且read1和read2的方向恰好是相反的（這裏排除mate-pair的情形）等，這些信息對於後面的變異檢測等分析來講都是很是有用的。

Pair-End 測序

另外，在read1在fq1文件中位置和read2在fq2文件中的文件中的位置是相同的，並且read ID之間只在末尾有一個’/1’或者’/2’的差異。

read1 ID和read2 ID的差異

既然有雙末端測序，那麼與之對應的就有單末端測序（Single End Sequecing，簡稱SE測序），即只測序其中一端。所以，咱們在使用bwa比對的時候，實際上，in2.fq是非強制性的（因此用方括號括起來），只有是雙末端測序的數據時才須要添加。

回到上面咱們的例子，大夥能夠看到我這裏除了用法中提到的參數以外，還多了2個額外的參數，分別是：-t，線程數，咱們在這裏使用4個線程；-R 接的是 Read Group的字符串信息，這是一個很是重要的信息，以@RG開頭，它是用來將比對的read進行分組的。不一樣的組之間測序過程被認爲是相互獨立的，這個信息對於咱們後續對比對數據進行錯誤率分析和Mark duplicate時很是重要。在Read Group中，有以下幾個信息很是重要：

(1) ID，這是Read Group的分組ID，通常設置爲測序的lane ID（不一樣lane之間的測序過程認爲是獨立的），下機數據中咱們都能看到這個信息的，通常都是包含在fastq的文件名中；

(2) PL，指的是所用的測序平臺，這個信息不要隨便寫！特別是當咱們須要使用GATK進行後續分析的時候，更是如此！這是一個不少新手都容易忽視的一個地方，在GATK中，PL只容許被設置爲：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN這幾個信息。基本上就是目前市場上存在着的測序平臺，固然，若是實在不知道，那麼必須設置爲UNKNOWN，名字方面不區分大小寫。若是你在分析的時候這裏沒設置正確，那麼在後續使用GATK過程當中可能會碰到相似以下的錯誤：

 

 

1	ERROR MESSAGE: The platform (xx) associated with read group GATKSAMReadGroupRecord @RG:xx is not a recognized platform.

這個時候你須要對比對文件的header信息進行重寫，就會稍微比較麻煩。

咱們上面的例子用的是PL:illumina。若是你的數據是CG測序的那麼記得不要寫成CG！而要寫COMPLETE。

(3) SM，樣本ID，一樣很是重要，有時候咱們測序的數據比較多的時候，那麼可能會分紅多個不一樣的lane分佈測出來，這個時候SM名字就是能夠用於區分這些樣本。

(4) LB，測序文庫的名字，這個重要性稍微低一些，主要也是爲了協助區分不一樣的group而存在。文庫名字通常能夠在下機的fq文件名中找到，若是上面的lane ID足夠用於區分的話，也能夠不用設置LB；

除了以上這四個以外，還能夠自定義添加其餘的信息，不過如無特殊的須要，對於序列比對而言，這4個就足夠了。這些信息設置好以後，在RG字符串中要用製表符（\t）將它們分開。

最後在咱們的例子中，咱們將比對的輸出結果直接重定向到一份sample_name.sam文件中，這類文件是BWA比對的標準輸出文件，它的具體格式我會在下一篇文章中進行詳細說明。但SAM文件是文本文件，通常整個文件都很是巨大，所以，爲了有效節省磁盤空間，通常都會用samtools將它轉化爲BAM文件（SAM的特殊二進制格式），並且BAM會更加方便於後續的分析。因此咱們上面比對的命令能夠和samtools結合並改進爲：

 

 

1	$ bwa mem -t 4 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' /path/to/human.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b - > sample_name.bam

咱們經過管道(「|」)把比對的輸出如同引導水流同樣導流給samtools去處理，上面samtools view的-b參數指的就是輸出爲BAM文件，這裏須要注意的地方是-b後面的’-‘，它表明就是上面管道引流過來的數據，通過samtools轉換以後咱們再重定向爲sample_name.bam。

關於BWA的其餘參數，我這裏不打算對其進行一一解釋，在絕大多數狀況下，採用默認是合適的作法。

[Tips] BWA MEM比對模塊是有必定適用範圍的：它是專門爲長read比對設計的，目的是爲了解決，第三代測序技術這種可以產生長達幾十kb甚至幾Mbp的read狀況。通常只有當read長度≥70bp的時候，才推薦使用，若是比這個要小，建議使用BWA ALN模塊。

排序

以上，咱們就完成了read比對的步驟。接下來是排序：

排序這一步咱們也是經過使用samtools來完成的，命令很簡單：

 

 

1	Usage: samtools sort [options...] [in.bam]

但在執行以前，咱們有必要先搞明白爲何須要排序，爲何BWA比對後輸出的BAM文件是沒順序的！緣由就是FASTQ文件裏面這些被測序下來的read是隨機分佈於基因組上面的，第一步的比對是按照FASTQ文件的順序把read逐必定位到參考基因組上以後，隨即就輸出了，它不會也不可能在這一步裏面可以自動識別比對位置的前後位置重排比對結果。所以，比對後獲得的結果文件中，每一條記錄之間位置的前後順序是亂的，咱們後續去重複等步驟都須要在比對記錄按照順序從小到大排序下來才能進行，因此這纔是須要進行排序的緣由。對於咱們的例子來講，這個排序的命令以下：

 

 

1	$ time samtools sort -@ 4 -m 4G -O bam -o sample_name.sorted.bam sample_name.bam

其中，-@，用於設定排序時的線程數，咱們設爲4；-m，限制排序時最大的內存消耗，這裏設爲4GB；-O 指定輸出爲bam格式；-o 是輸出文件的名字，這裏叫sample_name.sorted.bam。我會比較建議大夥在作相似分析的時候在文件名字將所作的關鍵操做包含進去，由於這樣即便過了很長時間，當你再去看這個文件的時候也可以馬上知道當時對它作了什麼；最後就是輸入文件——sample_name.bam。

【注意】排序後若是發現新的BAM文件比原來的BAM文件稍微小一些，不用以爲驚訝，這是壓縮算法致使的結果，文件內容是沒有損失的。

去除重複序列（或者標記重複序列）

在排序完成以後咱們就能夠開始執行去除重複（準確來講是 去除PCR重複序列）的步驟了。

首先，咱們須要先理解什麼是重複序列，它是如何產生的，以及爲何須要去除掉？要回答這幾個問題，咱們須要再次理解在建庫和測序時到底發生了什麼。

咱們在第1節中已經知道，在NGS測序以前都須要先構建測序文庫：經過物理（超聲）打斷或者化學試劑（酶切）切斷原始的DNA序列，而後選擇特定長度範圍的序列去進行PCR擴增並上機測序。

所以，這裏重複序列的來源實際上就是由PCR過程當中所引入的。由於所謂的PCR擴增就是把原來的一段DNA序列複製屢次。但是爲何須要PCR擴增呢？若是沒有擴增不就能夠省略這一步了嗎？

狀況確實如此，可是不少時候咱們構建測序文庫時能用的細胞量並不會很是充足，並且在打斷的步驟中也會引發部分DNA的降解，這兩點會使總體或者局部的DNA濃度太低，這時若是直接從這個溶液中取樣去測序就極可能漏掉本來基因組上的一些DNA片斷，致使測序不全。而PCR擴增的做用就是爲了把這些微弱的DNA多複製幾倍乃至幾十倍，以便增大它們在溶液中分佈的密度，使得可以在取樣時被獲取到。因此這裏你們須要記住一個重點，PCR擴增本來的目的是爲了增大微弱DNA序列片斷的密度，但因爲整個反應都在一個試管中進行，所以其餘一些密度並不低的DNA片斷也會被同步放大，那麼這時在取樣去上機測序的時候，這些DNA片斷就極可能會被重複取到相同的幾條去進行測序（下圖爲PCR擴增示意圖）。

PCR擴增示意圖：PCR擴增是一個指數擴增的過程，圖中本來只有一段雙鏈DNA序列，在通過3輪PCR後就被擴增成了8段

看到這裏，你或許會以爲，那不必去除不也應該能夠嗎？由於即使擴增了多幾回，不也一樣仍是原來的那一段DNA嗎？直接用來分析對結果也不會有影響啊！難道不是嗎？

會有影響，並且有時影響會很大！最直接的後果就是同時增大了變異檢測結果的假陰和假陽率。主要有幾個緣由：

• DNA在打斷的那一步會發生一些損失，主要表現是會引起一些鹼基發生顛換變換（嘌呤-變嘧啶或者嘧啶變嘌呤），帶來假的變異。PCR過程會擴大這個信號，致使最後的檢測結果中混入了假的結果；
• PCR反應過程當中也會帶來新的鹼基錯誤。發生在前幾輪的PCR擴增發生的錯誤會在後續的PCR過程當中擴大，一樣帶來假的變異；
• 對於真實的變異，PCR反應可能會對包含某一個鹼基的DNA模版擴增更加重烈（這個現象稱爲PCR Bias）。若是反應體系是對含有reference allele的模板擴增偏向強烈，那麼變異鹼基的信息會變小，從而會致使假陰。

PCR對真實的變異檢測和個體的基因型判斷都有很差的影響。GATK、Samtools、Platpus等這種利用貝葉斯原理的變異檢測算法都是認爲所用的序列數據都不是重複序列（即將它們和其餘序列一視同仁地進行變異的判斷，因此帶來誤導），所以必需要進行標記（去除）或者使用PCR-Free的測序方案（這個方案目前正變得愈來愈流行，特別是對於RNA-Seq來講尤其重要，如今著名的基因組學研究所——Broad Institute，基本都是使用PCR-Free的測序方案）。

那麼具體是如何作到去除這些PCR重複序列的呢？咱們能夠拋開任何工具，仔細想一想，既然PCR擴增是把同一段DNA序列複製出不少份，那麼這些序列在通過比對以後它們必定會定位到基因組上相同的位置，比對的信息看起來也將是同樣的！因而，咱們就能夠根據這個特色找到這些重複序列了！

事實上，現有的工具包括Samtools和Picard中去除重複序列的算法也的確是這麼作的。不一樣的地方在於，samtools的rmdup是直接將這些重複序列從比對BAM文件中刪除掉，而Picard的MarkDuplicates默認狀況則只是在BAM的FLAG信息中標記出來，而不是刪除，所以這些重複序列依然會被留在文件中，只是咱們能夠在變異檢測的時候識別到它們，並進行忽略。

考慮到儘量和如今主流的作法一致（但我並非說主流的作法就必定是對的，要分狀況看待，只是主流的作法容易被作成生產流程而已），咱們這裏也用Picard來完成這個事情：

 

 

1 2 3 4

java -jar picard.jar MarkDuplicates \ I=sample_name.sorted.bam \ O=sample_name.sorted.markdup.bam \ M=sample_name.markdup_metrics.txt

這裏只把重複序列在輸出的新結果中標記出來，但不刪除。若是咱們非要把這些序列徹底刪除的話能夠這樣作：

 

 

1 2 3 4 5

java -jar picard.jar MarkDuplicates \ REMOVE_DUPLICATES=true \ I=sample_name.sorted.bam \ O=sample_name.sorted.markdup.bam \ M=sample_name.markdup_metrics.txt

把參數REMOVE_DUPLICATES設置爲ture，那麼重複序列就被刪除掉，不會在結果文件中留存。我比較建議使用第一種作法，只是標記出來，並留存這些序列，以便在你須要的時候還能夠對其作分析。

這一步完成以後，咱們須要爲sample_name.sorted.markdup.bam建立索引文件，它的做用可以讓咱們能夠隨機訪問這個文件中的任意位置，並且後面的「局部重比對」步驟也要求這個BAM文件必定要有索引，命令以下：

 

 

1	$ samtools index sample_name.sorted.markdup.bam

完成以後，會生成一份sample_name.sorted.markdup.bam.bai文件，這就是上面這份BAM的index。

局部重比對

接下來是局部區域重比對，一般也叫Indel局部區域重比對。有時在進行這一步驟以前還有一個merge的操做，將同個樣本的全部比對結果合併成惟一一個大的BAM文件【注】，merge的例子以下：

 

 

1	$ samtools merge <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]

 

【注意】之因此會有這種狀況，是由於有些樣本測得很是深，其測序結果須要通過屢次測序（或者分佈在多個不一樣的測序lane中）才所有得到，這個時候咱們通常會先分別進行比對並去除重複序列後再使用samtools進行合併。

局部重比對的目的是將BWA比對過程當中所發現有 潛在序列插入或者序列刪除（insertion和deletion，簡稱Indel）的區域進行從新校訂。這個過程每每還會把一些已知的Indel區域一併做爲重比對的區域，但爲何須要進行這個校訂呢？

其根本緣由來自於參考基因組的序列特色和BWA這類比對算法自己，注意這裏不是針對BWA，而是針對全部的這類比對算法，包括bowtie等。這類在全局搜索最優匹配的算法在存在Indel的區域及其附近的比對狀況每每不是很準確，特別是當一些存在長Indel、重複性序列的區域或者存在長串單一鹼基（好比，一長串的TTTT或者AAAAA等）的區域中更是如此。

另外一個重要的緣由是在這些比對算法中，對鹼基錯配和開gap的容忍度是不一樣的。具體體如今罰分矩陣的偏向上，例如，在read比對時，若是發現鹼基錯配和開gap均可以的話，它們會更偏向於錯配。可是這種偏向錯配的方式，有時候卻還會反過來引發錯誤的開gap！這就會致使基因組上本來應該是一個長度比較大的Indel的地方，被錯誤地切割成多個錯配和短indel的混合集，這必然會讓咱們檢測到不少錯誤的變異。並且，這種狀況還會隨着所比對的read長度的增加（好比三代測序的Read，一般都有幾十kbp）而變得越加嚴重。

所以，咱們須要有一種算法來對這些區域進行局部的序列重比對。這個算法一般就是大名鼎鼎的Smith-Waterman算法，它很是適合於這類場景，能夠極其有效地實現對全局比對結果的校訂和調整，最大程度低地下降由全局比對算法的不足而帶來的錯誤。並且GATK的局部重比對模塊，除了應用這個算法以外，還會對這個區域中的read進行一次局部組裝，把它們鏈接成爲長度更大的序列，這樣可以更進一步提升局部重比對的準確性。

下圖給你們展現一個序列重比對以前和以後的結果，其中灰色的橫條指的是read，空白黑線指的是deletion，有顏色的鹼基指的是錯配鹼基。

Indel局部重比對的先後的對比

相信你們均可以明顯地看到在序列重比對以前，在這個區域的比對數據是多麼的糟糕，若是就這樣進行變異檢測，那麼必定會獲得不少假的結果。而在通過局部重比對以後，這個區域就變得很是清晰而分明，它本來發生的就只是一個比較長的序列刪除（deletion）事件，但在原始的比對結果中卻被錯誤地用鹼基錯配和短的Indel所代替。

說到這裏，那麼具體該怎麼作呢？咱們的WGS分析流程從這個步驟開始就須要用到GATK (GenomeAnalysisTK.jar)了，咱們的執行命令以下：

 

 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

java -jar /path/to/GenomeAnalysisTK.jar \ -T RealignerTargetCreator \ -R /path/to/human.fasta \ -I sample_name.sorted.markdup.bam \ -known /path/to/gatk/bundle/1000G_phase1.indels.b37.vcf \ -known /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \ -o sample_name.IndelRealigner.intervals   java -jar /path/to/GenomeAnalysisTK.jar \ -T IndelRealigner \ -R /path/to/human.fasta \ -I sample_name.sorted.markdup.bam \ -known /path/to/gatk/bundle/1000G_phase1.indels.b37.vcf \ -known /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \ -o sample_name.sorted.markdup.realign.bam \ --targetIntervals sample_name.IndelRealigner.intervals

這裏包含了兩個步驟：

• 第一步，RealignerTargetCreator ，目的是定位出全部須要進行序列重比對的目標區域（以下圖）；
• 第二步，IndelRealigner，對全部在第一步中找到的目標區域運用算法進行序列重比對，最後獲得捋順了的新結果。

IndelRealigner.intervals文件內容示例

以上這兩個步驟是缺一不可的，順序也是固定的。並且，須要指出的是，這裏的-R參數輸入的human.fasta不是BWA比對中的索引文件前綴，而是參考基因組序列（FASTA格式）文件，下同。

另外，在重比對步驟中，咱們還看到了兩個陌生的VCF文件，分別是：1000G_phase1.indels.b37.vcf和Mills_and_1000G_gold_standard.indels.b37.vcf。這兩個文件來自於千人基因組和Mills項目，裏面記錄了那些項目中檢測到的人羣Indel區域。我上面其實也提到了，候選的重比對區除了要在樣本自身的比對結果中尋找以外，還應該把人羣中已知的Indel區域也包含進來，而這兩個是咱們在重比對過程當中最經常使用到的。這些文件你能夠很方便地在GATK bundle ftp中下載，注意必定要選擇和你的參考基因組對應的版本，咱們這裏用的是b37版本。

GATK bundle

那麼既然Indel局部重比對這麼好，這麼重要，彷佛看起來在任何狀況下都應該是必須的。然鵝，個人回答是否認的！驚訝嗎！

但否認是有前提的！那就是咱們後面的變異檢測必須是使用GATK，並且必須使用GATK的HaplotypeCaller模塊，僅當這個時候才能夠減小這個Indel局部重比對的步驟。緣由是GATK的HaplotypeCaller中，會對潛在的變異區域進行相同的局部重比對！可是其它的變異檢測工具或者GATK的其它模塊就沒有這麼幹了！因此切記！

從新校訂鹼基質量值（BQSR）

在WGS分析中，變異檢測是一個極度依賴測序鹼基質量值的步驟。由於這個質量值是衡量咱們測序出來的這個鹼基到底有多正確的重要（甚至是惟一）指標。它來自於測序圖像數據的base calling。所以，基本上是由測序儀和測序系統來決定的。但不幸的是，影響這個值準確性的系統性因素有不少，包括物理和化學等對測序反應的影響，甚至連儀器自己和周圍環境都是其重要的影響因素。當把全部這些東西綜合在一塊兒以後，每每會發現計算出來的鹼基質量值要麼高於真實結果，要麼低於真實結果。那麼，咱們到底該如何才能得到符合真實狀況的鹼基質量值？

BQSR（Base Quality Score Recalibration）這個步驟就是爲此而存在的，這一步一樣很是重要。它主要是經過機器學習的方法構建測序鹼基的錯誤率模型，而後對這些鹼基的質量值進行相應的調整。

BQSR質量校訂對比

圖中，橫軸（Reported quality score）是測序結果在Base calling以後報告出來的質量值，也就是咱們在FASTQ文件中看到的那些；縱軸（Empirical quality score）表明的是「真實狀況的質量值」。

可是且慢，這個「真實狀況的質量值」是怎麼來的？由於實際上咱們並無辦法直接測得它們啊！沒錯，確實沒辦法直接測量到，可是咱們能夠經過統計學的技巧得到極其接近的分佈結果（所以我加了引號）。試想一下，若是咱們在看到某一個鹼基報告的質量值是20時，那麼它的預期錯誤率是1%，反過來想，就等因而說若是有100個質量值都是20的鹼基，那麼從統計上講它們中將只有1個是錯的！作了這個等效變換以後，咱們的問題就能夠轉變成爲尋找錯誤鹼基的數量了。

這時問題就簡單多了。咱們知道人與人之間的差別實際上是很小的，那麼在一個羣體中發現的已知變異，在某我的身上也極可能是一樣存在的。所以，這個時候咱們能夠對比對結果進行直接分析，首先排除掉全部的已知變異位點，而後計算每一個（報告出來的）質量值下面有多少個鹼基在比對以後與參考基因組上的鹼基是不一樣的，這些不一樣鹼基就被咱們認爲是錯誤的鹼基，它們的數目比例反映的就是真實的鹼基錯誤率，換算成Phred score（Phred score的定義能夠參考第2節的相關內容）以後，就是縱軸的Empirical quality score了。

上面‘BQSR質量校訂對比’的圖中左邊是原始質量值與真實質量值的比較，在這個圖的例子中咱們能夠發現，base calling給出的質量值並無正確地反映真實的錯誤率狀況，測序報告出來的鹼基質量值大部分被高估了，換句話說，就是錯誤率被低估了。

在咱們的流程中，BQSR的具體執行命令以下：

 

 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

java -jar /path/to/GenomeAnalysisTK.jar \ -T BaseRecalibrator \ -R /path/to/human.fasta \ -I sample_name.sorted.markdup.realign.bam \ --knownSites /path/to/gatk/bundle/1000G_phase1.indels.b37.vcf \ --knownSites /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \ --knownSites /path/to/gatk/bundle/dbsnp_138.b37.vcf \ -o sample_name.recal_data.table   java -jar /path/to/GenomeAnalysisTK.jar \ -T PrintReads \ -R /path/to/human.fasta \ -I sample_name.sorted.markdup.realign.bam \ --BQSR sample_name.recal_data.table \ -o sample_name.sorted.markdup.realign.BQSR.bam

這裏一樣包含了兩個步驟：

• 第一步，BaseRecalibrator，這裏計算出了全部須要進行重校訂的read和特徵值，而後把這些信息輸出爲一份校準表文件（sample_name.recal_data.table）
• 第二步，PrintReads，這一步利用第一步獲得的校準表文件（sample_name.recal_data.table）從新調整原來BAM文件中的鹼基質量值，並使用這個新的質量值從新輸出一份新的BAM文件。

注意，由於BQSR其實是爲了（儘量）校訂測序過程當中的系統性錯誤，所以，在執行的時候是按照不一樣的測序lane或者測序文庫來進行的，這個時候@RG信息（BWA比對時所設置的）就顯得很重要了，算法就是經過@RG中的ID來識別各個獨立的測序過程，這也是我開始強調其重要性的緣由。

變異檢測

事實上，這是目前全部WGS數據分析流程的一個目標——得到樣本準確的變異集合。這裏變異檢測的內容通常會包括：SNP、Indel，CNV和SV等，這個流程中咱們只作其中最主要的兩個：SNP和Indel。咱們這裏使用GATK HaplotypeCaller模塊對樣本中的變異進行檢測，它也是目前最適合用於對二倍體基因組進行變異（SNP+Indel）檢測的算法。

HaplotypeCaller和那些直接應用貝葉斯推斷的算法有所不一樣，它會先推斷羣體的單倍體組合狀況，計算各個組合的概率，而後根據這些信息再反推每一個樣本的基因型組合。所以它不但特別適合應用到羣體的變異檢測中，並且還可以依據羣體的信息更好地計算每一個個體的變異數據和它們的基因型組合。

通常來講，在實際的WGS流程中對HaplotypeCaller的應用有兩種作法，差異只在於要不要在中間生成一個gVCF：

（1）直接進行HaplotypeCaller，這適合於單樣本，或者那種固定樣本數量的狀況，也就是執行一次HaplotypeCaller以後就老死不相往來了。不然你會碰到僅僅只是增長一個樣本就得從新運行這個HaplotypeCaller的坑爹狀況（即，N+1難題），而這個時候算法須要從新去讀取全部人的BAM文件，這將會是一個很費時間的痛苦過程；

（2）每一個樣本先各自生成gVCF，而後再進行羣體joint-genotype。這其實就是GATK團隊爲了解決（1）中的N+1難題而設計出來的模式。gVCF全稱是genome VCF，是每一個樣本用於變異檢測的中間文件，格式相似於VCF，它把joint-genotype過程當中所需的全部信息都記錄在這裏面，文件不管是大小仍是數據量都遠遠小於原來的BAM文件。這樣一旦新增長樣本也不須要再從新去讀取全部人的BAM文件了，只需爲新樣本生成一份gVCF，而後從新執行這個joint-genotype就好了。

咱們先以第一種（直接HaplotypeCaller）作法爲例子：

 

 

1 2 3 4 5 6 7 8 9 10 11 12 13 14

java -jar /path/to/GenomeAnalysisTK.jar \ -T HaplotypeCaller \ -R /path/to/human.fasta \ -I sample_name.sorted.markdup.realign.BQSR.bam \ -D /path/to/gatk/bundle/dbsnp_138.b37.vcf \ -stand_call_conf 50 \ -A QualByDepth \ -A RMSMappingQuality \ -A MappingQualityRankSumTest \ -A ReadPosRankSumTest \ -A FisherStrand \ -A StrandOddsRatio \ -A Coverage \ -o sample_name.HC.vcf

這裏我特別提一下-D參數輸入的dbSNP一樣能夠再GATK bundle目錄中找到，這份文件聚集的是目前幾乎全部的公開人羣變異數據集。另外，因爲咱們的例子只有一個樣本所以只輸入一個BAM文件就能夠了，若是有多個樣本那麼能夠繼續用-I參數輸入：

 

 

1 2 3 4 5

java -jar GenomeAnalysisTK.jar \ -T HaplotypeCaller \ -R reference.fasta \ -I sample1.bam [-I sample2.bam ...] \ ...

以上的命令是直接對全基因組作變異檢測，這個過程會消耗很長的時間，一般須要幾十個小時甚至幾天。

然而，基因組上各個不一樣的染色體之間實際上是能夠理解爲相互獨立的（結構性變異除外），也就是說，爲了提升效率咱們能夠按照染色體一條條來獨立執行這個步驟，最後再把結果合併起來就行了，這樣的話就可以節省不少的時間。下面我給出一個按照染色體區分的例子：

 

 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

java -jar /path/to/GenomeAnalysisTK.jar \ -T HaplotypeCaller \ -R /path/to/human.fasta \ -I sample_name.sorted.markdup.realign.BQSR.bam \ -D /path/to/gatk/bundle/dbsnp_138.b37.vcf \ -L 1 \ -stand_call_conf 50 \ -A QualByDepth \ -A RMSMappingQuality \ -A MappingQualityRankSumTest \ -A ReadPosRankSumTest \ -A FisherStrand \ -A StrandOddsRatio \ -A Coverage \ -o sample_name.HC.1.vcf

注意到了嗎？其它參數都沒任何改變，就只增長了一個 -L 參數，經過這個參數咱們能夠指定特定的染色體（或者基因組區域）！咱們這裏指定的是 1 號染色體，有些地方會寫成chr1，具體看human.fasta中如何命名，與其保持一致便可。其餘染色體的作法也是如此，就再也不舉例了。最後合併：

 

 

1 2 3 4 5 6 7 8 9

java -jar /path/to/GenomeAnalysisTK.jar \ -T CombineVariants \ -R /path/to/human.fasta \ --genotypemergeoption UNSORTED \ --variant sample_name.HC.1.vcf \ --variant sample_name.HC.2.vcf \ ... --variant sample_name.HC.MT.vcf \ -o sample_name.HC.vcf

第二種，先產生gVCF，最後再joint-genotype的作法：

 

 

1 2 3 4 5 6 7 8 9 10 11 12 13

java -jar /path/to/GenomeAnalysisTK.jar \ -T HaplotypeCaller \ -R /path/to/human.fasta \ -I sample_name.sorted.markdup.realign.BQSR.bam \ --emitRefConfidence GVCF \ -o sample_name.g.vcf   #調用GenotypeGVCFs完成變異calling java -jar /path/to/GenomeAnalysisTK.jar \ -T GenotypeGVCFs \ -R /path/to/human.fasta \ --variant sample_name.g.vcf \ -o sample_name.HC.vcf

其實，就是加了–emitRefConfidence GVCF的參數。並且，假如嫌慢，一樣能夠按照染色體或者區域去產生一個樣本的gVCF，而後在GenotypeGVCFs中把它們所有做爲輸入文件完成變異calling。也許你會擔憂同個樣本被分紅多份gVCF以後，是否會被看成不一樣的多個樣本？回答是不會！由於生成gVCF文件的過程當中，GATK會根據@RG信息中的SM（也就是sample name）來判斷這些gVCF是否來自同一個樣本，若是名字相同，那麼就會被認爲是同一個樣本，不會產生多樣本問題。

變異檢測質控和過濾（VQSR）

這是咱們這個流程中最後的一步了。在得到了原始的變異檢測結果以後，咱們還須要作的就是質控和過濾。這一步或多或少都有着一些個性化的要求，我暫時就不作太多解釋吧（一旦解釋恐怕一樣是一篇萬字長文）。只用一句話來歸納，VQSR是經過構建GMM模型對好和壞的變異進行區分，從而實現對變異的質控，具體的原理暫時不展開了。

下面就直接給出例子吧：

 

 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

## SNP Recalibrator java -jar /path/to/GenomeAnalysisTK.jar \ -T VariantRecalibrator \ -R reference.fasta \ -input sample_name.HC.vcf \ -resource:hapmap,known=false,training=true,truth=true,prior=15.0 /path/to/gatk/bundle/hapmap_3.3.b37.vcf \ -resource:omini,known=false,training=true,truth=false,prior=12.0 /path/to/gatk/bundle/1000G_omni2.5.b37.vcf \ -resource:1000G,known=false,training=true,truth=false,prior=10.0 /path/to/gatk/bundle/1000G_phase1.snps.high_confidence.b37.vcf \ -resource:dbsnp,known=true,training=false,truth=false,prior=6.0 /path/to/gatk/bundle/dbsnp_138.b37.vcf \ -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \ -mode SNP \ -recalFile sample_name.HC.snps.recal \ -tranchesFile sample_name.HC.snps.tranches \ -rscriptFile sample_name.HC.snps.plots.R   java -jar /path/to/GenomeAnalysisTK.jar -T ApplyRecalibration \ -R human_g1k_v37.fasta \ -input sample_name.HC.vcf \ --ts_filter_level 99.5 \ -tranchesFile sample_name.HC.snps.tranches \ -recalFile sample_name.HC.snps.recal \ -mode SNP \ -o sample_name.HC.snps.VQSR.vcf   ## Indel Recalibrator java -jar /path/to/GenomeAnalysisTK.jar -T VariantRecalibrator \ -R human_g1k_v37.fasta \ -input sample_name.HC.snps.VQSR.vcf \ -resource:mills,known=true,training=true,truth=true,prior=12.0 /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \ -an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum \ -mode INDEL \ -recalFile sample_name.HC.snps.indels.recal \ -tranchesFile sample_name.HC.snps.indels.tranches \ -rscriptFile sample_name.HC.snps.indels.plots.R   java -jar /path/to/GenomeAnalysisTK.jar -T ApplyRecalibration \ -R human_g1k_v37.fasta\ -input sample_name.HC.snps.VQSR.vcf \ --ts_filter_level 99.0 \ -tranchesFile sample_name.HC.snps.indels.tranches \ -recalFile sample_name.HC.snps.indels.recal \ -mode INDEL \ -o sample_name.HC.snps.indels.VQSR.vcf

 

最後，sample_name.HC.snps.indels.VQSR.vcf即是咱們最終的變異檢測結果。對於人類而言，通常來講，每一個人最後檢測到的變異數據大概在400萬左右（包括SNP和Indel）。

這篇文章已經很長了，在變異檢測的這個過程當中GATK應用了不少重要的算法，包括如何構建模型、如何進行局部組裝和比對、如何應用貝葉斯、PariHMM、GMM、參數訓練、特徵選擇等等，這些只能留在後面介紹GATK的專題文章中再進行展開了。

小結

在這裏，整篇文章就結束了。如你所見，文章很是長，這裏基本包含了WGS最佳實踐中的全部內容，但其實我想說的還遠不止如此（包括CNV和SV的檢測），只是暫時只能做罷了，不然恐怕就沒人願意看下去了，呵呵。在這個WGS主流程的構建過程當中，我並不是只是硬邦邦地告訴你們幾條簡單的命令就了事了，由於我認爲那種作法要麼是極其不負責任的，要麼就是做者並不是真的懂。並且若是都以爲只要懂得幾條命令就能夠了的話，那麼咱們就活該被機器和人工智能所取替，它們必定會操做得更好更高效。我想掌握工具和技術的目的是爲了可以更好地發現並解決問題（包括科研和生產），全部的數據分析流程本質上是要服務於咱們所要解決的問題的。

畢竟工具是死的，人是活的，需求老是會變的。理解咱們所要處理的問題的本質，選擇合適的工具，而不是反過來被工具所束縛，這一點很重要。我的的能力不能只是會跑一個流程，或者只是會建立流程，由於那都是一個「術」的問題，我以爲咱們真正要去掌握的應該是如何分析數據的能力，如何發現問題和解決數據問題等的能力。