-
GATK
使用方法详解
一、使用
GATK
前须知事项：
<
/p>
（
1
）对
GAT
K
的测试主要使用的是人类全基因组和外显子组的测序数据，而且
全部是基于
illumina
数据格式，目前还没有提供其他
格式文件（如
Ion Torrent
）
或者实验设计（
RNA-
Seq
）的分析方法。
（
2
）
GATK
是一个应用
于前沿科学研究的软件，不断在更新和修正，因此，在
使用
GA
TK
进行变异检测时，最好是下载最新的版本，目前的版本是
2
.8.1
（
2014-02-25
）。
下载网站：
/gatk/download
。
< br>
（
3
）在
GATK
使用过程中（见下面图），有些步骤需要用到已知变异信息，对
于这些已知变
异，
GA
TK
只提供了人类的已知变异信息，可以在
GATK
的
FTP
站点下载（
GA
TK resource bundle
）。如果要研究的不是人类基因组，需要自行<
/p>
构建已知变异，
GATK
提供了详细的构
建方法。
（
4
）
GATK
在进行
BQSR
和
VQSR
的过程中会使用到
< br>R
软件绘制一些图，
因此，
在运
行
GATK
之前最好先检查一下是否
正确安装了
R
和所需要的包，所需要的
包大概包括
ggplot2
、
gplo
ts
、
bitops
、
caTools
、
colo
rspace
、
gdata
、
gsalib
、
reshape
、
RColorBrewer
等。如果画图时出现错误
，会提示需要安装的包的名称。
二、
GATK
的使用流程
GATK
最佳使用方案：
共
3
大步骤，即
:
原始数据的处理
-->
变异检测
-->
初步分析
。
原始数据的处理
1.
对原始下机
fastq
文件进行过滤和比对（
mapping
）
对于
Illumina
下机数据推荐使用
bwa
进行
mapping
。
Bwa
比对步骤大致如下：
（
1
）对参考基因组构建索引：
例子：
bwa index -a bwtsw
。
构建索引时需要注意的问题：
bwa
构建索引有两种算法，
两种算法都是
基于
BWT
的，
这两种算法通过参数<
/p>
-a is
和
-a
bwtsw
进行选择。
其中
-a bw
tsw
对于短的参考
序列是不工作的，必须要大于等于
10Mb
；
-a is
是默认参数，这个参数不适用于
大的参考序列，必须要小于等于
2G
。
（
2
）寻找输入
reads
文件的
SA
坐标。
对于
pair end
数据，每个
reads
文件单独做运算，
single
end
数据就不用说了，只
有一个文件。
pair end
：
bwa aln -t 4 -I >
bwa aln -t 4 -I >
single end
：
bwa
aln
-l
30
-k
2
-t
4
-I
>
主要参数说明：
-o int
：允许出现的最大
gap
数。
-e int
：每个
gap
允许的最大长度。
-d int
：不允许在
3
’端出现大于多少
bp
的
deletion
。
-i int
：不允许在
reads
两端出现大于多少
bp
的
indel
。
-l
int
：
Read
前多少个碱基作为
seed
，
如果设置的
seed
大于
read
长度，将无法继续，最好设置在
25-3
5
，与
-k 2
配合使用。
-k
< br>int
：在
seed
中的最大编
辑距离，使用默认
2
，与
-l
配合使用。
-t
int
：要使用的线程数。
-R
int
：此参数只应用于
pair end
中，当没有出现大于此值的最
佳比对结果时，
将会降低标准
再次进行比对。
增加这个值可以提高
配对比对的准确率，但是同
时会消耗更长的时间，默认是
32
。
-I
int
：表示输入的文件格式为
Illumina
1.3+
数据格式。
-B int<
/p>
：设置标记序列。从
5
’端开始多少个碱
基作为标记序列，
当
-B
为正值时，在
比对之前会将每个
read
的标记序列剪切，并将
此标记序列表示在
BC SAM
标签里，对于
pair end
数据，
两端的
标记序列会被连接。
-b <
/p>
：指定输入格式为
bam
格式。这是一个
很奇怪的功能，就是对
其它软件的
bam
文件进行重新比对的意思
bwa aln >
（
3
）生成
sam
格式的比对文件。如果一条
re
ad
比对到多个位置，会随机选择一
种。
例子：
single
end
：
bwa samse
>
参数：
-n int
：如果
reads
比对次数超过多少次，就不在
XA
标签显示。
-r str
：定义头文件。‘
@RG
tID:footSM:bar
’，如果在此步骤
不进行头文件
定义，
在后续
GATK
分析中还是
需要重新增加头文件。
pair
end
：
bwa sampe -a
500 >
参数：
-a
int
：最大插入片段大小。
-o
int
：
pair end
两
reads
中其中之一所允许配对的最大次数，
超过该次数，将被视为
single end
。降低这个参数，可以加快运算速度，对于少于
30bp
的
read
，建议降低
-o
值。
-r
str
：定义头文件。同
single
end
。
-n int
：每对
reads
输出到结果中的最多比对数。
对于最后得到的
sam
< br>文件，将比对上的结果提取出来（
awk
即可处理），即
可直
接用于
GATK
的分析。
注意：由于
GATK
在下游的
snp-
calling
时，是按染色体进行
call-snp
的。因此，
在准备原始
sam
< br>文件时，可以
先按染色体将文件分开，这样会提高运行速度
。
但是当数据量不足时，可能会影响后续的
VQSR
分析，这是需要注意的。
2.
对
sam
文件进行进行重新排序（
< br>reorder
）
由
BWA
生成的
sam
文
件时按字典式排序法进行的排序
（
lexicographic
ally
）
进行
排序的（
chr10
，
chr11…chr19
，
chr1
，
chr20…
chr22
，
chr2
，
chr3…chrM
，
chrX
，
chrY
），但是
GATK
在进行
callsnp
的时候是按照染色体
组型（
karyotypic
）进行
的
（
chrM
，
chr1
，
chr2…chr22
，<
/p>
chrX
，
chrY
）
，
因此要对原始
sam
文件进行
reorder
。
可以使用
picard-tools
中的
< br>ReorderSam
完成。
eg.
java -jar picard-
tools-1.96/
I=
O=r_
REFERENCE=
注意：
1)
这一步的头文件可以人工加上，同时要确保头文件中有的序号在下面序列中
也有对应的
。
虽然在
GATK
网站上的说明
chrM
可以在最前也可以在最后，
但是
当把
chrM
放在最后时可能会出错。
2)
在进行排序之前，要先构建参考序
列的索引
。
e.g.
samtools faidx
。最后生成的索引文件：
。
3)
如果在
上一步想把大文件切分成小文件的时候，头文件可以自己手工加上，
之后运行这一步就好
了。
3.
将
sam
文件转换成
bam
文件（
bam
是二进制文件，运算速度快）<
/p>
这一步可使用
samtools
view
完成。
e.g.
samtools view -bS r_ -o _
4.
对
bam
文件进行
sort
排序处理
这一步是将
sam
文件中
同一染色体对应的条目按照坐标顺序从小到大进
行排序
。
可以使用
picard-to
ols
中
SortSam
完成。
e.g.
java -jar
picard-tools-1.96/
INPUT=_
OUTPUT=_
SORT_ORDER=coordinate
5.
对
bam
文件进行加头（
head
）处理
GATK
2.0
以上版本将不再支持无头文件的变异检测。加头这一步可以在
BWA
比对的时候进行，通过
-r
参数的选择可以完成。如果
在
BWA
比对期间没有选择
-r
参数，可以增加这一步骤。可使用
picard-tools
中
AddOrReplaceReadGroups
完
成。
e.g.
java -jar picard-tools-1.96/
I=_
O=d_
ID=hg19ID
LB=hg19ID
PL=illumine
PU=hg19PU
SM=hg19
ID str
：
输入
reads
集
ID
号；
LB
：
read
集文库名；
PL
：
测序平台
（
illun
ima
或
solid
）
；
PU
：测序平台下级单位名称（
run
的名称）；
SM
：样本名称
。
注意：
这一步尽量不要手动加头，
本人尝试过多次手工加头，
虽然看起来与软件
< br>加的头是一样的，但是程序却无法运行。
6. Merge
如果一个样本分为多个
lane
进行测序，那么在进行下一步之前可以将每个
la
ne
的
bam
文件合并。
e.g.
java -jar picard-
tools-1.70/
INPUT=
INPUT=
INPUT=
INPUT=
……
INPUT=
OUTPUT=
7. Duplicates Marking
在制备文库的
过程中，由于
P
CR
扩增过程中会存在
一些偏差，也就是说有的序
列会被过量扩增
。
< br>这样，
在比对的时候，
这些过量扩增出来的完全相同的序
列
就
会比对到基因组的相同位置。而
这些过量扩增的
reads
并不是基因组自身固有序
列，不能作为变异检测的证据，因此，要尽量去除这些由
PCR
扩增所形成的
duplicates
，这一步可以使用
picard-tools
来完
成。去重复的过程是给这些序列设
置一个
flag
以标志它们，方便
GATK
的识别。还可以设置
REMOVE_DUPLICATES=true
来丢弃
duplicated
序列。对于是否选择标
记或者
删除，
对结果应该没有什么影响，
GATK
官方流程里面给出的例子是仅做标记不
删除。这里定
义的重复序列是这样的：如果两条
reads
具有相同的长度而
且比对
到了基因组的同一位置，那么就认为这样的
reads<
/p>
是由
PCR
扩增而来，就会被
GATK
标记。
e.g.
java -jar picard-
tools-1.96/
REMOVE_DUPLICATES= false
MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000
INPUT=d_
OUTPUT=_
METRICS_FILE=_s
注意：
dedup
< br>这一步只要在
library
层面上进行就可以了，例如
一个
sample
如果
建了多个库的话
，对每个库进行
dedup
即可，不需要把所有库合成一
个
sample
再
进行
dedup
操作。
其实并不能准确
的定义被
mask
的
reads
到底是不是
duplicates
，
重复序列的程度与测序深度和文库类型
都有关
系。
最主要目的就是尽量减小文
库构建时引入文库的
PCR bias
。
8.
要对上一步得到的结果生成索引文件
可以用
samtools
完成，生成的索引后缀是
bai
。
e.g.
samtools index _
realignment around indels
这一步的目的就是将比对到<
/p>
indel
附近的
reads
进行局部重新比对，将比对的错
误率降到最低。一般来说，绝大部分需要进
行重新比对的基因组区
域，都是因
为
插入
/
缺失的存在，因为在
indel
附近的比对会出现大量的碱基错配，这些碱基
的错配很容易被误
认为
SNP
。还有，在比对过程中，比对算法对
于每一条
read
的处理都
是独立的，不可能同时把多条
reads
与参考基因组比对来排
错。因此，
即使有一些
reads
能够
正确的比对到
indel
，
但那些
恰恰比对到
indel
< br>开始或者结
束位置的
read
也
会有很高的比对错误率，这都是需要重新比对的。
Local
realignment
就是将由
indel
< br>导致错配的区域进行重新比对，
将
indel
附近的比对错
误率降到最低。
主要分为两步：
第一步，通过运行<
/p>
RealignerTargetCreator
来确定要进行重
新比对的区域。
e.g.
java
-jar
-R
-T RealignerTargetCreator
-I _
-o
_als
-known
Mills_and_1000G_gold_
-known
1000G_
参数说明：
-R
：
参考基因组；
-T
：
选择
的
GATK
工具；
-I
：
输入
上一步所得
bam
文件；
-o
：
输出的需要重新比对的基因组区域结果；
-maxInterval
：
允许进行重新比对的基因组区域的最大值，不能
太大，太大耗费会很长
时间，默认值
500
；
-known
：
已知的可靠的
indel
位点，指定已知的可靠的
indel
位
点，重比对将主要围绕这些位点进
行，对于人类基因组数据而言，
可以直接指定
GATK
resource bundle
里面的
in
del
文件（必须是
vcf
文件）。<
/p>
对于
known sites
的选择很重要，
GATK
中
每一个用到
known
sites
的工具对于
known sites
的使用都是不一样的，
但是所有的都有一个共同目的，
那就是分辨真
实的变异位点和不可信的变异位点。如果不提供这些
known sites
的话，这些统
计工具就会产生偏差，
最后会严重影响结果的可信度。在这些需要知道
known
s
ites
的工具里面，只有
UnifiedGenotyper
和
HaplotypeCaller
对
known
sites
没有
太严格的要求。
如果你所研究的对象是人类基因组的话，那就简单多了，因为
GAT
K
网站上对
如何使用人类基因组的
kn
own sites
做出了详细的说明，
具体的选择方法如下表
，
这些文件都可以在
GATK
resource bundle
中下载。
Mills
dbSNP 129
dbSNP >132
indels
RealignerTargetCreator
X
IndelRealigner
X
Tool
1KG
indels
HapMap
Omni
X
X
-
-
-
-
-
-
-
-
-
上一篇：GPS术语
下一篇：计量经济学(重要名词解释).(精选)