【原创】GATK使用方法详解（包含bwa使用）第二部分
2014-03-03 11:16阅读：
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来生成校正后的数据文
