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

前言

在前面的一系列WGS文章中，我講述了很多基因數據分析的來龍去脈，雖然許多同學覺得很有幫助，但是卻缺了一個重要的環節——沒有提供實際可用的數據來實戰完成具體的流程，不能得到直觀的體會。許多讀者也紛紛在後台留言反饋這個問題，特別是在我寫了如何入門生物信息學的文章之後情況尤甚。所以我決定要再寫一篇文章來解決這個問題，但碰巧今年年末事情稍多，在寫作的過程中也曾多次中斷，直至今天才完成，各位久等了。本次文章分為上下兩篇，這是第一篇，都是WGS系列第四節的延伸和補充，本意是結合具體的數據，讓更多的人能夠更好地理解整個WGS數據的分析和處理過程，我也結合自身的工作經驗給出一些做項目過程中的建議，以作參考，希望能夠對你有更多的幫助。另外，接下來我將系統寫一個關於全基因組關聯分析（GWAS）的文章，同時還會有更多全面而且緊扣前沿的技術文章分享出來。

那麼，事不宜遲我們馬上開始。考慮到實際的數據分析環境，我們只在Linux命令行終端（Terminal）中進行，執行步驟都寫在shell，因此不會有窗口式的操作（使用Mac OS的同學可以使用Mac自帶的Terminal，與Linux操作一致），在這篇文章中我們會用到以下幾個工具：

• sratoolkit
• bgzip
• tabix
• bwa
• samtools
• GATK 4.0

這些軟體都可以在github上找到（包括GATK），需要各位自行安裝。這裡補充一句，目前GATK4.0的正式版本已經發布，它的使用方式與之前相比有著一些差異（變得更加簡單，功能也更加豐富了），增加了結構性變異檢測和很多Spark、Cloud-Only的功能，並集成了Mutect2和picard的所有功能（以及其他很多有用的工具），這為我們減少了許多額外的工具，更加有利於流程的構建和維護，4.0之後的GATK是一個新的篇章，大家最好是掌握這一個版本！另外，3.x的版本貌似也已經不提供下載通道了，如果你還想使用3.x的話可以在公眾號後台回復「GATK3」，我為你準備了一個GATK官方3.7的版本。我們這裡則使用最新的4.0版本。

項目目錄結構

清晰的目錄結構是管理眾多項目的有效途徑，經久不忘，隨時可查。雖然看起來有些原始，但在Linux終端下面，我目前還沒有發現更好的文件管理辦法。這個項目的目錄結構，我的建議是按照時間+項目的規則來命名，下面是我的目錄結構：

./201802_wgs_practice/├── bin├── input└── output

頂層的項目名就是20180203_wgs_practice，下面有三個主目錄：

• input：存儲所有輸入數據
• output：存儲所有輸出數據
• bin：存放所有執行程序和代碼

output只存放結果數據，它是由input和bin中的數據和程序流程生成的。這樣做的好處是層次分明，流程邏輯清楚，數據互不干擾。

使用E.coli K12完成比對和變異檢測

人類基因組數據很大，參考序列長度是3Gb。而一個人的高深度測序數據往往是這個數字的30倍——100Gb。如果直接用這樣的數據來完成本文的分析，那麼許多同學需要下載大量的原始數據。除了下載時間很長之外，如果沒有合適的集群，只是在自己的桌面電腦上干這樣的事情，那麼硬碟空間也將很快不夠用。而且，要在單機電腦上完成這樣一個高深度WGS數據的分析，除了對機器性能有要求之後，跑起來也需要花上差不多140個小時——相信大家都等不起呀。

因此，為了解決這個問題，我找了E.coli K12（一種實驗用的大腸桿菌）的數據作為代替，用來演示數據比對和變異檢測這兩個最消耗計算資源和存儲空間的步驟。E.coli K12的特點是數據很小，它的基因組長度只有4.6Mb，很適合大家用來快速學習WGS的數據分析，遇到人類的數據時，再做替換就行了。

下載E.coli K12的參考基因組序列

