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

GATK使用方法详解一、使用GATK前须知事项：（1）对GATK的测试主要使用的是人类全基因组和外显子组的测序数据，而且全部是基于illumina数据格式，目前还没有提供其他格式文件（如Ion Torrent）或者实验设计（RNA-Seq）的分析方法。

（2）GATK是一个应用于前沿科学研究的软件，不断在更新和修正，因此，在使用GATK进行变异检测时，最好是下载最新的版本，目前的版本是2.8.1（2014-02-25）。

下载网站：/gatk/download。

（3）在GATK使用过程中（见下面图），有些步骤需要用到已知变异信息，对于这些已知变异，GATK只提供了人类的已知变异信息，可以在GATK的FTP 站点下载（GATK resource bundle）。

如果要研究的不是人类基因组，需要自行构建已知变异，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）对参考基因组构建索引：例子：bwa index -a bwtsw hg19.fa。

构建索引时需要注意的问题：bwa构建索引有两种算法，两种算法都是基于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 -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：在seed中的最大编辑距离，使用默认2，与-l配合使用。

-t int：要使用的线程数。

-R int：此参数只应用于pair end中，当没有出现大于此值的最佳比对结果时，将会降低标准再次进行比对。

增加这个值可以提高配对比对的准确率，但是同时会消耗更长的时间，默认是32。

-I int：表示输入的文件格式为Illumina 1.3+数据格式。

-B int：设置标记序列。

从5’端开始多少个碱基作为标记序列，当-B为正值时，在比对之前会将每个read的标记序列剪切，并将此标记序列表示在BC SAM 标签里，对于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：定义头文件。

‘@RG\tID:foo\tSM: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：每对reads输出到结果中的最多比对数。

对于最后得到的sam文件，将比对上的结果提取出来（awk即可处理），即可直接用于GATK的分析。

注意：由于GATK在下游的snp-calling时，是按染色体进行call-snp的。

因此，在准备原始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完成。

eg.java -jar picard-tools-1.96/ReorderSam.jarI=hg19.samO=hg19.reorder_00.samREFERENCE=hg19.fa注意：1) 这一步的头文件可以人工加上，同时要确保头文件中有的序号在下面序列中也有对应的。

虽然在GATK网站上的说明chrM可以在最前也可以在最后，但是当把chrM放在最后时可能会出错。

2) 在进行排序之前，要先构建参考序列的索引。

e.g. samtools faidx hg19.fa。

最后生成的索引文件：hg19.fa.fai。

3) 如果在上一步想把大文件切分成小文件的时候，头文件可以自己手工加上，之后运行这一步就好了。

3. 将sam文件转换成bam文件（bam是二进制文件，运算速度快）这一步可使用samtools view完成。

e.g. samtools view -bS hg19.reorder_00.sam -o hg19.sam_01.bam4. 对bam文件进行sort排序处理这一步是将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序。

可以使用picard-tools中SortSam完成。

e.g.java -jar picard-tools-1.96/SortSam.jarINPUT=hg19.sam_01.bamOUTPUT=hg19.sam.sort_02.bamSORT_ORDER=coordinate5. 对bam文件进行加头（head）处理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=hg19IDLB=hg19IDPL=illuminePU=hg19PUSM=hg19ID str：输入reads集ID号；LB：read集文库名；PL：测序平台（illunima或solid）；PU：测序平台下级单位名称（run的名称）；SM：样本名称。

注意：这一步尽量不要手动加头，本人尝试过多次手工加头，虽然看起来与软件加的头是一样的，但是程序却无法运行。

6. Merge如果一个样本分为多个lane进行测序，那么在进行下一步之前可以将每个lane 的bam文件合并。

e.g.java -jar picard-tools-1.70/MergeSamFiles.jarINPUT=lane1.bamINPUT=lane2.bamINPUT=lane3.bamINPUT=lane4.bam……INPUT=lane8.bamOUTPUT=sample.bam7. Duplicates Marking在制备文库的过程中，由于P CR扩增过程中会存在一些偏差，也就是说有的序列会被过量扩增。

这样，在比对的时候，这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。

而这些过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由PCR扩增所形成的duplicates，这一步可以使用picard-tools来完成。

去重复的过程是给这些序列设置一个flag以标志它们，方便GATK的识别。

还可以设置REMOVE_DUPLICATES=true 来丢弃duplicated序列。

对于是否选择标记或者删除，对结果应该没有什么影响，GATK官方流程里面给出的例子是仅做标记不删除。

这里定义的重复序列是这样的：如果两条reads具有相同的长度而且比对到了基因组的同一位置，那么就认为这样的reads是由PCR扩增而来，就会被GATK标记。

e.g.java -jar picard-tools-1.96/MarkDuplicates.jarREMOVE_DUPLICATES= falseMAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000INPUT=hg19.reorder.sort.addhead_03.bamOUTPUT=hg19.reorder.sort.addhead.dedup_04.bamMETRICS_FILE=hg19.reorder.sort.addhead.dedup_04.metrics注意：dedup这一步只要在library层面上进行就可以了，例如一个sample如果建了多个库的话，对每个库进行dedup即可，不需要把所有库合成一个sample 再进行dedup操作。