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, 5 months ago) by gpertea
File size: 1967 byte(s)
Log Message:
wip fqtrim

Line User Rev File contents
1 gpertea 2 #include "gdna.h"
2     #include <string.h>
3    
4 gpertea 171 const char* IUPAC_2BIT ="AACCTTGGTTAAAAAACCCCGGAAAAAACCAAAAAA";
5     const char* IUPAC_2BITN ="001133223300000011112200000011000000";
6     const char* IUPAC_DEFS ="AaCcTtGgUuMmRrWwSsYyKkVvHhDdBbNnXx-*";
7     const char* IUPAC_COMP ="TtGgAaCcAaKkYyWwSsRrMmBbDdHhVvNnXx-*";
8 gpertea 2
9 gpertea 171 #define A_2BIT 0 // 00
10     #define C_2BIT 1 // 01
11     #define G_2BIT 2 // 10
12     #define T_2BIT 3 // 11
13 gpertea 2
14 gpertea 171 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 gpertea 2
18 gpertea 171 //----------------------
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 gpertea 2 char ntComplement(char c) {
38     return ntCompTable[(int)c];
39     }
40    
41 gpertea 171 char g2bit2base(byte v2bit) {
42     return v_2bit2nt[v2bit & 0x03 ];
43     }
44    
45 gpertea 2 //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 gpertea 144 seq[r]=c; //this was: Gswap(str[l],str[r]);
55 gpertea 2 l++;r--;
56     }
57     for (int i=0;i<slen;i++) seq[i]=ntComplement(seq[i]);
58     return seq;
59     }
60    
61 gpertea 171 bool gDnaInit() {
62     if (gdna_Ready) return true;
63 gpertea 2 int l=strlen(IUPAC_DEFS);
64     ntCompTable[0]=0;
65 gpertea 171 nt2bit[0]=0;
66 gpertea 2 for (int ch=1;ch<256;ch++) {
67     ntCompTable[ch]=0;
68 gpertea 171 nt2bit[ch]=0;
69 gpertea 2 for (int i=0;i<l;i++)
70 gpertea 171 if (ch==IUPAC_DEFS[i]) {
71     ntCompTable[ch]=IUPAC_COMP[i];
72     nt2bit[ch] = IUPAC_2BITN[i]-'0';
73 gpertea 2 break;
74     }
75 gpertea 171 if (ntCompTable[ch]==0) {
76 gpertea 2 ntCompTable[ch]='N';
77 gpertea 171 }
78 gpertea 2 }
79 gpertea 171 gdna_Ready=true;
80 gpertea 2 return true;
81     }