变异检测（BWA+SAMtools+picard+GATK）

一、简介

首先说说GATK可以做什么。它主要用于从sequencing 数据中进行variant calling，包括SNP、INDEL。比如现在风行的exome sequencing找variant，一般通过BWA+GATK的pipeline进行数据分析。
要run GATK，首先得了解它的网站(http://www.broadinstitute.org/gatk/)。

使用GATK前须知事项：
（1）对GATK的测试主要使用的是人类全基因组和外显子组的测序数据，而且全部是基于illumina数据格式，目前还没有提供其他格式文件（如Ion Torrent）或者实验设计（RNA-Seq）的分析方法。
（2）GATK是一个应用于前沿科学研究的软件，不断在更新和修正，因此，在使用GATK进行变异检测时，最好是下载最新的版本，目前的版本是3.5（2016-05）。下载网站：https://www.broadinstitute.org/gatk/download/auth?package=GATK，下载可能会需要注册！

二、GATK的使用流程

GATK最佳使用方案：共3大步骤。原始数据的处理—变异检测—初步分析。
方案一：使用GATK检测变异，官方推荐的流程：

方案二：在这里我们使用的流程可能更复杂一点，这也是目前很多科研团队使用的流程，即在官方流程上结合SAMtools工具的一套流程。

三、准备工作

3.1 软件下载安装

bwa：http://iweb.dl.sourceforge.net/project/bio-bwa/bwa-0.7.12.tar.bz2
picard：http://broadinstitute.github.io/picard/
samtools：https://github.com/samtools/samtools/archive/0.1.19.tar.gz（这里使用0.1.19版本）
GATK：https://www.broadinstitute.org/gatk/download/auth?package=GATK（需要注册登录）
R语言；

java；

3.2 软件安装

其中picard tools和GATK，解压到某个目录即可，无需安装。打开解压缩后的目录，会发现GenomeAnalysisTK.jar这个文件。我们要用的各种工具都在这个包里。务必下载最新版（当然只是建议）。
BWA、SAMtools、R、java的安装我就不再这里说了，我前面的博文已经讲过BWA、和SAMtools可以参考，还有R和java参考网上，也许我以后会加进来。

3.3 参考基因组和annotation下载

这里主要是针对需要用到已知变异信息的情况。对于这些已知变异，GATK只提供了人类的已知变异信息，可以在GATK的FTP站点下载（GATK resource bundle）。
GATK resource bundle介绍：这里
GATK resource bundle FTP地址：这里
例如呢，以hg19(human genome build 19.此外，b37也很常用)为参考序列,
输入如下命令： lftp ftp.broadinstitute.org -u gsapubftp-anonymous
回车后键入空格，便可进入resource bundle。进入其中名为bundle的目录，找到最新版的hg19目录。
要下载的文件包括：

ucsc.hg19.fasta
ucsc.hg19.fasta.fai
1000G_omni2.5.hg19.vcf
1000G_omni2.5.hg19.vcf.idx
1000G_phase1.indels.hg19.vcf
1000G_phase1.indels.hg19.vcf.idx
dbsnp_137.hg19.vcf
dbsnp_137.hg19.vcf.idx
hapmap_3.3.hg19.vcf
hapmap_3.3.hg19.vcf.idx
Mills_and_1000G_gold_standard.indels.hg19.vcf
Mills_and_1000G_gold_standard.indels.hg19.vcf.idx

当然，如果要研究的不是人类基因组，在这里我们要研究的物种，并没有vcf注释信息，这时可以自己去找一些该物种已有的vcf注释文件，如果没有的话，需要自行构建已知变异，GATK提供了详细的构建方法。

3.4 R设置

GATK在进行BQSR和VQSR的过程中会使用到R软件绘制一些图，因此，在运行GATK之前最好先检查一下是否正确安装了R和所需要的包，所需要的包大概包括ggplot2、gplots、bitops、caTools、colorspace、gdata、gsalib、reshape、RColorBrewer等。如果画图时出现错误，会提示需要安装的包的名称。

如何check这些packages有没有安装？以ggplot2为例，进入R，输入library(ggplot2)，如果没有error信息，就表明安装了ggplot2；如果没有安装，输入命令install.packages(“ggplot2”)，按照提示操作便可安装ggplot2。如果package安装不成功，很可能是权限问题：默认使用的R是系统管理员安装的，而你在/usr等目录下没有write权限。这是可以选择直接在自己目录下安装了R，然后install各种packages。也可以在 install.packages时指定将packages安装到自己的目录。

将R packages安装到自己目录的方法：
install.packages(“ggplot2”, lib=”/your/own/R-packages/“)
然后将指定的路径添加到R_LIBS变量即可。

gsalib的安装有点不一样。
从该网站下载gsalib：https://github.com/broadgsa/gatk
解压缩后，进入build.xml文件所在的目录，输入命令ant gsalib。这样就可以安装gsalib 了。如果该gsalib仍然无法load，可能需要将gsllib的path告诉R。
这里有较详细的步骤。（奇怪的是我没有~/.Rprofle文件，ant完后也没有add path，仍然work了）
此外，gaslib安装要求java的版本为1.6。

设置路径：
export R_HOME=/your/R/install/path
export R_LIBS=/your/R-packages/path

如何查看当前使用的R的路径？输入 which R即可。

四、GATK使用

4.1 原始数据处理

(1). 对原始下机fastq文件进行过滤和比对（mapping)

