GATK4.0全基因组和全外显子组分析实战

前言

GATK是目前业内最权威、使用最广的基因数据变异检测工具。相比samtools + bcftools call SNP/Indel，GATK更加精确，当然代价是流程复杂且耗时长。目前已经更新到GATK4，GATK3在官网已没有下载通道但其使用仍然很广泛。GATK4在核心算法层面并没太多的修改，但参数设置还是有些改变的，并且取消了RealignerTargetCreator、IndelRealigner，应该是HaplotypeCaller继承了这部分功能。GATK4使用了新的设计模式，做了很多功能的整合，已经把picard完全整合。本文使用的是GATK4.0.2.0，参考了其他人编写的GATK3x和GATK4x教程，对全基因组和全外显子组的每个步骤都进行了讲解并放上脚本，但读者还是要略作修改才能本地运行成功，所以具有基本的生信知识和编程、Linux技能才能更好的学习这篇教程，下面开始正文。

软件

需要fastqc、fastp、BWA、samtools、GATK、annovar（需要学术邮箱才能注册下载）

数据质控

这部分不赘述，公司交给你的数据肯定质量很高，否则不会交付给你，但是为了放心还是要检查下的。先用fastqc

wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip

unzip fastqc_v0.11.7.zip

cd FastQC

chmod755 fastqc

使用：

fastqc wes_1.fq.gz -o fastqc_out_dir/ && fastqcwes_2.fq.gz -o fastqc_out_dir/

数据过滤使用fastp

安装：wget http://opengene.org/fastp/fastp

chmod 755 fastp

使用：./fastp -i in.R1.fq.gz -o out.R1.fq.gz -I in.R2.fq.gz -O out.R2.fq.gz

比对

第一步下载参考基因组

for i in $(seq 1 22) X Y M;

do

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr$.fa.gz &

Done

下载完成后解压合并

gunzip *.gz

for i in $(seq 1 22) X Y M;

do cat chr$.fa >> hg19.fa;

Done

构建索引

bwa index -a bwtsw -p hg19 hg19.fa &

程序运行时间较长，建议使用nohup 或者 screen

最终生成文件 hg19.amb hg19.ann hg19.bwt hg19.pac hg19.sa

bwa比对

-t，线程数；

-M , -M 将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件;

-R 接的是 Read Group的字符串信息，它是用来将比对的read进行分组的，这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要。不设置-R参数不会报错，但使用GATK时是会报错的。

(1)ID，这是Read Group的分组ID，一般设置为测序的lane ID

(2) PL，指的是所用的测序平台

(3) SM，样本ID

(4) LB，测序文库的名字

这些信息设置好之后，在RG字符串中要用制表符（\t）将它们分开。

上一步生成的SAM文件是文本文件，一般整个文件都非常巨大，因此，为了有效节省磁盘空间，一般都会用samtools将它转化为BAM文件（SAM的特殊二进制格式），而且BAM会更加方便于后续的分析。

samtools view -b -S abc.sam > abc.bam ##由sam生成bam

samtools view -h abc.bam>abc.sam##由bam生成sam，-h代表是否带上header

想看bam文件使用命令samtools view -h abc.bam | less

bam文件是一个重要的文件，每一列的意义都需要理解。

前五列分别为：reads名flag染色体位置比对质量。第十列时read序列，11列是碱基质量值。第1，10，11列可以提取出来还原成我们的测序数据fastq格式。有人说第五列比对质量值为0时代表该read可以比对到基因组多个位置，说法是否正确我不确定。

第二列的flag包含很多信息，flag解释网站 http://broadinstitute.github.io/picard/explain-flags.html ，打开此网站

输入flag值便可以看到包含的信息，其中read unmapped是4，mate unmapped是8，如果双端reads都不能map到参考基因组上，flag值就包括4和8，和为12。可以用命令samtools view –f 12 wes.bam > unmapped.sam 提取双端reads都没有比对成功的reads，samtools view –F 12 wes.bam > mapped.sam 代表过滤双端都比对不上的reads。

sort 排序

BWA比对后输出的bam文件是没顺序的，比对后得到的结果文件中，每一条记录之间位置的先后顺序是乱的，我们后续去重复等步骤都需要在比对记录按照顺序从小到大排序下来才能进行。

