GATK这个软件在做snp-calling的时候使用率非常高，因为之前一直是简单粗略的看看snp情况而已，所以没有具体研究它。
这些天做一些外显子项目以找snp为重点，所以想了想还是用起它，报错非常多，调试了好久才成功。
所以记录一些注意事项!
GATK软件本身是受版权保护的，所以需要申请才能下载使用，大家自己去broad institute申请即可。
下载软件就可以直接使用，java软件不需要安装，但是需要你的机器上面有java，当然软件只是个开始，重点是你还得下载很多配套数据，https://software.broadinstitute.org/gatk/download/bundle（ps:这个链接可能会失效，下面的文件，请自己谷歌找到地址哈。），而且这个时候要明确你的参考基因组版本了！！！ b36/b37/hg18/hg19/hg38，记住b37和hg19并不是完全一样的，有些微区别哦！！！
比如我选择了hg19
第一点是hg19的下载：这个下载地址非常多，常用的就是NCBI，ensembl和UCSC了，但是这里推荐用这个脚本下载
for i in $(seq 1 22) X Y M;
do echo $i;
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg19.fasta;
done
rm -fr chr*.fasta
看得懂shell脚本的应该知道这是一个个的下载hg19的染色体，再用cat按照染色体的顺序拼接起来，因为GATK后面的一些步骤对染色体顺序要求非常变态，如果下载整个hg19，很难保证染色体顺序是1-22，X,Y,M。如下
然后需要对下载的hg19进行索引（bwa和samtools）和建立dict文件（用picard）
bwa index -a bwtsw hg19.fasta
samtools faidx hg19.fasta
然后还要下载几个参考文件，这个是可以选择的.
对我的hg19来说，就应该是去，ftp://ftp.broadinstitute.org/bundle/hg19/ 下载咯。
最后，所有必备的文件如下：
231M Jul 2 05:14 1000G_phase1.indels.hg19.sites.vcf
1.2M Jul 2 10:45 1000G_phase1.indels.hg19.sites.vcf.idx
11G Jul 2 08:05 dbsnp_138.hg19.vcf
2.5K Jul 1 04:31 hg19.dict
3.0G Jun 30 21:29 hg19.fasta
6.6K Jun 30 22:54 hg19.fasta.amb
944 Jun 30 22:54 hg19.fasta.ann
2.9G Jun 30 22:54 hg19.fasta.bwt
788 Jul 2 01:53 hg19.fasta.fai
739M Jun 30 22:54 hg19.fasta.pac
1.5G Jun 30 23:23 hg19.fasta.sa
87M Jul 2 05:37 Mills_and_1000G_gold_standard.indels.hg19.sites.vcf
2.3M Jul 2 10:45 Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.idx
接下来开始跑程序
第一步就是生成sam文件啦bwa mem -t 12 -M hg19.fasta tmp*fq >tmp.sam
第二步是sort，我用的是picard工具java -Xmx100g -jar AddOrReplaceReadGroups.jar I=tmp.sam O=tmp.sorted.bam
SORT_ORDER=coordinate
CREATE_INDEX=true
RGID=tmp
RGLB="pe"
RGPU="HiSeq-2000"
RGSM=PC3-2
RGCN="Human Genetics of Infectious Disease"
RGDS=hg19 RGPL=illumina
VALIDATION_STRINGENCY=SILENT
第三步是去除PCR重复，我还是选择用picard工具
java -Xmx100g -jar MarkDuplicates.jar
CREATE_INDEX=true REMOVE_DUPLICATES=True
ASSUME_SORTED=True VALIDATION_STRINGENCY=LENIENT
I=tmp.sorted.bam OUTPUT=tmp.dedup.bam METRICS_FILE=tmp.metrics
第四步是终于要开始用GATK啦，主要是确定要进行重新比对的区域，这个步骤分成三个小步骤：
首先用RealignerTargetCreator找到需要重新比对的区域，输出文件intervals
java -Xmx200g -jar ~/apps/gatk/GenomeAnalysisTK.jar
-R hg19.fasta #这里需要用这个参考基因组，所以参考基因组特别重要，DICT也要按照流程生成
-T RealignerTargetCreator
-I tmp.dedup.bam -o tmp.intervals
-known /home/ldzeng/EXON/ref/1000G_phase1.indels.hg19.sites.vcf
这一步骤好像非常耗时
可以看到，我总共就测试了5014个reads，结果就花了近半个小时才搞定，只有947个reads被过滤了。
输出的tmp.intervals 文件是一个1404946行的文件
chr1:13957-13958
chr1:46402-46403
chr1:47190-47191
chr1:52185-52188
chr1:53234-53236
chr1:55249-55250
chr1:63735-63738
人的外显子只有二三十万，所以我暂时也不确定这个文件是什么！
然后用输出的 tmp.intervals 做输入文件来进行重新比对，也就是用IndelRealigner在这些区域内进行重新比对
java -Xmx150g -jar ~/apps/gatk/GenomeAnalysisTK.jar \
-R hg19.fasta \
-T IndelRealigner \
-targetIntervals tmp.intervals \
-I tmp.dedup.bam -o tmp.dedup.realgn.bam \
-known /home/ldzeng/EXON/ref/1000G_phase1.indels.hg19.sites.vcf
我只需要它的重新比对，所以后面的一些功能没有怎么用，一个是call snp，一个是算比对质量值
java -Xmx200g -jar ~apps/gatk/GenomeAnalysisTK.jar
-nct 20 -T HaplotypeCaller -R hg19.fasta
-I tmp.dedup.realgn.bam
-o tmp.gatk.vcf
最后输出的文件如下
639K Jul 5 10:17 tmp1.fq
639K Jul 5 10:19 tmp2.fq
1.5M Jul 5 10:26 tmp.dedup.bai
403K Jul 5 10:26 tmp.dedup.bam
12K Jul 5 12:02 tmp.gatk.vcf
3.4K Jul 5 12:02 tmp.gatk.vcf.idx
32M Jul 5 11:24 tmp.intervals
950 Jul 5 10:26 tmp.metrics
1.5M Jul 5 11:31 tmp.realgn.bai
409K Jul 5 11:31 tmp.realgn.bam
1.6M Jul 5 10:20 tmp.sam
1.5M Jul 5 10:23 tmp.sorted.bai
399K Jul 5 10:23 tmp.sorted.bam
备注：GATK对基因组要求一个字典文件
使用picard工具包的CreateSequenceDictionary.jar生成。以hg19.fa为例，生成的命令为：