Reject reads with low quality bases from a Bam file through pysam

Question

I have a code below:

def FilterReads(in_file, out_file):

def read_ok(read):
    """
    read_ok - reject reads with a low quality (<5) base call
    read - a PySam AlignedRead object
    returns: True if the read is ok
    """
    if any([ord(c)-33 < _BASE_QUAL_CUTOFF for c in list(read.qual)]):
        return False
    else:
        return True

_BASE_QUAL_CUTOFF = 30

bam_in = pysam.Samfile(in_file, 'rb')
bam_out = pysam.Samfile(out_file, 'wb', template=bam_in)


for read in bam_in.fetch():
    if read_ok(read):
        bam_out.write(read)

This code works fine by

Rejecting reads having a base with phred quality score below five
But it first takes a BAM file
And then creates a filtered BAM file for further analysis but this takes a lot of time.

So is there a way to reject these reads using pileup in pysam so that I may not have to create a file then sort it and again read it? Or can I modify this code to perform the same function

Why are you sorting the file prior to iterating over it? Your function doesn't require it to be sorted and since you are not using the random indexing pysam doesn't require it either. — Bioathlete, Jul 31 '17 at 21:11
if …: return False else: return True is an anti-pattern. Just write return not …. Or, better yet, logically invert the condition; in your case, you could write return all([ord(c) - 33 >= _BASE_QUAL_CUTOFF for c in list(read.qual)]) (this follows from straightforward application of de Morgan’s laws). — Konrad Rudolph, Aug 24 '17 at 12:20
This is probably not the cause of the slowdown, but you can avoid two list creations per read by refactoring line 9 as any(ord(c)-33 < _BASE_QUAL_CUTOFF for c in read.qual): — winni2k, Nov 23 '18 at 09:26
I don't think pileup is going to save you any time here, as it has to load the same data that you are loading with the fetch method. — winni2k, Nov 23 '18 at 09:32

Scott Gigante · Answer 1 · 2018-11-24T18:25:05.413

For can do this by accessing the basecall qualities from PileupRead.alignment. For example:

_BASE_QUAL_CUTOFF = 30
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            qual = read.alignment.query_qualities
            if all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual]):
                # do some stuff

This involves checking every read at every position. You could memoise this task:

def check_quality(read):
    try:
        return passed_quality[read.query_name]
    except KeyError:
        qual = read.alignment.query_qualities 
        passed = all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])
        passed_quality[read.query_name] = passed
        return passed_quality[read.query_name]

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

Or, using functools.lru_cache (thanks @winni2k):

from functools import lru_cache

@lru_cache(maxsize=1000)
def check_quality(read):
    qual = read.alignment.query_qualities 
    return all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

FYI: There's a decorator in the library that does this and more! https://docs.python.org/3/library/functools.html#functools.lru_cache — winni2k, Nov 23 '18 at 09:30

Reject reads with low quality bases from a Bam file through pysam

1 Answers1