GATK的FilterMutectCalls如何才能成功呢

因为有粉丝求助，他学习前面我分享的GATK的Mutect2流程都快奔溃了，总是各种报错。为了证明我教程没有错，所以我赶紧检查了代码，自己走了一遍，重新写了教程，了：最新最全的mutect2教程，提到了因为GATK的Mutect2流程更新太频繁，导致这个软件出现了一些无法解决的报错。走完了体细胞突变（somatic mutation）检测流程（Mutect2命令），这个时候拿到的文件仍然是需要过滤（走FilterMutectCalls命令）的，但是很多人就卡在了这一步。

比如我运行这个软件的FilterMutectCalls命令，测试了下面的几个情况：

reference=$HOME/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta  
GATK=$HOME/biosoft/GATK/gatk-4.1.8.1/gatk 
GATK=$HOME/biosoft/GATK/gatk-4.0.3.0/gatk
GATK=$HOME/biosoft/GATK/gatk-4.0.2.1/gatk
ls $reference $GATK 
ls  *_mutect2.vcf  |while read id
do
sample=$(basename "$id" _mutect2.vcf)
$GATK  FilterMutectCalls -R $reference  --java-options  -DGATK_STACKTRACE_ON_USER_EXCEPTION=true  
   -V  $id 
   -O ${sample}_filtered.vcf 
done

如果是gatk-4.1.8.1，就会报错如下：

A USER ERROR has occurred: Mutect stats table  _mutect2.vcf.stats not found.  
When Mutect2 outputs a file calls.vcf it also creates a calls.vcf.stats file.  
Perhaps this file was not moved along with the vcf, 
or perhaps it was not delocalized from a virtual machine while running in the cloud.

gatk官方论坛的意思是，在集群运行的过程中，会丢失后缀为.vcf.stats的文件，所以FilterMutectCalls 命令失败。

如果是Gatk-4.0.3.0，就会报错如下：

java.lang.IllegalStateException: Key P_CONTAM found in VariantContext field INFO at chr1:14932 but this key isn't defined in the VCFHeader.  We require all VCFs to have complete VCF headers by default.

但是，我记得我以前写这个软件教程的时候，明明没有出现问题啊，所以就去检查了我的脚本，发现居然是 gatk-4.0.2.1 版本。

如果是是 gatk-4.0.2.1 版本

报错就更诡异了，运行到一半后戛然而止。仔细检查了vcf文件停止的地方，发现它对

chr2	112391072	.	GAAA	G,GA,GAA
chr2	131598742	.	CT	C,CTT,CTTT

所以我干脆仅仅是保留SNP吧：

reference=$HOME/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta   
GATK=$HOME/biosoft/GATK/gatk-4.0.2.1/gatk
ls $reference $GATK 
ls  *_mutect2.vcf  |while read id
do
sample=$(basename "$id" _mutect2.vcf)
cat $id  | perl -alne '{if(/^#/){print}else{ next if $F[0] =~ "_";print if (length($F[3])+length($F[4])) eq 2   } }'  >  ${sample}_snp.vcf 
$GATK  FilterMutectCalls -R $reference  --java-options  -DGATK_STACKTRACE_ON_USER_EXCEPTION=true  
   -V ${sample}_snp.vcf   
   -O ${sample}_filtered.vcf 
cat  ${sample}_filtered.vcf |perl -alne '{if(/^#/){print}else{next unless $F[6] eq "PASS";next if $F[0] =~/_/;print } }' >  ${sample}_pass.vcf 
done

讽刺的是，居然就看到了成功的log日志：

18:10:54.132 INFO  ProgressMeter - Starting traversal
18:10:54.132 INFO  ProgressMeter -        Current Locus  Elapsed Minutes    Variants Processed  Variants/Minute
18:10:54.504 INFO  ProgressMeter -             unmapped              0.0                   792         128086.3
18:10:54.504 INFO  ProgressMeter - Traversal complete. Processed 792 total variants in 0.0 minutes.
18:10:54.657 INFO  FilterMutectCalls - Shutting down engine
[September 29, 2020 6:10:54 PM CST] org.broadinstitute.hellbender.tools.walkers.mutect.FilterMutectCalls done. Elapsed time: 0.04 minutes.
Runtime.totalMemory()=357564416
Tool returned:
SUCCESS

接下来这些后缀为_pass.vcf 的文件，就需要走vcf2maf流程啦！

vcf2maf流程我前些天在生信技能树已经分享过了，见：mskcc的vcf2maf极简解决方案代码分享