搜索
查看: 16749|回复: 5

GATK之HaplotypeCaller

[复制链接]

18

主题

56

帖子

403

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
403
发表于 2017-5-3 18:48:34 | 显示全部楼层 |阅读模式
本帖最后由 hoptop 于 2017-5-3 18:57 编辑

GATK的主要功能其实就是识别变异位点,其他功能都是锦上添花。所以这一次学习GATK寻找变异位点的工具。
在GATK的文档中,与变异位点识别相关的有9个工具,分别是:

NameSummary
ApplyRecalibrationApply a score cutoff to filter variants based on a recalibration table
CalculateGenotypePosteriorsCalculate genotype posterior likelihoods given panel data
GATKPaperGenotyperSimple Bayesian genotyper used in the original GATK paper
GenotypeGVCFsPerform joint genotyping on gVCF files produced by HaplotypeCaller
HaplotypeCallerCall germline SNPs and indels via local re-assembly of haplotypes
MuTect2Call somatic SNPs and indels via local re-assembly of haplotypes
RegenotypeVariantsRegenotypes the variants from a VCF containing PLs or GLs.
UnifiedGenotyperCall SNPs and indels on a per-locus basis
VariantRecalibratorBuild 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:
  
[Bash shell] 纯文本查看 复制代码
 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.vcf


RNA-seq variant-calling:

[Bash shell] 纯文本查看 复制代码
   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_modeDISCOVERY如何处理识别的等位基因
—min_base_quality_score/ -mbq10最低碱基质量
—sample_ploidy/-ploidy2样本是几倍体?
—standard_min_confidence_threshold_for_calling/-stand_call_conf10.0确定variant的最低phred-scaled 置信阈值
—annotation/-A-variant注释内容
—pcr_indel_model/-pcrModelCONSERVATIVEPCR indel模式

注:对于我做mapping-by-sequencing而言,需要结果有ref和alt碱基的支持数,所以选项-A一定要跟上StrandAlleleCountsBySample。

这是我用于MBS的流程的GATK部分

[Bash shell] 纯文本查看 复制代码
############################ 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和生信技能树,开拓自己的眼界。也希望那些老司机们不要吝啬自己的知识,能够带带我们。
最后,让我们“一起学习,共同进步”









上一篇:在R里面操作SQLite
下一篇:reStructuredText(.rst)语法规则快速入门
回复

使用道具 举报

2

主题

21

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-8-11 22:33:17 | 显示全部楼层
是否可以将DNA-seq variant-calling流程精简为: 质控→bwa组装→Picard添加read groups,对bam文件排序,标记重复,创建索引→GATK的HaplotypeCaller进行Variant calling→过滤
回复 支持 反对

使用道具 举报

18

主题

56

帖子

403

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
403
 楼主| 发表于 2017-8-12 20:17:20 | 显示全部楼层
AdaWon 发表于 2017-8-11 22:33
是否可以将DNA-seq variant-calling流程精简为: 质控→bwa组装→Picard添加read groups,对bam文件排序,标 ...

我建议是用瓶中基因组计划跑不同的流程,比较一下结果,如果差异不大,那么在算力不够的情况下,就按照你的方法精简。
回复 支持 反对

使用道具 举报

2

主题

21

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-8-13 20:07:49 | 显示全部楼层
hoptop 发表于 2017-8-12 20:17
我建议是用瓶中基因组计划跑不同的流程,比较一下结果,如果差异不大,那么在算力不够的情况下,就按照你 ...

THANKS ^_^!!!
回复 支持 反对

使用道具 举报

0

主题

11

帖子

259

积分

中级会员

Rank: 3Rank: 3

积分
259
发表于 2018-3-1 17:46:19 | 显示全部楼层
但是在用HC call snp时,想用IGV结合bam文件来检查snp的可靠信时,这样是不是就无法检查了呢~~~
回复 支持 反对

使用道具 举报

0

主题

23

帖子

91

积分

注册会员

Rank: 2

积分
91
发表于 2018-7-24 15:21:34 | 显示全部楼层
我看 GATK 官方的最佳实践里面,BQSR 还是推荐的步骤之一
专注于 Spark 分布式快速处理基因数据。提问前请先搜索《提问的智慧》
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-19 12:47 , Processed in 0.035463 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.