1 #include "gdna.h"
2 #include <string.h>
3 
4 const char* IUPAC_2BIT  ="AACCTTGGTTAAAAAACCCCGGAAAAAACCAAAAAA";
5 const char* IUPAC_2BITN ="001133223300000011112200000011000000";
6 const char* IUPAC_DEFS  ="AaCcTtGgUuMmRrWwSsYyKkVvHhDdBbNnXx-*";
7 const char* IUPAC_COMP  ="TtGgAaCcAaKkYyWwSsRrMmBbDdHhVvNnXx-*";
8 
9 #define A_2BIT 0 // 00
10 #define C_2BIT 1 // 01
11 #define G_2BIT 2 // 10
12 #define T_2BIT 3 // 11
13 
14 static byte ntCompTable[256];
15 static byte nt2bit[256]; //maps any character to a 2bit base value (with N = A)
16 static char v_2bit2nt[4] = {'A','C','G','T'};
17 
18 //----------------------
19 
20 static bool gdna_Ready=gDnaInit();
21 
22 //----------------------
23 
gdna2bit(char * & nt,int n)24 byte gdna2bit(char* &nt, int n) {
25 // Pack n bases into a byte (n can be 1..4)
26 byte out = 0;
27 while (n && *nt) {
28     n--;
29     out <<= 2;
30     out += nt2bit[(int)*nt];
31     nt++;
32     }
33 #ifdef GDEBUG
34 if (n) {
35      GError("Error: attempt to read 6-mer beyond the end of the string!\n");
36      }
37 #endif
38 return out;
39 }
40 
41 
ntComplement(char c)42 char ntComplement(char c) {
43  return ntCompTable[(int)c];
44  }
45 
g2bit2base(byte v2bit)46 char g2bit2base(byte v2bit) {
47  return v_2bit2nt[v2bit & 0x03 ];
48 }
49 
50 //in place reverse complement of nucleotide (sub)sequence
reverseComplement(char * seq,int slen)51 char* reverseComplement(char* seq, int slen) {
52    if (slen==0) slen=strlen(seq);
53    //reverseChars(seq,len);
54    int l=0;
55    int r=slen-1;
56    register char c;
57    while (l<r) {
58       c=seq[l];seq[l]=seq[r];
59       seq[r]=c;   //this was: Gswap(str[l],str[r]);
60       l++;r--;
61       }
62    for (int i=0;i<slen;i++) seq[i]=ntComplement(seq[i]);
63    return seq;
64  }
65 
gDnaInit()66 bool gDnaInit() {
67        if (gdna_Ready) return true;
68        int l=strlen(IUPAC_DEFS);
69        ntCompTable[0]=0;
70        nt2bit[0]=0;
71        for (int ch=1;ch<256;ch++) {
72           ntCompTable[ch]=0;
73           nt2bit[ch]=0;
74           for (int i=0;i<l;i++)
75                 if (ch==IUPAC_DEFS[i]) {
76                   ntCompTable[ch]=IUPAC_COMP[i];
77                   nt2bit[ch] = IUPAC_2BITN[i]-'0';
78                   break;
79                   }
80           if (ntCompTable[ch]==0) {
81               ntCompTable[ch]='N';
82               }
83           }
84       gdna_Ready=true;
85       return true;
86      }
87 
88 
89 
90 
91