<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>

CC BY-SA 4.0
2018-12-21T03:04:57.140Z
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.
conchoecia
CC BY-SA 4.0
2018-12-21T02:50:45.937Z
2018-12-21T02:59:45.817Z
2018-12-21T02:59:45.817Z
nomad
<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>

CC BY-SA 4.0
2018-12-20T19:15:03.650Z
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!
conchoecia
<fastq><c++><multithreading>
Multithread fastq processing with kseq.h in C++11?
# 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 as my fastq/fasta parser for a C++ project. My actual implementation uses something like this, 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:
//compiled with g++ 7.3 using command g++ this_prog.cpp -o this_prog -lz
// execute with ./this_prog my_file.fastq.gz
#include //cout
#include 
#include 
#include 
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char **argv)
{
  //initialize prefix
  std::map 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.  

---

