GATK使用方法详解plob最详尽说明指导书.doc

GATK使用方法详解
一、使用GATK前须知事项：
（1）对GATK测试关键使用是人类全基因组和外显子组测序数据，而且全部是基于illumina数据格式，现在还没有提供其它格式文件（如Ion Torrent）或试验设计（RNA-Seq）分析方法。
（2）GATK是一个应用于前沿科学研究软件，不停在更新和修正，所以，在使用GATK进行变异检测时，最好是下载最新版本，现在版本是2.8.1（-02-25）。下载网站：。
（3）在GATK使用过程中（见下面图），有些步骤需要用到已知变异信息，对于这些已知变 异，GATK只提供了人类已知变异信息，能够在GATKFTP站点下载（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。
构建索引时需要注意问题：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 -t 4 -I > read1.fq.gz.sai
bwa aln hg19.fa read2.fq.gz -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’端出现大于多少bpdeletion。
-i int：不许可在reads两端出现大于多少bpindel。
-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格式。这是一个很奇怪功效，就是对其它软件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。降低这个参数，能够加紧运算速度，对于少于30bpread，提议降低-o值。
-r str：定义头文件。同single end。
-n int：每对reads输出到结果中最多比对数。
对于最终得到sam文件，将比对上结果提取出来（awk即可处理），即可直接用于GATK分析。
注意：因为GATK在下游snp-calling时，是按染色体进行call-snp。所以，在准备原始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完成。
eg.
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进行测序，那么在进行下一步之前能够将每个lanebam文件合并。
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操作。其实并不能正确定义被maskreads到底是不是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

不过假如你要研究不是人类基因组话，那就有点麻烦了，，这个网站上是做非人类基因组时，大家分享经验，能够参考一下。这个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比对中质量值为0read和在CIGAR信息中存在连续indelreads。
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时候，支持一个indelreads数量占比对到该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出来变异位点进行过滤，去掉不可信位点。这一步能够有两种方法，一个是经过GATKVariantFiltration，另一个是经过GATKVQSR（变异位点质量值重新校正）进行过滤。
经过GATK网站上提供最好方案能够看出，GATK是推荐使用VASR，但使用VQSR数据量一定要达成要求，数据量太小无法使用高斯模型。还有，在使用VAQR时，indel和snp要分别进行。
VQSR原理介绍：
这个模型是依据已经有真实变异位点（人类基因组通常使用HapMap3中位点，和这些位点在Omni 2.5M SNP芯片中出现多态位点）来训练，最终得到一个训练好能够很好评定真伪错误评定模型，能够叫她适应性错误评定模型。这个适应性错误评定模型可 以应用到call出来原始变异集合中已知变异位点和新发觉变异位点，进而去评定每一个变异位点发生错误概率，最终会给出一个得分。这个得分最终会 被写入vcf文件INFO信息里，叫做VQSLOD，就是在训练好混合高斯模型下，一个位点是真实概率比上这个位点可能是假阳性概率log odds ratio（对数差异比），所以，能够定性认为，这个值越大就越好。
VQSR关键分两个步骤，这两个步骤会使用两个不一样工具：VariantRecalibrator和ApplyRecalibration。
i) VariantRecalibrator
VariantRecalibrator：经过大量高质量已知变异集合各个注释（包含很多个，后面介绍）值来创建一个高斯混合模型，然后用于评定全部变异位点。这个文件最终将生成一个recalibration文件。
原理简单介绍： 这个模型首先要拿到真实变异数据集和上一步骤中得到原始变异数据集交集，然后对这些SNP值相对于具体注释信息分布情况进行模拟，将这些变异位点进 行聚类，最终依据聚类结果给予全部变异位点对应VQSLOD值。越靠近聚类关键变异位点得到VQSLOD值越高。

ii) ApplyRecalibration
ApplyRecalibration：这一步将模型各个参数应用于原始vcf文件中每一个变异位点，这时，每一个变异位点注释信息列中全部会出现一个VQSLOD值，然后模型会依据这个值对变异位点进行过滤，过滤后信息会写在vcf文件filter一列中。
原理简单介绍： 在VariantRecalibrator这一步中，每个变异位点已经得到了一个VQSLOD值了，同时，这些LOD值在训练集里也进行了排序。当你在这 一步中设置一个tranche sensitivity 阈值（这个阈值通常是一个百分数，如设置成99%），那么，假如LOD值从大到小排序话，这个程序就会认为在这个训练集中，LOD值在前99%是可 信，当这个值低于这个阈值，就认为是错误。最终，程序就会用这个标准来过滤上一步call出来原始变异集合。假如LOD值超出这个阈值，在 filter那一列就会显示PASS，假如低于这个值就会被过滤掉，不过这些位点仍然会显示在结果里面，只不过会在filter那一列标示出她所属于 tranche sensitivity 名称。在设置tranche sensitivity 阈值时，要兼顾敏感度和质量值。
初步分析
这一步关键是对上面所得到最终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选定模块。