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为例,生成的命令为: