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

-.z.GATK使用方法详解一、使用GATK前须知事项：〔1〕对GATK的测试主要使用的是人类全基因组和外显子组的测序数据，而且全部是基于illumina数据格式，目前还没有提供其他格式文件〔如IonTorrent〕或者实验设计〔RNA-Seq〕的分析方法。〔2〕GATK是一个应用于前沿科学研究的软件，不断在更新和修正，因此，在使用GATK进展变异检测时，最好是下载最新的版本，目前的版本是2.8.1〔2014-02-25〕。下载：./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〕对参考基因组构建索引：例子：bwainde*-abwtswhg19.fa。构建索引时需要注意的问题：bwa构建索引有两种算法，两种算法都是基于BWT的，这两种算法通过参数-ais和-abwtsw进展选择。其中-abwtsw对于短的参考序列是不工作的，必须要大于等于10Mb；-ais是默认参数，这个参数不适用于大的参考序列，必须要小于等于2G。〔2〕寻找输入reads文件的SA坐标。对于pairend数据，每个reads文件单独做运算，singleend数据就不用说了，只有一个文件。pairend：singleend：主要参数说明：-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格式。这是一个很奇怪的功能，就是对其它软件的bam文件进展重新比对的意思〔3〕生成sam格式的比对文件。如果一条read比对到多个位置，会随机选择一种。例子：singleend：参数：-nint：如果reads比对次数超过多少次，就不在*A标签显示。-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在下游的snp-calling时，是按染色体进展call-snp的。因此，在准备原始sam文件时，可以先按染色体将文件分开，这样会提高运行速度。但是当数据量缺乏时，可能会影响后续的VQSR分析，这是需要注意的。2.对sam文件进展进展重新排序〔reorder〕由BWA生成的sam文件时按字典式排序法进展的排序〔le*icographically〕进展排序的〔chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chr*，chrY〕，但是GATK在进展callsnp的时候是按照染色体组型〔karyotypic〕进展的〔chrM，chr1，chr2…chr22，chr*，chrY〕，因此要对原始sam文件进展reorder。可以使用picard-tools中的ReorderSam完成。eg.java-jarpicard-tools-1.96/ReorderSam.jarI=hg19.samREFERENCE=hg19.fa注意：1)这一步的头文件可以人工加上，同时要确保头文件中有的序号在下面序列中也有对应的。虽然在GATK上的说明chrM可以在最前也可以在最后，但是当把chrM放在最后时可能会出错。2)在进展排序之前，要先构建参考序列的索引。e.g.samtoolsfaid*hg19.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.jarSORT_ORDER=coordinate5.对bam文件进展加头〔head〕处理GATK2.0以上版本将不再支持无头文件的变异检测。加头这一步可以在BWA比对的时候进展，通过-r参数的选择可以完成。如果在BWA比对期间没有选择-r参数，可以增加这一步骤。可使用picard-tools中AddOrReplaceReadGroups完成。e.g.java-jarpicard-tools-1.96/AddOrReplaceReadGroups.jarID=hg19IDLB=hg19IDPL=illuminePU=hg19PUSM=hg19IDstr：输入reads集ID号；LB：read集文库名；PL：测序平台〔illunima或solid〕；PU：测序平台下级单位名称〔run的名称〕；SM：样本名称。注意：这一步尽量不要手动加头，本人尝试过屡次手工加头，虽然看起来与软件加的头是一样的，但是程序却无法运行。6.Merge如果一个样本分为多个lane进展测序，则在进展下一步之前可以将每个lane的bam文件合并。e.g.java-jarpicard-tools-1.70/MergeSamFiles.jarINPUT=lane1.bamINPUT=lane2.bamINPUT=lane3.bamINPUT=lane4.bam……INPUT=lane8.bamOUTPUT=sample.bam7.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.jarREMOVE_DUPLICATES=falseMA*_FILE_HANDLES_FOR_READ_ENDS_MAP=8000注意：dedup这一步只要在library层面上进展就可以了，例如一个sample如果建了多个库的话，对每个库进展dedup即可，不需要把所有库合成一个sample再进展dedup操作。其实并不能准确的定义被mask的reads到底是不是duplicates，重复序列的程度与测序深度和文库类型都有关系。最主要目的就是尽量减小文库构建时引入文库的PCRbias。8.要对上一步得到的结果生成索引文件可以用samtools完成，生成的索引后缀是bai。e.g.9.Localrealignmentaroundindels这一步的目的就是将比对到indel附近的reads进展局部重新比对，将比对的错误率降到最低。一般来说，绝大局部需要进展重新比对的基因组区域，都是因为插入/缺失的存在，因为在indel附近的比对会出现大量的碱基错配，这些碱基的错配很容易被误认为SNP。还有，在比对过程中，比对算法对于每一条read的处理都是独立的，不可能同时把多条reads与参考基因组比对来排错。因此，即使有一些reads能够正确的比对到indel，但那些恰恰比对到indel开场或者完毕位置的read也会有很高的比对错误率，这都是需要重新比对的。Localrealignment就是将由indel导致错配的区域进展重新比对，将indel附近的比对错误率降到最低。主要分为两步：第一步，通过运行RealignerTargetCreator来确定要进展重新比对的区域。e.g.java-jarGenomeAnalysisTK.jar-Rhg19.fa-TRealignerTargetCreator-known-known参数说明：-R：参考基因组；-T：选择的GATK工具；-I：输入上一步所得bam文件；-o：输出的需要重新比对的基因组区域结果；-ma*Interval：允许进展重新比对的基因组区域的最大值，不能太大，太大消耗会很长时间，默认值500；-known：的可靠的indel位点，指定的可靠的indel位点，重比对将主要围绕这些位点进展，对于人类基因组数据而言，可以直接指定GATKresourcebundle里面的indel文件〔必须是vcf文件〕。对于knownsites的选择很重要，GATK中每一个用到knownsites的工具对于knownsites的使用都是不一样的，但是所有的都有一个共同目的，那就是分辨真实的变异位点和不可信的变异位点。如果不提供这些knownsites的话，这些统计工具就会产生偏差，最后会严重影响结果的可信度。在这些需要知道knownsites的工具里面，只有UnifiedGenotyper和HaplotypeCaller对knownsites没有太严格的要求。如果你所研究的对象是人类基因组的话，那就简单多了，因为GATK上对如何使用人类基因组的knownsites做出了详细的说明，具体的选择方法如下表，这些文件都可以在GATKresourcebundle中下载。TooldbSNP129dbSNP>132Millsindels1KGindelsHapMapOmniRealignerTargetCreator**IndelRealigner**BaseRecalibrator***(UnifiedGenotyper/HaplotypeCaller)*VariantRecalibrator****VariantEval*但是如果你要研究的不是人类基因组的话，那就有点麻烦了，./gatk/guide/article"id=1243，这个上是做非人类基因组时，大家分享的经历，可以参考一下。这个knownsites如果实在没有的话，也是可以自己构建的：首先，先使用没有经过矫正的数据进展一轮SNPcalling；然后，挑选最可信的SNP位点进展BQSR分析；最后，在使用这些经过BQSR的数据进展一次真正的SNPcalling。这几步可能要重复好屡次才能得到可靠的结果。第二步，通过运行IndelRealigner在这些区域进展重新比对。e.g.java-jarGenomeAnalysisTK.jar-Rhg19.fa-TIndelRealigner运行完毕后，生成的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第二步：利用第一步生成的ChrALL.100.sam.recal.08-1.grp来生成校正后的数据文件，也是以“.grp〞命名，这一步主要是为了与校正之前的数据进展比拟，最后生成碱基质量值校正前后的比拟图，如果不想生成最后BQSR比拟图，这一步可以省略。e.g.java-jarGenomeAnalysisTK.jar-TBaseRecalibrator-Rhg19.fa-oGATK第三步：利用工具PrintReads将经过质量值校正的数据输出到新的bam文件中，用于后续的变异检测。e.g.java-jarGenomeAnalysisTK.jar-TPrintReads-Rhg19.fa主要参数说明：-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-TAnalyzeCovariates-Rhg19.fa参数解释：-before：基于原始比对结果生成的第一次校对表格。-after：基于第一次校对表格生成的第二次校对表格。-plots：评估BQSR结果的报告文件。-csv：生成报告中图标所需要的所有数据。

