GATK使用方法详解-plob最详尽说明书(2)

2018-12-09 23:57

注意:这一步尽量不要手动加头,本人尝试过多次手工加头,虽然看起来与软件加的头是一样的,但是程序却无法运行。

6. Merge

如果一个样本分为多个lane进行测序,那么在进行下一步之前可以将每个lane的bam文件合并。 e.g.

java -jar picard-tools-1.70/MergeSamFiles.jar INPUT=lane1.bam INPUT=lane2.bam INPUT=lane3.bam INPUT=lane4.bam ……

INPUT=lane8.bam OUTPUT=sample.bam

7. Duplicates Marking

在制备文库的过程中,由于PCR扩增过程中会存在一些偏差,也就是说有的序列会被过量扩增。这样,在比对的时候,这些过量扩增出来的完全相同的序列 就会比对到基因组的相同位置。而这些过量扩增的reads并不是基因组自身固有序列,不能作为变异检测的证据,因此,要尽量去除这些由PCR扩增所形成的 duplicates,这一步可以使用picard-tools来完成。去重复的过程是给这些序列设置一个flag以标志它们,方便GATK的识别。还可以设置

REMOVE_DUPLICATES=true 来丢弃duplicated序列。对于是否选择标记或者删除,对结果应该没有什么影响,GATK官方流程里面给出的例子是仅做标记不删除。这里定义的重复序列是这样的:如果两条reads具有相同的长度而且比对到了基因组的同一位置,那么就认为这样的reads是由PCR扩增而来,就会被GATK标记。 e.g.

java -jar picard-tools-1.96/MarkDuplicates.jar REMOVE_DUPLICATES= false

MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 INPUT=hg19.reorder.sort.addhead_03.bam

OUTPUT=hg19.reorder.sort.addhead.dedup_04.bam

METRICS_FILE=hg19.reorder.sort.addhead.dedup_04.metrics

注意: dedup这一步只要在library层面上进行就可以了,例如一个sample如果建了多个库的话,对每个库进行dedup即可,不需要把所有库合成一 个sample再进行dedup操作。其实并不能准确的定义被mask的reads到底是不是duplicates,重复序列的程度与测序深度和文库类型 都有关系。最主要目的就是尽量减小文库构建时引入文库的PCR bias。

8. 要对上一步得到的结果生成索引文件 可以用samtools完成,生成的索引后缀是bai。 e.g.

samtools index hg19.reorder.sort.addhead.dedup_04.bam 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中下载。 Mills dbSNP 129 dbSNP >132 indels RealignerTargetCreator X IndelRealigner X Tool 1KG indels HapMap Omni X X BaseRecalibrator (UnifiedGenotyper/ HaplotypeCaller) VariantRecalibrator VariantEval

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


GATK使用方法详解-plob最详尽说明书(2).doc 将本文的Word文档下载到电脑 下载失败或者文档不完整,请联系客服人员解决!

下一篇:理力典型习题

相关阅读
本类排行
× 注册会员免费下载(下载后可以自由复制和排版)

马上注册会员

注:下载文档有可能“只有目录或者内容不全”等情况,请下载之前注意辨别,如果您已付费且无法下载或内容有问题,请联系我们协助你处理。
微信: QQ: