From what I understand the information about multi mapping is provided in BWA MEM output in two ways:

1. XA tag that reports, by default, all additional hits if there are 5 hits with score >80% of the max score
2. MAPQ (column 5), which is the lower the larger number of highly scoring supplementary alignments.

I'm confused because I can see reads with XA tags with high MAPQ, for example:

illumina-read:100   163 chr1    892223  60  151M    =   892300  228 TAGGAGGCCAAGGCAGGAGTATCATGAGACCAGCCTGAGCAAAATAGTGAGATGCTAACTCTACAAAAAAAATAAAAATTAGCTGAGCATGGTGGTGCATGCCTGTAATTACAGCTACTTGGGAGGCTGAGGTGGGNGGATCCCTTGAGGG >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>+>>>>>>>>>>>>>>>>>>>>4>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>4>>>>>+>>>>>>!4>44+>>4>>>4>> XA:Z:chr3,+75802872,110M1D41M,6;    MC:Z:151M   MD:Z:136A14

Do you think this may lead to artefacts in variant calling? In particular, when using hg19, where many of the repetitive regions like immunoglobulins or HLA are incomplete, so there will be no XA tag.

If so, are there any best practices to call/interpret variants from repetitive regions of the genome? Does it make sense to mask ambiguous single loci that can be explained by multi-mapping?

Based on the XA tag for the supplementary alignment (110M1D41M), I see that there is 1 deletion and 6 mismatches, so I can imagine that the quality of this alignment is inferior to the primary alignment. But sometimes regions are more ambiguous, having reads with still high MAPQ of 40ish and XA tags with 2-3 hits with 4-3 mismatched. In addition, we know that in WES/targeted sequencing even 10% of mismatch between a probe and a DNA fragment can lead to a successful pull down. This means that the sequencing reads in the region above can be a mixture of molecules coming from the primary and supplementary location producing potential artefacts in variant calling.

FYI my BWA command is as in GATK best practices: bwa mem -K 100000000 -p -v 3 -t 16 -Y and I was using hg38 for the example above.