一、使用GATK前須知事項：
（1）對GATK的測試主要使用的是人類全基因組和外顯子組的測序數據，而且全部是基於illumina數據格式，目前還沒有提供其他格式文件（如Ion Torrent）或者實驗設計（RNA-Seq）的分析方法。
（2）GATK是一個應用於前沿科學研究的軟件，不斷在更新和修正，因此，在使用GATK進行變異檢測時，最好是下載最新的版本，目前的版本是2.8.1（2014-02-25）。下載網站：http://www.broadinstitute.org/gatk/download。
（3）在GATK使用過程中（見下面圖），有些步驟需要用到已知變異信息，對於這些已知變異，GATK只提供了人類的已知變異信息，可以在GATK的FTP站點下載（GATK resource bundle）。如果要研究的不是人類基因組，需要自行構建已知變異，GATK提供了詳細的構建方法。
（4）GATK在進行BQSR和VQSR的過程中會使用到R軟件繪制一些圖，因此，在運行GATK之前最好先檢查一下是否正確安裝了R和所需要的包，所需要的包大概包括ggplot2、gplots、bitops、caTools、colorspace、gdata、gsalib、reshape、RColorBrewer等。如果畫圖時出現錯誤，會提示需要安裝的包的名稱。
二、GATK的使用流程
GATK最佳使用方案：共3大步驟。原始數據的處理—變異檢測—初步分析。
第一大步：原始數據的處理
1. 對原始下機fastq文件進行過濾和比對（mapping）
對於Illumina下機數據推薦使用bwa進行mapping。
Bwa比對步驟大致如下：
（1）對參考基因組構建索引：
例子：bwa index -a bwtsw hg19.fa。最后生成文件：hg19.fa.amb、hg19.fa.ann、hg19.fa.bwt、hg19.fa.pac和hg19.fa.sa。
構建索引時需要注意的問題：bwa構建索引有兩種算法，兩種算法都是基於BWT的，這兩種算法通過參數-a is 和-a bwtsw進行選擇。其中-a bwtsw對於短的參考序列是不工作的，必須要大於等於10Mb；-a is是默認參數，這個參數不適用於大的參考序列，必須要小於等於2G。
（2）尋找輸入reads文件的SA坐標。
對於pair end數據，每個reads文件單獨做運算，single end數據就不用說了，只有一個文件。
例子：pair end：
bwa aln hg19.fa read1.fq.gz -l 30 -k 2 -t 4 -I > read1.fq.gz.sai
bwa aln hg19.fa read2.fq.gz -l 30 -k 2 -t 4 -I > read2.fq.gz.sai
single end：
bwa aln hg19.fa read.fq.gz -l 30 -k 2 -t 4 -I > read.fq.gz.sai
主要參數說明：
-o int：允許出現的最大gap數。
-e int：每個gap允許的最大長度。
-d int：不允許在3’端出現大於多少bp的deletion。
-i int：不允許在reads兩端出現大於多少bp的indel。
-l int：Read前多少個鹼基作為seed，如果設置的seed大於read長度，將無法繼續，最
好設置在25-35，與-k 2 配合使用。
-k int：在seed中的最大編輯距離，使用默認2，與-l配合使用。
-t int：要使用的線程數。
-R int：此參數只應用於pair end中，當沒有出現大於此值的最佳比對結果時，將會降低標
准再次進行比對。增加這個值可以提高配對比對的准確率，但是同時會消耗更長的
時間，默認是32。
-I int：表示輸入的文件格式為Illumina 1.3+數據格式。
-B int：設置標記序列。從5’端開始多少個鹼基作為標記序列，當-B為正值時，在比對之
前會將每個read的標記序列剪切，並將此標記序列表示在BC SAM 標簽里，對於
pair end數據，兩端的標記序列會被連接。
-b ：指定輸入格式為bam格式。bwa aln hg19.fa read.bam > read.fq.gz.sai
（3）生成sam格式的比對文件。如果一條read比對到多個位置，會隨機選擇一種。
例子：single end：bwa samse hg19.fa read.fq.gz.sai read.fq.gz > read.fq.gz.sam
參數：-n int：如果reads比對次數超過多少次，就不在XA標簽顯示。
-r str：定義頭文件。‘@RG\tID:foo\tSM:bar’，如果在此步驟不進行頭文件定義，在
后續GATK分析中還是需要重新增加頭文件。
pair end：bwa sampe -a 500 read1.fq.gz.sai read2.fq.gz.sai read1.fq.gz read2.fq.gz > read.sam
參數：-a int：最大插入片段大小。
-o int：pair end兩reads中其中之一所允許配對的最大次數，超過該次數，將被視為
single end。降低這個參數，可以加快運算速度，對於少於30bp的read，建
議降低-o值。
-r str：定義頭文件。同single end。
-n int：每對reads輸出到結果中的最多比對數。
對於最后得到的sam文件，將比對上的結果提取出來（awk即可處理），即可直接用於GATK的分析。
注意：由於GATK在下游的snpcalling時，是按染色體進行callsnp的。因此，在准備原始sam文件時，可以先按染色體將文件分開，這樣會提高運行速度。但是當數據量不足時，可能會影響后續的VQSR分析，這是需要注意的。
2. 對sam文件進行進行重新排序（reorder）
由BWA生成的sam文件時按字典式排序法進行的排序（lexicographically）進行排序的（chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chrX，chrY），但是GATK在進行callsnp的時候是按照染色體組型（karyotypic）進行的（chrM，chr1，chr2…chr22，chrX，chrY），因此要對原始sam文件進行reorder。可以使用picard-tools中的ReorderSam完成。
e.g．
java -jar picard-tools-1.96/ReorderSam.jar
I=hg19.sam
O=hg19.reorder_00.sam
REFERENCE=hg19.fa
注意：
1. 這一步的頭文件可以人工加上，同時要確保頭文件中有的序號在下面序列中也有對應的。雖然在GATK網站上的說明chrM可以在最前也可以在最后，但是當把chrM放在最后時可能會出錯。
2. 在進行排序之前，要先構建參考序列的索引。
e.g. samtools faidx hg19.fa。最后生成的索引文件：hg19.fa.fai。
3. 如果在上一步想把大文件切分成小文件的時候，頭文件可以自己手工加上，之后運行這一步就好了。
3. 將sam文件轉換成bam文件（bam是二進制文件，運算速度快）
這一步可使用samtools view完成。
e.g. samtools view -bS hg19.reorder_00.sam -o hg19.sam_01.bam
4. 對bam文件進行sort排序處理
這一步是將sam文件中同一染色體對應的條目按照坐標順序從小到大進行排序。可以使用picard-tools中SortSam完成。
e.g.
java -jar picard-tools-1.96/SortSam.jar
INPUT=hg19.sam_01.bam
OUTPUT=hg19.sam.sort_02.bam
SORT_ORDER=coordinate
5. 對bam文件進行加頭（head）處理
GATK2.0以上版本將不再支持無頭文件的變異檢測。加頭這一步可以在BWA比對的時候進行，通過-r參數的選擇可以完成。如果在BWA比對期間沒有選擇-r參數，可以增加這一步驟。可使用picard-tools中AddOrReplaceReadGroups完成。
e.g.
java -jar picard-tools-1.96/AddOrReplaceReadGroups.jar
I=hg19.sam.sort_02.bam
O=hg19.reorder.sort.addhead_03.bam
ID=hg19ID
LB=hg19ID
PL=illumine
PU=hg19PU
SM=hg19
ID str：輸入reads集ID號；LB：read集文庫名；PL：測序平台（illunima或solid）；PU：測序平台下級單位名稱（run的名稱）；SM：樣本名稱。
注意：這一步盡量不要手動加頭，本人嘗試過多次手工加頭，雖然看起來與軟件加的頭是一樣的，但是程序卻無法運行。
6. Merge
如果一個樣本分為多個lane進行測序，那么在進行下一步之前可以將每個lane的bam文件合並。
e.g.
java -jar picard-tools-1.70/MergeSamFiles.jar
INPUT=lane1.bam
INPUT=lane2.bam
INPUT=lane3.bam
INPUT=lane4.bam
……
INPUT=lane8.bam
OUTPUT=sample.bam
7. Duplicates Marking
在制備文庫的過程中，由於PCR擴增過程中會存在一些偏差，也就是說有的序列會被過量擴增。這樣，在比對的時候，這些過量擴增出來的完全相同的序列就會比對到基因組的相同位置。而這些過量擴增的reads並不是基因組自身固有序列，不能作為變異檢測的證據，因此，要盡量去除這些由PCR擴增所形成的duplicates，這一步可以使用picard-tools來完成。去重復的過程是給這些序列設置一個flag以標志它們，方便GATK的識別。還可以設置 REMOVE_DUPLICATES=true 來丟棄duplicated序列。對於是否選擇標記或者刪除，對結果應該沒有什么影響，GATK官方流程里面給出的例子是僅做標記不刪除。這里定義的重復序列是這樣的：如果兩條reads具有相同的長度而且比對到了基因組的同一位置，那么就認為這樣的reads是由PCR擴增而來，就會被GATK標記。
e.g.
java -jar picard-tools-1.96/MarkDuplicates.jar
REMOVE_DUPLICATES= false
MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000
INPUT=hg19.reorder.sort.addhead_03.bam
OUTPUT=hg19.reorder.sort.addhead.dedup_04.bam METRICS_FILE=hg19.reorder.sort.addhead.dedup_04.metrics
注意： dedup這一步只要在library層面上進行就可以了，例如一個sample如果建了多個庫的話，對每個庫進行dedup即可，不需要把所有庫合成一個sample再進行dedup操作。其實並不能准確的定義被mask的reads到底是不是duplicates，重復序列的程度與測序深度和文庫類型都有關系。最主要目的就是盡量減小文庫構建時引入文庫的PCR bias。
8. 要對上一步得到的結果生成索引文件
可以用samtools完成，生成的索引后綴是bai。
e.g.
samtools index hg19.reorder.sort.addhead.dedup_04.bam
9.Local realignment around indels
這一步的目的就是將比對到indel附近的reads進行局部重新比對，將比對的錯誤率降到最低。一般來說，絕大部分需要進行重新比對的基因組區域，都是因為插入/缺失的存在，因為在indel附近的比對會出現大量的鹼基錯配，這些鹼基的錯配很容易被誤認為SNP。還有，在比對過程中，比對算法對於每一條read的處理都是獨立的，不可能同時把多條reads與參考基因組比對來排錯。因此，即使有一些reads能夠正確的比對到indel，但那些恰恰比對到indel開始或者結束位置的read也會有很高的比對錯誤率，這都是需要重新比對的。Local realignment就是將由indel導致錯配的區域進行重新比對，將indel附近的比對錯誤率降到最低。
主要分為兩步：
第一步，通過運行RealignerTargetCreator來確定要進行重新比對的區域。
e.g.
java -jar GenomeAnalysisTK.jar
-R hg19.fa
-T RealignerTargetCreator
-I hg19.reorder.sort.addhead.dedup_04.bam
-o hg19.dedup.realn_06.intervals
-known Mills_and_1000G_gold_standard.indels.hg19.vcf
-known 1000G_phase1.indels.hg19.vcf
參數說明：
-R： 參考基因組；
-T： 選擇的GATK工具；
-I： 輸入上一步所得bam文件；
-o： 輸出的需要重新比對的基因組區域結果；
-maxInterval： 允許進行重新比對的基因組區域的最大值，不能太大，太大耗費會很長時間，默認值500；
-known： 已知的可靠的indel位點，指定已知的可靠的indel位點，重比對將主要圍繞這些位點進行，對於人類基因組數據而言，可以直接指定GATK resource bundle里面的indel文件（必須是vcf文件）。
對於known sites的選擇很重要，GATK中每一個用到known sites的工具對於known sites的使用都是不一樣的，但是所有的都有一個共同目的，那就是分辨真實的變異位點和不可信的變異位點。如果不提供這些known sites的話，這些統計工具就會產生偏差，最后會嚴重影響結果的可信度。在這些需要知道known sites的工具里面，只有UnifiedGenotyper和HaplotypeCaller對known sites沒有太嚴格的要求。
如果你所研究的對象是人類基因組的話，那就簡單多了，因為GATK網站上對如何使用人類基因組的known sites做出了詳細的說明，具體的選擇方法如下表，這些文件都可以在GATK resource bundle中下載。
| Tool |
dbSNP 129 |
dbSNP >132 |
Mills indels |
1KG indels |
HapMap |
Omni |
| RealignerTargetCreator |
|
|
X |
X |
|
|
| IndelRealigner |
|
|
X |
X |
|
|
| BaseRecalibrator |
|
X |
X |
X |
|
|
| (UnifiedGenotyper/ HaplotypeCaller) |
|
X |
|
|
|
|
| VariantRecalibrator |
|
X |
X |
|
X |
X |
| VariantEval |
X |
|
|
|
|
|
但是如果你要研究的不是人類基因組的話，那就有點麻煩了，http://www.broadinstitute.org/gatk/guide/article?id=1243，這個網站上是做非人類基因組時，大家分享的經驗，可以參考一下。這個known sites如果實在沒有的話，也是可以自己構建的：首先，先使用沒有經過矯正的數據進行一輪SNP calling；然后，挑選最可信的SNP位點進行BQSR分析；最后，在使用這些經過BQSR的數據進行一次真正的SNP calling。這幾步可能要重復好多次才能得到可靠的結果。
第二步，通過運行IndelRealigner在這些區域內進行重新比對。
e.g.
java -jar GenomeAnalysisTK.jar
-R hg19.fa
-T IndelRealigner
-targetIntervals hg19.dedup.realn_06.intervals
-I hg19.reorder.sort.addhead.dedup_04.bam
-o hg19.dedup.realn_07.bam
-known Mills_and_1000G_gold_standard.indels.hg19.vcf
-known 1000G_phase1.indels.hg19.vcf
運行結束后，生成的hg19.dedup.realn_07.bam即為最后重比對后的文件。
注意：1. 第一步和第二步中使用的輸入文件（bam文件）、參考基因組和已知indel文件必須是相同的文件。
2. 當在相同的基因組區域發現多個indel存在時，這個工具會從其中選擇一個最有可能存在比對錯誤的indel進行重新比對，剩余的其他indel不予考慮。
3. 對於454下機數據，本工具不支持。此外，這一步還會忽略bwa比對中質量值為0的read以及在CIGAR信息中存在連續indel的reads。
10.Base quality score recalibration
這一步是對bam文件里reads的鹼基質量值進行重新校正，使最后輸出的bam文件中reads中鹼基的質量值能夠更加接近真實的與參考基因組之間錯配的概率。這一步適用於多種數據類型，包括illunima、solid、454、CG等數據格式。在GATK2.0以上版本中還可以對indel的質量值進行校正，這一步對indel calling非常有幫助
舉例說明，在reads鹼基質量值被校正之前，我們要保留質量值在Q25以上的鹼基，但是實際上質量值在Q25的這些鹼基的錯誤率在1%，也就是說質量值只有Q20，這樣就會對后續的變異檢測的可信度造成影響。還有，在邊合成邊測序的測序過程中，在reads末端鹼基的錯誤率往往要比起始部位更高。另外，AC的質量值往往要低於TG。BQSR的就是要對這些質量值進行校正。
BQSR主要有三步：
第一步：利用工具BaseRecalibrator，根據一些known sites，生成一個校正質量值所需要的數據文件，GATK網站以“.grp”為后綴命名。
e.g.
java -jar GenomeAnalysisTK.jar
-T BaseRecalibrator
-R hg19.fa
-I ChrALL.100.sam.dedup.realn.07.bam
-knownSites dbsnp_137.hg19.vcf
-knownSites Mills_and_1000G_gold_standard.indels.hg19.vcf
-knownSites 1000G_phase1.indels.hg19.vcf
-o ChrALL.100.sam.recal.08-1.grp
第二步：利用第一步生成的ChrALL.100.sam.recal.08-1.grp來生成校正后的數據文件，也是以“.grp”命名，這一步主要是為了與校正之前的數據進行比較，最后生成鹼基質量值校正前后的比較圖，如果不想生成最后BQSR比較圖，這一步可以省略。
e.g.
java -jar GenomeAnalysisTK.jar
-T BaseRecalibrator
-R hg19.fa
-I ChrALL.100.sam.dedup.realn.07.bam
-BQSR ChrALL.100.sam.recal.08-1.grp
-o GATK/hg19.recal.08-2.grp
-knownSites dbsnp_137.hg19.vcf
-knownSites Mills_and_1000G_gold_standard.indels.hg19.vcf
-knownSites 1000G_phase1.indels.hg19.vcf
第三步：利用工具PrintReads將經過質量值校正的數據輸出到新的bam文件中，用於后續的變異檢測。
e.g.
java -jar GenomeAnalysisTK.jar
-T PrintReads
-R hg19.fa
-I ChrALL.100.sam.dedup.realn.07.bam
-BQSR ChrALL.100.sam.recal.08-1.grp
-o ChrALL.100.sam.recal.08-3.grp.bam
主要參數說明：
-bqsrBAQGOP：BQSR BAQ gap open 罰值，默認值是40，如果是對全基因組數據進行BQSR分析，設置為30會更好。
-lqt： 在計算過程中，該算法所能考慮的reads兩端的最小質量值。如果質量值小於該值，計算過程中將不予考慮，默認值是2。
注意：（1）當bam文件中的reads數量過少時，BQSR可能不會正常工作，GATK網站建議base數量要大於100M才能得到比較好的結果。
（2）除非你所研究的樣本所得到的reads數實在太少，或者比對結果中的mismatch基本上都是實際存在的變異，否則必須要進行BQSR這一步。對於人類基因組，即使有了dbSNP和千人基因組的數據，還有很多mismatch是錯誤的。因此，這一步能做一定要做。
11. 分析和評估BQSR結果
這一步會生成評估前后鹼基質量值的比較結果，可以選擇使用圖片和表格的形式展示。
e.g.
java -jar GenomeAnalysisTK.jar
-T AnalyzeCovariates
-R hg19.fa
-before ChrALL.100.sam.recal.08-1.grp
-after ChrALL.100.sam.recal.08-2.grp
-csv ChrALL.100.sam.recal.grp.09.csv
-plots ChrALL.100.sam.recal.grp.09.pdf
參數解釋：
-before： 基於原始比對結果生成的第一次校對表格。
-after： 基於第一次校對表格生成的第二次校對表格。
-plots： 評估BQSR結果的報告文件。
-csv： 生成報告中圖標所需要的所有數據。
12.Reduce bam file
這一步是使用ReduceReads這個工具將bam文件進行壓縮，生成新的bam文件，新的bam文件仍然保持bam文件的格式和所有進行變異檢測所需要的信息。這樣不僅能夠節省存儲空間，也方便后續變異檢測過程中對數據的處理。
e.g.
java -jar GenomeAnalysisTK.jar
-T ReduceReads
-R hg19.fa
-I ChrALL.100.sam.recal.08-3.grp.bam
-o ChrALL.100.sam.recal.08-3.grp.reduce.bam
到此為止，GATK流程中的第一大步驟就結束了，完成了variants calling所需要的所有准備工作，生成了用於下一步變異檢測的bam文件。
第二大步：變異檢測
1. Variant Calling
GATK在這一步里面提供了兩個工具進行變異檢測——UnifiedGenotyper和HaplotypeCaller。其中HaplotypeCaller一直還在開發之中，包括生成的結果以及計算模型和命令行參數一直在變動，因此，目前使用比較多的還是UnifiedGenotyper。此外，HaplotypeCaller不支持Reduce之后的bam文件，因此，當選擇使用HaplotypeCaller進行變異檢測時，不需要進行Reduce reads。
UnifiedGenotyper是集合多種變異檢測方法而成的一種Variants Caller，既可以用於單個樣本的變異檢測，也可以用於群體的變異檢測。UnifiedGenotyper使用貝葉斯最大似然模型，同時估計基因型和基因頻率，最后對每一個樣本的每一個變異位點和基因型都會給出一個精確的后驗概率。
e.g.
java -jar GenomeAnalysisTK.jar
-glm BOTH
-l INFO
-R hg19.fa
-T UnifiedGenotyper
-I ChrALL.100.sam.recal.08-3.grp.reduce.bam
-D dbsnp_137.hg19.vcf
-o ChrALL.100.sam.recal.10.vcf
-metrics ChrALL.100.sam.recal.10.metrics
-stand_call_conf 10
-stand_emit_conf 30
上述命令將對輸入的bam文件中的所有樣本進行變異檢測，最后生成一個vcf文件，vcf文件中會包含所有樣本的變異位點和基因型信息。但是現在所得到的結果是最原始的、沒有經過任何過濾和校正的Variants集合。這一步產生的變異位點會有很高的假陽性，尤其是indel，因此，必須要進行進一步的篩選過濾。這一步還可以指定對基因組的某一區域進行變異檢測，只需要增加一個參數 -L：target_interval.list，格式是bed格式文件。
主要參數解釋：
-A： 指定一個或者多個注釋信息，最后輸出到vcf文件中。
-XA： 指定不做哪些注釋，最后不會輸出到vcf文件中。
-D： 已知的snp文件。
-glm： 選擇檢測變異的類型。SNP表示只進行snp檢測；INDEL表示只對indel進行檢測；BOTH表示同時檢測snp和indel。默認值是SNP。
-hets： 雜合度的值，用於計算先驗概率。默認值是0.001。
-maxAltAlleles： 容許存在的最大alt allele的數目，默認值是6。這個參數要特別注意，不要輕易修改默認值，程序設置的默認值幾乎可以滿足所有的分析，如果修改了可能會導致程序無法運行。
-mbq： 變異檢測時，鹼基的最小質量值。如果小於這個值，將不會對其進行變異檢測。這個參數不適用於indel檢測，默認值是17。
-minIndelCnt： 在做indel calling的時候，支持一個indel的最少read數量。也就是說，如果同時有多少條reads同時支持一個候選indel時，軟件才開始進行indel calling。降低這個值可以增加indel calling的敏感度，但是會增加耗費的時間和假陽性。
-minIndelFrac： 在做indel calling的時候，支持一個indel的reads數量占比對到該indel位置的所有reads數量的百分比。也就是說，只有同時滿足-minIndelCnt和-minIndelFrac兩個參數條件時，才會進行indel calling。
-onlyEmitSamples：當指定這個參數時，只有指定的樣本的變異檢測結果會輸出到vcf文件中。
-stand_emit_conf：在變異檢測過程中，所容許的最小質量值。只有大於等於這個設定值的變異位點會被輸出到結果中。
-stand_call_conf：在變異檢測過程中，用於區分低質量變異位點和高質量變異位點的閾值。只有質量值高於這個閾值的位點才會被視為高質量的。低於這個質量值的變異位點會在輸出結果中標注LowQual。在千人基因組計划第二階段的變異檢測時，利用35x的數據進行snp calling的時候，當設置成50時，有大概10%的假陽性。
-dcov： 這個參數用於控制檢測變異數據的coverage(X)，4X的數據可以設置為40，大於30X的數據可以設置為200。
注意：GATK進行變異檢測的時候，是按照染色體排序順序進行的（先call chr1，然后chr2，然后chr3…最后chrY），並非多條染色體並行檢測的，因此，如果數據量比較大的話，建議分染色體分別進行，對性染色體的變異檢測可以同常染色體方法。
大多數參數的默認值可以滿足大多數研究的需求，因此，在做變異檢測過程中，如果對參數意義不是很明確，不建議修改。
2. 對原始變異檢測結果進行過濾（hard filter and VQSR）
這一步的目的就是對上一步call出來的變異位點進行過濾，去掉不可信的位點。這一步可以有兩種方法，一種是通過GATK的VariantFiltration，另一種是通過GATK的VQSR（變異位點質量值重新校正）進行過濾。
通過GATK網站上提供的最佳方案可以看出，GATK是推薦使用VASR的，但使用VQSR數據量一定要達到要求，數據量太小無法使用高斯模型。還有，在使用VAQR時，indel和snp要分別進行。
VQSR原理介紹：
這個模型是根據已有的真實變異位點（人類基因組一般使用HapMap3中的位點，以及這些位點在Omni 2.5M SNP芯片中出現的多態位點）來訓練，最后得到一個訓練好的能夠很好的評估真偽的錯誤評估模型，可以叫他適應性錯誤評估模型。這個適應性的錯誤評估模型可以應用到call出來的原始變異集合中已知的變異位點和新發現的變異位點，進而去評估每一個變異位點發生錯誤的概率，最終會給出一個得分。這個得分最后會被寫入vcf文件的INFO信息里，叫做VQSLOD，就是在訓練好的混合高斯模型下，一個位點是真實的概率比上這個位點可能是假陽性的概率的log odds ratio（對數差異比），因此，可以定性的認為，這個值越大就越好。
VQSR主要分兩個步驟，這兩個步驟會使用兩個不同的工具：VariantRecalibrator和ApplyRecalibration。
VariantRecalibrator：通過大量的高質量的已知變異集合的各個注釋（包括很多種，后面介紹）的值來創建一個高斯混合模型，然后用於評估所有的變異位點。這個文件最后將生成一個recalibration文件。
原理簡單介紹： 這個模型首先要拿到真實變異數據集和上一步驟中得到的原始變異數據集的交集，然后對這些SNP值相對於具體注釋信息的分布情況進行模擬，將這些變異位點進行聚類，最后根據聚類結果賦予所有變異位點相應的VQSLOD值。越接近聚類核心的變異位點得到的VQSLOD值越高。
ApplyRecalibration：這一步將模型的各個參數應用於原始vcf文件中的每一個變異位點，這時，每一個變異位點的注釋信息列中都會出現一個VQSLOD值，然后模型會根據這個值對變異位點進行過濾，過濾后的信息會寫在vcf文件的filter一列中。
原理簡單介紹： 在VariantRecalibrator這一步中，每個變異位點已經得到了一個VQSLOD值了，同時，這些LOD值在訓練集里也進行了排序。當你在這一步中設置一個tranche sensitivity 的閾值（這個閾值一般是一個百分數，如設置成99%），那么，如果LOD值從大到小排序的話，這個程序就會認為在這個訓練集中，LOD值在前99%的是可信的，當這個值低於這個閾值，就認為是錯誤的。最后，程序就會用這個標准來過濾上一步call出來的原始變異集合。如果LOD值超過這個閾值，在filter那一列就會顯示PASS，如果低於這個值就會被過濾掉，但是這些位點仍然會顯示在結果里面，只不過會在filter那一列標示出他所屬於的tranche sensitivity 的名稱。在設置tranche sensitivity 的閾值時，要兼顧敏感度和質量值。
對高斯混合模型生成圖片的解釋：
在VariantRecalibrator這一步，程序會通過已知位點來訓練概率模型，訓練完成后會生成一組圖片，而且每對注釋信息都對應一組圖片（上圖），這組圖片能夠幫助我們理解一個概率模型是否與我們的數據相匹配，也就是說這個模型能不能很好的區分假陽性和真實位點。
上圖是第一步完成后生成的一個報告的一部分，圖中只表示了一對注釋所對應的圖。左上角的圖表示的是適合當前數據的概率密度圖，綠色區域表示高質量變異位點所在位置，紅色區域表示低質量概率分布區域。如果變異位點分布在紅色區域，則會被過濾掉。右上角圖中紅色的點表示在經過VQSR之后被過濾掉的變異位點，黑色的表示的是留下來的。紅色的表示的都是沒有達到所設定的tranche sensitivity 閾值的點。左下角的圖表示的是用來訓練模型的點，綠色的點表示通過訓練進入到ApplyRecalibration的變異位點，紫色的點則表示質量值很低的，沒有達到質量要求的點。右下角的圖表示的是已知的和新發現的變異位點的分布，紅色的點表示新發現的變異位點，而藍色的點表示的是已知的變異位點，看這幅圖就是看這兩個注釋信息能不能很好的區分已知的點（大部分是真實的）和未知的點（大部分是假陽性）。
從圖中可以看出，這個模擬結果可以很好的將真實的變異位點和假陽性變異位點分開（左下圖），形成了明顯的界限，也就是說，如果一個變異位點的這兩個注釋值，只要有一個落在了界限之外，就會被過濾掉。最主要的是要看右邊兩個圖片，只要能很好的區分開novel和known以及filtered和retained就可以。其實在如何選擇注釋值存在一定得主觀性，因此，在做VariantRecalibrator時可以做兩次，第一次盡可能的多的選擇這些注釋值，第一遍跑完之后，選擇幾個區分好的，再做一次VariantRecalibrator，然后再做ApplyRecalibration。具體每個注釋值得意義可以參考：
http://www.broadinstitute.org/gatk/guide/tagged?tag=annotation這個網址中的內容，有每個注釋的詳細信息的鏈接。
tranche值的設定
前面提到了，這個值得設定是用來在后續的ApplyRecalibration中如何根據這個閾值來過濾變異位點的，也就是說，如果這個值設定的比較高的話，那么最后留下來的變異位點就會多，但同時假陽性的位點也會相應增加；如果設定的低的話，雖然假陽性會減少，但是會丟失很多真實的位點。因此，跟選擇注釋時一樣，可以run兩遍VariantRecalibrator，第一遍的時候多寫幾個閾值，第一遍跑完之后看結果，看那個閾值好，選擇一個最好的閾值，再run一遍VariantRecalibrator。至於說怎么區分好壞，有幾個標准：
1. 看結果中已知變異位點與新發現變異位點之間的比例，這個比例不要太大，因為大多數新發現的變異都是假陽性，如果太多的話，可能假陽性的比例就比較大；
2. 看保留的變異數目，這個就要根據具體的需求進行選擇了。
3. 看TI/TV值，對於人類全基因組，這個值應該在2.15左右，對於外顯子組，這個值應該在3.2左右，不要太小或太大，越接近這個數值越好，這個值如果太小，說明可能存在比較多的假陽性。
千人中所選擇的tranche值是99，僅供參考。
注意：Indel不支持tranche值的選擇，另外，一部分注釋類型在做indel的校正時也不支持，具體信息可以詳查GATK網站。
當數據量太小時，可能高斯模型不會運行，因為變異位點數滿足不了模型的統計需求。這時候可以通過降低--maxGaussian的值，讓程序運行。這個值表示的是程序將變異位點分成的最大的組數，降低這個值讓程序把變異位點聚類到更少的組里面，使每個組中的變異位點數增加來滿足統計需求，但是這樣做降低程序分辨真偽的能力。因此，在運行程序的時候，要對各方面進行權衡。
e.g.
:對SNP結果進行校正
第一步：
java -jar GenomeAnalysisTK.jar
-R hg19.fa
--maxGaussians 4
-numBad 10000 (這個參數在最新的GATK版本里面已經沒有了，用的時候注意版本，2.8.1里面不用自己設置
這個參數)
-T VariantRecalibrator
-mode SNP
-input ChrALL.100.sam.recal.10.vcf
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.hg19.vcf
-resource:omni,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.hg19.vcf
-resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg19.vcf
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_137.hg19.vcf
-an QD
-an HaplotypeScore
-an MQRankSum
-an ReadPosRankSum
-an FS
-an MQ
-an InbreedingCoeff
-recalFile hg19.vcf.snp_11_Q10.recal
-tranchesFile hg19.vcf.snp_11_Q10.tranches
-rscriptFile hg19.vcf.snp_11.plot_Q10.R
-nt 4
--TStranche 90.0
--TStranche 93.0
--TStranche 95.0
--TStranche 97.0
--TStranche 99.0
--TStranche 99.9
先run一下上面的代碼，這一步可以盡可能多的設置注釋類型和tranche的值，然后根據這次跑出來的結果選擇出最好的注釋類型和tranche值之后，再次運行VariantRecalibrator。
第二步：
java -jar GenomeAnalysisTK.jar
-R hg19.fa
--maxGaussians 4
-numBad 10000
-T VariantRecalibrator
-mode SNP
-input ChrALL.100.sam.recal.10.vcf
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.hg19.vcf
-resource:omni,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.hg19.vcf
-resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg19.vcf
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_137.hg19.vcf
-an HaplotypeScore
-an MQRankSum
--TStranche 97.0
-recalFile hg19.vcf.snp_11_Q10.recal
-tranchesFile hg19.vcf.snp_11_Q10.tranches
-rscriptFile hg19.vcf.snp_11.plot_Q10.R
-nt 4
這一步run出來的結果可以直接用於下一步的ApplyRecalibration。
第三步
java -jar GenomeAnalysisTK.jar
-R hg19.fa
-T ApplyRecalibration
-mode SNP
-input hg19.recal_10_Q10.vcf
-tranchesFile hg19.vcf.snp_12_Q10-2.tranches
-recalFile hg19.vcf.snp_12_Q10-2.recal
-o hg19.snp.filter.t97.Q10_13.snp.vcf
--ts_filter_level 97
最終生成的hg19.snp.filter.t97.Q10_13.snp.vcf這個文件中的SNP位點已經全部經過校正過濾，INDEL位點還是原始數據，需要對INDEL再進行一次校正過濾。
對INDEL結果進行校正，與SNP基本一致，只不過INDEL需要使用的known resource不一樣
第一步：
同SNP 多選擇一些注釋類型，但是不用選擇tranche值,tranche值是專門為SNP設定的，即使設定
了這個值（2.4版本是可以計算這個的，以后就不計算了），計算出來也都是錯的，這個在indel
里不需要考慮。
第二步：
java -jar GenomeAnalysisTK.jar
-T VariantRecalibrator
-R hg19.fa
-mode INDEL
--maxGaussians 4
-std 10.0
-percentBad 0.12
-input ChrALL.100.sam.recal.10.vcf
-resource:mills,known=true,training=true,truth=true,prior=12.0 Mills_and_1000G_gold_standard.indels.hg19
-an MQ
-an FS
-an InbreedingCoeff
-recalFile ChrALL.100.sam.recal.10.indel.recal
-tranchesFile ChrALL.100.sam.recal.10.indel.tranche
-rscriptFile ChrALL.100.sam.recal.10.indel.R
第三步：
java -jar GenomeAnalysisTK.jar
-T ApplyRecalibration
-R hg19.fa
-mode INDEL
-input hg19.snp.filter.t97.Q10_13.snp.vcf
-recalFile ChrALL.100.sam.recal.11.indel.recal
-tranchesFile ChrALL.100.sam.recal.11.indel.tranche
-o hg19.snp.filter.t97.Q10_13.both.vcf
最后得到的hg19.snp.filter.t97.Q10_13.both.vcf文件，就是我們最終想得到的過濾好的變異集合。
主要參數解釋：
VariantRecalibrator
-badLodCutoff 當LOD得分低於這個值的時候，就用於構建高斯混合模型的bad variants。默認值是-5。
-maxNumTrainingData 構建高斯模型過程中，用於訓練的最大位點數目。如果超過這個數目，將被隨機刪除。默認值是2500000。
-minNumBad 構建高斯模型的bad variants時的最少低質量值得位點數。
-recalFile 用於ApplyRecalibration的輸出文件。
-resource 已知的變異信息。
-rscriptFile 結果中生成圖片的腳本。
-tranchesFile 用於ApplyRecalibration的tranche結果輸出文件。
-tranche 設置tranche閾值。
-an 選擇填加注釋信息。
更多其他參數參考：
ApplyRecalibration
-ef 輸出結果中不顯示被過濾掉的位點。
-lodCutoff VQSLOD值低於這個值就過濾掉。
-recalFile 上一步生成的recalFile。
-tranchesFile 上一步生成的tranchesFile。
-ts_filter_level 上一步中確定的tranche值。
另外，關於如何選擇resource data可以參考：
http://www.broadinstitute.org/gatk/guide/article?id=1259
如果要分析的數據集不符合進行VQSR的標准，可以進行hard filter，這一步將使用GATK中的VariantFiltration工具來完成。具體使用方法參考：
最后生成的vcf文件的格式說明，即每一列所代表的的內容，可參考下面的網站，有詳細的說明：
http://www.broadinstitute.org/gatk/guide/article?id=1268
其實到這里就已經完成了變異檢測的所有步驟，最后生成的hg19.snp.filter.t97.Q10_13.both.vcf文件就可以用於你的下游分析了。
第三大步：初步分析
這一步主要是對上面所得到的最終vcf中的結果進行一些初步的分析，比如計算這些變異位點在dbsnp中的比例、Ti/Tv的比例、每個樣本中的snp數量……。此外，還可以對變異位點的同義/非同義突變進行統計，識別是否為CpG位點以及氨基酸的簡並信息等。這一步主要是利用GATK中的VariantEval來完成。
需要注意的是，有些計算內容不能同時進行，例如AlleleCount和VariantSummary或者Sample和VariantSummary。如果選擇了這樣的組合方式，程序就會報錯。但是GATK並沒有告訴我們到底哪些不能同時運行，所以當選擇計算內容的時候可以先做一下測試。
e.g.
java -jar GenomeAnalysisTK.jar
-R hg19.fa
-T VariantEval
--eval hg19.snp.filter.t97.Q10_13.both.vcf
-D dbsnp_137.hg19.vcf
-o hg19.PASS.Eval_15_Final.gatkreport
主要參數解釋：
--eval 輸入要進行summary的文件，也就是hg19.snp.filter.t97.Q10_13.both.vcf。
-EV 選擇模塊計算相應的分析內容，。
--list 列出可供選擇的計算模塊。
-noEV 不是用默認的模塊，只計算用-EV選定的模塊。
更多其他參數請參考：
以上就是GATK整個流程的詳細介紹了，我敢保證，如果你是第一次使用GATK，就算按照這個說明一步一步去做，也肯定會遇到很多問題，因為GATK需要注意的細節真的很多，我在這里面還有很多沒有提到，遇到問題的時候可以仔細看一下報錯信息，一般情況下是可以通過報錯信息知道錯在哪里。遇到問題的時候可以多瀏覽GATK網站，里面的FAQ基本上可以包括所有出現過的問題的解決方法了，可以耐心查一下。要是不想查可以在論壇上直接發起提問，管理員真的會很快給你回復的。不像同類的XXX網站只有問題沒有回答，沒人用也是可以理解的……。祝好運！！！
此時版本 GATK 2.8.1
原文鏈接：http://blog.sina.com.cn/s/blog_a0f8a3060102wapk.html