SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BWA MEM mapping to divergent reference ellenell Bioinformatics 5 07-08-2015 03:19 AM
Discrepancy between mapping qualities of bowtie & bwa mem jmartin Bioinformatics 0 02-12-2014 10:02 AM
bwa mem stringent mapping lorendarith Bioinformatics 0 11-13-2013 04:37 AM
how to set parameters for mapping reads to a bacterial genome using "bwa-mem/bowtie"? joyjane88 RNA Sequencing 1 11-13-2013 03:54 AM
how many mismatches does BWA-MEM tolerate during mapping Pengfei Liu Bioinformatics 6 10-04-2013 05:20 AM

Reply
 
Thread Tools
Old 10-28-2015, 05:38 AM   #1
yinshe
Member
 
Location: Norway

Join Date: Nov 2010
Posts: 20
Default bwa mem mapping quality on ambiguous mapping reads

Hello,

I used bwa mem to align 125bp single end reads to human decoy reference genome. I know bwa will assign mapping quality as zero when one read mapped to two or more locations in the genome. However, I noticed some reads which are mapped equally well to different genomic locations, e.g. one read is mapped to equally well to autosome chromosome (chr16) and one of the patches (GL000192). CIGAR for both alignments are 125M. However, mapping quality for the alignment on chr16 is 23, while the alignment mapped to GL000192 got mapping quality of zero. I thought both of them should have mapping quality as zero? Is this right or not?

thanks!
yinshe is offline   Reply With Quote
Old 10-28-2015, 04:43 PM   #2
N311V
Member
 
Location: Australia

Join Date: Jul 2013
Posts: 34
Default

It is also my understanding that mapping quality in that case would be zero for both.

Is there an option to randomly keep one of the multiple mappings rather than discard all of them in bwa mem?
N311V is offline   Reply With Quote
Old 10-29-2015, 03:20 AM   #3
yinshe
Member
 
Location: Norway

Join Date: Nov 2010
Posts: 20
Default

I just put some more detail about this question:

The fastq file used in the alignment is not a fastq file from sequencer. I sliced HYDIN2 sequence into small pieces, each is 125 bp long. I assigned base quality as 30 ("I") for all bases. So all bases have a high base quality. When I did alignment, I asked bwa to output also secondary alignment (using -a option). The record I mentioned here are as following:

b38_1:146691684-146691808 16 16 71053369 23 125M * 0 0 AGCTGAAA.... IIIIIIIIIIII.... NM:i:1 MD:Z:88T36 AS:i:120 XS:i:110
b38_1:146691684-146691808 272 GL000192.1 263206 0 125M * 0 0 * * NM:i:3 MD:Z:5G31G50T36 AS:i:110
yinshe is offline   Reply With Quote
Old 10-29-2015, 03:50 AM   #4
blancha
Senior Member
 
Location: Montreal

Join Date: May 2013
Posts: 367
Default

I can't find anywhere a formal definition for the meaning of MAPQ set to 0 by BWA.
There are only forum posts saying that a MAPQ set to 0 means that a read has multiple hits.

In your example, the second alignment has the NM tag set to 3, meaning the edit distance to the reference (number of nucleotide differences) is 3.
The NM tag is set to 1 in the first alignment.

One could surmise that the 1st alignment is unique in the sense that the second alignment is of such poor quality that it doesn't count.

Admittedly, this is just wild speculation.
There should be a formal definition of MAPQ set to 0 to which aligners should adhere, to make the interpretation of the mapping quality less arduous.

It is certain that the second alignment is of far lesser quality than the first, so it does make sense that the mapping quality is much lower.
blancha is offline   Reply With Quote
Old 10-29-2015, 04:18 AM   #5
yinshe
Member
 
Location: Norway

Join Date: Nov 2010
Posts: 20
Default

Hi blancha,

Thanks for the explanation! But both alignments says 125 base pair matching (CIGAR), so there is no base differences. It seems the SAM record gives different information? Or something I understand wrong?
yinshe is offline   Reply With Quote
Old 10-29-2015, 11:55 AM   #6
blancha
Senior Member
 
Location: Montreal

Join Date: May 2013
Posts: 367
Default

Quote:
But both alignments says 125 base pair matching (CIGAR), so there is no base differences. It seems the SAM record gives different information? Or something I understand wrong?
If you check the official SAM format specification, you'll see that M is for alignment match, and "can be a sequence match or mismatch". 125 bases aligned, but there still can be mismatches, in this case 3.
https://samtools.github.io/hts-specs/SAMv1.pdf

At least, that is my understanding of the convoluted SAM format.

Last edited by blancha; 10-29-2015 at 12:06 PM.
blancha is offline   Reply With Quote
Old 10-30-2015, 12:35 AM   #7
yinshe
Member
 
Location: Norway

Join Date: Nov 2010
Posts: 20
Default

Thanks! This is more clear! It seems 'M' and 'X','=' giving some redundant information.
yinshe is offline   Reply With Quote
Old 10-30-2015, 12:55 AM   #8
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Quote:
Originally Posted by blancha View Post
If you check the official SAM format specification, you'll see that M is for alignment match, and "can be a sequence match or mismatch". 125 bases aligned, but there still can be mismatches, in this case 3.

At least, that is my understanding of the convoluted SAM format.
Yep, that's correct. But the most recent SAM specification reports mismatches in the cigar string, as well. You can see this by mapping with BBMap, which uses the 'X' and '=' symbols.
Brian Bushnell is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 11:59 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO