I'm trying to align reads that belong to high similarity regions (pseudogenes)
I'm aligning reads with BWA (last version) using tha -a option.
When observing the BAM file, I notice that the reads are aligned with MAPQ0. Which is normal.
The problem is, every read is present once in the BAM file, i don't find all alignments (even with -a option).
I tried to extract one of these reads with MAPQ 0, and create a FASTQ (r1 & r2) file with only 1 read inside, and then align this read with BWA (using exactly same command as for the complet fq file), and I obtain this time a BAM file with all the possible alignments !!
How is it possible ?
thanks
I think the reason has shown in the flags. When aligning only one pair of reads (let's say R1 and R2), bwa didn't think your R1 and R2 are in proper pair, so the -a option (which will output all alignments for unpaired paired-end reads) led to the 4 alignments you see
But when aligning the whole FQ file, bwa marked R1 and R2 as proper pair, so the -a option won't report other alignments since it thought there was no need to report other alignments.
What made this difference I think is due to the criteria used by bwa on proper pair reads, it requires R1 and R2 have the right orientation and the correct insert size, but threshold for insert is determined by the insert size distribution from all data, that's why bwa mem has different behavior in your situation.
BWA MEM manual: http://bio-bwa.sourceforge.net/bwa.shtml Similar problem: Bwa-Mem And Sam-Flag "Read Paired In Proper Pair"
Thanks for this explanation !
Then why the read have MAPQ 0 when aligning the whole FQ, while in the BAM file I have an uniq alignement ?
Though bwa mem marked them as proper pair, the MAPQ still shows they are not unique map reads. In both situations, R1 and R2 are non-unique mapped reads. It seems bwa mem chose one alignment per read to be a proper pair, but I don't know how this selection is done.
So to resume, BWA gave a MAPQ 0 to reads because they map to multiple regions (without taking into account the pair information). But In the BAM file, I found only one occurrence per read-pair because the others combinations are not Proper Pairs ? (so -a would be useful only in case of multiple proper pairs ?)
How can avoid I remove correct these MAPQ 0 in order to use these reads for variant calling later ?
Your understanding of MAPQ 0 is right.
Maybe I didn't explain clearly, the default behavior of bwa is to report one random position for each read if having multiple alignments, so you will have one occurrence per read-pair when aligning the whole FQ or single read pair without -a option.
What makes difference is, when aligning single read pair, bwa doesn't think this read pair is a proper pair, so the -a option lets it report the other two alignments it found. But when aligning the whole FQ, -a won't let bwa report other alignments because it thought this read pair is a proper pair
-a option: Output all found alignments for single-end or unpaired paired-end reads.
I didn't get the meaning of your last question, can you be more specific?
So their are is no way to get multiple Proper Pairs (even with -a) ? BWA will take a random choice ?
My second question: since the MAPQ 0 will impact the variant calling process, can I replace the MAPQ0 by the true MAPQ ?
i will try Bowtie2. Thanks alot !
Log in to answer this question.
Could you please show the command you used and the corresponding output?
when aligning only one read:
When Algning the whole FQ file:
The length of the quality string is shorter that the sequence, I think part of the sam record is being left out - part that could help explain this behaviour, as it could contains the optional tags describing additional features from these reads.