Mutect2 doesn't call SNP despite high coverage and strong evidence unless forced with --alleles

[PeterResutik]

Hi there,

I'm using Mutect2 to call variants on amplicon sequencing data of human mitochondrial DNA. In one specific amplicon with a very clean signal and a read depth of nearly 2000x, I observe two SNPs. However, only one of them is called by Mutect2, while the other is consistently missed.

The top panel shows the aligned reads (.bam file) after BWA-MEM mapping. The bottom panel shows the bamout.bam produced by Mutect2. You can clearly see the variants at positions 596 and 750 in both views.

Here’s a snippet of my pipeline:
gatk --java-options "-Xmx16G"
Mutect2
-R ${reference}
-L '${detected_contig}'
--min-base-quality-score ${params.baseQ}
--callable-depth 2
--linked-de-bruijn-graph true
--dont-use-soft-clipped-bases true
--initial-tumor-lod -10
--tumor-lod-to-emit -10
--max-reads-per-alignment-start 0
--bam-output ${bam_file.baseName}.bamout.bam
--tmp-dir tmp_${sample_id}
-I ${bam_file}
-O raw.vcf.gz \

When I force the variant at position 596 using the --alleles flag (--alleles $params.force_sites), the SNP is called with very strong support:
chrM 596 . T C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=1,0|1834,0;DP=1836;ECNT=3;ECNTH=1;GERMQ=93;MBQ=39,38;MFRL=0,0;MMQ=60,60;MPOS=26;POPAF=7.3;TLOD=7752.76 GT:AD:AF:DP:F1R2:F2R1:FAD:SB 0/1:1,1834:0.999:1835:0,0:1,1819:1,1834:1,0,1834,0

This is how force_sites.vcf.gz looks like:
##fileformat=VCFv4.2
#CHROM POS ID REF ALT QUAL FILTER INFO
chrM 596 . T C . . .

Is there a reason why this site wouldn't be called by default, even though the evidence seems pretty strong?

Thanks in advance for your help!

Best,
Peter

[phylyc]

From my experience, while the linked de Bruijn graph implementation theoretically gives higher sensitivity, it has trouble in complex or high depth regions. I find it pretty necessary to --recover-all-dangling-branches true. This should already emit some more calls.

However, even then, some SNPs will be missed. My solution is to pre-call with Mutect1 and force-call the passing variants with Mutect2.

[apraga]

Hi,
Thanks for the tip. This flag does not help catching our missing variants. I've tried to re-open this post https://gatk.broadinstitute.org/hc/en-us/community/posts/9037989302171-Mutect2-low-allele-fraction-mutation-calling.
It would be great to have official guidance for that.