熟悉的同學應該第一時間能夠知道，這些物種的基因組參考序列都可以在NCBI上獲取，我們這裡也是一樣，可以在NCBI網站上直接搜索這個序列，為了簡化步驟，我直接給出E.coli K12參考序列的ftp地址給大家下載之用：

ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz

你可以在Linux（或者Mac OSX）命令行上直接使用wget，將這個fasta下載下來，由於它很小，所以幾秒之後我們就可以得到這個fasta序列。

$ wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz

為了接下來表達上的清晰和操作上的方便，我們使用bgzip將這個序列文件進行解壓並把名字重命名為E.coli_K12_MG1655.fa，這樣就一目了然了。

$ gzip -dc GCF_000005845.2_ASM584v2_genomic.fna.gz > E.coli_K12_MG1655.fa

E.coli K12隻有一條完整的染色體，你打開文件後將會看到和我一樣的內容：

>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genomeAGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGAC

接著，我們用samtools為它創建一個索引，這是為方便其他數據分析工具（比如GATK）能夠快速地獲取fasta上的任何序列做準備。

$ /Tools/common/bin/samtools faidx E.coli_K12_MG1655.fa

這時會生成一份E.coli_K12_MG1655.fa.fai文件。除了方便其他工具之外，我們可以通過這樣的索引來獲取fasta文件中任意位置的序列或者任意完整的染色體序列。可以很方便地完成對參考序列（或者任意fasta文件）特定區域序列的提取。舉個例子：

$ samtools faidx E.coli_K12_MG1655.fa NC_000913.3:1000000-1000200

我們就獲得了E.coli K12參考序列上的這一段序列：

>NC_000913.3:1000000-1000200GTGTCAGCTTTCGTGGTGTGCAGCTGGCGTCAGATGACAACATGCTGCCAGACAGCCTGAAAGGGTTTGCGCCTGTGGTGCGTGGTATCGCCAAAAGCAATGCCCAGATAACGATTAAGCAAAATGGTTACACCATTTACCAAACTTATGTATCGCCTGGTGCTTTTGAAATTAGTGATCTCTATTCCACGTCGTCGAGCG

這個小技巧在特定的時候非常實用。

下載E.coli K12的測序數據

基因組參考序列準備好之後，接下來我們需要下載它的測序數據。E.coli K12作為一種供研究使用的模式生物，自然已經有許多的測序數據在NCBI上了，在這裡我們選擇了其中的1個數據——SRR1770413。這個數據來自Illumina MiSeq測序平台（不用擔心平台的事情），read長度是300bp，測序類型Pair-End（沒了解過PE read同學可以參考我前面WGS系列第四節文章）。你可以在NCBI上直接搜到，這裡。

在NCBI給出的信息頁面中，我們可以清楚地看到這個數據的大小（如下圖）——差不多200MB，一般家庭網速也能夠較快下載完成。

從NCBI上下載下來的測序數據，不是我們熟悉的fastq格式，而是SRA（一種NCBI自己設計的測序數據存儲格式，具較高的壓縮率），我們需要對其進行轉換，下文詳述。現在我們先下載，有兩個下載方式（我在這裡告訴大家的方法同樣適用於其他類型的數據），第一個是如上面所說搜索到SRR1770413這個數據的ftq地址，然後直接在命令行中執行wget進行下載，如下：

$ wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR177/SRR1770413/SRR1770413.sra

注意，下載下來的這個SRA文件雖然只有一份，但是裡面其實存了read1和read2的測序數據，我們要把它解出來，轉換成為我們所需的fastq格式。這個時候，我們就需要用到NCBI的官方工具包sratoolkit—— https://github.com/ncbi/sra-tools/wiki/Downloads ，大家下載對應系統的版本，直接解壓之後就可以使用了。

sratoolkit是一個工具包，所有的執行程序都在它解壓後的bin文件夾下，我們要把SRA轉換為fastq，只需直接通過工具包中的fastq-dump即可完成。

$ /Tools/sratoolkit/2.8.2/bin/fastq-dump --split-files SRR1770413.sra

然後我們就會得到這個E.coli K12數據的read1和read2了：

SRR1770413_1.fastqSRR1770413_2.fastq

另一個數據下載的方法就是用上面說到的sratoolkit，我也比較推薦這個方法，操作很簡單。同樣是使用fastq-dump（沒錯，與上面一樣，區別在於下載的時候，輸入的是數據的SRA編號），它可以在我們下載的過程中就直接將SRA轉換為兩個fastq。

$ /Tools/sratoolkit/2.8.2/bin/fastq-dump --split-files SRR1770413

下載完成後，我們最好用bgzip（不推薦gzip）將其壓縮為.gz文件，這樣可以節省空間，而且也不會對接下來的數據分析產生影響。

$ /Tools/common/bin/bgzip -f SRR1770413_1.fastq$ /Tools/common/bin/bgzip -f SRR1770413_1.fastq

至此，E.coli K12相關的數據我們就都準備好了，先看一眼我現在的目錄結構。

./201802_wgs_practice/├── bin├── input│ ├── E.coli│ ├── fasta│ │ ├── E.coli_K12_MG1655.fa│ │ ├── E.coli_K12_MG1655.fa.fai│ │ └── work.log.sh│ └── fastq│ ├── SRR1770413_1.fastq.gz│ ├── SRR1770413_2.fastq.gz│ └── work.log.sh ├── output└── work.log.sh

其中各個目錄下的work.log.sh記錄了我在該目錄下的所有重要操作——這是我的個人習慣，目的是方便以後反查數據的需要。

數據準備完畢之後，接下來就可以進行具體的分析了。

質控

質控是必須做的，我們需要完整認識原始的測序數據質量到底如何，該步驟不能省略。我專門為此單獨寫了一篇文章（WGS系列第三節），在正式的數據分析過程中，大家可以參考它來完成數據的質控，然後再進行接下來的分析。本篇文章為了控制篇幅和儘可能扣住核心內容，就不再對此深入展開，大家如果碰到問題，可以到在後台留言，或者加入我的交流圈（知識星球：解螺旋技術交流圈）和更多有經驗的人一起交流。

比對

首先是比對。所謂比對就是把測序數據定位到參考基因組上，確定每一個read在基因組中的位置。這裡，我們依然用目前使用最廣的BWA來完成這個工作。在正式比對之前，需要先為參考序列構建BWA比對所需的FM-index（比對索引）。

$ /Tools/common/bin/bwa index E.coli_K12_MG1655.fa[bwa_index] Pack FASTA... 0.06 sec[bwa_index] Construct BWT for the packed sequence...[bwa_index] 1.78 seconds elapse.[bwa_index] Update BWT... 0.04 sec[bwa_index] Pack forward-only FASTA... 0.03 sec[bwa_index] Construct SA from BWT and Occ... 0.63 sec[main] Version: 0.7.16a-r1181[main] CMD: /Tools/common/bin/bwa index E.coli_K12_MG1655.fa[main] Real time: 3.266 sec; CPU: 2.550 sec

由於這個序列很短，只需幾秒就可以完成這個索引文件的構建（對於人類基因組則需要3個小時的時間）。創建完畢之後，將多出5份以E.coli_K12_MG1655.fa為前綴的序列索引文件。

E.coli_K12_MG1655.fa.ambE.coli_K12_MG1655.fa.annE.coli_K12_MG1655.fa.bwtE.coli_K12_MG1655.fa.pacE.coli_K12_MG1655.fa.sa

現在我們使用bwa完成比對，用samtools完成BAM格式轉換、排序並標記PCR重複序列。步驟分解如下：

#1 比對time /Tools/common/bin/bwa mem -t 4 -R @RG ID:foo PL:illumina SM:E.coli_K12 /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_1.fastq.gz /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_2.fastq.gz | /Tools/common/bin/samtools view -Sb - > /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo "** bwa mapping done **"#2 排序time /Tools/common/bin/samtools sort -@ 4 -m 4G -O bam -o /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo "** BAM sort done"rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam#3 標記PCR重複time /Tools/common/bin/gatk/4.0.1.2/gatk MarkDuplicates -I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam -M /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup_metrics.txt && echo "** markdup done **"#4 刪除不必要文件(可選)rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bamrm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam#5 創建比對索引文件time /Tools/common/bin/samtools index /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam && echo "** index done **"