对于Illumina下机数据推荐使用bwa进行mapping。
Bwa比对步骤大致如下：

(a) 对参考基因建立索引
eg:

1	$ bwa index ref.fa

构建索引时需要注意的问题：bwa构建索引有两种算法，两种算法都是基于BWT的，这两种算法通过参数-a is 和-a bwtsw进行选择。其中-a bwtsw对于短的参考序列是不工作的，必须要大于等于10Mb；-a is是默认参数，这个参数不适用于大的参考序列，必须要小于等于2G。

(b) BWA Alignment
详细的命令参数，参考我的博文：BWA命令详解其中bwa以及后面使用的一系列tools都能直接读取gzip压缩文件。为了节省硬盘空间，以及减少压缩/解压缩过程的read/write时间，最好直接用*.gz文件作为input。output也最好先在管道里使用gzip压缩一下再输出。
eg:

1	$ bwa mem -t 20 ref.fa read1.fq read2.fq > sample.sam

(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完成。
Picard的SortSam需指定一个tmp目录，用于存放中间文件，中间文件会很大，above 10G.注意指定目录的空间大小。
eg:

$ java  -jar $PICARD/SortSam.jar \
        REFERENCE=ref.fa \
        INPUT=sample.sam \
        OUTPUT=sample.reorder.sam \

注意：

这一步的头文件可以人工加上，同时要确保头文件中有的序号在下面序列中也有对应的。虽然在GATK网站上的说明chrM可以在最前也可以在最后，但是当把chrM放在最后时可能会出错。

在进行排序之前，要先构建参考序列的索引。
e.g：samtools faidx hg19.fa ;
最后生成的索引文件：hg19.fa.fai。

如果在上一步想把大文件切分成小文件的时候，头文件可以自己手工加上，之后运行这一步就好了。

(3). 将sam文件转换成bam文件（bam是二进制文件，运算速度快）

这一步可使用samtools view完成。如果上一步直接使用的bam后缀，则不需要再转换。
e.g:

1	$ samtools view -bS sample.reorder.sam -o sample.reorder.bam

(4). 对bam文件进行sort排序处理

这一步是将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序。可以使用picard-tools中SortSam完成。
e.g:

$ java  -jar picard-tools/SortSam.jar \
        INPUT=sample.reorder.bam \ #这里也可以直接用sample.sort.sam,则输出为bam格式，可以省去前面samtools的转换步骤
        OUTPUT=sample.sort.bam \ 
        SORT_ORDER=coordinate

注意：
SAMtools也提供了工具进行sam file的sort和bam file的生成，并且不会生成大的中间文件。但是不建议在GATK做SNP calling等步骤的时候使用samtools的排序方式进行sort。

(5). 对bam文件进行加头（head）处理

GATK2.0以上版本将不再支持无头文件的变异检测。加头这一步可以在BWA比对的时候进行，通过-R参数的选择可以完成。如果在BWA比对期间没有选择-R参数，可以增加这一步骤。可使用picard-tools中AddOrReplaceReadGroups完成。
e.g:

$ java  -jar picard-tools/AddOrReplaceReadGroups.jar \
        I=sample.sort.bam \
        O=sample.sort.addhead.bam \
        ID=sampleID \
        LB=sampleID \
        PL=illumina \
        PU=samplePU \
        SM=sample

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_all.bam

(7). Duplicates Marking

在制备文库的过程中，由于PCR扩增过程中会存在一些偏差，也就是说有的序列会被过量扩增。这样，在比对的时候，这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。而这些过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由PCR扩增所形成的duplicates，这一步可以使用picard-tools来完成。去重复的过程是给这些序列设置一个flag以标志它们，方便GATK的识别。还可以设置 REMOVE_DUPLICATES=true 来丢弃duplicated序列。对于是否选择标记或者删除，对结果应该没有什么影响，GATK官方流程里面给出的例子是仅做标记不删除。这里定义的重复序列是这样的：如果两条reads具有相同的长度而且比对到了基因组的同一位置，那么就认为这样的reads是由PCR扩增而来，就会被GATK标记。
e.g:

$ java  -Xmx10g -jar $PICARD/MarkDuplicates.jar \
        I=sample.sort.addhead.bam \
        O=sample.rmdup.bam \
        VALIDATION_STRINGENCY=LENIENT \
        MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000 \
        REMOVE_DUPLICATES= false \
        M=sample.sort.addhead.rmdup.metric

注意：
dedup这一步只要在library层面上进行就可以了，例如一个sample如果建了多个库的话，对每个库进行dedup即可，不需要把所有库合成一个sample再进行dedup操作。其实并不能准确的定义被mask的reads到底是不是duplicates，重复序列的程度与测序深度和文库类型都有关系。最主要目的就是尽量减小文库构建时引入文库的PCR bias。

(8). 要对上一步得到的结果生成索引文件

可以用samtools完成，生成的索引后缀是bai。
e.g:

1	$ samtools index sample.rmdup.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 ref.fa \
        -T RealignerTargetCreator \
        -I sample.rmdup.bam \
        -o sample.realign.intervals \
        -known known.indels.vcf #如果没有已知变异位点的信息，暂时可以不用，那么这一步的结果用来后面自己构建know sites

参数说明：

-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

第二步，通过运行IndelRealigner在这些区域内进行重新比对。
e.g:

$ java  -jar GenomeAnalysisTK.jar \
        -R ref.fa \
        -T IndelRealigner \
        -targetIntervals sample.realign.intervals \
        -I sample.rmdup.bam \
        -o sample.realign.bam \
        -known known.indels.vcf #如果没有已知变异位点的信息，暂时可以不用，那么这一步的结果用来后面自己构建know sites

运行结束后，生成的sample.realign.bam即为最后重比对后的文件。

注意：

第一步和第二步中使用的输入文件（bam文件）、参考基因组和已知indel文件必须是相同的文件。

当在相同的基因组区域发现多个indel存在时，这个工具会从其中选择一个最有可能存在比对错误的indel进行重新比对，剩余的其他indel不予考虑。

对于454下机数据，本工具不支持。此外，这一步还会忽略bwa比对中质量值为0的read以及在CIGAR信息中存在连续indel的reads。

GATK+samtools构建know sites
但是如果你要研究的不是人类基因组的话，那就有点麻烦了，这个网站上是做非人类基因组时，大家分享的经验，可以参考一下。这个known sites如果实在没有的话，也是可以自己构建的：首先，先使用没有经过矫正的数据进行一轮SNP calling；然后，挑选最可信的SNP位点进行BQSR分析；最后，在使用这些经过BQSR的数据进行一次真正的SNP calling。这几步可能要重复好多次才能得到可靠的结果。
第一步：
使用GATK检测变异：
eg:

$ java  -jar GenomeAnalysisTK.jar \
        -R ref.fa \
        -T UnifiedGenotyper \
        -I sample.realign.bam \
        -o sample.gatk.raw.vcf \
        -dcov 200 -A QualByDepth -A FisherStrand -A AlleleBalance -A Coverage -A MappingQualityZero -A TandemRepeatAnnotator -baq CALCULATE_AS_NECESSARY -stand_call_conf 30.0 -stand_emit_conf 10.0 \
        -glm BOTH -rf BadCigar \
        -log ./log/sample.gatk.raw.vcf.log

第二步:
使用Samtools检测变异：

1 2	$ samtools index sample.realign.bam $ samtools mpileup -DSugf $SOYBEAN_REF sample.realign.bam \| bcftools view -cvg - > sample.samtools.raw.vcf

第三步:
结合GATK和SAMtools的vcf结果：

$ java  -jar GenomeAnalysisTK.jar \
        -R ref.fa \
        -T SelectVariants \
        --variant sample.gatk.raw.vcf \
        --concordance sample.samtools.raw.vcf \
        -o sample.concordance.raw.vcf \
        -log ./log/sample.concordance.raw.vcf.log

第四步：
过滤，挑选最可信的SNP位点：

$ MEANQUAL=$(awk '{ if ($1 !~ /#/) { total += $6; count++  }  } END { print total/count  }' sample.concordance.raw.vcf)
$ java  -jar GenomeAnalysisTK.jar \
        -R ref.fa \
        -T VariantFiltration \
        --filterExpression "QD < 20.0 || ReadPosRankSum < -8.0 || FS > 10.0 || QUAL < $MEANQUAL" --filterName LowQualFilter \
        --missingValuesInExpressionsShouldEvaluateAsFailing \
        --variant sample.concordance.raw.vcf \
        --logging_level ERROR -o sample.concordance.flt.vcf \
        -log ./log/sample.concordance.flt.vcf.log
$ grep -v "Filter" sample.concordance.flt.vcf > sample.concordance.flt.knowsites.vcf

这样我们就得到了可信度比较高的know sites。

(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”为后缀命名。
eg:

$ java  -jar GenomeAnalysisTK.jar \
        -T BaseRecalibrator \
        -R ref.fa \
        -I sample.realign.bam \
        -knownSites sample.concordance.flt.knowsites.vcf \
        -o sample.recal_data.grp \
        -log ./log/sample.recal_data_1.grp.log

第二步：
利用第一步生成的ChrALL.100.sam.recal.08-1.grp来生成校正后的数据文件，也是以“.grp”命名，这一步主要是为了与校正之前的数据进行比较，最后生成碱基质量值校正前后的比较图，如果不想生成最后BQSR比较图，这一步可以省略。
eg:

$ java  -jar GenomeAnalysisTK.jar \
        -T BaseRecalibrator \
        -R ref.fa \
        -I sample.realign.bam \
        -BQSR sample.recal_data.grp \
        -o sample.recal_data_2.grp \
        -knownSites sample.concordance.flt.knowsites.vcf

第三步：
利用工具PrintReads将经过质量值校正的数据输出到新的bam文件中，用于后续的变异检测。
eg:

$ java  -jar GenomeAnalysisTK.jar \
        -R ref.fa \
        -T PrintReads \
        -I sample.realign.bam \
        -o sample.recal.bam \
        -BQSR sample.recal_data_1.grp \
        -log ./log/sample.recal.bam.log

主要参数说明：

1
2

-bqsrBAQGOP:    BQSR BAQ gap open 罚值，默认值是40，如果是对全基因组数据进行BQSR分析，设置为30会更好。
-lqt:           在计算过程中，该算法所能考虑的reads两端的最小质量值。如果质量值小于该值，计算过程中将不予考虑，默认值是2。

注意：

当bam文件中的reads数量过少时，BQSR可能不会正常工作，GATK网站建议base数量要大于100M才能得到比较好的结果。

除非你所研究的样本所得到的reads数实在太少，或者比对结果中的mismatch基本上都是实际存在的变异，否则必须要进行BQSR这一步。对于人类基因组，即使有了dbSNP和千人基因组的数据，还有很多mismatch是错误的。因此，这一步能做一定要做。

(11). 分析和评估BQSR结果

这一步会生成评估前后碱基质量值的比较结果，可以选择使用图片和表格的形式展示。
eg:

$ java  -jar GenomeAnalysisTK.jar \
        -T AnalyzeCovariates \
        -R ref.fa \
        -before sample.recal_data_1.grp \
        -after sample.recal_data_2.grp \
        -csv sample.recal.grp.csv \
        -plots sample.recal.grp.pdf

参数解释：

-before：   基于原始比对结果生成的第一次校对表格。
-after：    基于第一次校对表格生成的第二次校对表格。
-plots：    评估BQSR结果的报告文件。
-csv：      生成报告中图标所需要的所有数据。

(12).Reduce bam file

这一步是使用ReduceReads这个工具将bam文件进行压缩，生成新的bam文件，新的bam文件仍然保持bam文件的格式和所有进行变异检测所需要的信息。这样不仅能够节省存储空间，也方便后续变异检测过程中对数据的处理。
eg:

$ java  -jar GenomeAnalysisTK.jar \
        -T ReduceReads \
        -R hg19.fa \
        -I sample.recal.bam \
        -o sample.recal.reduce.bam

到此为止，GATK流程中的第一大步骤就结束了，完成了variants calling所需要的所有准备工作，生成了用于下一步变异检测的bam文件。

4.2 变异检测

(1). SNP Calling

e.g:

$ java  -jar GenomeAnalysisTK.jar \
        -R ref.fa \
        -T UnifiedGenotyper \
        -I sample.recal.bam \
        -o sample.snp.raw.vcf \
        -dcov 200 \
        -A QualByDepth -A FisherStrand -A AlleleBalance -A Coverage -A MappingQualityZero -A TandemRepeatAnnotator \
        -baq CALCULATE_AS_NECESSARY \
        -stand_call_conf 30.0 \
        -stand_emit_conf 10.0 \
        -glm SNP \
        -log ./log/sample.snp.raw.vcf.log

主要参数解释：

-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。

(2). InDel calling

eg:

$ java  -jar GenomeAnalysisTK.jar \
        -R ref.fa \
        -T UnifiedGenotyper \
        -I sample.recal.bam \
        -o sample.indel.raw.vcf \
        -dcov 200 \
        -A QualByDepth -A FisherStrand -A AlleleBalance -A Coverage -A MappingQualityZero -A TandemRepeatAnnotator \
        -baq CALCULATE_AS_NECESSARY \
        -stand_call_conf 30.0 \
        -stand_emit_conf 10.0 \
        -glm INDEL \
        -log ./log/sample.indel.raw.vcf.log

###附：###
以下均参考葡萄糖的博客，自己没有测试：
转接上面的4.2变异检测
变异检测
（1）Variant Calling
GATK在这一步里面提供了两个工具进行变异检测——UnifiedGenotyper和HaplotypeCaller。其中HaplotypeCaller一直还在开发之中，包括生成的结果以及计算模型和命令行参数一直在变动，因此，目前使用比较多的还是UnifiedGenotyper。此外，HaplotypeCaller不支持Reduce之后的bam文件，因此，当选择使用HaplotypeCaller进行变异检测时，不需要进行Reduce reads。

UnifiedGenotyper是集合多种变异检测方法而成的一种Variants Caller，既可以用于单个样本的变异检测，也可以用于群体的变异检测。UnifiedGenotyper使用贝叶斯最大似然模型，同时估计基因型和基因频率，最后对每一个样本的每一个变异位点和基因型都会给出一个精确的后验概率。

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格式文件。

注意：
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。具体每个注释值得意义可以参考：这个网址中的内容，有每个注释的详细信息的链接。

tranche值的设定
前面提到了，这个值得设定是用来在后续的ApplyRecalibration中如何根据这个阈值来过滤变异位点的，也就是说，如果这个值设定的比较高的话，那么最后留下来的变异位点就会多，但同时假阳性的位点也会相应增加；如果设定的低的话，虽然假阳性会减少，但是会丢失很多真实的位点。因此，跟选择注释时一样，可以run两遍VariantRecalibrator，第一遍的时候多写几个阈值，第一遍跑完之后看结果，看那个阈值好，选择一个最好的阈值，再run一遍VariantRecalibrator。至于说怎么区分好坏，有几个标准：

1. 看结果中已知变异位点与新发现变异位点之间的比例，这个比例不要太大，因为大多数新发现的变异都是假阳性，如果太多的话，可能假阳性的比例就比较大；
2. 看保留的变异数目，这个就要根据具体的需求进行选择了。
3. 看TI/TV值，对于人类全基因组，这个值应该在2.15左右，对于外显子组，这个值应该在3.2左右，不要太小或太大，越接近这个数值越好，这个值如果太小，说明可能存在比较多的假阳性。

千人中所选择的tranche值是99，仅供参考。

注意：
Indel不支持tranche值的选择，另外，一部分注释类型在做indel的校正时也不支持，具体信息可以详查GATK网站。
当数据量太小时，可能高斯模型不会运行，因为变异位点数满足不了模型的统计需求。这时候可以通过降低–maxGaussian的值，让程序运行。这个值表示的是程序将变异位点分成的最大的组数，降低这个值让程序把变异位点聚类到更少的组里面，使每个组中的变异位点数增加来满足统计需求，但是这样做降低程序分辨真伪的能力。因此，在运行程序的时候，要对各方面进行权衡。

eg:
对SNP结果进行校正
第一步：

$ java  -jar GenomeAnalysisTK.jar \
        -R ref.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可以参考：这里

如果要分析的数据集不符合进行VQSR的标准，可以进行hard filter，这一步将使用GATK中的VariantFiltration工具来完成。具体使用方法参考：这里

最后生成的vcf文件的格式说明，即每一列所代表的的内容，可参考下面的网站，有详细的说明：这里

其实到这里就已经完成了变异检测的所有步骤，最后生成的hg19.snp.filter.t97.Q10_13.both.vcf文件就可以用于你的下游分析了。

4.4 初步分析

这一步主要是对上面所得到的最终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选定的模块。

更多其他参数请参考：这里

参考：
糗世界：http://blog.qiubio.com:8080/archives/3207
葡萄糖的博客：http://blog.sina.com.cn/s/blog_12d5e3d3c0101qu6e.html
imitosis的博客：http://blog.sina.com.cn/s/blog_681aaa550101bhdc.html

一、 简介

二、GATK的使用流程

三、 准备工作

3.1 软件下载安装

3.2 软件安装

3.3 参考基因组和annotation下载

3.4 R设置

四、GATK使用

4.1 原始数据处理

(1). 对原始下机fastq文件进行过滤和比对（mapping)

(2). 对sam文件进行进行重新排序（reorder）

(3). 将sam文件转换成bam文件（bam是二进制文件，运算速度快）

(4). 对bam文件进行sort排序处理

(5). 对bam文件进行加头（head）处理

(6). Merge

(7). Duplicates Marking

(8). 要对上一步得到的结果生成索引文件

(9). Local realignment around indels

(10). Base quality score recalibration

(11). 分析和评估BQSR结果

(12).Reduce bam file

4.2 变异检测

(1). SNP Calling

(2). InDel calling

4.4 初步分析

一、简介

三、准备工作