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
•
written
by
Anupriya Dalmia •
try
bcftools mpileup -q -1 -f ucsc.hg19.fasta file.bam -r chr17:43000000-45000000 | bcftools call -mv > file.vcf
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?
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.
Log in to answer this question.
