BWA mem align reads with MQ = 0 to only one region (reference genome has two similar sequences).

Question

[captainlabman]

So, i use WBA mem to align my reads (from NGS panels) to reference genome. My reads align to two different, but very similar regions. I want to prevent align these reads to non target region when i use BWA mem. I hope that will make MQ of my reads > 0 because they will be aligned only to target region. Maybe i can use bed-file contained regions that BWA should ignore to align reads or something like that... I did not found something like that in BWA docs.

Answer 1

[Raony Guimarães]

You could do something like this:

bwa mem genome.fa reads.fastq | samtools view -b -L manifest.bed -U non-overlapping.bam - > overlapping.bam

[captainlabman]

Perhaps I can do that, but after alignment I would like to do a variant calling. As I know variant call tools (GATK) just ignore reads with MQ - 0. And even if i do that u wrote, I anyway get bam file with those reads with MQ = 0.

[Raony Guimarães]

I see, so the best way would be to mask the region, prior to performing the alignment, you could use bedtools maskfasta for that:

bedtools maskfasta -fi hg18.fa -bed probes.500bp.complement.bed \
-fo hg18.probecomplement.masked.fa

https://bedtools.readthedocs.io/en/latest/content/tools/maskfasta.html