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

25、ling；然后，挑选最可信的 SNP 位点进行 BQSR 分析；最后，在使用这些经 过 BQSR 的数据进行一次真正的 SNP calling。这几步可能要重复好多次才能得 到可靠的结果。第二步，通过运行 IndelRealigner 在这些区域内进行重新比对。e.g.java -jar GenomeAnalysisTK.jar-R hg19.fa-T IndelRealigner-targetIntervals hg19.dedup.realn_06.intervals-I hg19.reorder.sort.addhead.dedup_04.bam-o hg19.dedup.realn_0

26、7.bam-known Mills_and_1000G_gold_standard.indels.hg19.vcf-known 1000G_phase1.indels.hg19.vcf运行结束后，生成的 hg19.dedup.realn_07.bam 即为最后重比对后的文件。注意：注意：1. 第一步和第二步中使用的输入文件（bam 文件）、参考基因组和已知 indel 文件必须是相同的文件。2. 当在相同的基因组区域发现多个 indel 存在时，这个工具会从其中选择一个 最有可能存在比对错误的 indel 进行重新比对，剩余的其他 indel 不予考虑。3. 对于 454 下机数据，本工具不支

27、持。此外，这一步还会忽略 bwa 比对中质量 值为 0 的 read 以及在 CIGAR 信息中存在连续 indel 的 reads。10.Base10.Base qualityquality scorescore recalibrationrecalibration这一步是对 bam 文件里 reads 的碱基质量值进行重新校正，使最后输出的 bam 文件中 reads 中碱基的质量值能够更加接近真实的与参考基因组之间错配的概 率。这一步适用于多种数据类型，包括 illunima、solid、454、CG 等数据格式。 在 GATK2.0 以上版本中还可以对 indel 的质量值进行校正，这

28、一步对 indel calling 非常有帮助举例说明，在 reads 碱基质量值被校正之前，我们要保留质量值在 Q25 以上的 碱基，但是实际上质量值在 Q25 的这些碱基的错误率在 1%，也就是说 质量值 只有 Q20，这样就会对后续的变异检测的可信度造成影响。还有，在边合成边 测序的测序过程中，在 reads 末端碱基的错误率往往要比起始部位更高。 另外， AC 的质量值往往要低于 TG。BQSR 的就是要对这些质量值进行校正。BQSR 主要有三步：第一步：利用工具 BaseRecalibrator，根据一些 known sites，生成一个校正质量 值所需要的数据文件，GATK 网站以

29、“.grp”为后缀命名。e.g.java -jar GenomeAnalysisTK.jar-T BaseRecalibrator-R hg19.fa-I ChrALL.100.sam.dedup.realn.07.bam-knownSites dbsnp_137.hg19.vcf-knownSites Mills_and_1000G_gold_standard.indels.hg19.vcf-knownSites 1000G_phase1.indels.hg19.vcf-o ChrALL.100.sam.recal.08-1.grp第二步：利用第一步生成的 ChrALL.100.sam.re

30、cal.08-1.grp 来生成校正后的数据 文件，也是以“.grp”命名，这一步主要是为了与校正之前的数据进行比较，最后 生成碱基质量值校正前后的比较图，如果不想生成最后 BQSR 比较图，这一步 可以省略。e.g.java -jar GenomeAnalysisTK.jar-T BaseRecalibrator-R hg19.fa-I ChrALL.100.sam.dedup.realn.07.bam-BQSR ChrALL.100.sam.recal.08-1.grp-o GATK/hg19.recal.08-2.grp-knownSites dbsnp_137.hg19.vcf-kno

31、wnSites Mills_and_1000G_gold_standard.indels.hg19.vcf-knownSites 1000G_phase1.indels.hg19.vcf第三步：利用工具 PrintReads 将经过质量值校正的数据输出到新的 bam 文件中， 用于后续的变异检测。e.g.java -jar GenomeAnalysisTK.jar-T PrintReads-R hg19.fa-I ChrALL.100.sam.dedup.realn.07.bam-BQSR ChrALL.100.sam.recal.08-1.grp-o ChrALL.100.sam.recal

