Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

C++ - Output all possible DNA kmers of a given length

A "kmer" is a DNA sequence of length K. A valid DNA sequence (for my purposes) can only contain the 4 following bases: A,C,T,G. I am looking for a C++ algorithm that simply outputs all possible combinations of these bases in alphabetical order into a string array. For example, if K = 2, The program should generate the following array:

kmers[0]  = AA
kmers[1]  = AC
kmers[2]  = AG
kmers[3]  = AT
kmers[4]  = CA
kmers[5]  = CC
kmers[6]  = CG
kmers[7]  = CT
kmers[8]  = GA
kmers[9]  = GC
kmers[10] = GG
kmers[11] = GT
kmers[12] = TA
kmers[13] = TC
kmers[14] = TG
kmers[15] = TT

If I'm thinking about this correctly, the problem really breaks down to converting a decimal integer to base 4 then substituting the appropriate bases. I thought I could use itoa for this, but itoa is not C standard, and my compiler did not support it. I welcome any clever ideas. Here is my sample code:

#include <iostream>
#include <string>
#include <math.h>

#define K 3

using namespace std;

int main() {

  int num_kmers = pow(4,K);
  string* kmers = NULL;

  /* Allocate memory for kmers array */
  kmers = new string[num_kmers];

  /* Populate kmers array */
  for (int i=0; i< pow(4,K); i++) {

    // POPULATE THE kmers ARRAY HERE                                                                                                                                                         

  }

  /* Display all possible kmers */
  for (int i=0; i< pow(4,K); i++)
    cout << kmers[i] << "\n";

  delete [] kmers;
}
like image 797
nedblorf Avatar asked Apr 06 '11 16:04

nedblorf


1 Answers

You would need to use recursion to be flexible (i.e. so that you could change K easily).

void populate(int depth, string base, string* kmers, int* kmers_offset)
{
    if(depth == K)
    {
        kmers[*kmers_offset].assign(base);
        (*kmers_offset)++;
    }
    else
    {
        static char bases[] = { 'A', 'C', 'G', 'T' };
        for(int i = 0; i < 4; ++i)
            populate(depth + 1, base + bases[i], kmers, kmers_offset);
    }
}

And then call it like this:

int kmers_offset = 0;
populate(0, "", kmers, &kmers_offset);

Cheers.

like image 139
Chris Eberle Avatar answered Oct 01 '22 00:10

Chris Eberle