GATK之HaplotypeCaller
GATK的主要功能其实就是识别变异位点，其他功能都是锦上添花。所以这一次学习GATK寻找变异位点的工具。
在GATK的文档中，与变异位点识别相关的有9个工具，分别是：
| Name | Summary |
|---|---|
| Apply a score cutoff to filter variants based on a recalibration table | |
| Calculate genotype posterior likelihoods given panel data | |
| Simple Bayesian genotyper used in the original GATK paper | |
| Perform joint genotyping on gVCF files produced by HaplotypeCaller | |
| Call germline SNPs and indels via local re-assembly of haplotypes | |
| Call somatic SNPs and indels via local re-assembly of haplotypes | |
| Regenotypes the variants from a VCF containing PLs or GLs. | |
| Call SNPs and indels on a per-locus basis | |
| Build a recalibration model to score variant quality for filtering purposes |
稍微看了一下GATK最佳实践文档，发现目前用的最多的是HaplotypeCaller，准确率高，但是比较吃配置，所以运行时间会比较久。不过由于HaplotypeCaller的工作原理，直接省去了BQSR和indel realignment步骤，所以对于一个variant calling流程而言，可以直接比对，去重复后运行HaplotypeCaller。
简介
功能： 通过局部单倍型重组装寻找体细胞的SNP和indel
分类： 变异位点探索功能
HaplotypeCaller，简称HC，能过通过对活跃区域（也就是与参考基因组不同处较多的区域）局部重组装，同时寻找SNP和INDEL。换句话说，当HC看到一个地方好活跃呀，他就不管之前的比对结果了，直接对这个地方进行重新组装，所以就比传统的基于位置（position-based)的工具，如前任UnifiedGenotyper好用多了。
HC可以处理非二倍体物种和混池实验数据，但是不推荐用于癌症或者体细胞。
HC还可以正确处理可变剪切，所以也能用于RNA-Seq数据。
HC有一个GVCF模式，产生各个样本的中间基因组文件(gVCF),然后用于多样本的联合基因定型(joint genotyping)，效率很高呢。
工作原理
HaplotypeCaller的核心操作就是四步：
寻找活跃区域，就是和参考基因组不同部分较多的区域
通过对该区域进行局部重组装，确定单倍型（haplotypes）。就是这一步可以省去indel realignment
在给定的read数据下，计算单倍型的可能性。
分配样本的基因型
使用说明
输入：BAM文件
输出：原始的VCF文件或者gVCF文件，未过滤SNP和INDEL。一般而言，在进行下游分析，需要对这些数据要么手动要么VQSR过滤。
案例:
DNA-seq variant-calling
java -jar GenomeAnalysisTK.jar \
-R reference.fasta \
-T HaplotypeCaller \
-I sample1.bam [-I sample2.bam ...] \
[--dbsnp dbSNP.vcf] \
[-stand_call_conf 30] \
[-L targets.interval_list] \
-o output.raw.snps.indels.vcfRNA-seq variant-calling
java -jar GenomeAnalysisTK.jar \
-R reference.fasta \
-T HaplotypeCaller \
-I sample1.bam \
[--dbsnp dbSNP.vcf] \
-stand_call_conf 20 \
-o output.raw.snps.indels.vcf其他感觉比较使用的参数：
| 参数名 | 默认值 | 概要 |
|---|---|---|
| —genotyping_mode/-gt_mode | DISCOVERY | 如何处理识别的等位基因 |
| —min_base_quality_score/ -mbq | 10 | 最低碱基质量 |
| —sample_ploidy/-ploidy | 2 | 样本是几倍体？ |
| —standard_min_confidence_threshold_for_calling/-stand_call_conf | 10.0 | 确定variant的最低phred-scaled 置信阈值 |
| —annotation/-A | - | variant注释内容 |
| —pcr_indel_model/-pcrModel | CONSERVATIVE | PCR indel模式 |
注：对于我做mapping-by-sequencing而言，需要结果有ref和alt碱基的支持数，所以选项-A一定要跟上StrandAlleleCountsBySample。
这是我用于MBS的流程的GATK部分
##############################
# Vaiant Discovery with GATK #
##############################mkdir -p ${wd_dir}/analysis/variant_gatk
recal_files=$(ls *recal_reads.bam)
for file in $recal_files
do
java -Xmx16g -jar $gatk -T HaplotypeCaller \
-R $reference -I $file --genotyping_mode DISCOVERY \
--standard_min_confidence_threshold_for_calling 10 \
-A StrandAlleleCountsBySample \
-o ../variant_gatk/${file%%.*}_raw_variants.vcf
done更多参数详细说明见，可以阅读原文。
写在最后
GATK best practice写在2015年，现在是2017年，GATK的版本到3.7，下一个版本就是4.0了，所以很多工具都不能想当然的用了。比如说之前写的BQSR，有一篇文献就说已经没有必要了，文献名见下：
《 impact of post-alignment processing in variant discovery from whole exome data>
而这个也是我发完BQSR，留言的大神告诉我的。后来我去biostar搜索，也发现有人说BQSR和VQSR都不是必须的（我正在在问回答者原文出处）
技术的发展越来越快，对于一个刚入门的bioinformatician，很多东西都要从头开始学习，压力非常大。所以我非常感谢一个论坛—技能树，非常感谢一本书—biostar handbook.还有我们生信媛的小伙伴们。
没事多去逛biostar，seqanswers和生信技能树，开拓自己的眼界。也希望那些老司机们不要吝啬自己的知识，能够带带我们。
最后，让我们“一起学习，共同进步”