Answers.Body: <p>You might take a look at my FQFeeder library (<a href="https://github.com/rob-p/FQFeeder" rel="noreferrer">https://github.com/rob-p/FQFeeder</a>), which pairs the kseq parser with a multi-producer / multi-consumer queue to allow efficient parallel processing of the parsed reads.  The library contains an example program to demonstrate how to use it with multiple threads.  It's worth noting that kseq is already <em>very</em> fast, so multi-threading only really makes sense if worker threads are doing non-trivial work on each read so that you can then make the parsing and processing concurrent (have the processing keep pace with kseq).</p>
  
  <p>I'll note that my library relies on a header-only MPMC implementation and so technically isn't "only using the standard library", but you're already relying on <code>kseq.h</code>, so omitting third-party libraries must not be a hard constraint.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-21T03:04:57.140Z
Answers.Text: Thanks for posting this - I will take a look. Yes, the work on each pair is nontrivial so I've been wondering how to best store reads in chunks/process them or have several threads of kseq.h read different portions of a single file.
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-21T02:50:45.937Z
Answers.LastActivityDate: 2018-12-21T02:59:45.817Z
Answers.LastEditDate: 2018-12-21T02:59:45.817Z
Answers.UserDisplayName: nomad
Body: <h2>Background</h2>
  
  <p>I am using the wonderful <a href="http://lh3lh3.users.sourceforge.net/parsefastq.shtml" rel="noreferrer"><code>kseq.h</code></a> as my fastq/fasta parser for a C++ project. <a href="https://bioinformatics.stackexchange.com/a/5471/2085">My actual implementation uses something like this</a>, where I am counting kmers of each sequence and performing a bunch of other processing on every read pair.</p>
  
  <p>Thanks to @user1782818's code my program is very fast, but I want it to be faster. The only way I know how to do this is with multithreading. However, I have never worked with multithreading in C++.</p>
  
  <p>For each read pair I am accessing a vector or map of structs. Within each struct I am incrementing counters or pushing values onto vectors. Essentially I am doing something like this:</p>
  
  <pre class="lang-cpp prettyprint-override"><code>//compiled with g++ 7.3 using command g++ this_prog.cpp -o this_prog -lz
  // execute with ./this_prog my_file.fastq.gz
  #include &lt;iostream&gt;//cout
  #include &lt;map&gt;
  #include &lt;zlib.h&gt;
  #include &lt;stdio.h&gt;
  #include "kseq.h"
  KSEQ_INIT(gzFile, gzread)
  
  int main(int argc, char **argv)
  {
    //initialize prefix
    std::map&lt;std::string, int&gt; prefix_counter;
    std::string prefix;
    prefix.resize(2);
  
    gzFile fp;
    kseq_t *seq;
    int n = 0, slen = 0, qlen = 0;
    fp = gzopen(argv[1], "r");
    seq = kseq_init(fp);
    while (kseq_read(seq) &gt;= 0){
      ++n;
      //COUNT PREFIXES
      std::copy(seq-&gt;seq.s, seq-&gt;seq.s + 2, begin(prefix));
      if (prefix_counter.count(prefix) == 0){
        prefix_counter[prefix] = 1;
      } else prefix_counter[prefix]++;
    }
    std::cout &lt;&lt; "Total n. of sequences: " &lt;&lt; n &lt;&lt; '\n';
    kseq_destroy(seq);
    gzclose(fp);
  
    //print counts
    for( const auto&amp; sm_pair : prefix_counter ){
      std::cout &lt;&lt; sm_pair.first &lt;&lt; ": " &lt;&lt; sm_pair.second &lt;&lt; '\n';
    }
    return 0;
  }
  </code></pre>
  
  <p>Example output from the above program is something like this for a 100,000-read test fastq file.</p>
  
  <pre><code>Total n. of sequences: 100000
  AA: 8561
  AC: 5178
  AG: 7634
  AT: 8523
  CA: 7526
  </code></pre>
  
  <h2>Question</h2>
  
  <p><strong>Q1:</strong> Is it possible to perform a task like the one shown above using multithreading and <code>kseq.h</code> to read in a fastq file and process the reads in chunks? I only want to use standard libraries. C++11 is OK - I know that multithreading support was added in C++11.</p>
  
  <p><strong>Q2:</strong> Can you please provide a C++ example of how to convert the above or similar program to use multithreading using <code>kseq.h</code>?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-12-20T19:15:03.650Z
Comments.Text: Note: I'm not posting this on SO to avoid "read the manual"-type responses and I believe that others in the bioinformatics SE community will benefit from an answer to this question. Thanks!
Comments.UserDisplayName: conchoecia
Tags: <fastq><c++><multithreading>
Title: Multithread fastq processing with kseq.h in C++11?
ThreadText: # Multithread fastq processing with kseq.h in C++11?
  
  (Asked by: username_0 on 2018-12-20)
  
  ## Background
  
  I am using the wonderful [`kseq.h`](http://lh3lh3.users.sourceforge.net/parsefastq.shtml) as my fastq/fasta parser for a C++ project. [My actual implementation uses something like this](https://bioinformatics.stackexchange.com/a/5471/2085), where I am counting kmers of each sequence and performing a bunch of other processing on every read pair.
  
  Thanks to @user1782818's code my program is very fast, but I want it to be faster. The only way I know how to do this is with multithreading. However, I have never worked with multithreading in C++.
  
  For each read pair I am accessing a vector or map of structs. Within each struct I am incrementing counters or pushing values onto vectors. Essentially I am doing something like this:
  
  ```cpp
  //compiled with g++ 7.3 using command g++ this_prog.cpp -o this_prog -lz
  // execute with ./this_prog my_file.fastq.gz
  #include <iostream>//cout
  #include <map>
  #include <zlib.h>
  #include <stdio.h>
  #include "kseq.h"
  KSEQ_INIT(gzFile, gzread)
  
  int main(int argc, char **argv)
  {
    //initialize prefix
    std::map<std::string, int> prefix_counter;
    std::string prefix;
    prefix.resize(2);
  
    gzFile fp;
    kseq_t *seq;
    int n = 0, slen = 0, qlen = 0;
    fp = gzopen(argv[1], "r");
    seq = kseq_init(fp);
    while (kseq_read(seq) >= 0){
      ++n;
      //COUNT PREFIXES
      std::copy(seq->seq.s, seq->seq.s + 2, begin(prefix));
      if (prefix_counter.count(prefix) == 0){
        prefix_counter[prefix] = 1;
      } else prefix_counter[prefix]++;
    }
    std::cout << "Total n. of sequences: " << n << '\n';
    kseq_destroy(seq);
    gzclose(fp);
  
    //print counts
    for( const auto& sm_pair : prefix_counter ){
      std::cout << sm_pair.first << ": " << sm_pair.second << '\n';
    }
    return 0;
  }
  
  ```
  
  Example output from the above program is something like this for a 100,000-read test fastq file.
  
  ```
  Total n. of sequences: 100000
  AA: 8561
  AC: 5178
  AG: 7634
  AT: 8523
  CA: 7526
  
  ```
  
  ## Question
  
  **Q1:** Is it possible to perform a task like the one shown above using multithreading and `kseq.h` to read in a fastq file and process the reads in chunks? I only want to use standard libraries. C++11 is OK - I know that multithreading support was added in C++11.
  
  **Q2:** Can you please provide a C++ example of how to convert the above or similar program to use multithreading using `kseq.h`?
  
  **Tags:** fastq, c++, multithreading
  
  > **Comments:**
  > **username_0**: Note: I'm not posting this on SO to avoid "read the manual"-type responses and I believe that others in the bioinformatics SE community will benefit from an answer to this question. Thanks!  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2018-12-21. Score: 7 | Accepted answer
  
  You might take a look at my FQFeeder library (<https://github.com/rob-p/FQFeeder>), which pairs the kseq parser with a multi-producer / multi-consumer queue to allow efficient parallel processing of the parsed reads. The library contains an example program to demonstrate how to use it with multiple threads. It's worth noting that kseq is already *very* fast, so multi-threading only really makes sense if worker threads are doing non-trivial work on each read so that you can then make the parsing and processing concurrent (have the processing keep pace with kseq).
  
  I'll note that my library relies on a header-only MPMC implementation and so technically isn't "only using the standard library", but you're already relying on `kseq.h`, so omitting third-party libraries must not be a hard constraint.
  
  > **Comments:**
  > **username_0**: Thanks for posting this - I will take a look. Yes, the work on each pair is nontrivial so I've been wondering how to best store reads in chunks/process them or have several threads of kseq.h read different portions of a single file.  
  
  ---
  
  