GATK best practice每个步骤都是必须的吗？

GATK4这个新的版本已经发布了，我前面写的一系列教程都是基于GATK3.X的，当然，整体数据处理流程其实并没有变化，但是时间消耗，步骤选择全部需要重新评价了。

一个全基因组重测序分析实战

GATK best practice每个步骤耗时如何？

最后这个步骤的必要性必须得马上发出来了，就是因为有些步骤也的确太耗时了，尤其是recal，不仅仅耗时还特别占硬盘存储！ 那么我们就从最后找到的变异文件的比较情况来说明这些步骤的必要性吧！ 我们再回顾一下GATK best practice流程吧！

流程介绍

bwa(MEM alignment)

picard(SortSam)

picard(MarkDuplicates)

picard(FixMateInfo)

GATK(RealignerTargetCreator)

GATK(IndelRealigner)

GATK(BaseRecalibrator)

GATK(PrintReads)

GATK(HaplotypeCaller)

GATK(GenotypeGVCFs)

我上一讲说到过，我对realign和recal以及原始的bam都用了HaplotypeCaller找变异，得到的vcf文件如下：

1. 1.1G Jun 21 02:29 jmzeng_merge_raw.snps.indels.vcf

2. 1.1G Jun 21 02:15 jmzeng_realigned_raw.snps.indels.vcf

3. 1.1G Jun 20 21:20 jmzeng_recal_raw.snps.indels.vcf

这是最原始的变异文件，我们需要进行过滤，拆分SNP和INDEL，这样才能更好的对它们两两比较。(这样就是比较查看realign和recal步骤对最后找变异的影响)

拆分SNP和INDEL并过滤低质量位点

代码如下：

1. module load java/1.8.0_91

2. GENOME=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta

3. GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar

4. TMPDIR=/home/jianmingzeng/tmp/software

6. sample=$1

8. ## for SNP

9. : '

10. '

11. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T SelectVariants  -R $GENOME  \

12. -selectType SNP \

13. -V ${sample}_raw.snps.indels.vcf -o ${sample}_raw_snps.vcf

16. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T VariantFiltration -R $GENOME  \

17. --filterExpression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"  \

18. --filterName "my_snp_filter" \

19. -V ${sample}_raw_snps.vcf  -o ${sample}_filtered_snps.vcf  

21. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T SelectVariants -R $GENOME  \

22. --excludeFiltered \

23. -V ${sample}_filtered_snps.vcf  -o  ${sample}_filtered_PASS.snps.vcf

26. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T VariantEval -R $GENOME  \

27. -eval ${sample}_filtered_PASS.snps.vcf -o  ${sample}_filtered_PASS.snps.vcf.eval

29. ## for  INDEL

30. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T SelectVariants  -R $GENOME \

31. -selectType INDEL  \

32. -V ${sample}_raw.snps.indels.vcf   -o ${sample}_raw_indels.vcf

35. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T VariantFiltration -R $GENOME  \

36. --filterExpression "QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0"  \

37. --filterName "my_indel_filter" \

38. -V ${sample}_raw_indels.vcf  -o ${sample}_filtered_indels.vcf  

40. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T SelectVariants -R $GENOME  \

41. --excludeFiltered \

42. -V ${sample}_filtered_indels.vcf  -o  ${sample}_filtered_PASS.indels.vcf

44. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T VariantEval -R $GENOME  \

45. -eval ${sample}_filtered_PASS.indels.vcf -o  ${sample}_filtered_PASS.indels.vcf.eval

要深度理解这个代码请点击自行阅读文档。 其实我本人不大喜欢用这个工具的各种参数，我比较喜欢自行写脚本来过滤vcf文件。

这样得到的文件如下：

1.     801242 jmzeng_merge_filtered_indels.vcf

2.     760890 jmzeng_merge_filtered_PASS.indels.vcf

3.    3812094 jmzeng_merge_filtered_PASS.snps.vcf

4.    4034168 jmzeng_merge_filtered_snps.vcf

5.     801240 jmzeng_merge_raw_indels.vcf

6.    4837892 jmzeng_merge_raw.snps.indels.vcf

7.    4034166 jmzeng_merge_raw_snps.vcf

8.     801963 jmzeng_realigned_filtered_indels.vcf

9.     761510 jmzeng_realigned_filtered_PASS.indels.vcf

10.    3812442 jmzeng_realigned_filtered_PASS.snps.vcf

11.    4034797 jmzeng_realigned_filtered_snps.vcf

12.     801961 jmzeng_realigned_raw_indels.vcf

13.    4839256 jmzeng_realigned_raw.snps.indels.vcf

14.    4034795 jmzeng_realigned_raw_snps.vcf

15.     793611 jmzeng_recal_filtered_indels.vcf

16.     754755 jmzeng_recal_filtered_PASS.indels.vcf

17.    3784343 jmzeng_recal_filtered_PASS.snps.vcf

18.    4010670 jmzeng_recal_filtered_snps.vcf

19.     793609 jmzeng_recal_raw_indels.vcf

20.    4806696 jmzeng_recal_raw.snps.indels.vcf

21.    4010668 jmzeng_recal_raw_snps.vcf

很容易理解：

• 对recal的bam得到的变异vcf文件来说，总共是480万位点，拆分成401万的SNP和79万的INDEL，然后经过过滤后剩下378万的SNP和75万的INDEL。

