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

Comments are closed.