Read sorted and indexed BAM files faster in C++? I have some sorted and indexed alignment large BAM files. I'd like to read all the reads, and apply custom operation on each of them. The orders aren't important. I'm using htslib in C++ to do the reading. I have a single powerful machine.

The problem is ... reading a large BAM file is very slow even in C++.

Q: Is there any trick to speed up the reading? Reading BAM files multithreaded (I have many cores)?
In this answer I want to show you some benchmarks that compares three different serial ways of reading data in C++, the third one being the fastest. 

In the first method I use an std::istreambuf_iterator, in the second I read the file line by line into a std::vector and in the third one I use C flavoured operations.

/*
compile with g++ -std=c++11 -O2 reading_files.cpp -o reading_files
test the program with the following
dd if=/dev/random of=./dummy5GiB bs=1024 count=$[1024*1024*5]
 for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
*/

#include &lt;fstream&gt;
#include &lt;iostream&gt;
#include &lt;vector&gt;
#include &lt;cstring&gt; //for strcmp
#include &lt;string&gt;
#include &lt;algorithm&gt; //for min

using namespace std;

int main(int argc, char * argv[])
{
    if(argc != 3) {
        cerr &lt;&lt; "usage: reading_files &lt;path&gt; &lt;number 0 or 1 or 2&gt;\n";
        return 1;
    }

    ifstream file(argv[1], ios::in|ios::binary); 
    if(!file.is_open()) {
        cerr &lt;&lt; argv[1] &lt;&lt; " file not opened!\n";
        return 1;
    }
    cout &lt;&lt; "---------------------------\n";
    if(!strcmp(argv[2],"0")) {
        cout &lt;&lt; "method 0\n";
        vector&lt;char&gt; content((std::istreambuf_iterator&lt;char&gt;(file)), std::istreambuf_iterator&lt;char&gt;());
        cout &lt;&lt; content.size() &lt;&lt; " bytes read\n";
    } else if(!strcmp(argv[2],"1")) {
        cout &lt;&lt; "method 1\n";
        string line;
        vector&lt;string&gt; content;
        while(getline(file,line)) {
            content.push_back(line);
        }
        cout &lt;&lt; content.size() &lt;&lt; " lines read\n";
    } else if(!strcmp(argv[2],"2")) {
        cout &lt;&lt; "method 2\n";
        file.seekg(0, std::ios::end);
        unsigned long long file_size = file.tellg();
        cout &lt;&lt; "file size is " &lt;&lt; file_size &lt;&lt; "\n";
        //maximum 5 GiB of memory, so to fit my RAM
        unsigned long long max_block_size = 1024ULL * 1024ULL * 1024ULL * 5; 
        file.clear();
        file.seekg(0, std::ios::beg);

        unsigned long long block_size = min(max_block_size,file_size);
        char * buffer = (char*)malloc(block_size*sizeof(char));
        unsigned long long processed = 0;

        while(!file.eof()) {
            file.read(buffer, block_size);
            cout &lt;&lt; "read a block of " &lt;&lt; min(file_size - processed, block_size) &lt;&lt; " bytes\n";
            //do work here with the block in memory (if the block size is &gt; 0)
            processed += block_size;
        }

        free(buffer);
    }
    file.close();
    return 0;
}


In the third method I read the file in blocks of 5GiB, so to fit the free RAM of my machine.

Here is the benchmark of the program on a 5 GiB randomly generated file.

bash-3.2$ for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
---------------------------
method 0
5368709120 bytes read

real    0m35.194s
user    0m19.665s
sys 0m12.999s
---------------------------
method 1
20973135 lines read

real    0m40.478s
user    0m34.046s
sys 0m5.870s
---------------------------
method 2
file size is 5368709120
read a block of 5368709120 bytes
read a block of 0 bytes

real    0m4.413s
user    0m1.757s
sys 0m2.617s


Here is the benchmark on a real 32GiB .bam file using only the third method, here the splitting of the file read plays a fundamental role. I can't tell you how big because my hard disk memory gets saturated before the end of the execution and so I have to kill the program.

time ./reading_files ./HG00252.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam 2
---------------------------
method 2
file size is 34316054058
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 2103799338 bytes

real    0m50.656s
user    0m9.135s
sys 0m35.728s


The key point is to minimise the file operations reading large block of data at once (as long as they fit the memory) and operate on them.

Using mmap could lead to an improvement, but I don't think it would be as large as using the third method instead of the first two.