做类似分析的时候在文件名字将所做的关键操作包含进去，因为这样即使过了很长时间，当你再去看这个文件的时候也能够立刻知道当时对它做了什么。

bam信息统计

做到这一步需要对序列比对情况进行统计，如果比对情况很差需要查找原因。

覆盖度，深度等信息统计

覆盖度和深度是我们关心的重要参数，如果是全外显子组可以用picard（已经整合到GATK4中）进行统计。

生成参考基因组的dict文件

$GATK CreateSequenceDictionary -R hg19.fa -O hg19.dict

生成interval

外显子组是用试剂盒捕获再进行测序，不同试剂盒捕获的区域不同，要下载相应的包含捕获区域bed文件，本文用的是安捷伦的捕获区域的bed文件。

覆盖度，深度等信息统计

On targeted bases相对总bases的比例（PCT_USABLE_BASES_ON_BAIT）

On and near targeted bases相对总bases的比例（PCT_SELECTED_BASES）

MEAN_TARGET_COVERAGE平均覆盖度

如果是全基因组wgs，运行以下命令，

生成的wgs.metrics包含多种信息，可自行查阅研究。

去除由于PCR扩增引起的重复reads

在NGS测序之前都需要先构建测序文库：通过物理（超声）打断或者化学试剂（酶切）切断原始的DNA序列，然后选择特定长度范围的序列去进行PCR扩增并上机测序。这个过程中产生的重复reads，增大了变异检测结果的假阴率和假阳率！！原因如下：

1.PCR反应过程中也会带来新的碱基错误。发生在前几轮的PCR扩增发生的错误会在后续的PCR过程中扩大，同样带来假的变异；

2. PCR反应可能会对包含某一个碱基的DNA模版扩增更加剧烈（这个现象称为PCR Bias）

3.如果某个变异位点的变异碱基都是来自于PCR重复，而我们却认为它深度足够判断是真的变异位点，这个结论其实有很大可能是假阳性。

利用picard标记重复序列

查看被标记重复的reads

samtools view –f 1024wes.sorted.MarkDuplicates.bam| less

为何是1024去查阅flag解释网站。

下一步为wes.sorted.MarkDuplicates.bam创建索引文件，它的作用能够让我们可以随机访问这个文件中的任意位置，而且后面的步骤也要求这个BAM文件一定要有索引.

samtools index wes.sorted.MarkDuplicates.bam

变异检测

开始检测变异前还是要做一些准备工作，首先是重新校正碱基质量值（BQSR）。

变异检测是一个极度依赖测序碱基质量值，因为这个质量值是衡量我们测序出来的这个碱基到底有多正确的重要指标。它来自于测序图像数据的base calling，因此，基本上是由测序仪和测序系统来决定的，计算出来的碱基质量值要么高于真实结果，要么低于真实结果。

BQSR（Base Quality Score Recalibration）这个步骤就是为此而存在的，这一步非常重要。它主要是通过机器学习的方法构建测序碱基的错误率模型，然后对这些碱基的质量值进行相应的调整。

这里包含了两个步骤：

第一步，BaseRecalibrator，这里计算出了所有需要进行重校正的read和特征值，然后把这些信息输出为一份校准表文件（wes.recal_data.table）

第二步，ApplyBQSR，这一步利用第一步得到的校准表文件（wes.recal_data.table）重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新的BAM文件。

-R $GENOME -I wes.sorted.MarkDuplicates.bam \

--known-sites $hg19_VCF/dbsnp_138.hg19.vcf \

-O wes.recal_data.table

（注意：此步骤以及后面几个步骤中外显子数据要加上外显子捕获区域的bed文件，并把-ip设为reads长，全基因组数据则不需要加-L 和 -ip）

$GATK --java-options "-Xmx8G -Djava.io.tmpdir=./" ApplyBQSR \

-R $GENOME -I wes.sorted.MarkDuplicates.bam \

-bqsr wes.recal_data.table \

-O wes.sorted.MarkDuplicates.BQSR.bam

变异检测前确定bam文件是否符合GATK要求，运行

$GATK ValidateSamFile -I wes.sorted.MarkDuplicates.BQSR.bam，如果显示 no error，则可以用HaplotypeCaller call SNP/Indel。

