Some of variants detected by the Mutect2 are contradictory
Hi, I used GATK Mutect2 to detect the mitochondrial variant of the amplicon data, but I found that some adjacent Indel/SNP results were contradictory, as follows:version: Version:GATK4.1.8.1
GATK commands:
java -jar gatk-package-4.1.8.1-local.jar Mutect2 -R hs37d5.fa -L MT --max-reads-per-alignment-start 0 --mitochondria-mode -I test.sort.bam -O test.vcf --max-mnp-distance 0
java -jar gatk-package-4.1.8.1-local.jar FilterMutectCalls --mitochondria-mode -R hs37d5.fa -V test.vcf -O test.filter.vcf
Part of the output test.filter.vcf is:
MT 547 . A G . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,4058|0,13;DP=4072;ECNT=22;MBQ=37,38;MFRL=0,0;MMQ=60,60;MPOS=37;OCM=0;POPAF=2.40;TLOD=14.20 GT:AD:AF:DP:F1R2:F2R1:SB
0/1:4058,13:3.388e-03:4071:4045,13:0,0:0,4058,0,13
MT 555 . A G . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,4056|0,15;DP=4072;ECNT=22;MBQ=37,38;MFRL=0,0;MMQ=60,60;MPOS=45;OCM=0;POPAF=2.40;TLOD=11.08 GT:AD:AF:DP:F1R2:F2R1:SB
0/1:4056,15:3.095e-03:4071:4031,13:0,0:0,4056,0,15
MT 556 . A G . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,4047|0,18;DP=4072;ECNT=22;MBQ=37,37;MFRL=0,0;MMQ=60,60;MPOS=46;OCM=0;POPAF=2.40;TLOD=17.57 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:4047,18:4.034e-03:4065:4003,17:0,0:0,4047,0,18
MT 567 . A G,AC . PASS AS_FilterStatus=SITE|SITE;AS_SB_TABLE=0,4004|0,21|0,46;DP=4072;ECNT=22;MBQ=37,38,37;MFRL=0,0,0;MMQ=60,60,60;MPOS=48,50;OCM=0;POPAF=2.40,2.40;TLOD=23.85,12.25 GT:AD:AF:DP:F1R2:F2R1:SB 0/1/2:4004,21,46:5.073e-03,7.694e-03:4071:3950,21,46:0,0,0:0,4004,0,67
MT 572 . CCA C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,4032|0,34;DP=4072;ECNT=22;MBQ=37,37;MFRL=0,0;MMQ=60,60;MPOS=42;OCM=0;POPAF=2.40;RPA=2,1;RU=CA;STR;TLOD=10.35 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:4032,34:5.733e-03:4066:3894,34:0,0:0,4032,0,34
MT 573 . CA C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,437|0,3634;DP=4072;ECNT=22;MBQ=37,37;MFRL=0,0;MMQ=60,60;MPOS=42;OCM=0;POPAF=2.40;TLOD=10425.69 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:437,3634:0.894:4071:426,3448:0,0:0,437,0,3634
MT 574 . A C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,57|0,3962;DP=4072;ECNT=22;MBQ=37,37;MFRL=0,0;MMQ=60,60;MPOS=41;OCM=0;POPAF=2.40;TLOD=610.85 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:57,3962:0.986:4019:48,364:0,0:0,57,0,3962
MT 593 . T C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,4044|0,19;DP=4072;ECNT=22;MBQ=37,34;MFRL=0,0;MMQ=60,60;MPOS=23;OCM=0;POPAF=2.40;TLOD=8.20 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:4044,19:3.165e-03:4063:3917,13:0,0:0,4044,0,19
MT 596 . T C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,4052|0,15;DP=4072;ECNT=22;MBQ=37,37;MFRL=0,0;MMQ=60,60;MPOS=20;OCM=0;POPAF=2.40;TLOD=9.06 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:4052,15:2.886e-03:4067:3934,11:0,0:0,4052,0,15
MT 608 . A G . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,4057|0,14;DP=4072;ECNT=22;MBQ=37,38;MFRL=0,0;MMQ=60,60;MPOS=8;OCM=0;POPAF=2.40;TLOD=16.31 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:4057,14:3.637e-03:4071:4034,14:0,0:0,4057,0,14
Obviously, on the one hand, MT:573 is del (CA->C), which means that base A of MT:574 is del.But MT:574 is a SNP site (A->C) in the vcf. This result is obviously contradictory. Why does this happen? Is there a problem with my parameters? Or is the software vulnerable?
-
Hello rime,
My first suggestion is probably unrelated to this issue but we always recommend you submit GATK commands with the wrapper script. You can find more information about the best way to run the GATK command line at this link: https://gatk.broadinstitute.org/hc/en-us/articles/360035531892-GATK4-command-line-syntax
We have an article outlining how to troubleshoot HaplotypeCaller and Mutect2 results. In your case, I would definitely recommend looking at the bam out file in IGV to determine why these variants are being called. Here is that article: https://gatk.broadinstitute.org/hc/en-us/articles/360043491652-When-HaplotypeCaller-and-Mutect2-do-not-call-an-expected-variant
Hope this helps to figure out what is going on!
Best,
Genevieve
-
Thanks very much. I've checked the bamout file by add --bam-output argument. Please look at the following IGV screenshots. Obviously, both variants do exist (MT:574 del and MT:574 A->C), most of the reads support MT:574 DEL, and the rest of reads support SNP mutation. But another problem is that the DP of MT:574 not match the original bamout in IGV: the DP of MT:573 and MT:574 are more than 4000 in the VCF, while the DP of MT:574 in IGV is only up to 532. In the article of troubleshoot, it is mentioned that realign may lead to depth difference. Is this the reason for my situation? Is there any way to avoid it?
“Reads in GATK are realigned to their best haplotype before genotyping. This can cause reads that appear to overlap a site to move (which either lowers or raises the DP).”
MT 573 . CA C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,437|0,3634;DP=4072;ECNT=22;MBQ=37,37;MFRL=0,0;MMQ=60,60;MPOS=42;OCM=0;POPAF=2.40;TLOD=10425.69 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:437,3634:0.894:4071:426,3448:0,0:0,437,0,3634
MT 574 . A C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=0,57|0,3962;DP=4072;ECNT=22;MBQ=37,37;MFRL=0,0;MMQ=60,60;MPOS=41;OCM=0;POPAF=2.40;TLOD=610.85 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:57,3962:0.986:4019:48,364:0,0:0,57,0,3962 -
Hi rime,
We have some documentation explaining possible reasons for these differences. Please see these resources:
- https://gatk.broadinstitute.org/hc/en-us/articles/360035532252-Allele-Depth-AD-is-lower-than-expected
- https://gatk.broadinstitute.org/hc/en-us/community/posts/360057830232-Wrong-Calculation-of-DP-AD-AF
Genevieve
Please sign in to leave a comment.
3 comments