GATK使用方法.docx

GATK使用方法.docx

《GATK使用方法.docx》由会员分享，可在线阅读，更多相关《GATK使用方法.docx（32页珍藏版）》请在冰豆网上搜索。

GATK使用方法

原创】GATK使用方法详解（包含bwa使用）第一部分

（2014-03-0311:

07:

29）

转载▼

标签：

gatk

bwa

snp

indel

分类：

 生物信息

由于新浪博客规定，每篇文章不可超过2万字符，因此分4篇发布。

一、使用GATK前须知事项：

（1）对GATK的测试主要使用的是人类全基因组和外显子组的测序数据，而且全部是基于illumina数据格式，目前还没有提供其他格式文件（如IonTorrent）或者实验设计（RNA-Seq）的分析方法。

（2）GATK是一个应用于前沿科学研究的软件，不断在更新和修正，因此，在使用GATK进行变异检测时，最好是下载最新的版本，目前的版本是2.8.1（2014-02-25）。

下载网站：

http:

//www.broadinstitute.org/gatk/download。

（3）在GATK使用过程中（见下面图），有些步骤需要用到已知变异信息，对于这些已知变异，GATK只提供了人类的已知变异信息，可以在GATK的FTP站点下载（GATKresourcebundle）。

如果要研究的不是人类基因组，需要自行构建已知变异，GATK提供了详细的构建方法。

（4）GATK在进行BQSR和VQSR的过程中会使用到R软件绘制一些图，因此，在运行GATK之前最好先检查一下是否正确安装了R和所需要的包，所需要的包大概包括ggplot2、gplots、bitops、caTools、colorspace、gdata、gsalib、reshape、RColorBrewer等。

如果画图时出现错误，会提示需要安装的包的名称。

二、GATK的使用流程

GATK最佳使用方案：

共3大步骤。

原始数据的处理—变异检测—初步分析。

第一大步：

原始数据的处理

1. 对原始下机fastq文件进行过滤和比对（mapping）

对于Illumina下机数据推荐使用bwa进行mapping。

Bwa比对步骤大致如下：

（1）对参考基因组构建索引：

     例子：

bwaindex-abwtswhg19.fa。

最后生成文件：

hg19.fa.amb、hg19.fa.ann、hg19.fa.bwt、hg19.fa.pac和hg19.fa.sa。

     构建索引时需要注意的问题：

bwa构建索引有两种算法，两种算法都是基于BWT的，这两种算法通过参数-ais 和-abwtsw进行选择。

其中-abwtsw对于短的参考序列是不工作的，必须要大于等于10Mb；-ais是默认参数，这个参数不适用于大的参考序列，必须要小于等于2G。

（2）寻找输入reads文件的SA坐标。

     对于pairend数据，每个reads文件单独做运算，singleend数据就不用说了，只有一个文件。

     例子：

pairend：

bwa  aln  hg19.fa  read1.fq.gz  -l30  -k2  -t4  -I  >read1.fq.gz.sai

bwa  aln  hg19.fa  read2.fq.gz  -l30  -k2  -t4  -I  >read2.fq.gz.sai

singleend：

bwa  aln  hg19.fa  read.fq.gz  -l30  -k2  -t4  -I  >read.fq.gz.sai

主要参数说明：

-oint：

允许出现的最大gap数。

-eint：

每个gap允许的最大长度。

-dint：

不允许在3’端出现大于多少bp的deletion。

-iint：

不允许在reads两端出现大于多少bp的indel。

-lint：

Read前多少个碱基作为seed，如果设置的seed大于read长度，将无法继续，最

好设置在25-35，与-k2 配合使用。

-kint：

在seed中的最大编辑距离，使用默认2，与-l配合使用。

-tint：

要使用的线程数。

-Rint：

此参数只应用于pairend中，当没有出现大于此值的最佳比对结果时，将会降低标

准再次进行比对。

增加这个值可以提高配对比对的准确率，但是同时会消耗更长的

时间，默认是32。

-Iint：

表示输入的文件格式为Illumina1.3+数据格式。

-Bint：

设置标记序列。

从5’端开始多少个碱基作为标记序列，当-B为正值时，在比对之

