title: Fast counting of nucleotide types in a large number of sequences
body: <p>First, a bit of background about my question.<br>
  I work as a bioinformatician, which means that I do informatics treatment to try to answer a biological question. In my problem, I have to manipulate a file called a FASTA file which looks like this :</p>
  
  <pre><code>&gt;Header 1  
  ATGACTGATCGNTGACTGACTGTAGCTAGC  
  &gt;Header 2  
  ATGCATGCTAGCTGACTGATCGTAGCTAGC  
  ATCGATCGTAGCT
  </code></pre>
  
  <p>So a FASTA file is basically just a header, preceded by a '>' character, then a sequence on one or multiple lines that is composed of nucleotides. Nucleotides are characters that can take 5 possible values : A, T, C, G or N.</p>
  
  <p>The thing I would like to do is count the number of times each nucleotide type appears so if we consider this dummy FASTA file :</p>
  
  <pre><code>&gt;Header 1  
  ATTCGN
  </code></pre>
  
  <p>I should have, as a result :<br>
  <code>A:1 T:2 C:1 G:1 N:1</code></p>
  
  <p>Here is what I got so far :</p>
  
  <pre><code>ifstream sequence_file(input_file.c_str());
  string line;
  string sequence = "";
  map&lt;char, double&gt; nucleotide_counts;
  
  while(getline(sequence_file, line)) {
      if(line[0] != '&gt;') {
          sequence += line;
      }
      else {
          nucleotide_counts['A'] = boost::count(sequence, 'A');
          nucleotide_counts['T'] = boost::count(sequence, 'T');
          nucleotide_counts['C'] = boost::count(sequence, 'C');
          nucleotide_counts['G'] = boost::count(sequence, 'G');
          nucleotide_counts['N'] = boost::count(sequence, 'N');
          sequence = "";
      }
  }
  </code></pre>
  
  <p>So it reads the file line by line, if it encounters a '>' as the first character of the line, it knows that the sequence is complete and starts to count. Now the problem I'm facing is that I have millions of sequences with several billions of nucleotides in total. I can see that my method is not optimized because I call <code>boost::count</code> five times on the same sequence.</p>
  
  <p>Other things I have tried :</p>
  
  <ul>
  <li>Parsing the sequence to increment a counter for each nucleotide types. I tried using a <code>map&lt;char, double&gt;</code> to map each nucleotide to a value but this was slower than the boost solution.</li>
  <li>Using the <code>std::count</code> of the algorithm library but this was too slow too.</li>
  </ul>
  
  <p>I searched the internet for solutions but every solution I found was good if the number of sequences was low, which is not my case. Would you have any idea that could help me speed things up ?</p>
  
  <p><strong>EDIT 1</strong> :
  I also tried this version but it was 2 times slower than the boost one :</p>
  
  <pre><code>ifstream sequence_file(input_file.c_str());
  string line;
  string sequence = "";
  map&lt;char, double&gt; nucleotide_counts;
  
  while(getline(sequence_file, line)) {
      if(line[0] != '&gt;') {
          sequence += line;
      }
      else {
          for(int i = 0; i &lt; sequence.size(); i++) {
             nucleotide_counts[sequence[i]]++;
          }
          sequence = "";
      }
  }
  </code></pre>
  
  <p><strong>EDIT 2</strong> : Thanks to everyone in this thread, I was able to obtain a speed up of about 30 times compared to the boost original solution. Here is the code :</p>
  
  <pre><code>#include &lt;map&gt; // std::array
  #include &lt;fstream&gt; // std::ifstream
  #include &lt;string&gt; // std::string  
  
  void count_nucleotides(std::array&lt;double, 26&gt; &amp;nucleotide_counts, std::string sequence) {
      for(unsigned int i = 0; i &lt; sequence.size(); i++) {
          ++nucleotide_counts[sequence[i] - 'A'];
      }
  }  
  
  std::ifstream sequence_file(input_file.c_str());
  std::string line;
  std::string sequence = "";
  std::array&lt;double, 26&gt; nucleotide_counts;
  
  while(getline(sequence_file, line)) {
      if(line[0] != '&gt;') {
          sequence += line;
      }
      else {
          count_nucleotides(nucleotide_counts, sequence);
          sequence = "";
      }
  }
  </code></pre>