• 对原始的bam得到的变异vcf文件来说， 总共是483万位点，拆分成403万的SNP和80万的INDEL，然后经过过滤后剩下381万的SNP和76万的INDEL。

• 对重排的bam得到的变异vcf文件来说， 总共是483万位点，拆分成403万的SNP和80万的INDEL，然后经过过滤后剩下381万的SNP和76万的INDEL。

从位点数量级来说，原始的bam和重排的bam得到的变异vcf文件没区别，所以需要用专业的工具来具体比较它们里面的每一个位点。 我这里还是选择SnpEff套件里面的SnpSift工具。

首先比较SNP位点

代码如下：

1. java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar  concordance -v ../jmzeng_merge_filtered_PASS.snps.vcf ../jmzeng_realigned_filtered_PASS.snps.vcf 1>concordance.txt 2>SnpSift_Concordance.log

2. java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar  concordance -v ../jmzeng_recal_filtered_PASS.snps.vcf  ../jmzeng_realigned_filtered_PASS.snps.vcf 1>concordance.txt 2>SnpSift_Concordance.log

3. java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar  concordance -v ../jmzeng_recal_filtered_PASS.snps.vcf  ../jmzeng_merge_filtered_PASS.snps.vcf 1>concordance.txt 2>SnpSift_Concordance.log

比较的结果如下：

1. Number of samples:

2.    1   File ../jmzeng_merge_filtered_PASS.snps.vcf

3.    1   File ../jmzeng_realigned_filtered_PASS.snps.vcf

4.    1   Both files

6. # Errors:

7.    ALT field does not match    31

9. Number of samples:

10.    1   File ../jmzeng_recal_filtered_PASS.snps.vcf

11.    1   File ../jmzeng_realigned_filtered_PASS.snps.vcf

12.    1   Both files

14. # Errors:

15.    ALT field does not match    204

17. Number of samples:

18.    1   File ../jmzeng_recal_filtered_PASS.snps.vcf

19.    1   File ../jmzeng_merge_filtered_PASS.snps.vcf

20.    1   Both files

22. # Errors:

23.    ALT field does not match    208

可以看到对高质量的SNP位点来说，3种bam文件得到SNP信息差别很微弱，在可接受范围内。但是我们不能忽视原始的bam和重排的bam得到的变异vcf文件要比recal后的少了近两万这个事实！！！

接着比较INDEL位点

代码如下：

1. java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar  concordance -v ../jmzeng_merge_filtered_PASS.indels.vcf ../jmzeng_realigned_filtered_PASS.indels.vcf 1>concordance.txt 2>SnpSift_Concordance.log

2. java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar  concordance -v ../jmzeng_recal_filtered_PASS.indels.vcf  ../jmzeng_realigned_filtered_PASS.indels.vcf 1>concordance.txt 2>SnpSift_Concordance.log

3. java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar  concordance -v ../jmzeng_recal_filtered_PASS.indels.vcf  ../jmzeng_merge_filtered_PASS.indels.vcf 1>concordance.txt 2>SnpSift_Concordance.log

比较的结果如下：

1. Number of samples:

2.    1   File ../jmzeng_merge_filtered_PASS.indels.vcf

3.    1   File ../jmzeng_realigned_filtered_PASS.indels.vcf

4.    1   Both files

6. # Errors:

7.    REF fields does not match   28

8.    ALT field does not match    45

11. Number of samples:

12.    1   File ../jmzeng_recal_filtered_PASS.indels.vcf

13.    1   File ../jmzeng_merge_filtered_PASS.indels.vcf

14.    1   Both files

16. # Errors:

17.    REF fields does not match   1295

18.    ALT field does not match    964

21. Number of samples:

22.    1   File ../jmzeng_recal_filtered_PASS.indels.vcf

23.    1   File ../jmzeng_realigned_filtered_PASS.indels.vcf

24.    1   Both files

26. # Errors:

27.    REF fields does not match   1276

28.    ALT field does not match    947

INDEL本身对参数非常敏感，这时候不仅仅是数量的差异，而且本身找到的位点突变情况的差异也不少。

所以我的结论是，GATK的BEST PRACTISE的每个步骤都是必须的！虽然他们非常耗时，但是对结果准确性的影响的确非常显著。 如果要想把测序数据在临床上面应用，那么每一个步骤都是非常有意义的。

比如，如果我们来分析realign之后的高质量SNP为什么会有31个的ALT改变了呢？

1. 21    10716592

2. 21    10701512

3. 21    10700605

4. 20    26317761

5. 19    15787221

6. 18    18515822

7. 17    25266293

8. 16    35230302

9. 16    33975649

10. 16    33972478

11. 10    42400348

12. 10    42393437

13. 9    66455306

14. 4    49111623

15. 1    142537167

16. Y    58977742

17. Y    13478115

18. X    61909282

19. X    61730552

20. X    61721098

简单看了几眼， 发现都是在染色体的中心粒的地方！

虽然GATK best practice中的BSQR步骤的确非常耗时(WGS数据约40小时)，但是对最后的结果影响还是蛮大的，所以是否需要省略这个步骤取决于你自己对找变异的精度要求。毕竟这个步骤说明书上面可说的是用了机器学习的方法来矫正测序误差呀！

最后，上两张GATK4的PPT吧~

又要开始学新的软件，写新的教程啦~~~

GATK4性能提升如此显著，不得不换！

尤其是在生物信息学这个日新月异的领域。