java queue GATK_GATK使用注意事项

GATK这个软件在做snp-calling的时候使用率非常高，因为之前一直是简单粗略的看看snp情况而已，所以没有具体研究它。

这些天做一些外显子项目以找snp为重点，所以想了想还是用起它，报错非常多，调试了好久才成功。

所以记录一些注意事项!

GATK软件本身是受版权保护的，所以需要申请才能下载使用，大家自己去broad institute申请即可。

下载软件就可以直接使用，java软件不需要安装，但是需要你的机器上面有java，当然软件只是个开始，重点是你还得下载很多配套数据，https://software.broadinstitute.org/gatk/download/bundle(ps:这个链接可能会失效，下面的文件，请自己谷歌找到地址哈。)，而且这个时候要明确你的参考基因组版本了！！！ b36/b37/hg18/hg19/hg38，记住b37和hg19并不是完全一样的，有些微区别哦！！！

比如我选择了hg19

第一点是hg19的下载：这个下载地址非常多，常用的就是NCBI，ensembl和UCSC了，但是这里推荐用这个脚本下载

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

do echo $i;

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;

done

gunzip *.gz

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

do cat chr${i}.fa >> hg19.fasta;

done

rm -fr chr*.fasta

看得懂shell脚本的应该知道这是一个个的下载hg19的染色体，再用cat按照染色体的顺序拼接起来，因为GATK后面的一些步骤对染色体顺序要求非常变态，如果下载整个hg19，很难保证染色体顺序是1-22，X,Y,M。如下

然后需要对下载的hg19进行索引(bwa和samtools)和建立dict文件(用picard)

bwa index -a bwtsw hg19.fasta

samtools faidx hg19.fasta

然后还要下载几个参考文件，这个是可以选择的.

对我的hg19来说，就应该是去，ftp://ftp.broadinstitute.org/bundle/hg19/ 下载咯。

最后，所有必备的文件如下：

231M Jul 2 05:14 1000G_phase1.indels.hg19.sites.vcf

1.2M Jul 2 10:45 1000G_phase1.indels.hg19.sites.vcf.idx

11G Jul 2 08:05 dbsnp_138.hg19.vcf

2.5K Jul 1 04:31 hg19.dict

3.0G Jun 30 21:29 hg19.fasta

6.6K Jun 30 22:54 hg19.fasta.amb

944 Jun 30 22:54 hg19.fasta.ann

2.9G Jun 30 22:54 hg19.fasta.bwt

788 Jul 2 01:53 hg19.fasta.fai

739M Jun 30 22:54 hg19.fasta.pac

1.5G Jun 30 23:23 hg19.fasta.sa

87M Jul 2 05:37 Mills_and_1000G_gold_standard.indels.hg19.sites.vcf

2.3M Jul 2 10:45 Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.idx

接下来开始跑程序

第一步就是生成sam文件啦bwa mem -t 12 -M  hg19.fasta tmp*fq >tmp.sam

第二步是sort，我用的是picard工具java  -Xmx100g -jar AddOrReplaceReadGroups.jar I=tmp.sam  O=tmp.sorted.bam

SORT_ORDER=coordinate

CREATE_INDEX=true

RGID=tmp

RGLB="pe"

RGPU="HiSeq-2000"

RGSM=PC3-2

RGCN="Human Genetics of Infectious Disease"

RGDS=hg19 RGPL=illumina

VALIDATION_STRINGENCY=SILENT

第三步是去除PCR重复，我还是选择用picard工具

java  -Xmx100g  -jar MarkDuplicates.jar

CREATE_INDEX=true REMOVE_DUPLICATES=True

ASSUME_SORTED=True VALIDATION_STRINGENCY=LENIENT

I=tmp.sorted.bam OUTPUT=tmp.dedup.bam METRICS_FILE=tmp.metrics

第四步是终于要开始用GATK啦，主要是确定要进行重新比对的区域，这个步骤分成三个小步骤：

首先用RealignerTargetCreator找到需要重新比对的区域，输出文件intervals

java -Xmx200g -jar ~/apps/gatk/GenomeAnalysisTK.jar

-R hg19.fasta  #这里需要用这个参考基因组，所以参考基因组特别重要，DICT也要按照流程生成

-T RealignerTargetCreator

-I tmp.dedup.bam -o tmp.intervals

-known /home/ldzeng/EXON/ref/1000G_phase1.indels.hg19.sites.vcf

这一步骤好像非常耗时

可以看到，我总共就测试了5014个reads，结果就花了近半个小时才搞定，只有947个reads被过滤了。

输出的tmp.intervals 文件是一个1404946行的文件

chr1:13957-13958

chr1:46402-46403

chr1:47190-47191

chr1:52185-52188

chr1:53234-53236

chr1:55249-55250

chr1:63735-63738

人的外显子只有二三十万，所以我暂时也不确定这个文件是什么！

然后用输出的 tmp.intervals 做输入文件来进行重新比对，也就是用IndelRealigner在这些区域内进行重新比对

java -Xmx150g -jar ~/apps/gatk/GenomeAnalysisTK.jar \

-R hg19.fasta \

-T IndelRealigner \

-targetIntervals tmp.intervals \

-I tmp.dedup.bam -o tmp.dedup.realgn.bam \

-known /home/ldzeng/EXON/ref/1000G_phase1.indels.hg19.sites.vcf

我只需要它的重新比对，所以后面的一些功能没有怎么用，一个是call snp，一个是算比对质量值

java -Xmx200g -jar ~apps/gatk/GenomeAnalysisTK.jar

-nct 20 -T HaplotypeCaller -R hg19.fasta

-I tmp.dedup.realgn.bam

-o tmp.gatk.vcf

最后输出的文件如下

639K Jul 5 10:17 tmp1.fq

639K Jul 5 10:19 tmp2.fq

1.5M Jul 5 10:26 tmp.dedup.bai

403K Jul 5 10:26 tmp.dedup.bam

12K Jul 5 12:02 tmp.gatk.vcf

3.4K Jul 5 12:02 tmp.gatk.vcf.idx

32M Jul 5 11:24 tmp.intervals

950 Jul 5 10:26 tmp.metrics

1.5M Jul 5 11:31 tmp.realgn.bai

409K Jul 5 11:31 tmp.realgn.bam

1.6M Jul 5 10:20 tmp.sam

1.5M Jul 5 10:23 tmp.sorted.bai

399K Jul 5 10:23 tmp.sorted.bam

备注：GATK对基因组要求一个字典文件

使用picard工具包的CreateSequenceDictionary.jar生成。以hg19.fa为例，生成的命令为：

java -Xmx2g -jar /path_to_picard/CreateSequenceDictionary.jar R=hg19.fa O=hg19.dict