Biostar Beta. Not for public use.
How To Get Number Of Reads In Bam File Efficiently In Python?
15
Entering edit mode
6.1 years ago
User 9996 • 800

how can I get the number of reads total in a BAM file, using Pysam (or something else Pythonic) in Python? I want to avoid iterating through the file and counting each read since it takes a very long time and is inefficient.

thank you.

ADD COMMENTlink
0
Entering edit mode

I am not sure about doing this on Python but guess we can get the nuber of reads in BAM file using FASTQC.

ADD REPLYlink
7
Entering edit mode
6.1 years ago
Jts ♦ 1.2k

I don't think you can get the number of reads from pysam without iterating through the file. See this post to the samtools mailing list about doing the same with the java api.

I would make an external call to samtools flagstat as suggested.

ADD COMMENTlink
4
Entering edit mode

Perhaps this is obvious, but if you're going to be getting these read counts many times, it makes sense to save the counts so that you don't have to iterate through multiple times. This can be as simple as counting the reads once and dumping those counts to a file in the same directory as the BAM files.

ADD REPLYlink
0
Entering edit mode

I would use samtools idxstats as it's the fastest solution...

ADD REPLYlink
7
Entering edit mode
6.1 years ago
Marvin • 850

Straight forward without iterating through the whole file:

samtools idxstats file.bam

Naturally, that requires the file to be indexed. I don't know about a Python API to do that, but a pipe should work.

ADD COMMENTlink
1
Entering edit mode

This is the preferred solution when the BAM is coordinate sorted.

ADD REPLYlink
1
Entering edit mode

+1, this is the fastest solution when your bam is sorted and indexed.

ADD REPLYlink
0
Entering edit mode

and aggregate to obtain mapped + unmapped reads

samtools idxstats file.bam | awk -F '\t' '{s+=$3+$4}END{print s}'
ADD REPLYlink
0
Entering edit mode

I think different versions of idxstats produces different outputs

ADD REPLYlink
6
Entering edit mode
6.1 years ago
Chris Penkett • 480
Cambridge, UK

Short python/pysam answer:

import pysam
filename = 'test.bam'
print reduce(lambda x, y: x + y, [ eval('+'.join(l.rstrip('\n').split('\t')[2:]) ) for l in pysam.idxstats(filename) ])
ADD COMMENTlink
0
Entering edit mode

Cool! That's a fast way for pysam.

If you want mapped reads, you can switch the third line as:

print reduce(lambda x, y: x + y, [ int(l.rstrip('\n').split('\t')[2]) for l in pysam.idxstats(filename) ])

and for unmapped read number:

print reduce(lambda x, y: x + y, [ int(l.rstrip('\n').split('\t')[3]) for l in pysam.idxstats(filename) ])
ADD REPLYlink
3
Entering edit mode
6.1 years ago
Marcelo • 30

The straightforward way:

samtools view -c file.bam
ADD COMMENTlink
1
Entering edit mode

Unfortunately, that prints the number of alignments in the file, which can be (and usually is) distinct from the number of reads (because each read can align more than once).

ADD REPLYlink
1
Entering edit mode
6.1 years ago
Ying • 10

Multiple different ways (including flagstat suggestion) here

ADD COMMENTlink
1
Entering edit mode
12 months ago
Cambridge, MA

I like Chris's answer, but if you don't have pysam:

https://gist.github.com/mdshw5/73c6591237c7a9f88518

This function will do the same thing and takes the same amount of time regardless of the number of reads in your file as long as it's indexed.

ADD COMMENTlink
0
Entering edit mode
15 months ago
John 12k
Germany

There is unfortunately no efficient way to read a whole file without reading a whole file - there are however shortcuts with a time/accuracy tradeoff. For the SeQC project, we calculate a read count estimate with:

import os
import subprocess
import pysam
def bamCheck(fileName):
    xamfile = open(fileName, 'rb')
    try:
        for i in range(10):
            section = xamfile.read(2048)
            if '\0' in section: return True
            if len(section) < 2048: break
    finally: xamfile.close()
    return False

def guessReads(fileName,samtools):
    DEVNULL = open(os.devnull, 'wb')
    bytesUsed = 10000000
    sizeInBytes = int(subprocess.check_output("ls -ln '" + str(fileName) + "' | awk '{print $5}'", stderr=subprocess.STDOUT, shell=True))
    if bamCheck(fileName):
        if bytesUsed*2 < sizeInBytes:
            readsIn100MB = int(subprocess.check_output('head -c ' + str(bytesUsed) + ' "' + str(fileName) + '" | "' + samtools + '" view - | wc -l', stderr=DEVNULL, shell=True))
            estimatedReads = (readsIn100MB/float(bytesUsed))*sizeInBytes
            estimatedReads = estimatedReads * 1.1 #basically, because of the BAM header (or worse, SAM header) we typically underestimate the number of reads by between 8 and 12 percent. So adding an extra 10% on to this number gets it closer to accurate :) Of course, this all depends on the bytesUsed to sample the file - the bigger the number, the more accurate our guess, but the slower it will take to estimate read count.
        else:
            estimatedReads = int(subprocess.check_output('"' + samtools + '" view -c "' + str(fileName) + '"', stderr=DEVNULL, shell=True)) # because if we can read the whole file in twice the time as it took to sample it, we might as well just read the whole file!
    else:
        bytesIn10000Lines = int(subprocess.check_output('head -11000 "' + str(fileName) + '" | tail -10000 | wc -c', stderr=DEVNULL, shell=True))
        estimatedReads = ( sizeInBytes/float(bytesIn10000Lines) )*10000
        estimatedReads = int(estimatedReads * 1.05) # because we skipped the header by head/tailing the file, we have a much more accurate value here for the total number of reads, so we only need to increase by about 5% to make sure our estimates are conservative/unreachable ;)
    return int(estimatedReads)

>>> guessReads('44_Mm01_WEAd_C23_H3K27me3_F.bam','/package/samtools/samtools')
237033936

Hope that helps someone - I originally wrote it to estimate read counts for use in a progress bar - so i actually much preferred a slight over estimate than an underestimate (causing 101% of file read..) - but you can tweak all that kind of stuff pretty obviously.

The only other possible answer to the OP with regards to 'efficient without reading whole file' would be to find the total read count as a byproduct of reading the whole file during some other process. People have given the example of using the index (which is created by reading the whole file), as a source for getting the total read count for little extra work - but actually any process which reads the whole file (say, MD5'ing the file) could also be used where chunked bytes are also piped to pysam and linecount++'d

Since you're IO bound anyway, using pysam in such a situation is, strictly speaking, even more efficient than using the index ;) hahah

ADD COMMENTlink

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.1