Mutect2 - AD does not match AF neither IGV

Follow

Fran

• March 25, 2025 09:03

Hi all,

I have a question regarding the AD and AF fields reported in Mutect2 (tumor-only mode). The reported AD does not match what I see in IGV neither the reported AF for a given variant.

• reported AD = 1042,1056
• reported AF = 0.018

what I see in IGV: AF aprox 1%

 

if AD is formatted as"REF,ALT", I would expect "2136,14" or similar numbers.

The mentioned variant is:
chr17    43082434    .    G    A    .    .    AS_SB_TABLE=5,1037|1048,8;DP=2129;ECNT=15;MBQ=20,10;MFRL=140,140;MMQ=60,60;MPOS=27;POPAF=7.30;TLOD=39.03    GT:AD:AF:DP:F1R2:F2R1:FAD:SB0/1:1042,1056:0.018:2098:0,0:10,8:802,31:5,1037,1048,8

I find this weird since it is a big difference between 1% and 50% and 14 and 1056 reads.

 

REQUIRED for all errors and issues:
a) GATK version used:

v4.4.0.0

b) Exact command used:

/nfs/home/software/packages/gatk-4.4.0.0/gatk --java-options -Xmx40g Mutect2 --dont-use-soft-clipped-bases true --genotype-germline-sites --min-base-quality-score 20 --max-reads-per-alignment-start 0 --native-pair-hmm-threads 16 -R /nfs/home/references/genomes/human/GRCh38/bwa_mem2/Homo_sapiens.GRCh38.fa -I /media/scratch1/20250314_UB485_runNS252_Oncocarta//analysis_BRC2_100x_sub50/mapping/BRC2-ONC1-1316-S-sub50.bam -L /media/scratch1/20250314_UB485_runNS252_Oncocarta//analysis_BRC2_100x_sub50/QC/mapping/BRC2-ONC1-1316-S-sub50_covered_regions_100x.bed -O /media/scratch1/20250314_UB485_runNS252_Oncocarta//tmp_BRC2_100x_sub50/variant_calling/BRC2-ONC1-1316-S-sub50_mutect2_tmp1.vcf

c) Entire program log: 

Using GATK jar /nfs/home/software/packages/gatk-4.4.0.0/gatk-package-4.4.0.0-local.jar

Running:

java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -Xmx40g -jar /nfs/home/software/packages/gatk-4.4.0.0/gatk-package-4.4.0.0-local.jar Mutect2 --dont-use-soft-clipped-bases true --genotype-germline-sites --min-base-quality-score 20 --max-reads-per-alignment-start 0 --native-pair-hmm-threads 16 -R /nfs/home/references/genomes/human/GRCh38/bwa_mem2/Homo_sapiens.GRCh38.fa -I /media/scratch1/20250314_UB485_runNS252_Oncocarta//analysis_BRC2_100x_sub50/mapping/BRC2-ONC1-1316-S-sub50.bam -L /media/scratch1/20250314_UB485_runNS252_Oncocarta//analysis_BRC2_100x_sub50/QC/mapping/BRC2-ONC1-1316-S-sub50_covered_regions_100x.bed -O /media/scratch1/20250314_UB485_runNS252_Oncocarta//tmp_BRC2_100x_sub50/variant_calling/BRC2-ONC1-1316-S-sub50_mutect2_tmp1.vcf

13:43:57.602 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/nfs/home/software/packages/gatk-4.4.0.0/gatk-package-4.4.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so

13:43:57.634 INFO Mutect2 - ------------------------------------------------------------

13:43:57.636 INFO Mutect2 - The Genome Analysis Toolkit (GATK) v4.4.0.0

13:43:57.637 INFO Mutect2 - For support and documentation go to https://software.broadinstitute.org/gatk/

13:43:57.637 INFO Mutect2 - Executing as fmartinez@vlinuxcervantes6srv on Linux v6.8.0-40-generic amd64

13:43:57.637 INFO Mutect2 - Java runtime: OpenJDK 64-Bit Server VM v18.0.2-ea+9-Ubuntu-222.04

13:43:57.637 INFO Mutect2 - Start Date/Time: March 18, 2025 at 1:43:57 PM CET

13:43:57.637 INFO Mutect2 - ------------------------------------------------------------

13:43:57.637 INFO Mutect2 - ------------------------------------------------------------

13:43:57.638 INFO Mutect2 - HTSJDK Version: 3.0.5

13:43:57.638 INFO Mutect2 - Picard Version: 3.0.0

13:43:57.638 INFO Mutect2 - Built for Spark Version: 3.3.1

13:43:57.639 INFO Mutect2 - HTSJDK Defaults.COMPRESSION_LEVEL : 2

13:43:57.639 INFO Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false

13:43:57.639 INFO Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true

13:43:57.639 INFO Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false

13:43:57.639 INFO Mutect2 - Deflater: IntelDeflater

13:43:57.639 INFO Mutect2 - Inflater: IntelInflater

13:43:57.640 INFO Mutect2 - GCS max retries/reopens: 20

13:43:57.640 INFO Mutect2 - Requester pays: disabled

13:43:57.640 INFO Mutect2 - Initializing engine

13:43:57.801 INFO FeatureManager - Using codec BEDCodec to read file file:///media/scratch1/20250314_UB485_runNS252_Oncocarta/analysis_BRC2_100x_sub50/QC/mapping/BRC2-ONC1-1316-S-sub50_covered_regions_100x.bed

13:43:57.808 INFO IntervalArgumentCollection - Processing 18870 bp from intervals

13:43:57.814 INFO Mutect2 - Done initializing engine

13:43:57.830 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/nfs/home/software/packages/gatk-4.4.0.0/gatk-package-4.4.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so

13:43:57.832 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/nfs/home/software/packages/gatk-4.4.0.0/gatk-package-4.4.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so

13:43:57.845 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM

13:43:57.846 INFO IntelPairHmm - Available threads: 20

13:43:57.846 INFO IntelPairHmm - Requested threads: 16

13:43:57.846 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation

13:43:57.881 INFO ProgressMeter - Starting traversal

13:43:57.882 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute

13:45:02.170 INFO ProgressMeter - chr13:32319323 1.1 10 9.3

13:56:16.365 INFO ProgressMeter - chr13:32326268 12.3 20 1.6

13:56:40.081 INFO ProgressMeter - chr13:32332470 12.7 30 2.4

13:57:35.106 INFO ProgressMeter - chr13:32336649 13.6 40 2.9

14:02:24.112 INFO ProgressMeter - chr13:32337858 18.4 50 2.7

14:11:54.795 INFO ProgressMeter - chr13:32338387 27.9 60 2.1

14:12:35.488 INFO ProgressMeter - chr13:32339573 28.6 70 2.4

14:13:02.528 INFO ProgressMeter - chr13:32340498 29.1 80 2.8

14:13:14.425 INFO ProgressMeter - chr13:32346942 29.3 100 3.4

14:13:46.562 INFO ProgressMeter - chr13:32363154 29.8 110 3.7

14:16:26.377 INFO ProgressMeter - chr13:32370539 32.5 120 3.7

14:20:55.928 INFO ProgressMeter - chr13:32379297 37.0 130 3.5

14:21:56.874 INFO ProgressMeter - chr13:32396932 38.0 140 3.7

14:58:12.923 INFO ProgressMeter - chr17:43047608 74.3 150 2.0

14:58:27.972 INFO ProgressMeter - chr17:43067670 74.5 170 2.3

16:22:33.639 INFO ProgressMeter - chr17:43091576 158.6 190 1.2

16:31:11.549 INFO ProgressMeter - chr17:43092804 167.2 200 1.2

16:31:25.543 INFO ProgressMeter - chr17:43094405 167.5 210 1.3

16:31:48.257 INFO ProgressMeter - chr17:43097284 167.8 220 1.3

16:32:01.839 INFO ProgressMeter - chr17:43104894 168.1 230 1.4

16:32:02.015 INFO Mutect2 - 0 read(s) filtered by: MappingQualityReadFilter

0 read(s) filtered by: MappingQualityAvailableReadFilter

0 read(s) filtered by: MappingQualityNotZeroReadFilter

0 read(s) filtered by: MappedReadFilter

46 read(s) filtered by: NotSecondaryAlignmentReadFilter

0 read(s) filtered by: NotDuplicateReadFilter

0 read(s) filtered by: PassesVendorQualityCheckReadFilter

0 read(s) filtered by: NonChimericOriginalAlignmentReadFilter

0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter

0 read(s) filtered by: ReadLengthReadFilter

43 read(s) filtered by: GoodCigarReadFilter

0 read(s) filtered by: WellformedReadFilter

89 total reads filtered out of 6958926 reads processed

16:32:02.016 INFO ProgressMeter - chr17:43104894 168.1 236 1.4

16:32:02.016 INFO ProgressMeter - Traversal complete. Processed 236 total regions in 168.1 minutes.

16:32:02.057 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 11.263194235

16:32:02.057 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 1027.6256414450002

16:32:02.057 INFO SmithWatermanAligner - Total compute time in java Smith-Waterman : 178.05 sec

16:32:02.058 INFO Mutect2 - Shutting down engine

[March 18, 2025 at 4:32:02 PM CET] org.broadinstitute.hellbender.tools.walkers.mutect.Mutect2 done. Elapsed time: 168.07 minutes.

Runtime.totalMemory()=1811939328

Tool returned:

SUCCESS

0

• 

  Gökalp Çelik

  • March 26, 2025 06:24

  Hi Fran

  This is all about informative reads. Mutect2/HaplotypeCaller uses local reassembly and realignment to capture variants in active regions therefore certain reads may not be used due to filters applied onto them. Once  filtered remaining reads will be used to generate kmers and reassembled haplotypes. So AD and DP values may not match exactly what you see in IGV. 

  https://gatk.broadinstitute.org/hc/en-us/articles/360035532252-Allele-Depth-AD-is-lower-than-expected 

  I hope this helps.

  Regards. 

  1

  Comment actions Permalink
• 

  Fran

  • March 26, 2025 12:41

  Thanks for your answer. I have other questions:

  • Then why I can't see the realigned reads that match the reported AD in the bam I got using the "--bam-output" option. I thought this bam contained the filtered and reassembled reads.

  

  • Why is AF reported as 0.018? Since AD is "1042,1056", I would expect that, following those values, the reported AF to be AF = 1042 / (1056 + 1042) = 0.49.​

   

  0

  Comment actions Permalink
• 

  Gökalp Çelik

  • March 31, 2025 15:22

  Hi again. 

  Can you try running your sample through a more recent version such as 4.6.1.0. There seems to be a bug fix somewhere in between 4.4.0.0 and 4.5.0.0 about the way AD is calculated. So this could be due to that particular bug present in 4.4.0.0. 

  0

  Comment actions Permalink
• 

  alejandra

  • June 12, 2025 14:59

  Hello, 

  I have a similar issue that the AF value does not match the AD values. I would expect that the informative values are used to estimate AF.

  e.g. 

  chr3 51311924 . G C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=111,0|10,0;DP=129;ECNT=1;GERMQ=93;MBQ=38,39;MFRL=231,233;MMQ=60,60;MPOS=33;POPAF=7.3;ROQ=51;TLOD=27.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB 0/1:111,10:0.089:121:58,6:53,4:111,10:111,0,10,0

  Based on this, AF should be 10/121=0.08264, rather than 0.089. 

  The difference is not that big as it's been reported in the comment above, but still there is a difference. Is this a bug? 

   

   

  0

  Comment actions Permalink

Please sign in to leave a comment.