首页 GATK使用方法详解

GATK使用方法详解

举报
开通vip

GATK使用方法详解GATK使用方法详解一、使用GATK前须知事项:对GATK的测试主要使用的是人类全基因组和外显子组的测序数据,而且全部是基于illumina数据格式,目前还没有提供其他格式文件(如IonTorrent)或者实验设计(RNA-Seq)的分析方法。GATK是一个应用于前沿科学研究的软件,不断在更新和修正,因此,在使用GATK进行变异检测时,最好是下载最新的版本,目前的版本是2.8.1(2014-02-25)。下载网站:。在GATK使用过程中(见下面图),有些步骤需要用到已知变异信息,对于这些已知变异,GATK只提供了人...

GATK使用方法详解
GATK使用方法详解一、使用GATK前须知事项:对GATK的测试主要使用的是人类全基因组和外显子组的测序数据,而且全部是基于illumina数据格式,目前还没有提供其他格式文件(如IonTorrent)或者实验设计(RNA-Seq)的分析方法。GATK是一个应用于前沿科学研究的软件,不断在更新和修正,因此,在使用GATK进行变异检测时,最好是下载最新的版本,目前的版本是2.8.1(2014-02-25)。下载网站:。在GATK使用过程中(见下面图),有些步骤需要用到已知变异信息,对于这些已知变异,GATK只提供了人类的已知变异信息,可以在GATK的FTP站点下载(GATKresourcebundle)。如果要研究的不是人类基因组,需要自行构建已知变异,GATK提供了详细的构建方法。GATK在进行BQSR和VQSR的过程中会使用到R软件绘制一些图,因此,在运行GATK之前最好先检查一下是否正确安装了R和所需要的包,所需要的包大概包括ggplot2、gplots、bitops、caTools、colorspace、gdata、gsalib、reshape、RColorBrewer等。如果画图时出现错误,会提示需要安装的包的名称。二、GATK的使用流程»Variant»Pr^lliTiJnaryAnalysesDataPfe-processlngi^-H31UDZIndAlRrtPkligr-iiin^nlSamRiRRCofnpnrasionAmql«ifl-RpodtfnIVar|Hotftg-GalitiraJlkjnFunchorulAnruMstlmAnai^tia-AaAdy[9NPifrwWejIwU-incinp 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 :共3大步骤,即:原始数据的处理-->变异检测-->初步分析原始数据的处理对原始下机fastq文件进行过滤和比对(mapping)对于Illumina下机数据推荐使用bwa进行mapping。Bwa比对步骤大致如下:(1)对参考基因组构建索引:例子:bwaindex-abwtswhgl9.fa。构建索引时需要注意的问题:bwa构建索引有两种算法,两种算法都是基于BWT的,这两种算法通过参数-ais和-abwtsw进行选择。其中-abwtsw对于短的参考序列是不工作的,必须要大于等于10Mb;-ais是默认参数,这个参数不适用于大的参考序列,必须要小于等于2G。(2)寻找输入reads文件的SA坐标。对于pairend数据,每个reads文件单独做运算,singleend数据就不用说了,只有一个文件。pairend:bwaalnhg19.faread1.fq.gz-t4-I>read1.fq.gz.saibwaalnhg19.faread2.fq.gz-t4-I>read2.fq.gz.saisingleend:bwaalnhg19.faread.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格式。这是一个很奇怪的功能,就是对其它软件的bam文件进行重新比对的意思bwaalnhg19.faread.bam>read.fq.gz.sai(3)生成sam格式的比对文件。如果一条read比对到多个位置,会随机选择一种。例子:singleend:bwasamsehg19.faread.fq.gz.sairead.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在下游的snp-calling时,是按染色体进行call-snp的。因此,在准备原始sam文件时,可以先按染色体将文件分开,这样会提高运行速度但是当数据量不足时,可能会影响后续的VQSR分析,这是需要注意的。对sam文件进行进行重新排序(reorder)由BWA生成的sam文件时按字典式排序法进行的排序(lexicographically)进行排序的(chr10,chrll...chrl9,chrl,chr20...chr22,chr2,chr3...chrM,chrX,chrY),但是GATK在进行callsnp的时候是按照染色体组型(karyotypic)进行的(chrM,chrl,chr2...chr22,chrX,chrY),因此要对原始sam文件进行reorder。可以使用picard-tools中的ReorderSam完成。eg.java-jarpicard-toolsT.96/ReorderSam.jarI=hg19.samO=hg19.reorder_00.samREFERENCE=hg19.fa1)这一步的头文件可以人工加上,同时要确保头文件中有的序号在下面序列中也有对应的。虽然在GATK网站上的说明chrM可以在最前也可以在最后,但是当把chrM放在最后时可能会出错。2)在进行排序之前,要先构建参考序列的索引。e.g.samtoolsfaidxhg19.fa。最后生成的索引文件:hg19.fa.fai。3)如果在上一步想把大文件切分成小文件的时候,头文件可以自己手工加上,之后运行这一步就好了。将sam文件转换成bam文件(bam是二进制文件,运算速度快)这一步可使用samtoolsview完成。e.g.samtoolsview-bShg19.reorder_00.sam-ohg19.sam_01.bam对bam文件进行sort排序处理这一步是将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序可以使用picard-tools中SortSam完成。e.g.java-jarpicard-tools-1.96/SortSam.jarINPUT=hg19.sam_01.bamOUTPUT=hg19.sam.sort_02.bamSORT_ORDER=coordinate对bam文件进行加头(head)处理GATK2.0以上版本将不再支持无头文件的变异检测。加头这一步可以在BWA比对的时候进行,通过-r参数的选择可以完成。如果在BWA比对期间没有选择-r参数,可以增加这一步骤。可使用picard-tools中AddOrReplaceReadGroups完成。e.g.java-jarpicard-tools-1.96/AddOrReplaceReadGroups.jarI=hg19.sam.sort_02.bamO=hg19.reorder.sort.addhead_03.bamID=hg19IDLB=hg19IDPL=illuminePU=hg19PUSM=hg19IDstr:输入reads集ID号;LB:read集文库名;PL:测序平台(illunima或solid);PU:测序平台下级单位名称(run的名称);SM:样本名称。注意:这一步尽量不要手动加头,本人尝试过多次手工加头,虽然看起来与软件加的头是一样的,但是程序却无法运行。Merge如果一个样本分为多个lane进行测序,那么在进行下一步之前可以将每个lane的bam文件合并。e.g.java-jarpicard-tools-1.70/MergeSamINPUT=lane1.bamINPUT=lane2.bamINPUT=lane3.bamINPUT=lane4.bamINPUT=lane8.bamOUTPUT=sample.bamDuplicatesMarking在制备文库的过程中,由于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=falseMAX_INPUT=hg19.reorder.sort.addhead_03.bamOUTPUT=hg19.reorder.sort.addhead.dedup_04.bamMETRICS_注意:dedup这一步只要在library层面上进行就可以了,例如一个sample如果建了多个库的话,对每个库进行dedup即可,不需要把所有库合成一个sample再进行dedup操作。其实并不能准确的定义被mask的reads到底是不是duplicates,重复序列的程度与测序深度和文库类型都有关系。最主要目的就是尽量减小文库构建时引入文库的PCRbias。要对上一步得到的结果生成索引文件可以用samtools完成,生成的索引后缀是bai。e.g.samtoolsindexhg19.reorder.sort.addhead.dedup_04.bamLocalrealignmentaroundindels这一步的目的就是将比对至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做出了详细的说明,具体的选择方法如下表,这些文件都可以在GATKresourcebundle中下载。TooldbSNP129dbSNP>132Millsindels1KGindelsHapMapOmniRealignerTargetCreatorXXIndelRealignerXXBaseRecalibratorXXX(UnifiedGenotyper/HaplotypeCaller)XVariantRecalibratorXXXXVariantEvalX但是如果你要研究的不是人类基因组的话,那就有点麻烦了,/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运行结束后,生成的hgl9.dedup.realn_07.bam即为最后重比对后的文件。注意:第一步和第二步中使用的输入文件(bam文件)、参考基因组和已知indel文件必须是相同的文件。当在相同的基因组区域发现多个indel存在时,这个工具会从其中选择一个最有可能存在比对错误的indel进行重新比对,剩余的其他indel不予考虑。对于454下机数据,本工具不支持。此外,这一步还会忽略bwa比对中质量值为0的read以及在CIGAR信息中存在连续indel的reads。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/hgl9.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。注意:当bam文件中的reads数量过少时,BQSR可能不会正常工作,GATK网站建议base数量要大于100M才能得到比较好的结果。除非你所研究的样本所得到的reads数实在太少,或者比对结果中的mismatch基本上都是实际存在的变异,否则必须要进行BQSR这一步。对于人类基因组,即使有了dbSNP和千人基因组的数据,还有很多mismatch是错误的。因此,这一步能做一定要做。分析和评估BQSR结果这一步会生成评估前后碱基质量值的比较结果,可以选择使用图片和表格的形式展示。e.g.java-jarGenomeAnalysisTK.jar-TAnalyzeCovariates-Rhg19.fa-beforeChrALL.100.sam.recal.08-1.grp-afterChrALL.100.sam.recal.08-2.grp-csvChrALL.100.sam.recal.grp.09.csv-plotsChrALL.100.sam.recal.grp.09.pdf参数解释:-before:基于原始比对结果生成的第一次校对表格。-after:基于第一次校对表格生成的第二次校对表格。-plots:评估BQSR结果的报告文件。-csv:生成报告中图标所需要的所有数据。Reducebamfile这一步是使用ReduceReads这个工具将bam文件进行压缩,生成新的bam文件,新的bam文件仍然保持bam文件的格式和所有进行变异检测所需要的信息。这样不仅能够节省存储空间,也方便后续变异检测过程中对数据的处理。e.g.java-jarGenomeAnalysisTK.jar-TReduceReads-Rhg19.fa-IChrALL.100.sam.recal.08-3.grp.bam-oChrALL.100.sam.recal.08-3.grp.reduce.bam到此为止,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-IChrALL.100.sam.recal.08-3.grp.reduce.bam-Ddbsnp_137.hg19.vcf-oChrALL.100.sam.recal.10.vcf-metricsChrALL.100.sam.recal.10.metrics-stand_call_conf10-stand_emit_conf30上述命令将对输入的bam文件中的所有样本进行变异检测,最后生成一个vcf文件,vcf文件中会包含所有样本的变异位点和基因型信息。但是现在所得到的结果是最原始的、没有经过任何过滤和校正的Variants集合。这一步产生的变异位点会有很高的假阳性,尤其是indel,因此,必须要进行进一步的筛选过滤。这一步还可以指定对基因组的某一区域进行变异检测,只需要增加一个参数-L:target_interval.list,格式是bed格式文件。主要参数解释:-A:指定一个或者多个注释信息,最后输出到vcf文件中。-XA:指定不做哪些注释,最后不会输出到vcf文件中。-D:已知的snp文件。-glm:选择检测变异的类型。SNP表示只进行snp检测;INDEL表示只对indel进行检测;BOTH表示同时检测snp和indel。默认值是SNP。-hets:杂合度的值,用于计算先验概率。默认值是0.001。-maxAltAlleles:容许存在的最大altallele的数目,默认值是6。这个参数要特别注意,不要轻易修改默认值,程序设置的默认值几乎可以满足所有的分析,如果修改了可能会导致程序无法运行。-mbq:变异检测时,碱基的最小质量值。如果小于这个值,将不会对其进行变异检测。这个参数不适用于indel检测,默认值是17。-minlndelCnt:在做indelcalling的时候,支持一个indel的最少read数量。也就是说,如果同时有多少条reads同时支持一个候选indel时,软件才开始进行indelcalling。降低这个值可以增加indelcalling的敏感度,但是会增加耗费的时间和假阳性。-minlndelFrac:在做indelcalling的时候,支持一个indel的reads数量占比对到该indel位置的所有reads数量的百分比。也就是说,只有同时满足-minlndelCnt和-minlndelFrac两个参数条件时,才会进行indelcalling。-onlyEmitSamples:当指定这个参数时,只有指定的样本的变异检测结果会输出到vcf文件中。-stand_emit_conf:在变异检测过程中,所容许的最小质量值。只有大于等于这个设定值的变异位点会被输出到结果中。-stand_call_conf:在变异检测过程中,用于区分低质量变异位点和高质量变异位点的阈值。只有质量值高于这个阈值的位点才会被视为高质量的。低于这个质量值的变异位点会在输出结果中标注LowQual。在千人基因组 计划 项目进度计划表范例计划下载计划下载计划下载课程教学计划下载 第二阶段的变异检测时,利用35x的数据进行snpcalling的时候,当设置成50时,有大概10%的假阳性。-dcov:这个参数用于控制检测变异数据的coverage(X),4X的数据可以设置为40,大于30X的数据可以设置为200。注意:GATK进行变异检测的时候,是按照染色体排序顺序进行的(先callchrl,然后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文件。原理简单介绍:这个模型首先要拿到真实变异数据集和上一步骤中得到的原始变异数据集的交集,然后对这些SNP值相对于具体注释信息的分布情况进行模拟,将这些变异位点进行聚类,最后根据聚类结果赋予所有变异位点相应的VQSLOD值。越接近聚类核心的变异位点得到的VQSLOD值越高。ii)ApplyRecalibrationApplyRecalibration:这一步将模型的各个参数应用于原始vcf文件中的每一个变异位点,这时,每一个变异位点的注释信息列中都会出现一个VQSLOD值,然后模型会根据这个值对变异位点进行过滤,过滤后的信息会写在vcf文件的filter—列中。原理简单介绍:在VariantRecalibrator这一步中,每个变异位点已经得到了一个VQSLOD值了,同时,这些LOD值在训练集里也进行了排序。当你在这一步中设置一个tranchesensitivity的阈值(这个阈值一般是一个百分数,如设置成99%),那么,如果LOD值从大到小排序的话,这个程序就会认为在这个训练集中,LOD值在前99%的是可信的,当这个值低于这个阈值,就认为是错误的。最后,程序就会用这个标准来过滤上一步call出来的原始变异集合。如果LOD值超过这个阈值,在filter那一列就会显示PASS,如果低于这个值就会被过滤掉,但是这些位点仍然会显示在结果里面,只不过会在filter那一列标示出他所属于的tranchesensitivity的名称。在设置tranchesensitivity的阈值时,要兼顾敏感度和质量值。初步分析这一步主要是对上面所得到的最终vcf中的结果进行一些初步的分析,比如计算这些变异位点在dbsnp中的比例、Ti/Tv的比例、每个样本中的snp数量。此外,还可以对变异位点的同义/非同义突变进行统计,识别是否为CpG位点以及氨基酸的简并信息等。这一步主要是利用GATK中的VariantEval来完成。需要注意的是,有些计算内容不能同时进行,例如AlleleCount和VariantSummary或者Sample和VariantSummary。如果选择了这样的组合方式,程序就会报错。但是GATK并没有告诉我们到底哪些不能同时运行,所以当选择计算内容的时候可以先做一下测试。e.g.java-jarGenomeAnalysisTK.jar-Rhg19.fa-TVariantEval--evalhg19.snp.filter.t97.Q10_13.both.vcf-Ddbsnp_137.hg19.vcf-ohg19.PASS.Eval_15_Final.gatkreport主要参数解释:--eval输入要进行summary的文件,也就是hg19.snp.filter.t97.Q10_13.both.vcf。-EV选择模块计算相应的分析内容,。--list列出可供选择的计算模块。-noEV不是用默认的模块,只计算用-EV选定的模块。
本文档为【GATK使用方法详解】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑, 图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
该文档来自用户分享,如有侵权行为请发邮件ishare@vip.sina.com联系网站客服,我们会及时删除。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。
下载需要: 免费 已有0 人下载
最新资料
资料动态
专题动态
is_072127
暂无简介~
格式:doc
大小:64KB
软件:Word
页数:18
分类:
上传时间:2019-11-18
浏览量:1