前会将每个read的标记序列剪切，并将此标记序列表示在BCSAM 标签里，对于

pairend数据，两端的标记序列会被连接。

-b ：

指定输入格式为bam格式。

bwa  aln  hg19.fa  read.bam  >read.fq.gz.sai

（3）生成sam格式的比对文件。

如果一条read比对到多个位置，会随机选择一种。

     例子：

singleend：

bwa  samse  hg19.fa  read.fq.gz.sai  read.fq.gz  >read.fq.gz.sam

           参数：

-nint：

如果reads比对次数超过多少次，就不在XA标签显示。

                 -rstr：

定义头文件。

‘@RG\tID:

foo\tSM:

bar’，如果在此步骤不进行头文件定义，在

后续GATK分析中还是需要重新增加头文件。

           pairend：

bwasampe-a500read1.fq.gz.sairead2.fq.gz.sairead1.fq.gzread2.fq.gz>read.sam

           参数：

-aint：

最大插入片段大小。

                 -oint：

pairend两reads中其中之一所允许配对的最大次数，超过该次数，将被视为

singleend。

降低这个参数，可以加快运算速度，对于少于30bp的read，建

议降低-o值。

                 -rstr：

定义头文件。

同singleend。

                 -nint：

每对reads输出到结果中的最多比对数。

对于最后得到的sam文件，将比对上的结果提取出来（awk即可处理），即可直接用于GATK的分析。

注意：

由于GATK在下游的snpcalling时，是按染色体进行callsnp的。

因此，在准备原始sam文件时，可以先按染色体将文件分开，这样会提高运行速度。

但是当数据量不足时，可能会影响后续的VQSR分析，这是需要注意的。

2. 对sam文件进行进行重新排序（reorder）

由BWA生成的sam文件时按字典式排序法进行的排序（lexicographically）进行排序的（chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chrX，chrY），但是GATK在进行callsnp的时候是按照染色体组型（karyotypic）进行的（chrM，chr1，chr2…chr22，chrX，chrY），因此要对原始sam文件进行reorder。

可以使用picard-tools中的ReorderSam完成。

e.g．

java-jarpicard-tools-1.96/ReorderSam.jar

I=hg19.sam

O=hg19.reorder_00.sam

REFERENCE=hg19.fa

注意：

1. 这一步的头文件可以人工加上，同时要确保头文件中有的序号在下面序列中也有对应的。

虽然在GATK网站上的说明chrM可以在最前也可以在最后，但是当把chrM放在最后时可能会出错。

2. 在进行排序之前，要先构建参考序列的索引。

   e.g.samtoolsfaidxhg19.fa。

最后生成的索引文件：

hg19.fa.fai。

3.如果在上一步想把大文件切分成小文件的时候，头文件可以自己手工加上，之后运行这一步就好了。

3. 将sam文件转换成bam文件（bam是二进制文件，运算速度快）

这一步可使用samtoolsview完成。

e.g.samtoolsview-bShg19.reorder_00.sam-ohg19.sam_01.bam

4. 对bam文件进行sort排序处理

这一步是将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序。

可以使用picard-tools中SortSam完成。

e.g.

java-jarpicard-tools-1.96/SortSam.jar

INPUT=hg19.sam_01.bam

OUTPUT=hg19.sam.sort_02.bam

SORT_ORDER=coordinate

5. 对bam文件进行加头（head）处理

GATK2.0以上版本将不再支持无头文件的变异检测。

加头这一步可以在BWA比对的时候进行，通过-r参数的选择可以完成。

如果在BWA比对期间没有选择-r参数，可以增加这一步骤。

可使用picard-tools中AddOrReplaceReadGroups完成。

e.g.

java-jarpicard-tools-1.96/AddOrReplaceReadGroups.jar

I=hg19.sam.sort_02.bam

O=hg19.reorder.sort.addhead_03.bam

ID=hg19ID

LB=hg19ID

PL=illumine

PU=hg19PU

SM=hg19

IDstr：

输入reads集ID号；LB：

read集文库名；PL：

测序平台（illunima或solid）；PU：

测序平台下级单位名称（run的名称）；SM：

样本名称。

注意：

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

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. DuplicatesMarking

