第9章

章节实例：第9章 变异检测


GATK（全称Genome Analysis Toolkit，https://software.broadinstitute.org/gatk/）是由美国哈佛-麻省理工的博德研究所（Broad Institute）开发的用于变异检测包括SNP和INDEL的检测强大工具包，主要针对深度测序数据，如人类全基因组和外显子组的测序数据。GATK一直在不断更新，目前的版本是3.6。

GATK是基于MapReduce框架开发的Java程序，可轻易的实现并行和分布式处理。无论是DNA还是RNA的测序数据，GATK的分析流程都包括三个部分：数据预处理、变异检测、结果评估。


Genome Analysis Toolkit安装

操作系统

GATK可运行在基于Unix的操作系统，具有MacOSX和Linux版本，不支持Microsoft Windows。至少有500M硬盘容量和2G内存。这里是在Ubuntu系统下运行为例。

所需软件

1）Java 8或者JRE1.8以上版本；

2）IGV 基因组浏览器2.3.60（或最新的版本）：http://software.broadinstitute.org/software/igv/；

3）Picard 2.4.1或更新的版本，这里我们选用2.5.0版本：http://broadinstitute.github.io/picard/；

4）Samtools：http://samtools.sourceforge.net/；

5）RStudio和软件包ggplot2、gplots、bitops、caTools、colorspace、gdata、gsalib、reshape、RColorBrewer等。如果画图时出现错误，会提示需要安装的包的名称：https://www.rstudio.com/；

6）RTG Tools 3.6.2：http://realtimegenomics.com/products/rtg-tools/。

以上软件在各自的主页都有详细的安装方法，可以根据操作系统选择相应的版本安装这里不再赘述。

GATK的安装

1）在GATK的下载页：http://www.broadinstitute.org/gatk/download下载最新的版本；

2）解压下载的文件GenomeAnalysisTk­3.6­0.tar.bz2，得到一个名为GenomeAnalysisTK­3.6­0的文件夹，包括了GenomeAnalysisTK.jar以及示例文件。

$ tar xvjf GenomeAnalysisTk­3.6­0.tar.bz2

3）测试GATK是否安装成功

$ java ­jar GenomeAnalysisTK.jar ­h/p>

如果GATK正常工作，可得到如下图所示的内容。

图9-1 GATK正常工作截图


Genome Analysis Toolkit使用

这里将主要介绍GATK提供的两个变异检测工具——UnifiedGenotyper和HaplotypeCaller。

1. UnifiedGenotyper

UnifiedGenotyper使用贝叶斯最大似然模型，同时估计基因型和基因频率，最后对每一个样本的每一个变异位点和基因型都会给出一个后验概率。它可以用于单个样本的变异检测，也可以用于群体的变异检测。

命令语法：

java -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -R reference.fasta -I input.bam -o raw_variants.vcf -glm variants type [-L targets.interval_list]

其中相关参数的含义如下

-T：变异检测的工具

-R：参考序列，fasta格式

-I：BAM格式的输入数据文件

-o：变异检测的输出文件，vcf格式

-glm：变异检测的类型，SNP、INDEL和both（代表前面两种变异都检测）。默认为SNP

-L：可选参数，如指定时表示BAM文件中需要预测的区域

示例：

这里选取GATK官网教程中提供的数据集：data ，包括了命令所需的各种数据。检测NA12878_wgs_20.bam文件中20:10,000,000-10,200,000区域的SNP和INDEL变异，命令如下：

java -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -R data/ref/human_g1k_b37_20.fasta -I data/bams/exp_design/NA12878_wgs_20.bam -o data/sandbox/NA12878_wgs_20_UG_calls.vcf -glm BOTH -L 20:10,000,000-10,200,000

2. HaplotypeCaller

HaplotypeCaller是GATK3.x主力推荐的变异检测工具，在-ERC GVCF模式下单独对每个样的BAM文件进行变异检测，生成每个样本的gVCF文件；然后将合并这些gVCF文件进行genotyping，优化了流程，减少了消耗时间。

命令语法：

java -jar GenomeAnalysisTK.jar -R reference.fasta -T HaplotypeCaller -I sample1.bam -ERC GVCF [-D dbSNP.vcf] [-L targets.interval_list] -o output.raw.snps.indels.g.vcf

其中相关参数的含义如下：

-R、-T、-I、-L和-o的参数含义和设置与UnifiedGenotyper一致。

-ERC：参考置信度分数模式

-D：dbSNP文件

示例：检测多个文件的变异

这里同样使用数据集data。

对每个样本的BAM文件进行变异检测，这里共对三个bam文件，即NA12878_wgs_20.bam、NA12877_wgs_20.bam和NA12882_wgs_20.bam分别变异检测。仅以NA12878_wgs_20.bam为例，命令如下：

java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R data/ref/human_g1k_b37_20.fasta -I data/bams/exp_design/NA12878_wgs_20.bam -o data/sandbox/NA12878_wgs_20.g.vcf -ERC GVCF -L 20:10,000,000-10,200,000

对其他两个BAM文件同样运行如上代码，得到NA12877_wgs_20.g.vcf和NA12882_wgs_20.g.vcf

整合每个gVCF文件，进行joint genotyping，生成最后的VCF文件，代码如下：

java -jar GenomeAnalysisTK.jar -T GenotypeGVCFs -R data/ref/human_g1k_b37_20.fasta -V data/sandbox/NA12878_wgs_20.g.vcf -V data/gvcfs/NA12877_wgs_20.g.vcf -V data/gvcfs/NA12882_wgs_20.g.vcf -o data/sandbox/CEUTrio_wgs_20_GGVCFs_jointcalls.vcf -L 20:10,000,000-10,200,000

用IGV查看变异（可视化）。在终端用如下命令打开IGV

java -Xmx750m -jar ~ /IGV_2.3.80/igv.jar

加载变异检测的最终结果CEUTrio_wgs_20_GGVCFs_jointcalls.vcf，更改视图到20:10,002,584­10,002,665区域。