5 * Created by Pat Schloss on 12/15/08.
6 * Copyright 2008 Patrick D. Schloss. All rights reserved.
17 /**************************************************************************************************/
19 Kmer::Kmer(int size) : kmerSize(size) {
21 int power4s[9] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536 };
22 maxKmer = power4s[kmerSize]+1;// (int)pow(4.,k)+1;
26 /**************************************************************************************************/
28 string Kmer::getKmerString(string sequence){
29 int length = sequence.length();
30 int nKmers = length - kmerSize + 1;
31 vector<int> counts(maxKmer, 0);
33 for(int i=0;i<nKmers;i++){
34 int kmerNumber = getKmerNumber(sequence, i);
38 string kmerString = "";
39 for(int i=0;i<maxKmer;i++){
40 kmerString += getASCII(counts[i]);
46 /**************************************************************************************************/
48 int Kmer::getKmerNumber(string sequence, int index){
49 int power4s[9] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536 };
51 for(int i=0;i<kmerSize;i++){
52 if(toupper(sequence[i+index]) == 'A') { kmer += (0 * power4s[kmerSize-i-1]); }
53 else if(toupper(sequence[i+index]) == 'C') { kmer += (1 * power4s[kmerSize-i-1]); }
54 else if(toupper(sequence[i+index]) == 'G') { kmer += (2 * power4s[kmerSize-i-1]); }
55 else if(toupper(sequence[i+index]) == 'U') { kmer += (3 * power4s[kmerSize-i-1]); }
56 else if(toupper(sequence[i+index]) == 'T') { kmer += (3 * power4s[kmerSize-i-1]); }
57 else if(toupper(sequence[i+index]) == 'N') { return (int)power4s[kmerSize]; }
62 /**************************************************************************************************/
64 string Kmer::getKmerBases(int kmerNumber){
65 int power4s[9] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536 };
69 if(kmerNumber == power4s[kmerSize]){//pow(4.,7)){
70 for(int i=0;i<kmerSize;i++){
75 for(int i=0;i<kmerSize;i++){
76 int nt = (int)(kmerNumber / (float)power4s[i]) % 4;
77 if(nt == 0) { kmer = 'A' + kmer; }
78 else if(nt == 1){ kmer = 'C' + kmer; }
79 else if(nt == 2){ kmer = 'G' + kmer; }
80 else if(nt == 3){ kmer = 'T' + kmer; }
86 /**************************************************************************************************/
88 char Kmer::getASCII(int number) { return (char)(33+number); }
90 /**************************************************************************************************/
92 int Kmer::getNumber(char character) { return ((int)(character-'!')); }
94 /**************************************************************************************************/