12.Reducebamfile这一步是使用ReduceReads这个工具将bam文件进展压缩，生成新的bam文件，新的bam文件仍然保持bam文件的格式和所有进展变异检测所需要的信息。这样不仅能够节省存储空间，也方便后续变异检测过程中对数据的处理。e.g.java-jarGenomeAnalysisTK.jar-TReduceReads-Rhg19.fa到此为止，GATK流程中的第一大步骤就完毕了，完成了variantscalling所需要的所有准备工作，生成了用于下一步变异检测的bam文件。变异检测1.VariantCallingGATK在这一步里面提供了两个工具进展变异检测——UnifiedGenotyper和HaplotypeCaller。其中HaplotypeCaller一直还在开发之中，包括生成的结果以及计算模型和命令行参数一直在变动，因此，目前使用比拟多的还是UnifiedGenotyper。此外，HaplotypeCaller不支持Reduce之后的bam文件，因此，中选择使用HaplotypeCaller进展变异检测时，不需要进展Reducereads。UnifiedGenotyper是集合多种变异检测方法而成的一种VariantsCaller，既可以用于单个样本的变异检测，也可以用于群体的变异检测。UnifiedGenotyper使用贝叶斯最大似然模型，同时估计基因型和基因频率，最后对每一个样本的每一个变异位点和基因型都会给出一个准确的后验概率。e.g.java-jarGenomeAnalysisTK.jar-glmBOTH-lINFO-Rhg19.fa-TUnifiedGenotyper-stand_call_conf10-stand_emit_conf30上述命令将对输入的bam文件中的所有样本进展变异检测，最后生成一个vcf文件，vcf文件中会包含所有样本的变异位点和基因型信息。但是现在所得到的结果是最原始的、没有经过任何过滤和校正的Variants集合。这一步产生的变异位点会有很高的假阳性，尤其是indel，因此，必须要进展进一步的筛选过滤。这一步还可以指定对基因组的*一区域进展变异检测，只需要增加一个参数-L：target_interval.list，格式是bed格式文件。主要参数解释：-A：指定一个或者多个注释信息，最后输出到vcf文件中。-*A：指定不做哪些注释，最后不会输出到vcf文件中。-D：的snp文件。-glm：选择检测变异的类型。SNP表示只进展snp检测；INDEL表示只对indel进展检测；BOTH表示同时检测snp和indel。默认值是SNP。-hets：杂合度的值，用于计算先验概率。默认值是0.001。-ma*AltAlleles：容许存在的最大altallele的数目，默认值是6。这个参数要特别注意，不要轻易修改默认值，程序设置的默认值几乎可以满足所有的分析，如果修改了可能会导致程序无法运行。-mbq：变异检测时，碱基的最小质量值。如果小于这个值，将不会对其进展变异检测。这个参数不适用于indel检测，默认值是17。-minIndelt：在做indelcalling的时候，支持一个indel的最少read数量。也就是说，如果同时有多少条reads同时支持一个候选indel时，软件才开场进展indelcalling。降低这个值可以增加indelcalling的敏感度，但是会增加消耗的时间和假阳性。-minIndelFrac：在做indelcalling的时候，支持一个indel的reads数量占比对到该indel位置的所有reads数量的百分比。也就是说，只有同时满足-minIndelt和-minIndelFrac两个参数条件时，才会进展indelcalling。-onlyEmitSamples：当指定这个参数时，只有指定的样本的变异检测结果会输出到vcf文件中。-stand_emit_conf：在变异检测过程中，所容许的最小质量值。只有大于等于这个设定值的变异位点会被输出到结果中。-stand_call_conf：在变异检测过程中，用于区分低质量变异位点和高质量变异位点的阈值。只有质量值高于这个阈值的位点才会被视为高质量的。低于这个质量值的变异位点会在输出结果中标注LowQual。在千人基因组方案第二阶段的变异检测时，利用35*的数据进展snpcalling的时候，当设置成50时，有大概10%的假阳性。-dcov：这个参数用于控制检测变异数据的coverage(*)，4*的数据可以设置为40，大于30*的数据可以设置为200。注意：GATK进展变异检测的时候，是按照染色体排序顺序进展的〔先callchr1，然后chr2，然后chr3…最后chrY〕，并非多条染色体并行检测的，因此，如果数据量比拟大的话，建议分染色体分别进展，对性染色体的变异检测可以同常染色体方法。大多数参数的默认值可以满足大多数研究的需求，因此，在做变异检测过程中，如果对参数意义不是很明确，不建议修改。2.对原始变异检测结果进展过滤〔hardfilterandVQSR〕这一步的目的就是对上一步call出来的变异位点进展过滤，去掉不可信的位点。这一步可以有两种方法，一种是通过GATK的VariantFiltration，另一种是通过GATK的VQSR〔变异位点质量值重新校正〕进展过滤。通过GATK上提供的最正确方案可以看出，GATK是推荐使用VASR的，但使用VQSR数据量一定要到达要求，数据量太小无法使用高斯模型。还有，在使用VAQR时，indel和snp要分别进展。VQSR原理介绍：这个模型是根据已有的真实变异位点〔人类基因组一般使用HapMap3中的位点，以及这些位点在Omni2.5MSNP芯片中出现的多态位点〕来训练，最后得到一个训练好的能够很好的评估真伪的错误评估模型，可以叫他适应性错误评估模型。这个适应性的错误评估模型可以应用到call出来的原始变异集合中的变异位点和新发现的变异位点，进而去评估每一个变异位点发生错误的概率，最终会给出一个得分。这个得分最后会被写入vcf文件的INFO信息里，叫做VQSLOD，就是在训练好的混合高斯模型下，一个位点是真实的概率比上这个位点可能是假阳性的概率的logoddsratio〔对数差异比〕，因此，可以定性的认为，这个值越大就越好。VQSR主要分两个步骤，这两个步骤会使用两个不同的工具：VariantRecalibrator和ApplyRecalibration。i)VariantRecalibratorVariantRecalibrator：通过大量的高质量的变异集合的各个注释〔包括很多种，后面介绍〕的值来创立一个高斯混合模型，然后用于评估所有的变异位点。这个文件最后将生成一个recalibration文件。原理简单