samtools filter PE reads on mapping quality of only one of the reads

[pguilha]

Hi,

I'm trying to filter a BAM file using samtools for mapping quality, but I want to keep both reads in a pair if at least one of them has MAPQ>=10.

My understanding is that using

-q 10

as an option would remove all reads with MAPQ<10 regardless of whether the other read in the pair is of high quality. Is there a way to get around that?

Thanks

[dpryan]

There's no way to do that directly with samtools, you'll need to script something. It'd be pretty straight forward with pysam, if you know a bit of python (I'm in the middle of a bunch of C code at the moment or I'd through together a quick script for you).

[pguilha]

ahh, I was afraid that would be the case...Thanks for your reply anyway!
I don't know python so I'm going to have to do it painfully in R.

[dpryan]

Too bad, I recall that the R method of dealing with SAM/BAM files is to read the whole thing into memory, which would really suck. Post the script you use when you have it working in case someone else ends up needing this in the future.

[pguilha]

Yep, that's exactly why it's going to be painful...Sure, I'll post it if/when it works!

[dpryan]

I ended up getting side-tracked by a dissection anyway, so here's a quick method using python and the pysam package. It assumes that you have everything name-sorted and that there aren't any random single-end reads in the file:

#!/usr/bin/env python
import pysam
import argparse
import sys

parser = argparse.ArgumentParser(description='Filter a paired-end BAM file such that both reads have MAPQs >= a threshold.')
parser.add_argument('threshold', type=int, help="MAPQ threshold")
parser.add_argument('input_file', help="input file name")
parser.add_argument('output_file', help="output file name")
args = parser.parse_args()
if(args.threshold == None or args.input_file == None or args.output_file == None) :
    parser.print_help()
    sys.exit()

#One could look at the file name too see if we're dealing with a SAM or BAM file, but this is simpler
ifile = pysam.Samfile(args.input_file, "rb")
ofile = pysam.Samfile(args.output_file, "wb", template=ifile)

#iterate
while(1) :
    try :
        r1 = ifile.next()
        r2 = ifile.next()
    except :
        break

    if(r1.mapq >= int(args.threshold) and r2.mapq >= int(args.threshold)) :
        ofile.write(r1)
        ofile.write(r2)

ifile.close()
ofile.close()

[pguilha]

Thanks a lot for doing this.
Quick question: what I'd like is to filter the file so that AT LEAST ONE of the reads has MAPQs>= threshold (as opposed to BOTH);

For that should I just change this section of your code:

if(r1.mapq >= int(args.threshold) and r2.mapq >= int(args.threshold)) :

to this:?

if(r1.mapq >= int(args.threshold) or r2.mapq >= int(args.threshold)) :

or would that not keep pairs where both are above the threshold?
thanks again!

[dpryan]

That should work. "or" isn't an "exclusive or" in python.