32、.08-3.grp.bam主要参数说明：-bqsrBAQGOP：BQSR BAQ gap open 罚值，默认值是 40，如果是 对全基因组数据进行 BQSR 分析，设置为 30 会更好。-lqt： 在计算过程中，该算法所能考虑的 reads 两端的最小质量 值。如果质量值小于该值，计算过程中将不予考虑，默认值是 2。注意：（1）当 bam 文件中的 reads 数量过少时，BQSR 可能不会正常工 作，GATK 网站建议 base 数量要大于 100M 才能得到比较好的结果。（2）除非你所研究的样本所得到的 reads 数实在太少，或者比对 结果中的 mismatch 基本上都是实际存在的变

33、异，否则必须要进行 BQSR 这一步。对于人类基因组，即使有了 dbSNP 和千人基因组的 数据，还有很多 mismatch 是错误的。因此，这一步能做一定要做。11.11. 分析和评估分析和评估 BQSRBQSR 结果结果这一步会生成评估前后碱基质量值的比较结果，可以选择使用图片和表格的形 式展示。e.g.java -jar GenomeAnalysisTK.jar-T AnalyzeCovariates-R hg19.fa-before ChrALL.100.sam.recal.08-1.grp-after ChrALL.100.sam.recal.08-2.grp-csv ChrALL.

34、100.sam.recal.grp.09.csv-plots ChrALL.100.sam.recal.grp.09.pdf参数解释：-before： 基于原始比对结果生成的第一次校对表格。-after： 基于第一次校对表格生成的第二次校对表格。-plots： 评估 BQSR 结果的报告文件。-csv： 生成报告中图标所需要的所有数据。12.Reduce12.Reduce bambam filefile这一步是使用 ReduceReads 这个工具将 bam 文件进行压缩，生成新的 bam 文件， 新的 bam 文件仍然保持 bam 文件的格式和所有进行变异检测所需要的信息。这 样不仅能够节省

35、存储空间，也方便后续变异检测过程中对数据的处理。e.g.java -jar GenomeAnalysisTK.jar-T ReduceReads-R hg19.fa-I ChrALL.100.sam.recal.08-3.grp.bam-o ChrALL.100.sam.recal.08-3.grp.reduce.bam到此为止，GATK 流程中的第一大步骤就结束了，完成了 variants calling 所需 要的所有准备工作，生成了用于下一步变异检测的 bam 文件。变异检测变异检测1.1. VariantVariant CallingCallingGATK 在 这一步里面提供了两个工具

36、进行变异检测UnifiedGenotyper 和 HaplotypeCaller。其中 HaplotypeCaller 一直 还在开发之中，包括生成的结果以 及计算模型和命令行参数一直在变动，因此，目前使用比较多的还是 UnifiedGenotyper。此 外，HaplotypeCaller 不支持不支持 Reduce 之后的之后的 bam 文件文件， 因此，当选择使用 HaplotypeCaller 进行变异检测时，不需要进行 Reduce reads。UnifiedGenotyper 是集合多种变异检测方法而成的一种 Variants Caller，既可以 用于单个样本的变异检测，也可以用

37、于群体的变异检测。UnifiedGenotyper 使用 贝叶斯最大似然模型，同时估计基因型和基 因频率，最后对每一个样本的每一 个变异位点和基因型都会给出一个精确的后验概率。e.g.java -jar GenomeAnalysisTK.jar-glm BOTH-l INFO-R hg19.fa-T UnifiedGenotyper-I ChrALL.100.sam.recal.08-3.grp.reduce.bam-D dbsnp_137.hg19.vcf-o ChrALL.100.sam.recal.10.vcf-metrics ChrALL.100.sam.recal.10.metric

38、s-stand_call_conf 10-stand_emit_conf 30上述命令将对输入的 bam 文件中的所有样本进行变异检测，最后生成一个 vcf 文件，vcf 文件中会包含所有样本的变异位点和基因型信息。但是现在所 得到 的结果是最原始的、没有经过任何过滤和校正的 Variants 集合。这一步产生的 变异位点会有很高的假阳性，尤其是 indel，因此，必须要进行进一 步的筛选 过滤。这一步还可以指定对基因组的某一区域进行变异检测，只需要增加一个 参数 -L：target_interval.list，格式是 bed 格式文件。主要参数解释：-A： 指定一个或者多个注释信息，最后输出

39、到 vcf 文件中。-XA： 指定不做哪些注释，最后不会输出到 vcf 文件中。-D： 已知的 snp 文件。-glm： 选择检测变异的类型。SNP 表示只进行 snp 检测；INDEL 表示只对indel 进行检测；BOTH 表示同时检测 snp 和 indel。默认值是 SNP。-hets： 杂合度的值，用于计算先验概率。默认值是 0.001。-maxAltAlleles： 容许存在的最大 alt allele 的数目，默认值是 6。这个参数要特别注意，不要轻易修改默认值，程序设置的默认值几乎可以满足所有的分析，如果修改了可能会导致程序无法运行。-mbq： 变异检测时，碱基的最小质量值。如

40、果小于这个值，将不会对其进行变异检测。这个参数不适用于 indel 检测，默认值是 17。-minIndelCnt： 在做 indel calling 的时候，支持一个 indel 的最少 read 数量。也就是说，如果同时有多少条 reads 同时支持一个候选 indel 时，软件才开始进行 indel calling。降低这个值可以增加 indel calling 的敏感度，但是会增加耗费的时间和假阳性。-minIndelFrac： 在做 indel calling 的时候，支持一个 indel 的 reads 数量占比对到该 indel 位置的所有 reads 数量的百分比。也就是说，只

41、有同时满足 -minIndelCnt 和-minIndelFrac 两个参数条件时，才会进行 indel calling。-onlyEmitSamples：当指定这个参数时，只有指定的样本的变异检测结果会输出到 vcf 文件中。-stand_emit_conf：在变异检测过程中，所容许的最小质量值。只有大于等于这个设定值的变异位点会被输出到结果中。-stand_call_conf：在变异检测过程中，用于区分低质量变异位点和高质量变异位点的阈值。只有质量值高于这个阈值的位点才会被视为高 质量的。低于这个质量值的变异位点会在输出结果中标注 LowQual。在千人基因组计划第二阶段的变异检测时，利用

42、 35x 的数据进行 snp calling 的时候，当设置成 50 时，有大概 10%的假阳性。-dcov： 这个参数用于控制检测变异数据的 coverage(X)，4X 的数据可以设置为 40，大于 30X 的数据可以设置为 200。注意：GATK 进行变异检测的时候，是按照染色体排序顺序进行的（先 call chr1，然后 chr2，然后 chr3最后 chrY），并非多条染色体并行检测的，因此， 如果数据量比较大的话，建议分染色体分别进行，对性染色体的变异检测可以 同常染色体方法。大多数参数的默认值可以满足大多数研究的需求，因此，在做变异检测过程中， 如果对参数意义不是很明确，不建议修

43、改。2.2. 对原始变异检测结果进行过滤（对原始变异检测结果进行过滤（hardhard filterfilter andand VQSRVQSR）这一步的目的就是对上一步 call 出来的变异位点进行过滤，去掉不可信的位点。 这一步可以有两种方法，一种是通过 GATK 的 VariantFiltration，另一种是通过 GATK 的 VQSR（变异位点质量值重新校正）进行过滤。通过 GATK 网站上提供的最佳方案可以看出，GATK 是推荐使用 VASR 的，但 使用 VQSR 数据量一定要达到要求，数据量太小无法使用高斯模型。还有，在 使用 VAQR 时，indel 和 snp 要分别进行。