利用samtools为hg19.fa创建一个索引 samtools faidx hg19.fa

利用这个索引可以查看参考基因组任何位置，如运行

利用HaplotypeCaller检测突变还需利用许多数据库，这些数据库包含以前研究过的突变位点，利用这些位点提高变异检测的准确率。

运行上面命令下载所需数据库。注意下载文件均为压缩文件，需要解压才能使用。

下面开始HaplotypeCaller突变检测，

HaplotypeCaller的应用有两种做法，区别在于是否生成中间文件gVCF：

（1）直接进行HaplotypeCaller，这适合于单样本，只执行一次HaplotypeCaller。如果多样本，你每增加一个样本数据都需要重新运行这个HaplotypeCaller，而这个时候算法需要重新去读取所有人的BAM文件，浪费大量时间；

（2）每个样本先各自生成gVCF，然后再进行群体joint-genotype。gVCF全称是genome VCF，是每个样本用于变异检测的中间文件，格式类似于VCF，它把joint-genotype过程中所需的所有信息都记录在这里面，文件无论是大小还是数据量都远远小于原来的BAM文件。这样一旦新增加样本也不需要再重新去读取所有人的BAM文件了，只需为新样本生成一份gVCF，然后重新执行这个joint-genotype就行了。

第一种方法

-R $GENOME -I wes.sorted.MarkDuplicates.BQSR.bam \

第二种方法

#1 生成中间文件gvcf

-R $GENOME --emit-ref-confidence GVCF \

-I wes.sorted.MarkDuplicates.BQSR.bam -D $hg19_VCF/dbsnp_138.hg19.vcf -L S31285117_Regions.bed -ip 90 -O wes.gvcf

#2 通过gvcf检测变异

若有多个样本的gvcf文件，运行$GATK CombineGVCFs -V 1.gvcf –V 2.gvcf …… -O final.gvcf ，再用final.gvcf运行下一步

变异质控和过滤

质控的含义和目的是指通过一定的标准，最大可能地剔除假阳性的结果，并尽可能地保留最多的正确数据。

第一种方法 GATK VQSR，它通过机器学习的方法利用多个不同的数据特征训练一个模型（高斯混合模型）对变异数据进行质控，使用VQSR需要具备以下两个条件：

第一，需要一个精心准备的已知变异集，它将作为训练质控模型的真集。比如，Hapmap、OMNI，1000G和dbsnp等这些国际性项目的数据，这些可以作为高质量的已知变异集。

第二，要求新检测的结果中有足够多的变异，不然VQSR在进行模型训练的时候会因为可用的变异位点数目不足而无法进行。适合全基因组分析。

第二种方法通过过滤指标过滤。

QualByDepth（QD）

FisherStrand (FS)

StrandOddsRatio (SOR)

RMSMappingQuality (MQ)

MappingQualityRankSumTest (MQRankSum)

ReadPosRankSumTest (ReadPosRankSum)

GATK VQSR

-resource hapmap,known=false,training=true,truth=true,prior=15.0:$hg19_VCF/hapmap_3.3.hg19.sites.vcf \

-resource omini,known=false,training=true,truth=false,prior=12.0:$hg19_VCF/1000G_omni2.5.hg19.sites.vcf \

-resource dbsnp,known=true,training=false,truth=false,prior=6.0:$hg19_VCF/dbsnp_138.hg19.vcf \

-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an SOR -mode SNP \

$GATK ApplyVQSR -V wes.raw.vcf -O wes.snps.VQSR.vcf \

此方法要求新检测的结果中有足够多的变异，不然VQSR在进行模型训练的时候会因为可用的变异位点数目不足而无法进行。可能很多非人的物种在完成变异检测之后没法使用GATK VQSR的方法进行质控，一些小panel、外显子测序，由于最后的变异位点不够，也无法使用VQSR。全基因组分析或多个样本的全外显子组分析适合用此方法。

通过过滤指标过滤

# 使用SelectVariants，选出SNP

$GATK SelectVariants -select-type SNP -V wes.GATK.vcf -O wes.GATK.snp.vcf

# 为SNP作过滤

# 使用SelectVariants，选出Indel

$GATK SelectVariants -select-type INDEL -V wes.GATK.vcf -O wes.GATK.indel.vcf

# 为Indel作过滤

VCF文件

