GATK使用方法.docx
- 文档编号:30284943
- 上传时间:2023-08-13
- 格式:DOCX
- 页数:32
- 大小:217.39KB
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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- GATK 使用方法
![提示](https://static.bdocx.com/images/bang_tan.gif)