ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gclib/gdna.cpp
(Generate patch)
# Line 1 | Line 1
1   #include "gdna.h"
2   #include <string.h>
3  
4 < #define IUPAC_DEFS "AaCcTtGgUuMmRrWwSsYyKkVvHhDdBbNnXx-*"
5 < #define IUPAC_COMP "TtGgAaCcAaKkYyWwSsRrMmBbDdHhVvNnXx-*"
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  
7 unsigned char ntCompTable[256];
8
9 static bool gdna_ntCompTableReady=ntCompTableInit();
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);
# Line 21 | Line 51
51     register char c;
52     while (l<r) {
53        c=seq[l];seq[l]=seq[r];
54 <      seq[r]=c;   //this was: swap(str[l],str[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 ntCompTableInit() {
62 <       //if (gdna_ntCompTableReady) return true;
33 <       char n[]=IUPAC_DEFS;
34 <       char c[]=IUPAC_COMP;
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==n[i]) {
71 <                  ntCompTable[ch]=c[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)
75 >          if (ntCompTable[ch]==0) {
76                ntCompTable[ch]='N';
77 +              }
78            }
79 <      //gdna_ntCompTableReady=true;
79 >      gdna_Ready=true;
80        return true;
81       }
82  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines