Question: Samtools: How to extract reads with MAPQ 0?
0
gravatar for stefano.meucci
5 months ago by
stefano.meucci0 wrote:

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

sequencing alignment assembly • 399 views
ADD COMMENTlink modified 3 months ago • written 5 months ago by stefano.meucci0

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 REPLYlink modified 5 months ago • written 5 months ago by Bastien Hervé4.4k

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

ADD REPLYlink written 3 months ago by stefano.meucci0

Many thanks to all of you for your help!

ADD REPLYlink written 3 months ago by stefano.meucci0

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 REPLYlink modified 3 months ago • written 3 months ago by finswimmer12k
0
gravatar for genomax
5 months ago by
genomax73k
United States
genomax73k wrote:

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 COMMENTlink modified 5 months ago • written 5 months ago by genomax73k
0
gravatar for Pierre Lindenbaum
5 months ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum123k wrote:

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

 java -jar dist/samjdk.jar -e 'return record.getMappingQuality()==0;'  input.bam
ADD COMMENTlink written 5 months ago by Pierre Lindenbaum123k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2350 users visited in the last hour