從上面的命令大家也可以看到，我嚴格按照上文提到的項目目錄規範來執行（步驟中涉及到的數據路徑也都儘可能使用全路徑），這個比對的shell存放在bin目錄下，名稱是bwa_and_markdup.sh，從名稱也能夠一眼可以看出是要做什麼的。

簡單解釋一下這個shell所做的事情：首先利用bwa mem比對模塊將E.coli K12質控後的測序數據定位到其參考基因組上（我們這裡設置了4個線程來完成比對，根據電腦性能可以適當調大），同時通過管道（| 操作符）將比對數據流引到samtools轉換為BAM格式（SAM的二進位壓縮格式），然後重定向(>操作符)輸出到文件中保存下來。

-R 設置Read Group信息，雖然我在以前的文章中已經反覆強調過它的重要性，但這裡還是再說一次，它是read數據的組別標識，並且其中的ID，PL和SM信息在正式的項目中是不能缺少的(如果樣本包含多個測序文庫的話，LB信息也不要省略)，另外由於考慮到與GATK的兼容關係，PL（測序平台）信息不能隨意指定，必須是：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN這12個中的一個。

接著用samtools對原始的比對結果按照參考序列位置從小到大進行排序（同樣是4個線程），只有這個步驟完成之後才可以繼續往下。

然後，我們使用GATK標記出排完序的數據中的PCR重複序列。這個步驟完成後，如無特殊需要，我們就可以直接刪除前面那兩個BAM文件了（原始比對結果和排序後的結果）——後續幾乎不會再用到那兩份文件了。關於標記PCR重複序列的操作比較簡單，不再細說（如果希望了解更多有關重複序列特徵的信息可以回看WGS系列第四節中的內容）。

最後，我們再用samtools為E_coli_K12.sorted.markdup.bam創建索引。我認為不論是否有後續分析，為BAM文件創建索引應該作為一個常規步驟，它可以讓我們快速地訪問基因組上任意位置的比對情況，這一點非常有助於我們隨時了解數據。

至於每個步驟最前面的time，則是用於記錄執行時間的，有助於我們清楚地知道每一個分析過程都花了多少時間，當需要優化流程的時候這個信息會很有用。

變異檢測

接下來是用GATK完成變異檢測。但在開始之前之前我們還需要先為E.coli K12的參考序列生成一個.dict文件，這可以通過調用CreateSequenceDictonary模塊來完成(這是原來picard的功能)。

$ /Tools/common/bin/gatk/4.0.1.2/gatk CreateSequenceDictionary -R E.coli_K12_MG1655.fa -O E.coli_K12_MG1655.dict && echo "** dict done **"

唯一需要注意的是.dict文件的名字前綴需要和fasta的一樣，並跟它在同一個路徑下，這樣GATK才能夠找到。

OK，現在我們就可以進行變異檢測了，同樣使用GATK 4.0的HaplotypeCaller模塊來完成。由於我們只有一個樣本，要完成這個工作其實很簡單，直接輸入比對文件和參考序列就行了，但是考慮到實際的情況，我想告訴大家一個更好的方式（雖然這會多花些時間），就是：先為每個樣本生成一個GVCF，然後再用GenotypeGVCFs對這些GVCF進行joint calling，如下 ，我把命令都寫在gatk.sh中，並執行。

#1 生成中間文件gvcftime /Tools/common/bin/gatk/4.0.1.2/gatk HaplotypeCaller -R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa --emit-ref-confidence GVCF -I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf && echo "** gvcf done **"#2 通過gvcf檢測變異time /Tools/common/bin/gatk/4.0.1.2/gatk GenotypeGVCFs -R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa -V /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf && echo "** vcf done **"