在制备文库的过程中，由于PCR扩增过程中会存在一些偏差，也就是说有的序列会被过量扩增。

这样，在比对的时候，这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。

而这些过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由PCR扩增所形成的duplicates，这一步可以使用picard-tools来完成。

去重复的过程是给这些序列设置一个flag以标志它们，方便GATK的识别。

还可以设置 REMOVE_DUPLICATES=true 来丢弃duplicated序列。

对于是否选择标记或者删除，对结果应该没有什么影响，GATK官方流程里面给出的例子是仅做标记不删除。

这里定义的重复序列是这样的：

如果两条reads具有相同的长度而且比对到了基因组的同一位置，那么就认为这样的reads是由PCR扩增而来，就会被GATK标记。

e.g.

java-jarpicard-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.bamMETRICS_FILE=hg19.reorder.sort.addhead.dedup_04.metrics

注意：

 dedup这一步只要在library层面上进行就可以了，例如一个sample如果建了多个库的话，对每个库进行dedup即可，不需要把所有库合成一个sample再进行dedup操作。

其实并不能准确的定义被mask的reads到底是不是duplicates，重复序列的程度与测序深度和文库类型都有关系。

最主要目的就是尽量减小文库构建时引入文库的PCRbias。

8. 要对上一步得到的结果生成索引文件

可以用samtools完成，生成的索引后缀是bai。

e.g.

samtoolsindexhg19.reorder.sort.addhead.dedup_04.bam

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

 （2014-03-0311:

16:

18）

转载▼

标签：

gatk

indel

snp

bwa

分类：

 生物信息

9.Localrealignmentaroundindels

这一步的目的就是将比对到indel附近的reads进行局部重新比对，将比对的错误率降到最低。

一般来说，绝大部分需要进行重新比对的基因组区域，都是因为插入/缺失的存在，因为在indel附近的比对会出现大量的碱基错配，这些碱基的错配很容易被误认为SNP。

还有，在比对过程中，比对算法对于每一条read的处理都是独立的，不可能同时把多条reads与参考基因组比对来排错。

因此，即使有一些reads能够正确的比对到indel，但那些恰恰比对到indel开始或者结束位置的read也会有很高的比对错误率，这都是需要重新比对的。

Localrealignment就是将由indel导致错配的区域进行重新比对，将indel附近的比对错误率降到最低。

主要分为两步：

第一步，通过运行RealignerTargetCreator来确定要进行重新比对的区域。

e.g.

java-jarGenomeAnalysisTK.jar

-Rhg19.fa

-TRealignerTargetCreator

-Ihg19.reorder.sort.addhead.dedup_04.bam

-ohg19.dedup.realn_06.intervals

-knownMills_and_1000G_gold_standard.indels.hg19.vcf

-known1000G_phase1.indels.hg19.vcf

参数说明：

-R：

               参考基因组；

-T：

               选择的GATK工具；

-I：

                输入上一步所得bam文件；

-o：

                输出的需要重新比对的基因组区域结果；

-maxInterval：

 允许进行重新比对的基因组区域的最大值，不能太大，太大耗费会很长时间，默认值500；

-known：

        已知的可靠的indel位点，指定已知的可靠的indel位点，重比对将主要围绕这些位点进行，对于人类基因组数据而言，可以直接指定GATKresourcebundle里面的indel文件（必须是vcf文件）。

对于knownsites的选择很重要，GATK中每一个用到knownsites的工具对于knownsites的使用都是不一样的，但是所有的都有一个共同目的，那就是分辨真实的变异位点和不可信的变异位点。

如果不提供这些knownsites的话，这些统计工具就会产生偏差，最后会严重影响结果的可信度。

在这些需要知道knownsites的工具里面，只有UnifiedGenotyper和HaplotypeCaller对knownsites没有太严格的要求。

如果你所研究的对象是人类基因组的话，那就简单多了，因为GATK网站上对如何使用人类基因组的knownsites做出了详细的说明，具体的选择方法如下表，这些文件都可以在GATK resourcebundle中下载。

Tool

dbSNP129 

dbSNP>132 

 Millsindels 

 1KGindels 

 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，这个网站上是做非人类基因组时，大家分享的经验，可以参考一下。

这个knownsites如果实在没有的话，也是可以自己构建的：

首先，先使用没有经过矫正的数据进行一轮SNPcalling；然后，挑选最可信的SNP位点进行BQSR分析；最后，在使用这些经过BQSR的数据进行一次真正的SNPcalling。

这几步可能要重复好多次才能得到可靠的结果。

第二步，通过运行IndelRealigner在这些区域内进行重新比对。

e.g.

java-jarGenomeAnalysisTK.jar

-Rhg19.fa

-TIndelRealigner

-targetIntervalshg19.dedup.realn_06.intervals

-Ihg19.reorder.sort.addhead.dedup_04.bam

-ohg19.dedup.realn_07.bam

-knownMills_and_1000G_gold_standard.indels.hg19.vcf

-known1000G_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.Basequalityscorerecalibration

这一步是对bam文件里reads的碱基质量值进行重新校正，使最后输出的bam文件中reads中碱基的质量值能够更加接近真实的与参考基因组之间错配的概率。

这一步适用于多种数据类型，包括illunima、solid、454、CG等数据格式。

在GATK2.0以上版本中还可以对indel的质量值进行校正，这一步对indelcalling非常有帮助

举例说明，在reads碱基质量值被校正之前，我们要保留质量值在Q25以上的碱基，但是实际上质量值在Q25的这些碱基的错误率在1%，也就是说质量值只有Q20，这样就会对后续的变异检测的可信度造成影响。

还有，在边合成边测序的测序过程中，在reads末端碱基的错误率往往要比起始部位更高。

另外，AC的质量值往往要低于TG。

BQSR的就是要对这些质量值进行校正。

BQSR主要有三步：

第一步：

利用工具BaseRecalibrator，根据一些knownsites，生成一个校正质量值所需要的数据文件，GATK网站以“.grp”为后缀命名。

e.g.

java-jarGenomeAnalysisTK.jar

-TBaseRecalibrator

-Rhg19.fa

-IChrALL.100.sam.dedup.realn.07.bam

-knownSitesdbsnp_137.hg19.vcf

-knownSitesMills_and_1000G_gold_standard.indels.hg19.vcf

-knownSites1000G_phase1.indels.hg19.vcf

-oChrALL.100.sam.recal.08-1.grp

第二步：

利用第一步生成的ChrALL.100.sam.recal.08-1.grp来生成校正后的数据文件，也是以“.grp”命名，这一步主要是为了与校正之前的数据进行比较，最后生成碱基质量值校正前后的比较图，如果不想生成最后BQSR比较图，这一步可以省略。

e.g.

java-jarGenomeAnalysisTK.jar

-TBaseRecalibrator

-Rhg19.fa

-IChrALL.100.sam.dedup.realn.07.bam

-BQSRChrALL.100.sam.recal.08-1.grp

-oGATK/hg19.recal.08-2.grp

-knownSitesdbsnp_137.hg19.vcf

-knownSitesMills_and_1000G_gold_standard.indels.hg19.vcf

-knownSites1000G_phase1.indels.hg19.vcf

第三步：

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

e.g.

java-jarGenomeAnalysisTK.jar

-TPrintReads

-Rhg19.fa

-IChrALL.100.sam.dedup.realn.07.bam

-BQSRChrALL.100.sam.recal.08-1.grp

-oChrALL.100.sam.recal.08-3.grp.bam

主要参数说明：

-bqsrBAQGOP：

BQSRBAQgapopen 罚值，默认值是40，如果是对全基因组数据进行BQSR分析，设置为30会更好。

-lqt：

                在计算过程中，该算法所能考虑的reads两端的最小质量值。

如果质量值小于该值，计算过程中将不予考虑，默认值是2。

    注意：

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

（2）除非你所研究的样本所得到的reads数实在太少，或者比对结果中的mismatch基本上都是实际存在的变异，否则必须要进行BQSR这一步。

对于人类基因组，即使有了dbSNP和千人基因组的数据，还有很多mismatch是错误的。

因此，这一步能做一定要做。

11. 分析和评估BQSR结果

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

e.g.

java-jarGenomeAnalysisTK.jar

-TAnaly