kmer.cpp

   1 /*
   2  *  kmer.cpp
   3  *
   4  *
   5  *  Created by Pat Schloss on 12/15/08.
   6  *  Copyright 2008 Patrick D. Schloss. All rights reserved.
   7  *
   8  */
   9
  10 #include "kmer.hpp"
  11
  12 /**************************************************************************************************/
  13
  14 Kmer::Kmer(int size) : kmerSize(size) { //      The constructor sets the size of the kmer
  15
  16         int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 };
  17                                                                                                                                                 //      No reason to waste the time of recalculating
  18         maxKmer = power4s[kmerSize]+1;// (int)pow(4.,k)+1;                                      //      powers of 4 everytime through.  We need an
  19                                                                                                                                                 //      extra kmer if we get a non-ATGCU base
  20 }
  21
  22 /**************************************************************************************************/
  23
  24 string Kmer::getKmerString(string sequence){    //      Calculate kmer for each position in the sequence, count the freq
  25         int length = sequence.length();                         //      of each kmer, and convert it to an ascii character with base '!'.
  26         int nKmers = length - kmerSize + 1;                     //      Export the string of characters as a string
  27         vector<int> counts(maxKmer, 0);
  28
  29         for(int i=0;i<nKmers;i++){                                      //      Go though sequence and get the number between 0 and maxKmer for that
  30                 int kmerNumber = getKmerNumber(sequence, i);//  kmer.  Increase the count for the kmer in the counts vector
  31                 counts[kmerNumber]++;
  32         }
  33
  34         string kmerString = "";
  35         for(int i=0;i<maxKmer;i++){                                     //      Scan through the vector of counts and convert each element of the
  36                 kmerString += getASCII(counts[i]);              //      vector to an ascii character that is equal to or greater than '!'
  37         }
  38
  39         return kmerString;
  40 }
  41 /**************************************************************************************************/
  42
  43 vector< map<int, int> > Kmer::getKmerCounts(string sequence){   //      Calculate kmer for each position in the sequence, save info in a map
  44         int length = sequence.length();                         //      so you know at each spot in the sequence what kmers were found
  45         int nKmers = length - kmerSize + 1;                     //
  46         vector< map<int, int> > counts; counts.resize(nKmers);  // a map kmer counts for each spot
  47          map<int, int>::iterator it;
  48
  49         for(int i=0;i<nKmers;i++){                                      //      Go though sequence and get the number between 0 and maxKmer for that
  50                 if (i == 0) {
  51                         int kmerNumber = getKmerNumber(sequence, i);//  kmer.
  52                         counts[i][kmerNumber] = 1;                      // add this kmer if not already there
  53                 }else {
  54                         //your count is everything that came before and whatever you find now
  55                         counts[i] = counts[i-1];
  56                         int kmerNumber = getKmerNumber(sequence, i);//  kmer.
  57
  58                         it = counts[i].find(kmerNumber);
  59                         if (it!= counts[i].end()) {   counts[i][kmerNumber]++;  }  //increment number of times you have seen this kmer
  60                         else {   counts[i][kmerNumber] = 1;   }  // add this kmer since not already there
  61                 }
  62         }
  63
  64         return counts;
  65 }
  66
  67
  68 /**************************************************************************************************/
  69
  70 int Kmer::getKmerNumber(string sequence, int index){
  71
  72 //      Here we convert a kmer to a number between 0 and maxKmer.  For example, AAAA would equal 0 and TTTT would equal 255.
  73 //      If there's an N in the kmer, it is set to 256 (if we are looking at 4mers).  The largest we can look at are 8mers,
  74 //      this could be easily increased.
  75 //
  76 //      Example:   ATGCAT (kSize = 6)
  77 //                i:   012345
  78 //
  79 //              Score   = 0*4^(6-0-1) + 3*4^(6-1-1) + 2*4^(6-2-1) + 1*4^(6-3-1) + 0*4^(6-4-1) + 3*4^(6-5-1)
  80 //                              = 0*4^5 +       3*4^4   +       2*4^3   +       1*4^2   +       0*4^1   +       3*4^0
  81 //                              = 0 + 3*256 + 2*64 + 1*16 + 0*4 + 3*1
  82 //                              = 0 + 768 + 128 + 16 + 0 + 3
  83 //                              = 915
  84
  85         int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 };
  86
  87         int kmer = 0;
  88
  89         for(int i=0;i<kmerSize;i++){
  90                 if(toupper(sequence[i+index]) == 'A')           {       kmer += (0 * power4s[kmerSize-i-1]);    }
  91                 else if(toupper(sequence[i+index]) == 'C')      {       kmer += (1 * power4s[kmerSize-i-1]);    }
  92                 else if(toupper(sequence[i+index]) == 'G')      {       kmer += (2 * power4s[kmerSize-i-1]);    }
  93                 else if(toupper(sequence[i+index]) == 'U')      {       kmer += (3 * power4s[kmerSize-i-1]);    }
  94                 else if(toupper(sequence[i+index]) == 'T')      {       kmer += (3 * power4s[kmerSize-i-1]);    }
  95                 else if(toupper(sequence[i+index]) == 'N')      {       return (int)power4s[kmerSize];                  }
  96         }
  97         return kmer;
  98 }
  99
 100 /**************************************************************************************************/
 101
 102 string Kmer::getKmerBases(int kmerNumber){
 103
 104 //      Here we convert the kmer number into the kmer in terms of bases.
 105 //
 106 //      Example:        Score = 915 (for a 6-mer)
 107 //                              Base6 = (915 / 4^0) % 4 = 915 % 4 = 3 => T      [T]
 108 //                              Base5 = (915 / 4^1) % 4 = 228 % 4 = 0 => A      [AT]
 109 //                              Base4 = (915 / 4^2) % 4 = 57 % 4 = 1 => C       [CAT]
 110 //                              Base3 = (915 / 4^3) % 4 = 14 % 4 = 2 => G       [GCAT]
 111 //                              Base2 = (915 / 4^4) % 4 = 3 % 4 = 3 => T        [TGCAT]
 112 //                              Base1 = (915 / 4^5) % 4 = 0 % 4 = 0 => A        [ATGCAT] -> this checks out with the previous method
 113
 114         int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 };
 115
 116         string kmer = "";
 117
 118         if(kmerNumber == power4s[kmerSize]){//pow(4.,7)){       //      if the kmer number is the same as the maxKmer then it must
 119                 for(int i=0;i<kmerSize;i++){                                    //      have had an N in it and so we'll just call it N x kmerSize
 120                         kmer += 'N';
 121                 }
 122         }
 123         else{
 124                 for(int i=0;i<kmerSize;i++){
 125                         int nt = (int)(kmerNumber / (float)power4s[i]) % 4;             //      the '%' operator returns the remainder
 126                         if(nt == 0)             {       kmer = 'A' + kmer;      }                               //      from int-based division ]
 127                         else if(nt == 1){       kmer = 'C' + kmer;      }
 128                         else if(nt == 2){       kmer = 'G' + kmer;      }
 129                         else if(nt == 3){       kmer = 'T' + kmer;      }
 130                 }
 131         }
 132         return kmer;
 133 }
 134 /**************************************************************************************************/
 135
 136 int Kmer::getReverseKmerNumber(int kmerNumber){
 137
 138         string kmerString = getKmerBases(kmerNumber);
 139
 140         //get Reverse
 141         string reverse = "";
 142         for(int i=kmerString.length()-1;i>=0;i--){
 143                 if(kmerString[i] == 'A')                {       reverse += 'T'; }
 144                 else if(kmerString[i] == 'T'){  reverse += 'A'; }
 145                 else if(kmerString[i] == 'G'){  reverse += 'C'; }
 146                 else if(kmerString[i] == 'C'){  reverse += 'G'; }
 147                 else                                            {       reverse += 'N'; }
 148         }
 149
 150         int reverseNumber = getKmerNumber(reverse, 0);
 151
 152         return reverseNumber;
 153
 154 }
 155
 156 /**************************************************************************************************/
 157
 158 char Kmer::getASCII(int number)         {       return (char)(33+number);                       }       // '!' is the first printable char and
 159                                                                                                                                                                 // has the int value of 33
 160 /**************************************************************************************************/
 161
 162 int Kmer::getNumber(char character)     {       return ((int)(character-'!'));          }       // '!' has the value of 33
 163
 164 /**************************************************************************************************/