GATK使用方法详解-plob最详尽介绍说明手册.doc
《GATK使用方法详解-plob最详尽介绍说明手册.doc》由会员分享,可在线阅读,更多相关《GATK使用方法详解-plob最详尽介绍说明手册.doc(18页珍藏版)》请在得力文库 - 分享文档赚钱的网站上搜索。
1、GATKGATK 使用方法详解使用方法详解一、使用一、使用 GATK 前须知事项:前须知事项:(1)对 GATK 的测试主要使用的是人类全基因组和外显子组的测序数据,而 且全部是基于 illumina 数据格式,目前还没有提供其他格式文件(如 Ion Torrent)或者实验设计(RNA-Seq)的分析方法。(2)GATK 是一个应用于前沿科学研究的软件,不断在更新和修正,因此,在 使用 GATK 进行变异检测时,最好是下载最新的版本,目前的版本是 2.8.1(2014-02-25)。下载网站:http:/www.broadinstitute.org/gatk/download。(3)在 GA
2、TK 使用过程中(见下面图),有些步骤需要用到已知变异信息, 对于这些已知变 异,GATK 只提供了人类的已知变异信息,可以在 GATK 的 FTP 站点下载(GATK resource bundle)。如果要研究的不是人类基因组,需要 自行构建已知变异,GATK 提供了详细的构建方法。(4)GATK 在进行 BQSR 和 VQSR 的过程中会使用到 R 软件绘制一些图,因 此,在运行 GATK 之前最好先检查一下是否正确安装了 R 和所需要的包,所 需要的包大概包括 ggplot2、gplots、bitops、caTools、 colorspace、gdata、gsalib、reshape、
3、RColorBrewer 等。如果画图时出现错误, 会提示需要安装的包的名称。二、二、GATK 的使用流程的使用流程GATKGATK 最佳使用方案:最佳使用方案:共 3 大步骤,即:原始数据的处理 变异检测 初步分析。原始数据的处理原始数据的处理1.1. 对原始下机对原始下机 fastqfastq 文件进行过滤和比对(文件进行过滤和比对(mappingmapping)对于 Illumina 下机数据推荐使用 bwa 进行 mapping。Bwa 比对步骤大致如下:(1)对参考基因组构建索引:例子:bwa index -a bwtsw hg19.fa。构建索引时需要注意的问题:bwa 构建索引有
4、两种算法,两种算法都是基于 BWT 的,这两种算法通过参数-a is 和-a bwtsw 进行选择。其中-a bwtsw 对于短 的参考序列是不工作的,必须要大于等于 10Mb;-a is 是默认参数,这个参数不 适用于大的参考序列,必须要小于等于 2G。(2)寻找输入 reads 文件的 SA 坐标。对于 pair end 数据,每个 reads 文件单独做运算,single end 数据就不用说了, 只有一个文件。pair end:bwa aln hg19.fa read1.fq.gz -t 4 -I read1.fq.gz.saibwa aln hg19.fa read2.fq.gz -
5、t 4 -I read2.fq.gz.saisingle end:bwa aln hg19.fa read.fq.gz -l 30 -k 2 -t 4 -I read.fq.gz.sai主要参数说明:-o int:允许出现的最大 gap 数。-e int:每个 gap 允许的最大长度。-d int:不允许在 3端出现大于多少 bp 的 deletion。-i int:不允许在 reads 两端出现大于多少 bp 的 indel。-l int:Read 前多少个碱基作为 seed,如果设置的 seed 大于 read 长度,将无法继续,最好设置在 25-35,与-k 2 配合使用。-k int:
6、在 seed 中的最大编辑距离,使用默认 2,与-l 配合使用。-t int:要使用的线程数。-R int:此参数只应用于 pair end 中,当没有出现大于此值的最 佳比对结果时,将会降低标准再次进行比对。增加这个值可以提 高配对比对的准确率,但是同时会消耗更长的时间,默认是 32。-I-I intint:表示输入的文件格式为:表示输入的文件格式为 IlluminaIllumina 1.3+1.3+数据格式。数据格式。-B int:设置标记序列。从 5端开始多少个碱基作为标记序列, 当-B 为正值时,在比对之前会将每个 read 的标记序列剪切,并 将此标记序列表示在 BC SAM 标签里
7、,对于 pair end 数据,两端 的标记序列会被连接。-b :指定输入格式为 bam 格式。这是一个很奇怪的功能,就是对 其它软件的 bam 文件进行重新比对的意思bwa aln hg19.fa read.bam read.fq.gz.sai(3)生成 sam 格式的比对文件。如果一条 read 比对到多个位置,会随机选择 一种。例子:single end:bwa samse hg19.fa read.fq.gz.sai read.fq.gz read.fq.gz.sam参数:-n int:如果 reads 比对次数超过多少次,就不在 XA 标签显示。-r str:定义头文件。RGtID:
8、footSM:bar,如果在此步骤 不进行头文件定义,在后续 GATK 分析中还是需要重新增加头文件。pair end:bwa sampe -a 500 read1.fq.gz.sai read2.fq.gz.sai read1.fq.gz read2.fq.gz read.sam参数:-a int:最大插入片段大小。-o int:pair end 两 reads 中其中之一所允许配对的最大次数, 超过该次数,将被视为single end。降低这个参数,可以加快运算速度,对于少于 30bp 的 read,建议降低-o 值。-r str:定义头文件。同 single end。-n int:每对
9、reads 输出到结果中的最多比对数。对于最后得到的 sam 文件,将比对上的结果提取出来(awk 即可处理),即可 直接用于 GATK 的分析。注意:由于 GATK 在下游的 snp-calling 时,是按染色体进行 call-snp 的。因此, 在准备原始 sam 文件时,可以先按染色体将文件分开,这样会提高运行速度。 但是当数据量不足时,可能会影响后续的 VQSR 分析,这是需要注意的。2.2. 对对 samsam 文件进行进行重新排序(文件进行进行重新排序(reorderreorder)由 BWA 生成的 sam 文件时按字典式排序法进行的排序(lexicographically)进
10、 行排序的 (chr10,chr11chr19,chr1,chr20chr22,chr2,chr3chrM,chrX,chrY ),但是 GATK 在进行 callsnp 的时候是按照染色体组型(karyotypic)进行的 (chrM,chr1,chr2chr22,chrX,chrY),因此要对原始 sam 文件进行 reorder。可以使用 picard-tools 中的 ReorderSam 完成。eg.java -jar picard-tools-1.96/ReorderSam.jarI=hg19.samO=hg19.reorder_00.samREFERENCE=hg19.fa注意:
11、1) 这一步的头文件可以人工加上,同时要确保头文件中有的序号在下面序列中 也有对应的。虽然在 GATK 网站上的说明 chrM 可以在最前也可以在最后,但 是当把 chrM 放在最后时可能会出错。2) 在进行排序之前,要先构建参考序列的索引。e.g. samtools faidx hg19.fa。最后生成的索引文件:hg19.fa.fai。3) 如果在上一步想把大文件切分成小文件的时候,头文件可以自己手工加上, 之后运行这一步就好了。3.3. 将将 samsam 文件转换成文件转换成 bambam 文件(文件(bambam 是二进制文件,运算速度快)是二进制文件,运算速度快)这一步可使用 sa
12、mtools view 完成。e.g. samtools view -bS hg19.reorder_00.sam -o hg19.sam_01.bam 4.4. 对对 bambam 文件进行文件进行 sortsort 排序处理排序处理这一步是将 sam 文件中同一染色体对应的条目按照坐标顺序从小到大进行排序。 可以使用 picard-tools 中 SortSam 完成。e.g.java -jar picard-tools-1.96/SortSam.jarINPUT=hg19.sam_01.bamOUTPUT=hg19.sam.sort_02.bamSORT_ORDER=coordinate
13、5.5. 对对 bambam 文件进行加头(文件进行加头(headhead)处理)处理GATK2.0 以上版本将不再支持无头文件的变异检测。加头这一步可以在 BWA 比对的时候进行,通过-r 参数的选择可以完成。如果在 BWA 比对期间没有选 择-r 参数,可以增加这一步骤。可使用 picard-tools 中 AddOrReplaceReadGroups 完成。e.g.java -jar picard-tools-1.96/AddOrReplaceReadGroups.jarI=hg19.sam.sort_02.bamO=hg19.reorder.sort.addhead_03.bamID=
14、hg19IDLB=hg19IDPL=illuminePU=hg19PUSM=hg19ID str:输入 reads 集 ID 号;LB:read 集文库名;PL:测序平台(illunima 或 solid);PU:测序平台下级单位名称(run 的名称);SM:样本名称。注意:这一步尽量不要手动加头,本人尝试过多次手工加头,虽然看起来与软 件加的头是一样的,但是程序却无法运行。6.6. MergeMerge如果一个样本分为多个 lane 进行测序,那么在进行下一步之前可以将每个 lane 的 bam 文件合并。e.g.java -jar picard-tools-1.70/MergeSamFil
15、es.jarINPUT=lane1.bamINPUT=lane2.bamINPUT=lane3.bamINPUT=lane4.bamINPUT=lane8.bamOUTPUT=sample.bam7.7. DuplicatesDuplicates MarkingMarking在制备文库的过程中,由于 PCR 扩增过程中会存在一些偏差,也就是说有的序 列会被过量扩增。这样,在比对的时候,这些过量扩增出来的完全相同的序列 就会比对到基因组的相同位置。而这些过量扩增的 reads 并不是基因组自身固 有序列,不能作为变异检测的证据,因此,要尽量去除这些由 PCR 扩增所形成 的 duplicates
16、,这一步可以使用 picard-tools 来完成。去重复的过程是给这些序 列设置一个 flag 以标志它们,方便 GATK 的识别。还可以设置 REMOVE_DUPLICATES=true 来丢弃 duplicated 序列。对于是否选择标记或者 删除,对结果应该没有什么影响,GATK 官方流程里面给出的例子是仅做标记 不删除。这里定义的重复序列是这样的:如果两条 reads 具有相同的长度而且 比对到了基因组的同一位置,那么就认为这样的 reads 是由 PCR 扩增而来,就 会被 GATK 标记。e.g.java -jar picard-tools-1.96/MarkDuplicates
17、.jarREMOVE_DUPLICATES= falseMAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000INPUT=hg19.reorder.sort.addhead_03.bamOUTPUT=hg19.reorder.sort.addhead.dedup_04.bam METRICS_FILE=hg19.reorder.sort.addhead.dedup_04.metrics注意: dedup 这一步只要在 library 层面上进行就可以了,例如一个 sample 如 果建了多个库的话,对每个库进行 dedup 即可,不需要把所有库合成一 个 sample
18、再进行 dedup 操作。其实并不能准确的定义被 mask 的 reads 到底是不是 duplicates,重复序列的程度与测序深度和文库类型 都有关系。最主要目的就是最主要目的就是 尽量减小文库构建时引入文库的尽量减小文库构建时引入文库的 PCR bias。8.8. 要对上一步得到的结果生成索引文件要对上一步得到的结果生成索引文件可以用 samtools 完成,生成的索引后缀是 bai。e.g.samtools index hg19.reorder.sort.addhead.dedup_04.bam9.Local9.Local realignmentrealignment aroundar
19、ound indelsindels这一步的目的就是将比对到 indel 附近的 reads 进行局部重新比对,将比对的错 误率降到最低。一般来说,绝大部分需要进行重新比对的基因组区 域,都是因 为插入/缺失的存在,因为在 indel 附近的比对会出现大量的碱基错配,这些碱 基的错配很容易被误认为 SNP。还有,在比对过程中,比对算法对 于每一条 read 的处理都是独立的,不可能同时把多条 reads 与参考基因组比对来排错。因 此,即使有一些 reads 能够正确的比对到 indel,但那些 恰恰比对到 indel 开始 或者结束位置的 read 也会有很高的比对错误率,这都是需要重新比对的
20、。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.in
21、dels.hg19.vcf-known 1000G_phase1.indels.hg19.vcf参数说明:-R: 参考基因组;-T: 选择的 GATK 工具;-I: 输入上一步所得 bam 文件;-o: 输出的需要重新比对的基因组区域结果;-maxInterval: 允许进行重新比对的基因组区域的最大值,不能 太大,太大耗费会很长时间,默认值 500;-known: 已知的可靠的 indel 位点,指定已知的可靠的 indel 位 点,重比对将主要围绕这些位点进行,对于人类基因组数据而言, 可以直接指定 GATK resource bundle 里面的 indel 文件(必须是 vcf 文件)
22、。对于 known sites 的选择很重要,GATK 中 每一个用到 known sites 的工具对于 known sites 的使用都是不一样的,但是所有的都有一个共同目的,那就是分辨 真实的变异位点和不可信的变异位点。如果不提供这些 known sites 的话,这些 统计工具就会产生偏差,最后会严重影响结果的可信度。在这些需要知道 known sites 的工具里面,只有 UnifiedGenotyper 和 HaplotypeCaller 对 known sites 没有太严格的要求。如果你所研究的对象是人类基因组的话,那就简单多了,因为 GATK 网站上对 如何使用人类基因组的
23、known sites 做出了详细的说明,具体的选择方法如下表, 这些文件都可以在 GATK resource bundle 中下载。ToolTooldbSNPdbSNP 129129 dbSNPdbSNP 132132MillsMills indelsindels1KG1KG indelsindels HapMapHapMapOmniOmniRealignerTargetCreatorRealignerTargetCreatorXXIndelRealignerIndelRealignerXX BaseRecalibratorBaseRecalibratorXXX (UnifiedGenoty
24、per/(UnifiedGenotyper/ HaplotypeCaller)HaplotypeCaller)XVariantRecalibratorVariantRecalibratorXXXX VariantEvalVariantEvalX但是如果你要研究的不是人类基因组的话,那就有点麻烦了, http:/www.broadinstitute.org/gatk/guide /article?id=1243,这个网站上是做非人类 基因组时,大家分享的经验,可以参考一下。这个 known sites 如果实在没有的 话,也是可以自己构建的:首先,先使用没有经过矫正的数据进行一轮 SNP cal
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- GATK 使用方法 详解 plob 详尽 介绍 说明 手册
限制150内