Calling variants on reads with MAPQ=0 on HaplotypeCaller or bcftools mpileup
3
0
Entering edit mode
2.1 years ago

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

whole-exome-sequencing variants GATK hg19 • 1.0k views
0
Entering edit mode
2.1 years ago

try

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

0
Entering edit mode

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?

0
Entering edit mode
19 months ago
YT • 0

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

0
Entering edit mode
9 months ago
Gama313 ▴ 110

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.