44、VQSR 原理介绍：这个模型是根据已有的真实变异位点（人类基因组一般使用 HapMap3 中的位点， 以及这些位点在 Omni 2.5M SNP 芯片中出现的多态位点）来训练，最后得到一 个训练好的能够很好的评估真伪的错误评估模型，可以叫他适应性错误评估模 型。这个适应性的错误评估模型可 以应用到 call 出来的原始变异集合中已知的 变异位点和新发现的变异位点，进而去评估每一个变异位点发生错误的概率， 最终会给出一个得分。这个得分最后会 被写入 vcf 文件的 INFO 信息里，叫做 VQSLOD，就是在训练好的混合高斯模型下，一个位点是真实的概率比上这个 位点可能是假阳性的概率的 log

45、odds ratio（对数差异比），因此，可以定性的 认为，这个值越大就越好。VQSR 主要分两个步骤，这两个步骤会使用两个不同的工具： VariantRecalibrator 和 ApplyRecalibration。i)i) VariantRecalibratorVariantRecalibratorVariantRecalibrator：通过大量的高质量的已知变异集合的各个 注释（包括很多种，后面介绍）的值来创建一个高斯混合模型， 然后用于评估所有的变异位点。这个文件最后将生成一个 recalibration 文件。原理简单介绍： 这个模型首先要拿到真实变异数据集和上一步骤 中得到的原始

46、变异数据集的交集，然后对这些 SNP 值相对于具体 注释信息的分布情况进行模拟，将这些变异位点进 行聚类，最后 根据聚类结果赋予所有变异位点相应的 VQSLOD 值。越接近聚类核 心的变异位点得到的 VQSLOD 值越高。ii)ii) ApplyRecalibrationApplyRecalibrationApplyRecalibration：这一步将模型的各个参数应用于原始 vcf 文件中的每一个变异位点，这时，每一个变异位点的注释信息列 中都会出现一个 VQSLOD 值，然后模型会根据这个值对变异位点进 行过滤，过滤后的信息会写在 vcf 文件的 filter 一列中。原理简单介绍： 在

47、VariantRecalibrator 这一步中，每个变异位 点已经得到了一个 VQSLOD 值了，同时，这些 LOD 值在训练集里也进行了排序。当你在这 一步中设置一个 tranche sensitivity 的阈值（这个阈值一般是一个百分数，如设置成 99%），那么， 如果 LOD 值从大到小排序的话，这个程序就会认为在这个训练集 中，LOD 值在前 99%的是可 信的，当这个值低于这个阈值，就认 为是错误的。最后，程序就会用这个标准来过滤上一步 call 出来 的原始变异集合。如果 LOD 值超过这个阈值，在 filter 那一列 就会显示 PASS，如果低于这个值就会被过滤掉，但是这些

48、位点仍 然会显示在结果里面，只不过会在 filter 那一列标示出他所属于 的 tranche sensitivity 的名称。在设置 tranche sensitivity 的阈值时，要兼顾敏感度和质量值。初步分析初步分析这一步主要是对上面所得到的最终 vcf 中的结果进行一些初步的分析，比如计 算这些变异位点在 dbsnp 中的比例、Ti/Tv 的比例、每个样本中的 snp 数 量。此外，还可以对变异位点的同义同义/非同义突变进行统计非同义突变进行统计，识别是否为 CpG 位点以及氨基酸的简并信息等。这一步主要是利用 GATK 中的 VariantEval 来完成。需要注意的是，有些计算内

49、容不能同时进行，例如 AlleleCount 和 VariantSummary 或者 Sample 和 VariantSummary。如果选择了这样的组合方式， 程序就会报错。但是 GATK 并没有告诉我们到底哪些不能同时运行，所以当选 择计算内容的时候可以先做一下测试。e.g.java -jar GenomeAnalysisTK.jar-R hg19.fa-T VariantEval-eval hg19.snp.filter.t97.Q10_13.both.vcf-D dbsnp_137.hg19.vcf-o hg19.PASS.Eval_15_Final.gatkreport主要参数解释：-eval 输入要进行 summary 的文件，也就是 hg19.snp.filter.t97.Q10_13.both.vcf。-EV 选择模块计算相应的分析内容，。-list 列出可供选择的计算模块。-noEV 不是用默认的模块，只计算用-EV 选定的模块。