加载中…
个人资料
葡萄糖
葡萄糖
  • 博客等级:
  • 博客积分:0
  • 博客访问:87,031
  • 关注人气:63
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

【原创】GATK使用方法详解(包含bwa使用)第二部分

(2014-03-03 11:16:18)
标签:

gatk

indel

snp

bwa

分类: 生物信息

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的工具里面,只有UnifiedGenotyperHaplotypeCallerknown 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比对中质量值为0read以及在CIGAR信息中存在连续indelreads

 

10.Base quality score recalibration

这一步是对bam文件里reads的碱基质量值进行重新校正,使最后输出的bam文件中reads中碱基的质量值能够更加接近真实的与参考基因组之间错配的概率。这一步适用于多种数据类型,包括illunimasolid454CG等数据格式。在GATK2.0以上版本中还可以对indel的质量值进行校正,这一步对indel calling非常有帮助

 

举例说明,在reads碱基质量值被校正之前,我们要保留质量值在Q25以上的碱基,但是实际上质量值在Q25的这些碱基的错误率在1%,也就是说质量值只有Q20,这样就会对后续的变异检测的可信度造成影响。还有,在边合成边测序的测序过程中,在reads末端碱基的错误率往往要比起始部位更高。另外,AC的质量值往往要低于TGBQSR的就是要对这些质量值进行校正。

 

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来生成校正后的数据文件,也是以“.grp”命名,这一步主要是为了与校正之前的数据进行比较,最后生成碱基质量值校正前后的比较图,如果不想生成最后BQSR比较图,这一步可以省略。

e.g.

java -jar GenomeAnalysisTK.jar

-T BaseRecalibrator

-R hg19.fa

-I ChrALL.100.sam.dedup.realn.07.bam

-BQSR ChrALL.100.sam.recal.08-1.grp

-o GATK/hg19.recal.08-2.grp

-knownSites dbsnp_137.hg19.vcf

-knownSites Mills_and_1000G_gold_standard.indels.hg19.vcf

-knownSites 1000G_phase1.indels.hg19.vcf

 

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

e.g.

java -jar GenomeAnalysisTK.jar

-T PrintReads

-R hg19.fa

-I ChrALL.100.sam.dedup.realn.07.bam

-BQSR ChrALL.100.sam.recal.08-1.grp

-o ChrALL.100.sam.recal.08-3.grp.bam

 

主要参数说明:

-bqsrBAQGOPBQSR BAQ gap open 罚值,默认值是40,如果是对全基因组数据进行BQSR分析,设置为30会更好。

-lqt                在计算过程中,该算法所能考虑的reads两端的最小质量值。如果质量值小于该值,计算过程中将不予考虑,默认值是2

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

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

 

11. 分析和评估BQSR结果

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

e.g.

java -jar GenomeAnalysisTK.jar

-T AnalyzeCovariates

-R hg19.fa

-before ChrALL.100.sam.recal.08-1.grp

-after ChrALL.100.sam.recal.08-2.grp

-csv ChrALL.100.sam.recal.grp.09.csv

-plots ChrALL.100.sam.recal.grp.09.pdf

 

参数解释:

-before 基于原始比对结果生成的第一次校对表格。

-after    基于第一次校对表格生成的第二次校对表格。

-plots    评估BQSR结果的报告文件。

-csv      生成报告中图标所需要的所有数据。

 

12.Reduce bam file

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

e.g.

java -jar GenomeAnalysisTK.jar

-T ReduceReads

-R hg19.fa

-I ChrALL.100.sam.recal.08-3.grp.bam

-o ChrALL.100.sam.recal.08-3.grp.reduce.bam

 

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

 

 

作者:葡萄糖

欢迎加入群:分子进化与基因组学

0

阅读 收藏 喜欢 打印举报/Report
  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有