Samtools: How to extract reads with MAPQ 0?
2
0
Entering edit mode
5.6 years ago

After performing BWA MEM alignment, I want to extract only the reads with MAPQ 0, how to do it?

The "-q" parameter only filter out reads with MAPQ less then a specific value...therefore it is not working for me.

If this is not possible, is there any other strategy?

Here i explain my case: I´m analysing the chloroplast genome. First, i need to filter out all reads with MAPQ less then 40 in order to improve the mapping. Using: samtools view -q 40 input.sam output.bam But...since the reads mapped on the inverted repeats (IRs) regions have MAPQ 0 (because they perfectly map on both IRs), they are removed as well. Therefore, in order to recover them, i would extract the reads with MAPQ 0 e merge them back.

Is there any other way to solve it? Thank you for your help

Assembly sequencing alignment • 6.5k views
ADD COMMENT
0
Entering edit mode

The "-q" parameter only filter out reads with MAPQ less then a specific value...therefore it is not working for me.

What do you mean ? Option -q will extract reads with mapping quality above or equal the threshold. You cannot recover your 0 mapping quality this way.

Even if you are able to get reads with 0 mapping quality you will not get exclusively reads from IRs regions but also reads with bad mapping quality from elsewhere.

I do not know a lot about IRs but if you have the positions of those IRs regions, you can extract reads from IRs falling in those regions using samtools view, then filter your original file on mapping quality equal to 40 and add IRs reads to the result.

ADD REPLY
0
Entering edit mode

Hey Bastien, thank you very much! I followed your suggestion and solved it!

ADD REPLY
0
Entering edit mode

Many thanks to all of you for your help!

ADD REPLY
0
Entering edit mode

Hello stefano.meucci ,

If an answer was helpful, you should upvote it; if the answer resolved your question, you should mark it as accepted. You can accept more than one if they work.

Upvote|Bookmark|Accept

ADD REPLY
0
Entering edit mode
5.6 years ago
GenoMax 147k

Using reformat.sh from BBMap suite. You will need to have samtools available in $PATH. There should be no reads with mapq=1 but if there are then you may need to adjust maxmapq=1.

Check Sam and Bam processing in-line help options for reformat.sh for full list.

reformat.sh in=your.bam out=filtered.bam maxmapq=2
ADD COMMENT
0
Entering edit mode
5.6 years ago

using samjdk http://lindenb.github.io/jvarkit/SamJdk.html

 java -jar dist/samjdk.jar -e 'return record.getMappingQuality()==0;'  input.bam
ADD COMMENT

Login before adding your answer.

Traffic: 2336 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6