Calling variants on reads with MAPQ=0 on HaplotypeCaller or bcftools mpileup

Question

[Anupriya Dalmia]

I am working with about 500 samples of human exome data. used hg19 to align my reads and ran a standard best-practices GATK workflow. Later only to realise that a small 1Mb loci has not mapped properly due to the presence of a haplotype specific contig in the hg19 FASTA.

How do I force HaplotypeCaller or BCFTools to call a variant with mapping quality 0? I have disabled the read filters but that does not work. I just want the variants from one particular gene and since I have 500 samples, re-aligning them to a modified reference is not an option at the moment due to time constraints. I am aware of false positives occurring but we will validate any variants seen.

This is the command I've tried on BCFtools:

bcftools mpileup -q 0 -f ucsc.hg19.fasta file.bam -r chr17:43000000-45000000 | bcftools call -mv > file.vcf

this only outputs variants with MQ =/= 0

Answer 1

[Pierre Lindenbaum]

try

bcftools mpileup -q -1  -f ucsc.hg19.fasta file.bam -r chr17:43000000-45000000 | bcftools call -mv > file.vcf

[Anupriya Dalmia]

Thank you for responding! Unfortunately, this command returns the same result as the one I have posted above.

I've checked it on IGV and the variants are definitely present.

Can you think of any other solution to this problem, perhaps?

Answer 2

[YT]

Hi,

-q 0 is the default value. Maybe try bcftools mpileup -A -f reference.fa Try with the -A option. Hope it works for you.

Best wishes, Yutang

Answer 3

[Gama313]

As far as I know, a MAPQ equal to 0 identifies reads mapping to multiple positions (i.e. pseudogenes). This means that you cannot call that mutation with enough confidence. If the region of interest is filled with mutations, it is definetly misaligned. Try blast sequences supporting your mutation and see where they come from. Let me know how it goes.