很快我們就獲得了E.coli K12這個樣本初步的變異結果——E_coli_K12.vcf。之所以非要分成兩個步驟，是因為我想藉此告訴大家，變異檢測不是一個樣本的事情，有越多的同類樣本放在一起joint calling結果將會越準確，而如果樣本足夠多的話，在低測序深度的情況下也同樣可以獲得完整並且準確的結果，而這樣的分步方式是應對多樣本的好方法。

最後，我們用bgzip對這個VCF進行壓縮，並用tabix為它構建索引，方便以後的分析。

#1 壓縮 time /Tools/common/bin/bgzip -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf#2 構建tabix索引time /Tools/common/bin/tabix -p vcf /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf.gz

bgzip壓縮完成之後，原來的VCF文件會被自動刪除。

為了保持一致，現在再看一下完成到這裡之後我們的目錄長什麼樣了，供大家對照。

./201802_wgs_practice/├── bin│ ├── bwa_and_markdup.sh│ └── gatk.sh├── input│ └── E.coli│ ├── fasta│ │ ├── E.coli_K12_MG1655.dict│ │ ├── E.coli_K12_MG1655.fa│ │ ├── E.coli_K12_MG1655.fa.amb│ │ ├── E.coli_K12_MG1655.fa.ann│ │ ├── E.coli_K12_MG1655.fa.bwt│ │ ├── E.coli_K12_MG1655.fa.fai│ │ ├── E.coli_K12_MG1655.fa.pac│ │ ├── E.coli_K12_MG1655.fa.sa│ │ └── work.log.sh│ └── fastq│ ├── SRR1770413_1.fastq.gz│ ├── SRR1770413_2.fastq.gz│ └── work.log.sh├── output│ └── E.coli│ ├── E_coli_K12.g.vcf│ ├── E_coli_K12.g.vcf.idx│ ├── E_coli_K12.sorted.markdup.bam│ ├── E_coli_K12.sorted.markdup.bam.bai│ ├── E_coli_K12.sorted.markdup_metrics.txt│ ├── E_coli_K12.vcf│ └── E_coli_K12.vcf.idx└── work.log.sh

如果大家仔細看過WGS系列第四節的話，會發現我這裡缺少了兩個步驟：重比對和BQSR。沒有執行BQSR是因為E.coli K12沒有那些必須的known變異集（或者有但我沒找到），所以無法進行；但沒有重比對，則是因為我在GATK 4.0中沒發現IndelRealigner這個功能，雖然我們使用GATK HaplotypeCaller或者Mutect2的話確實可以省略這個步驟，但如果是其他軟體來進行變異檢測那麼該步驟依然十分重要，我目前不太清楚為何GATK 4.0沒有將這個功能單獨分離出來。

後面要談到的就是變異的質控了。很遺憾我們這個E.coli K12的變異結果並不適合通過VQSR來進行過濾，原因上面也提到了一些，它不像人類的基因組數據，有著一套適合用來訓練過濾模型的已知變異集（dbSNP，1000G，Hapmap和omini等）。其實這種情況有時候我們在工作中也會碰到，比如有些捕獲測序（Panel測序數據，甚至外顯子測序）的數據，由於它的區域較小，獲得的變異也不多，導致最終沒法滿足VQSR進行模型訓練時所需的最低變異數要求，那時你也不能通過這個方式協助變異質控。那麼碰到這種情況的時候該怎麼辦？我將這部分的內容放在了下一篇文章中，在那裡我們再來討論這個問題。我也會告訴大家變異質控的基本邏輯，而不是簡單羅列一個命令，同時也會再用NA12878這個人的數據來進一步告訴大家如何比較和評估變異結果。

小結

至此，這個篇文章的上半部分就到此為止了。除了那些重要的內容之外，在上文中，你會看到我反覆提到了創建「索引」這個事情，比如為fasta，為BAM，為VCF。我為什麼非要反覆強調這個事情不可呢？因為我發現許多初學者並不知道索引的作用，當被問到如何從巨大的比對文件或者變異文件中提取某個信息時，總是要走彎路——努力寫程序去提取，既慢又費力，結果還不一定好，甚至有些有一定經驗的同學也不知道使用bgzip和tabix的好處，因此我才反覆在文章里提及。

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