What are "multireads" with BWA mem?

[Jane M]

Hello everybody,

I try to suppress "multireads", I mean reads with more than 1 best location.
From this forum and biostars, I found different answers, that is why I open a new topic concerning this recurrent question.
Maybe things change upon the versions. Here I am using BWA mem version 0.7.12.

From what I read, I saw that mutireads are reads with:

1. mapping quality=0
2. mapping quality=3 (if there are two best mappings)
3. one high mapping quality and one mapping quality=0
4. reads whith AS=XS
   where AS represents the best alignment and XS the second best alignment

I don't know finally what is right. And I am more confused by some reads:

K00103:19:H2MFLBBXX:4:1101:3953:1598	83	chr9	67316168	24	29S46M	=	67316086	-128	GAAACATCCTTGTGAGGTGTGCACTGAAGTCACAGAGTTGAAACTGTCT
TTTGATTCAGCAGTTTTGAATCTCTC	KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKFFAAA	NM:i:2	MD:Z:14A1A29	AS:i:36	XS:i:36	RG:Z:WES

Here AS=XS and mapping score=24. Not 0, not 3...

More strange for me:

K00103:19:H2MFLBBXX:4:1101:21675:1598	99	chr21	9540637	31	68M7S	=	9540815	245	GCCGCCTCTCCGCCGCGCCGCCGCTGGCCGCCTNTCCGCCGCGCCAACGCTGGCCGCCTCTCCGCCGCTGTCCGC	AAFFFKKKKKKKKKKKKKKKK
KKKKKKKKKKKK#KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK	NM:i:6	MD:Z:0T25T6C11G0C5T15	AS:i:45	XS:i:47	RG:Z:WES
K00103:19:H2MFLBBXX:4:1101:21675:1598	147	chr21	9540815	31	67M	=	9540637	-245	TGCGTTCTGCTCTGCACAGACCTTGGGGCACTGCCTCGCTTTGGGACAACTCGGGGCCGCATCGACG	KKKKKKKKKKKKKKKKKFFAKKKFFFAKK
FFFFFFKKKKKKKKAKKKKKKKKKKKKKKKKKKFFFAA	NM:i:1	MD:Z:0C66	AS:i:66	XS:i:67	RG:Z:WES

XS>AS
Does this mean that the best alignment maps less bases but with a higher score?
In the second read in the example, 67 bases are mapped over 67. How could it give a lower score than the alignment of 66/67 bases?

Thank you for some clarifications,
Jane

[dpryan]

The AS==XS thing is specific to bowtie2. For bwa just filter by MAPQ. Actually, just do that for everything and use some reasonable threshold (5 or so is sufficient for most purposes).

[Jane M]

Thank you for your answer.
I intended to use a mapping score of 20. I am interesting in pretty highly mapped reads.

If it is not with AS==XS, does someone know how to identify multireads?
How is it possible to get a XS score higher than a AS score?

Jane

[dpryan]

The meaning of XS (in fact, any tag starting with "X", "Y", or a lower case letter) is dependent on the aligner (I actually haven't a clue what it means for bwa). I think bwa just sets an appropriately low MAPQ. If you filtered out everything below a MAPQ of 5 you wouldn't have any alignments whose next best was equally good.

[Jane M]

Thank you again.

If you filtered out everything below a MAPQ of 5 you wouldn't have any alignments whose next best was equally good.

That is what I would like to be sure of. But I am afraid that reads with 2 best equal locations can have pretty good mapping score as in my first example (assumming I understood well AS and XS):

K00103:19:H2MFLBBXX:4:1101:3953:1598	83	chr9	67316168	[B]24[/B]	29S46M	=	67316086	-128	GAAACATCCTTGTGAGGTGTGCACTGAAGTCACAGAGTTGAAACTGTCT
TTTGATTCAGCAGTTTTGAATCTCTC	KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKFFAAA	NM:i:2	MD:Z:14A1A29	[B]AS:i:36[/B]	[B]XS:i:36[/B]	RG:Z:WES

[dpryan]

It's likely that at least one of the scores reported isn't what's actually used to calculate the MAPQ score. If it'll put your mind at ease, generate some fake data (e.g., with WGSIM) and align that. Then look at accuracy as a function of MAPQ score. You'll then be able to make an educated choice (and also not need to go through the code).

[HESmith]

I believe that the apparent discrepancies are a consequence of having paired-end reads. Pairing information is used as a criterion for identifying the primary alignment. Secondary alignments for the individual reads may occasionally have higher scores, but positions that are inconsistent with pairing. You can confirm by aligning as single end and reporting the multi-mapped alignments.

I second Devon's recommendation to use MAPQ for filtering.

[dpryan]

I believe that the apparent discrepancies are a consequence of having paired-end reads.

Good catch!

[HESmith]

Like the proverbial blind squirrel...

[Jane M]

Thank both of you for your answers.

The percentage of reads with MAPQ <20 is around 2.5-2.8% for my 10 samples. I will probably take this threshold but, to be sure, I am checking this percentage for lower MAPQ scores but.