VCF文件是记录突变信息的重要文件

前五列信息为

1. 染色体(Chromosome)

2. 起始位置(Start)

3. 结束位置(End)

4. 参考等位基因(Reference Allele)

5. 替代等位基因(Alternative Allele)

ANNOVAR注释时主要也是利用前五列信息对数据库进行比对，注释变异。Info和Format信息同样重要，比如DP代表测序深度，这些内容的含义再vcf文件的开头都有介绍，请仔细阅读并理解相应内容的意义。

突变注释

vcf文件中保存的突变位点不能直接使用，只有根据已有数据库进行注释，才能知道该位点的有何功能，是否与疾病相关以及其它信息。annovar是突变注释的常用软件，ANNOVAR是一个perl编写的命令行工具，能在安装了perl解释器的多种操作系统上执行。允许多种输入文件格式，包括最常被使用的VCF格式。输出文件也有多种格式，包括注释过的VCF文件、用tab或者逗号分隔的txt文件。ANNOVAR能快速注释遗传变异并预测其功能，这个软件需要edu邮箱注册才能下载。

http://www.openbioinformatics.org/annovar/annovar_download_form.php

软件下载后无需安装直接使用，但要下载数据库才能对突变位点进行注释，annotate_variation.pl可以方便快捷下载数据库。

perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar refGene humandb/

# -buildver 表示version

# -downdb 下载数据库的指令

# -webfrom annovar 从annovar提供的镜像下载，不加此参数将寻找数据库本身的源

# humandb/ 存放于humandb/目录下

其它数据库下载

perl annotate_variation.pl--buildver hg19 --downdb gwascatalog humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb ljb26_all --webfrom annovarhumandb/ &

perl annotate_variation.pl--buildver hg19 --downdb esp6500siv2_ea --webfromannovar humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb esp6500siv2_all --webfromannovar humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb 1000g2014oct humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb cytoBand humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb avsift -webfrom annovarhumandb/ &

perl annotate_variation.pl--buildver hg19 --downdb snp138 humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb genomicSuperDups humandb/ &

perl annotate_variation.pl--buildver hg19 --downdbphastConsElements46wayhumandb/ &

perl annotate_variation.pl--buildver hg19 --downdb tfbs humandb/

annovar有三种注释方式Gene-based Annotation(基于基因的注释)，Region-based Annotation（基于区域的注释），Filter-based Annotation（基于过滤的注释）。看起来很复杂实际做起来很简单，用table_annovar.pl进行注释，可一次性完成三种类型的注释。

注释的第一步是把要注释的文件转化成annovar需要的文件，

convert2annovar.pl #将多种格式转为.avinput的程序

avinput文件由tab分割，最重要的地方为前5列，分别是:

1. 染色体(Chromosome)

2. 起始位置(Start)

3. 结束位置(End)

4. 参考等位基因(Reference Allele)

5. 替代等位基因(Alternative Allele)

annovar主要也是利用前五列信息对数据库进行比对，注释变异。

SNP注释：

用table_annovar.pl进行注释，可一次性完成三种类型的注释。

table_annovar.pl snp.avinput $humandb -buildver hg19 -out snpanno -remove -protocol refGene,cytoBand,genomicSuperDups,esp6500siv2_all -operation g,r,r,f -nastring . –csvout

# -buildver hg19 表示使用hg19版本

# -out snpanno 表示输出文件的前缀为snpanno

# -remove 表示删除注释过程中的临时文件

# -protocol 表示注释使用的数据库，用逗号隔开，且要注意顺序

# -operation 表示对应顺序的数据库的类型（g代表gene-based、r代表region-based、f代表filter-based），用逗号隔开，注意顺序

# -nastring . 表示用点号替代缺省的值

# -csvout 表示最后输出.csv文件

Indel注释同上

最后生成的注释文件如下

总结

本教程基本按照GATK4 Best PracticesGermline short variant discovery (SNPs + Indels)，GATK3与GATK4的分析思路也几乎没有变化，除了GATK4取消了RealignerTargetCreator 和 IndelRealigner 这两步，另外一些命令发生了变化。在外显子分析步骤中记得加上外显子捕获区域的bed文件（-L参数），这样可以防止检测出捕获区域以外的突变，还可以显著提高运行速度。