GATK使用(一)

  首先说说GATK可以做什么。它主要用于从sequencing 数据中进行variant calling，包括SNP、INDEL。比如现在风行的exome sequencing找variant，一般通过BWA+GATK的pipeline进行数据分析。

  要run GATK，首先得了解它的网站(http://www./gatk/)。

一.准备工作

(1)程序下载

http://www./gatk/download

务必下载最新版，注意不是那个GATK-lite，而是GATK2.（目前是2，不知何时会更新到3.更新很快啊有木有，所以一定要用新版）解压到某个目录即可，无需安装。打开解压缩后的目录，会发现GenomeAnalysisTK.jar这个文件。我们要用的各种工具都在这个包里。

(2)Reference sequence和annotation下载

务必使用GATK resource bundle上各种序列和注释来run GATK。 

GATK resource bundle介绍：

http://gatkforums./discussion/1213/whats-in-the-resource-bundle-and-how-can-i-get-it

GATK resource bundle FTP地址：

http://gatkforums./discussion/1215/how-can-i-access-the-gsa-public-ftp-server

例如呢，以hg19(human genome build 19.此外，b37也很常用)为参考序列,

输入如下命令： lftp ftp. -u gsapubftp-anonymous

回车后键入空格，便可进入resource bundle。进入其中名为bundle的目录，找到最新版的hg19目录。

要下载的文件包括：

ucsc.hg19.fasta

ucsc.hg19.fasta.fai

1000G_omni2.5.hg19.vcf

1000G_omni2.5.hg19.vcf.idx

1000G_phase1.indels.hg19.vcf

1000G_phase1.indels.hg19.vcf.idx

dbsnp_137.hg19.vcf

dbsnp_137.hg19.vcf.idx

hapmap_3.3.hg19.vcf

hapmap_3.3.hg19.vcf.idx

Mills_and_1000G_gold_standard.indels.hg19.vcf

Mills_and_1000G_gold_standard.indels.hg19.vcf.idx

(3)R 设置

GATK某些步骤需要使用R画一些图，此时会调用一些R packages。如果你使用的R没有装这些packages，就无法生成pdf文件。因此得安装一些必要的R packages，包括（但不限于）如下：

• bitops
• caTools
• colorspace
• gdata
• ggplot2
• gplots
• gtools
• RColorBrewer
• reshape
• gsalib

如何check这些packages有没有安装？以ggplot2为例，进入R，输入library(ggplot2)，如果没有error信息，就表明安装了ggplot2；如果没有安装，输入命令install.packages("ggplot2")，按照提示操作便可安装ggplot2。如果package安装不成功，很可能是权限问题：默认使用的R是系统管理员安装的，而你在/usr等目录下没有write权限。这是可以选择直接在自己目录下安装了R，然后install各种packages。也可以在 install.packages时指定将packages安装到自己的目录。

将R packages安装到自己目录的方法：

 install.packages("ggplot2", lib="/your/own/R-packages/")

然后将指定的路径添加到R_LIBS变量即可。

gsalib的安装有点不一样。

从该网站下载gsalib：https://github.com/broadgsa/gatk

解压缩后，进入build.xml文件所在的目录，输入命令ant gsalib。这样就可以安装gsalib 了。如果该gsalib仍然无法load，可能需要将gsllib的path告诉R。

http://gatkforums./discussion/1244/what-is-a-gatkreport

这里有较详细的步骤。（奇怪的是我没有~/.Rprofle文件，ant完后也没有add path，仍然work了）

此外，gaslib安装要求java的版本为1.6。

设置路径：

export R_HOME=/your/R/install/path

export R_LIBS=/your/R-packages/path

如何查看当前使用的R的路径？输入 which R即可。

二.流程概览

图片摘自http://www./gatk/guide/topic?name=best-practices。话说这个Best practices，刚开始接触GATK时每个人都推荐我看，可是我就是看不懂...许多东西要亲自做过后，再回头看才能懂。不过既然那么元老级人物都说要看看，所以如果你是新手，还是看看吧。

好消息是现在该页面多了个BroadE Workshop的东东，里面的一些材料非常有用。http://www./gatk/guide/events?id=2038#materials

在工作之前，最好把这些work materials通览一遍。

seqanswers上的这个wiki也是很好的入门材料，可以让你看看每一步在做什么。

http:///wiki/How-to/exome_analysis

但是呢，它比较过时了，它使用的GATK版本应该是2.0以前的，现在的流程有了一些调整，一些参数设置也有变化。所以该pipeline仅供参考，切勿照抄代码。例如它关于reference sequence的那一部分，是自己手动generate hg19的genome sequence。但是，由于后面要用到dbSNP和indel的各种annotation，而这些file的index很可能和自己generate的hg19的index不一样(resource bundle上hg19除了chr1-22，X,Y，chrC，ChrM外还会多出一些没有组装的序列)，这样会导致程序出错，到时候还得重头来过。所以，务必使用GATK resource bundle上的东西！

我后面贴出的代码也一样，切勿照抄。鬼知道过几个月GATK又会有什么变化。It keeps changing all the time...时刻跟上GATK网站上的update才是王道。（虽然它的网站做得...实在算不上好...好多时候我都找不到自己想要的东西）

三.流程

(1)Mapping。

因我只处理过Miseq和Hiseq的pair-end数据，以此为例。Mapping多使用BWA，因为它能允许indels。另一常用mapping软件Bowtie则不能generate indels，多用于perfect match的场合如sRNA mapping。

首先建立reference的index。

bwa index -a bwtsw -p hg19 hg19.fasta

这个过程耗时良久...above 5h吧。

-p是给reference起的名字。完成后会生成如下几个文件：

 hg19.amb

 hg19.ann

 hg19.bwt

 hg19.dict

 hg19.pac

 hg19.sa

然后是mapping喽，开始贴代码...先对pair-end数据的pair1和pair2分别进行mapping，生成相应的sai file，然后再结合两者生成sam file。接着利用picard对sam file进行sort，并且生成sort后的bam file。

sample=$1 data_dir=*** work_dir=***  ref_dir=~/database/hg19 sai_dir=$work_dir/saifiles sam_dir=$work_dir/samfiles bam_dir=$work_dir/bamfiles log_dir=$work_dir/logs met_dir=$work_dir/metrics other_dir=$work_dir/otherfiles gatk=/software/sequencing/GenomeAnalysisTK-2.3-4-g57ea19f/GenomeAnalysisTK.jar picard=/software/sequencing/picard_latest bwa aln $ref_dir/hg19 -t 4 $data_dir/${sample}_1.fastq.gz > $sai_dir/${sample}_1.sai bwa aln $ref_dir/hg19 -t 4 $data_dir/${sample}_2.fastq.gz > $sai_dir/${sample}_2.sai bwa sampe -r "@RG\tID:$sample\tLB:$sample\tSM:$sample\tPL:ILLUMINA"   $ref_dir/hg19 $sai_dir/${sample}_1.sai $sai_dir/${sample}_2.sai $data_dir/${sample}_1.fastq.gz $data_dir/${sample}_2.fastq.gz | gzip > $sam_dir/$sample.sam.gz java -Xmx10g -Djava.io.tmpdir=/tmp \  -jar $picard/SortSam.jar \  INPUT=$sam_dir/$sample.sam.gz \  OUTPUT=$bam_dir/$sample.bam \  SORT_ORDER=coordinate\  VALIDATION_STRINGENCY=LENIENT \  CREATE_INDEX=true

Tips：

a. bwa以及后面使用的一系列tools都能直接读取gzip压缩文件。为了节省硬盘空间，以及减少压缩/解压缩过程的read/write时间，最好直接用*.gz文件作为input。output也最好先在管道里使用gzip压缩一下再输出。

b.bwa sample一步里的-r选项不可忽略。

c.Picard的SortSam需指定一个tmp目录，用于存放中间文件，中间文件会很大，above 10G.注意指定目录的空间大小。

d.对于压缩后大小约3G的fastq.gz文件，bwa aln需run 30-60min。将生成的两个sai file merge到一起生成sam file需60-90min。利用picard sort sam需30-60min。

e.Samtools也提供了工具进行sam file的sort和bam file的生成，并且不会生成大的中间文件。

(2) Duplicate marking。

这一步的目的是mark那些PCR的duplicate。由于PCR有biases，有的序列在会被overpresented, 它们有着相同的序列和一致的map位置，对GATK的work会造成影响。这一步给这些序列设置一个flag以标志它们，方便GATK的识别。还可以设置 REMOVE_DUPLICATES=true 来丢弃duplicated序列，不过还未探索过只标记不丢弃和丢弃对于后续分析的影响。官方流程里只用标记就好。

java -Xmx15g -Djava.io.tmpdir=/tmp \  -jar $picard/MarkDuplicates.jar \  INPUT= $bam_dir/$sample.bam \ OUTPUT= $bam_dir/$sample.dedupped.bam \  METRICS_FILE= $met_dir/$sample.dedupped.metrics \  CREATE_INDEX=true \ VALIDATION_STRINGENCY=LENIENT

(3)Local relignment

INDEL附近的alignment通常不准，需要利用已知的indel信息进行realignment，分两步，第一步进行 RealignerTargetCreator，输出一个包含着possible indels的文件。第一步IndelRealigner，利用此文件进行realign。

第一步需要用到一些已知的indel信息。那么该使用哪些文件呢？

参考如下的链接：http://gatkforums./discussion/1247/what-should-i-use-as-known-variantssites-for-running-tool-x 

其中这张图表尤其有用，可以看到RealignerTargetCreator和IndelRealigner步骤中推荐使用两个indel注释file，如代码中-known所示。

Tool	dbSNP 129 -	- dbSNP >132 -	- Mills indels -	- 1KG indels -	- HapMap -	- Omni
RealignerTargetCreator			X	X
IndelRealigner			X	X
BaseRecalibrator		X	X	X
(UnifiedGenotyper/ HaplotypeCaller)		X
VariantRecalibrator		X	X		X	X
VariantEval	X

java -Xmx15g -jar $gatk \ -T RealignerTargetCreator \ -R $ref_dir/hg19.fasta \ -o $other_dir/$sample.realigner.intervals \ -I $bam_dir/$sample.dedupped.bam \ -known $ref_dir/1000G_phase1.indels.hg19.vcf \ -known $ref_dir/Mills_and_1000G_gold_standard.indels.hg19.vcf java -Xmx15g -jar $gatk \ -T IndelRealigner \  -R $ref_dir/hg19.fasta \  -I $bam_dir/$sample.dedupped.bam \ -targetIntervals $other_dir/$sample.realigner.intervals \  -o $bam_dir/$sample.dedupped.realigned.bam \ -known $ref_dir/1000G_phase1.indels.hg19.vcf \ -known $ref_dir/Mills_and_1000G_gold_standard.indels.hg19.vcf

在一些老版本的pipeline中(如seqanswers上的那个)，这一步做完后还要对pair-end数据的mate information进行fix。不过新版本已经不需要这一步了，Indel realign可以自己fix mate。参考http://gatkforums./discussion/1562/need-to-run-a-step-with-fixmateinformation-after-realignment-step

(4) Base quality score recalibration

这一步简称BQSR，对base的quality score进行校正。同样要用到一些SNP和indel的已知文件(参考上表和代码中-known选项)。第一步BaseRecalibrator生成一个用于校正的recal.grp文件和一个校正前的plot；第一步BaseRecalibrator生成一个校正后的plot；第三步PringReads利用recal.grp进行校正，输出校正后的bam file。

java -Xmx15g -jar $gatk \  -T BaseRecalibrator \  -R $ref_dir/hg19.fasta \  -I $bam_dir/$sample.dedupped.realigned.bam \  -knownSites $ref_dir/dbsnp_137.hg19.vcf \  -knownSites $ref_dir/Mills_and_1000G_gold_standard.indels.hg19.vcf \  -knownSites $ref_dir/1000G_phase1.indels.hg19.vcf \  -o $other_dir/$sample.recal.grp \  -plots $other_dir/$sample.recal.grp.pdf java -Xmx15g -jar $gatk \  -T BaseRecalibrator \  -R $ref_dir/hg19.fasta \  -I $bam_dir/$sample.dedupped.realigned.bam \  -BQSR $other_dir/$sample.recal.grp \  -o $other_dir/$sample.post_recal.grp \  -plots $other_dir/$sample.post_recal.grp.pdf \  -knownSites $ref_dir/dbsnp_137.hg19.vcf \  -knownSites $ref_dir/Mills_and_1000G_gold_standard.indels.hg19.vcf \  -knownSites $ref_dir/1000G_phase1.indels.hg19.vcf java -Xmx15g -jar $gatk \  -T PrintReads \  -R $ref_dir/hg19.fasta \  -I $bam_dir/$sample.dedupped.realigned.bam \  -BQSR $other_dir/$sample.recal.grp \  -o $bam_dir/$sample.dedupped.realigned.recal.bam

TIPS：

这一步里如果产生不了PDF文件，很可能是Rscript的执行出了问题。检查R里面是否安装了该安装的packages。如果仍无法找出原因，可以在BaseRecalibrator那两步加上-l DEBUG再运行，此时的log信息中会详细写明为什么Rscript执行不了。