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)24byte 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)42char ntComplement(char c) { 43 return ntCompTable[(int)c]; 44 } 45 g2bit2base(byte v2bit)46char 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)51char* 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()66bool 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