Answers.Body: <p>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. </p>
  
  <p>In the first method I use an <code>std::istreambuf_iterator</code>, in the second I read the file line by line into a <code>std::vector</code> and in the third one I use C flavoured operations.</p>
  
  <pre><code>/*
  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;
  }
  </code></pre>
  
  <p>In the third method I read the file in blocks of 5GiB, so to fit the free RAM of my machine.</p>
  
  <p>Here is the benchmark of the program on a 5 GiB randomly generated file.</p>
  
  <pre><code>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
  </code></pre>
  
  <p>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.</p>
  
  <pre><code>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
  </code></pre>
  
  <p>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.</p>
  
  <p>Using <code>mmap</code> 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.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-02T14:51:11.020Z
Answers.Text: Why is your “method 2” using gratuitous manual (C-style!) memory management instead of a vector? — I’m also not sure how much this actually helps OP since the question is about parsing BAM, and as far as I know neither Staden IO nor htslib supports reading from a custom 
  buffer (but should perform its own buffering anyway).
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-30T16:08:15.067Z
Answers.LastActivityDate: 2018-03-30T20:06:50.703Z
Answers.LastEditDate: 2018-03-30T20:06:50.703Z
Answers.UserDisplayName: Nisba
Answers.Body: <p>If you need to iterate through the entire BAM file, decompression is usually the main bottleneck.  Three things you can do:</p>
  
  <ul>
  <li>Install libdeflate (<a href="https://github.com/ebiggers/libdeflate" rel="nofollow noreferrer">https://github.com/ebiggers/libdeflate</a> ) and then the development version of htslib, and make sure the latter is compiled to use libdeflate.  This roughly doubles decompression efficiency in my experience.</li>
  <li>Call htslib's bgzf_mt() function right after opening the file for reading, to request multiple decompressor threads.  I generally go with 2-6, depending on what I'm actually doing with the reads.</li>
  <li>Have a separate read/decompress-ahead thread perform the bam_read1() calls.  Ideally, when your main thread is processing read #X, background threads are already working on decompressing reads e.g. (X+100) through (X+200).</li>
  </ul>
  
  <p>Also, if your workload is "embarrassingly parallel", you may also benefit from launching multiple processes simultaneously, configuring each to handle a similarly-sized chunk of the BAM.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-28T20:21:56.647Z
Answers.Text: I stumbled upon this when looking for ways to improve my BAM processing and I was wondering if you have a code example that you could share?
Answers.UserDisplayName: Paghillect
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-30T16:47:21.160Z
Answers.LastActivityDate: 2018-03-30T16:47:21.160Z
Answers.UserDisplayName: Christopher Chang
Answers.Body: <p>In both <a href="https://github.com/fidelram/deepTools" rel="nofollow noreferrer">deepTools</a> and <a href="https://github.com/dpryan79/MethylDackel" rel="nofollow noreferrer">MethylDackel</a> we implement multithreading by having individual threads process separate genomic regions. At some point you hit an IO limit, but unless your processing is trivial then IO won't typically be limiting with single-threaded applications. Since you already have sorted/index BAM files, you can take advantage of the same strategy.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-10-13T06:39:20.667Z
Answers.LastActivityDate: 2017-10-13T06:39:20.667Z
Answers.UserDisplayName: Devon Ryan
Body: <p>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 <code>htslib</code> in C++ to do the reading. I have a single powerful machine.</p>
  
  <p>The problem is ... reading a large BAM file is very slow even in C++.</p>
  
  <p><strong>Q:</strong> Is there any trick to speed up the reading? Reading BAM files multithreaded (I have many cores)?</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-13T00:25:37.500Z
Comments.Text: If you're reading the entire file, indexing won't help. That's for quickly seeking to specific regions in a BAM file.
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-13T01:08:09.460Z
Comments.Text: @gringer Can I speed up by splitting the file into partitions, then each of those by multithreading? I have many cores to do that.
Comments.UserDisplayName: SmallChess
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-13T03:06:08.970Z
Comments.Text: Reading data is an I/O problem, not a CPU one. Do you have your data on an SSD?
Comments.UserDisplayName: Alex Reynolds
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-13T03:07:00.293Z
Comments.Text: @AlexReynolds Yes. But is there anything else we can do in programming?
Comments.UserDisplayName: SmallChess
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-13T04:55:12.657Z
Comments.Text: Using mmap routines to read data, maybe? It might give you a slight benefit over reading with fopen/fseek/fread etc. But there's no way to guarantee a benefit, short of benchmarking whatever it is you are trying to do. And you have to be careful about purging caches between tests. It's easy for people to make spurious claims about their tools being fast if they aren't purging the read cache.
Comments.UserDisplayName: Alex Reynolds
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-13T05:01:23.917Z
Comments.Text: You might investigate terms like "striping" or "parallel I/O" to gain benefit from multiple read operations via hardware. Having multiple processors or cores can be useful for algorithms, but I don't think it will benefit reading sequential data from storage.
Comments.UserDisplayName: Alex Reynolds
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-13T22:07:53.893Z
Comments.Text: The optimal solution is entirely determined by the "custom operation" you want to apply. What do you want to do exactly?
Comments.UserDisplayName: user172818
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-14T04:45:15.497Z
Comments.Text: What do you mean by "very slow"? How does it compare with `time cat <file.bam> > /dev/null`? How does it compare to the code when nothing is done with the records?
Comments.UserDisplayName: gringer
Tags: <sam><c++>
Title: Read sorted and indexed BAM files faster in C++?
ThreadText: # Read sorted and indexed BAM files faster in C++?
  
  (Asked by: username_0 on 2017-10-12)
  
  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)?
  
  **Tags:** sam, c++
  
  > **Comments:**
  > **username_1**: If you're reading the entire file, indexing won't help. That's for quickly seeking to specific regions in a BAM file.  
  > **username_0**: @username_1 Can I speed up by splitting the file into partitions, then each of those by multithreading? I have many cores to do that.  
  > **username_2**: Reading data is an I/O problem, not a CPU one. Do you have your data on an SSD?  
  > **username_0**: @username_2 Yes. But is there anything else we can do in programming?  
  > **username_2**: Using mmap routines to read data, maybe? It might give you a slight benefit over reading with fopen/fseek/fread etc. But there's no way to guarantee a benefit, short of benchmarking whatever it is you are trying to do. And you have to be careful about purging caches between tests. It's easy for people to make spurious claims about their tools being fast if they aren't purging the read cache.  
  > **username_2**: You might investigate terms like "striping" or "parallel I/O" to gain benefit from multiple read operations via hardware. Having multiple processors or cores can be useful for algorithms, but I don't think it will benefit reading sequential data from storage.  
  > **username_3**: The optimal solution is entirely determined by the "custom operation" you want to apply. What do you want to do exactly?  
  > **username_1**: What do you mean by "very slow"? How does it compare with `time cat <file.bam> > /dev/null`? How does it compare to the code when nothing is done with the records?  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_4 on 2018-03-30. Score: 4
  
  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 <fstream>
  #include <iostream>
  #include <vector>
  #include <cstring> //for strcmp
  #include <string>
  #include <algorithm> //for min
  
  using namespace std;
  
  int main(int argc, char * argv[])
  {
      if(argc != 3) {
          cerr << "usage: reading_files <path> <number 0 or 1 or 2>\n";
          return 1;
      }
  
      ifstream file(argv[1], ios::in|ios::binary); 
      if(!file.is_open()) {
          cerr << argv[1] << " file not opened!\n";
          return 1;
      }
      cout << "---------------------------\n";
      if(!strcmp(argv[2],"0")) {
          cout << "method 0\n";
          vector<char> content((std::istreambuf_iterator<char>(file)), std::istreambuf_iterator<char>());
          cout << content.size() << " bytes read\n";
      } else if(!strcmp(argv[2],"1")) {
          cout << "method 1\n";
          string line;
          vector<string> content;
          while(getline(file,line)) {
              content.push_back(line);
          }
          cout << content.size() << " lines read\n";
      } else if(!strcmp(argv[2],"2")) {
          cout << "method 2\n";
          file.seekg(0, std::ios::end);
          unsigned long long file_size = file.tellg();
          cout << "file size is " << file_size << "\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 << "read a block of " << min(file_size - processed, block_size) << " bytes\n";
              //do work here with the block in memory (if the block size is > 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.
  
  > **Comments:**
  > **username_5**: Why is your “method 2” using gratuitous manual (C-style!) memory management instead of a vector? — I’m also not sure how much this actually helps OP since the question is about parsing BAM, and as far as I know neither Staden IO nor htslib supports reading from a custom  buffer (but should perform its own buffering anyway).  
  
  ---
  
  ### Answer by username_6 on 2018-03-30. Score: 3
  
  If you need to iterate through the entire BAM file, decompression is usually the main bottleneck. Three things you can do:
  
  * Install libdeflate (<https://github.com/ebiggers/libdeflate> ) and then the development version of htslib, and make sure the latter is compiled to use libdeflate. This roughly doubles decompression efficiency in my experience.
  * Call htslib's bgzf\_mt() function right after opening the file for reading, to request multiple decompressor threads. I generally go with 2-6, depending on what I'm actually doing with the reads.
  * Have a separate read/decompress-ahead thread perform the bam\_read1() calls. Ideally, when your main thread is processing read #X, background threads are already working on decompressing reads e.g. (X+100) through (X+200).
  
  Also, if your workload is "embarrassingly parallel", you may also benefit from launching multiple processes simultaneously, configuring each to handle a similarly-sized chunk of the BAM.
  
  > **Comments:**
  > **username_7**: I stumbled upon this when looking for ways to improve my BAM processing and I was wondering if you have a code example that you could share?  
  
  ---
  
  ### Answer by username_8 on 2017-10-13. Score: 2
  
  In both [deepTools](https://github.com/fidelram/deepTools) and [MethylDackel](https://github.com/dpryan79/MethylDackel) we implement multithreading by having individual threads process separate genomic regions. At some point you hit an IO limit, but unless your processing is trivial then IO won't typically be limiting with single-threaded applications. Since you already have sorted/index BAM files, you can take advantage of the same strategy.
  
  
  ---
  
  