ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gclib/gdna.cpp
Revision: 171
Committed: Tue Feb 14 22:36:26 2012 UTC (7 years, 8 months ago) by gpertea
File size: 1967 byte(s)
Log Message:
wip fqtrim

Line File contents
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
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 return out;
34 }
35
36
37 char ntComplement(char c) {
38 return ntCompTable[(int)c];
39 }
40
41 char g2bit2base(byte v2bit) {
42 return v_2bit2nt[v2bit & 0x03 ];
43 }
44
45 //in place reverse complement of nucleotide (sub)sequence
46 char* reverseComplement(char* seq, int slen) {
47 if (slen==0) slen=strlen(seq);
48 //reverseChars(seq,len);
49 int l=0;
50 int r=slen-1;
51 register char c;
52 while (l<r) {
53 c=seq[l];seq[l]=seq[r];
54 seq[r]=c; //this was: Gswap(str[l],str[r]);
55 l++;r--;
56 }
57 for (int i=0;i<slen;i++) seq[i]=ntComplement(seq[i]);
58 return seq;
59 }
60
61 bool gDnaInit() {
62 if (gdna_Ready) return true;
63 int l=strlen(IUPAC_DEFS);
64 ntCompTable[0]=0;
65 nt2bit[0]=0;
66 for (int ch=1;ch<256;ch++) {
67 ntCompTable[ch]=0;
68 nt2bit[ch]=0;
69 for (int i=0;i<l;i++)
70 if (ch==IUPAC_DEFS[i]) {
71 ntCompTable[ch]=IUPAC_COMP[i];
72 nt2bit[ch] = IUPAC_2BITN[i]-'0';
73 break;
74 }
75 if (ntCompTable[ch]==0) {
76 ntCompTable[ch]='N';
77 }
78 }
79 gdna_Ready=true;
80 return true;
81 }