root/branches/stable_5.0/AWT/AWT_pro_a_nucs.cxx

Revision 6033, 15.1 KB (checked in by westram, 3 years ago)
  • DB browser now is only included in DEBUG mode
  • include wrappers for ad_spec.hxx + ad_trees.hxx
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1#include <stdio.h>
2#include <stdlib.h>
3#include <ctype.h>
4// #include <malloc.h>
5#include <string.h>
6
7#include <arbdb.h>
8#include <arbdbt.h>
9#include <awt_pro_a_nucs.hxx>
10#include <awt_codon_table.hxx>
11
12#include <inline.h>
13
14char *AP_create_dna_to_ap_bases(){
15    int i;
16    AP_BASES    val;
17    char *table = new char[256];
18    for (i=0;i<256;i++) {
19        switch ((char)i) {
20            case 'a':
21            case 'A':
22                val = AP_A;
23                break;
24            case 'g':
25            case 'G':
26                val = AP_G;
27                break;
28            case 'c':
29            case 'C':
30                val = AP_C;
31                break;
32            case 't':
33            case 'T':
34            case 'u':
35            case 'U':
36                val = AP_T;
37                break;
38            case 'n':
39            case 'N':
40                val = (AP_BASES)(AP_A + AP_G + AP_C + AP_T);
41                break;
42            case '?':
43            case '.':
44                val = (AP_BASES)(AP_A + AP_G + AP_C + AP_T + AP_S);
45                break;
46            case '-':
47                val = AP_S;
48                break;
49            case 'm':
50            case 'M':
51                val = (AP_BASES)(AP_A+AP_C);
52                break;
53            case 'r':
54            case 'R':
55                val = (AP_BASES)(AP_A+AP_G);
56                break;
57            case 'w':
58            case 'W':
59                val = (AP_BASES)(AP_A+AP_T);
60                break;
61            case 's':
62            case 'S':
63                val = (AP_BASES)(AP_C+AP_G);
64                break;
65            case 'y':
66            case 'Y':
67                val = (AP_BASES)(AP_C+AP_T);
68                break;
69            case 'k':
70            case 'K':
71                val = (AP_BASES)(AP_G+AP_T);
72                break;
73            case 'v':
74            case 'V':
75                val = (AP_BASES)(AP_A+AP_C+AP_G);
76                break;
77            case 'h':
78            case 'H':
79                val = (AP_BASES)(AP_A+AP_C+AP_T);
80                break;
81            case 'd':
82            case 'D':
83                val = (AP_BASES)(AP_A+AP_G+AP_T);
84                break;
85            case 'b':
86            case 'B':
87                val = (AP_BASES)(AP_C+AP_G+AP_T);
88                break;
89            default:
90                val = AP_N;
91                break;
92        }
93        table[i] = (char)val;
94    } /*for*/
95    return table;
96}
97
98
99long *AWT_translator::create_pro_to_bits() const {
100    int i;
101    int j;
102    long *table = (long *)GB_calloc(sizeof(long),256);
103    for ( i = 0 ; i < max_aa; i++) {
104        j = index_2_spro[i];
105        if (j== '.')  {
106            table[i] = -1;
107            continue;
108        }
109        j = s2str[j]->index;
110        table[i] = 1<<j;
111    }
112    return table;
113}
114
115void AWT_translator::build_table(unsigned char pbase, const char *tri_pro, const char *nuc) {
116    struct arb_r2a_pro_2_nuc  *str = s2str[pbase];
117   
118    // search for existing protein, else generate new
119    if (!str) {
120        str                               = new arb_r2a_pro_2_nuc;
121        s2str[pbase]          = str;
122        s2str[tolower(pbase)] = str;
123
124        str->index      = max_aa++;
125        str->single_pro = pbase;
126        str->tri_pro[0] = tri_pro[0];
127        str->tri_pro[1] = tri_pro[1];
128        str->tri_pro[2] = tri_pro[2];
129       
130        index_2_spro[str->index] = pbase;
131    }
132    // fast hash table
133    GBS_write_hash(t2i_hash, nuc, pbase);
134
135    int n0 = nuc_2_bitset[safeCharIndex(nuc[0])];
136    int n1 = nuc_2_bitset[safeCharIndex(nuc[1])];
137    int n2 = nuc_2_bitset[safeCharIndex(nuc[2])];
138
139    struct arb_r2a_pro_2_nucs *nucs;
140    for (nucs = str->nucs; nucs; nucs = nucs->next) {
141        if ((!(nucs->nucbits[0] & ~n0)) &&      // search superset
142            (!(nucs->nucbits[1] & ~n1)) &&
143            (!(nucs->nucbits[2] & ~n2)) ) break;
144       
145        int c = 0;
146        if ( nucs->nucbits[0] != n0 ) c++;
147        if ( nucs->nucbits[1] != n1 ) c++;
148        if ( nucs->nucbits[2] != n2 ) c++;
149        if (c <= 1) break;
150    }
151    if (!nucs) {
152        //         nucs       = (struct arb_r2a_pro_2_nucs *)GB_calloc(sizeof(struct arb_r2a_pro_2_nucs),1);
153        nucs       = new arb_r2a_pro_2_nucs;
154        nucs->next = str->nucs;
155        str->nucs  = nucs;
156    }
157    nucs->nucbits[0] |= n0;
158    nucs->nucbits[1] |= n1;
159    nucs->nucbits[2] |= n2;
160}
161
162// ---------------------------
163//      arb_r2a_pro_2_nucs
164// ---------------------------
165
166arb_r2a_pro_2_nucs::arb_r2a_pro_2_nucs()
167    : next(0)
168{
169    memset(nucbits, 0, sizeof(nucbits));
170}
171arb_r2a_pro_2_nucs::~arb_r2a_pro_2_nucs() {
172    delete next;
173}
174
175// --------------------------
176//      arb_r2a_pro_2_nuc
177// --------------------------
178
179arb_r2a_pro_2_nuc::arb_r2a_pro_2_nuc()
180    : single_pro(0)
181    , index(0)
182    , nucs(0)
183{
184    memset(tri_pro, 0, sizeof(tri_pro));
185}
186arb_r2a_pro_2_nuc::~arb_r2a_pro_2_nuc() {
187    delete nucs;
188}
189
190// -----------------------
191//      AWT_translator
192// -----------------------
193
194static int codon_defined_in(const char *codon, const char *codons) {
195    for (int off=0; codons[off]; off+=3) {
196        if (codon[0]==codons[off] && codon[1]==codons[off+1] && codon[2]==codons[off+2]) {
197            return 1;
198        }
199    }
200    return 0;
201}
202
203// order of tables generated with build_table() by AWT_translator ctor is important:
204// must be compatible with DIST/PH_protdist.cxx !!
205// except that this table has an 's' insertion !!!
206
207AWT_translator::AWT_translator(int arb_protein_code_nr) :
208    distance_meter(0), 
209    code_nr(arb_protein_code_nr), 
210    pro_2_bitset(0), 
211    realmax_aa(0), 
212    max_aa(0)
213{
214    memset(s2str, 0, sizeof(s2str));
215    memset(index_2_spro, 0, sizeof(index_2_spro));
216
217    nuc_2_bitset = AP_create_dna_to_ap_bases();
218    t2i_hash     = GBS_create_hash(1024, GB_IGNORE_CASE); // case insensitive
219
220    AWT_initialize_codon_tables();
221
222    {
223        char *D_codons = strdup(AWT_get_codons('D', code_nr));
224        char *N_codons = strdup(AWT_get_codons('N', code_nr));
225        char *E_codons = strdup(AWT_get_codons('E', code_nr));
226        char *Q_codons = strdup(AWT_get_codons('Q', code_nr));
227
228        char protein;
229        for (protein='*'; protein<='Z'; protein = (protein=='*' ? 'A' : protein+1)) {
230            if (protein!='J' && protein!='O' && protein!='U') { // JOU are no aminos
231                const char *codons;
232                if (protein=='D')   codons = D_codons;
233                else if (protein=='N')  codons = N_codons;
234                else if (protein=='E')  codons = E_codons;
235                else if (protein=='Q')  codons = Q_codons;
236                else            codons = AWT_get_codons(protein, code_nr);
237                // codons now contains a 0-terminated-string containing all possible codons for protein
238
239                const char *protein_name = AWT_get_protein_name(protein);
240
241                for (int off=0; codons[off]; off+=3) {
242                    char codon[4];
243                    memcpy(codon, codons+off, 3);
244                    codon[3] = 0;
245
246                    if (protein=='B') {
247                        if (!codon_defined_in(codon, D_codons) && !codon_defined_in(codon, N_codons)) {
248                            build_table(protein, protein_name, codon);
249                        }
250                    }
251                    else if (protein=='Z') {
252                        if (!codon_defined_in(codon, E_codons) && !codon_defined_in(codon, Q_codons)) {
253                            build_table(protein, protein_name, codon);
254                        }
255                    }
256                    else {
257                        build_table(protein, protein_name, codon);
258                    }
259                }
260            }
261        }
262
263        free(Q_codons);
264        free(E_codons);
265        free(N_codons);
266        free(D_codons);
267    }
268
269    realmax_aa = max_aa;
270
271    build_table('-', "---", "---");
272    build_table('.', "...", "...");
273    build_table('.', "???", "???");
274    build_table('X', "NNN", "NNN");
275
276    pro_2_bitset = create_pro_to_bits();
277}
278
279AWT_translator::~AWT_translator() {
280    free(pro_2_bitset);
281
282    delete [] nuc_2_bitset;
283    GBS_free_hash(t2i_hash);
284    for (int i=0;i<256;i++){
285        if (i!= tolower(i)) continue; // do not delete duplicated entries (a-z == A-Z!)
286        delete s2str[i];
287    }
288
289    delete distance_meter;
290}
291
292const AWT_distance_meter *AWT_translator::getDistanceMeter() const {
293    if (!distance_meter) distance_meter = new AWT_distance_meter(this);
294    return distance_meter;
295}
296
297
298// ---------------------------
299//      Distance functions
300// ---------------------------
301
302static int nuc_dist(const AWT_translator *translator, unsigned char p1, unsigned char p2) {  // nucleotid changes between two aa
303    const struct arb_r2a_pro_2_nuc *s1,*s2;
304    s1 = translator->S2str(p1);
305    s2 = translator->S2str(p2);
306    if ( (!s1) || (!s2) ) return -1;
307    struct arb_r2a_pro_2_nucs *n1,*n2;
308    long mindist = 3;
309    // Check all combinations, if any combination is valid -> zero distance
310    for (   n1 = s1->nucs; n1; n1=n1->next){
311        for (   n2 = s2->nucs; n2; n2=n2->next){
312            int dist = 0;
313            int i;
314            for (i=0;i<3;i++) {
315                if ( n1->nucbits[i] & n2->nucbits[i] ) continue;
316                dist++;
317            }
318            if (dist< mindist) mindist = dist;
319        }
320    }
321    return mindist;
322}
323
324#if defined(DEBUG)
325static void awt_pro_a_nucs_debug(const AWT_translator *translator, const AWT_distance_meter *distmeter) {
326    int max_aa     = translator->MaxAA();
327    int realmax_aa = translator->RealmaxAA();
328
329    for (int s = 0; s< max_aa; s++){
330        const AWT_PDP *dist = distmeter->getDistance(s);
331
332        // check bits should not be present in distpad
333        if (s<realmax_aa) {
334            for (int i=0; i<2; i++) {
335                // assertion: bit is set in patd[i] -> bit is clear in patd[i+1]
336                assert_or_exit((dist->patd[i] & ~dist->patd[i+1]) == 0);
337            }
338        }
339        printf("Base %c[%i]: Dist to ", translator->Index2Spro(s), s);
340        for (int d = 0; d< max_aa; d++){
341            int i;
342            for (i=0; i<3; i++) {
343                if (dist->patd[i] & (1<<d)) break;
344            }
345            printf ("%c%i ",translator->Index2Spro(d),i);
346        }
347        printf ("\n");
348    }
349}
350#endif // DEBUG
351
352// ---------------------------
353//      AWT_distance_meter
354// ---------------------------
355
356AWT_distance_meter::AWT_distance_meter(const AWT_translator *translator) {
357    memset(dist_, 0, sizeof(dist_));
358    memset(transform07, 0, sizeof(transform07));
359    memset(transform815, 0, sizeof(transform815));
360    memset(transform1623, 0, sizeof(transform1623));
361
362    int s;
363    int i;
364    int max_aa     = translator->MaxAA();
365    int realmax_aa = translator->RealmaxAA();
366
367    for (s = 0; s< max_aa; s++) {
368        dist_[s] = (AWT_PDP *)calloc(sizeof(AWT_PDP), max_aa);
369
370        const arb_r2a_pro_2_nuc *s2str = translator->S2str(translator->Index2Spro(s));
371        for (i=0;i<3;i++) {
372            dist_[s]->nucbits[i] = s2str->nucs->nucbits[i];
373        }
374    }
375
376    for (s = 0; s< max_aa; s++){
377        for (int d = 0; d< realmax_aa; d++){
378            int dist = nuc_dist(translator, translator->Index2Spro(s), translator->Index2Spro(d) );
379
380            if (dist==0) dist_[s]->patd[0] |= 1<<d; // distance == 0
381            if (dist<=1) dist_[s]->patd[1] |= 1<<d; // distance <= 1
382        }
383        dist_[s]->patd[2] |= dist_[s]->patd[1]; // (distance <= 1) => (distance <= 2)
384        dist_[s]->patd[0] |= 1<<s; // set 'distance to self' == 0
385    }
386
387    for (s = 0; s< max_aa; s++) {
388        long sum = 0;
389        for (int d = 0; d< realmax_aa; d++) {
390            if ( (1 << d) & dist_[s]->patd[1] ) { // if distance(s, d) <= 1
391                sum |= dist_[d]->patd[1]; // collect all proteins which have 'distance <= 1' to 'd'
392            }
393        }
394        dist_[s]->patd[2] |= sum; // and store them in 'distance <= 2'
395    }
396
397    for (i=0; i<256; i++) {
398        for (s = 0; s<8; s++) {
399            if (i & (1<<s) ) {
400                transform07[i]   |= dist_[s]->patd[1];
401                transform815[i]  |= dist_[s+8]->patd[1];
402                transform1623[i] |= dist_[s+16]->patd[1];
403            }
404        }
405    }
406#ifdef DEBUG
407    awt_pro_a_nucs_debug(translator, this);
408#endif
409}
410
411AWT_distance_meter::~AWT_distance_meter() {
412    for (int i=0;i<64; i++){
413        delete dist_[i];
414    }
415   
416}
417
418// --------------------------------------------------------------------------------
419// Translator factory:
420
421static int current_user_code_nr = -1;                // always contain same value as AWAR_PROTEIN_TYPE (after calling AWT_default_protein_type once)
422
423static void user_code_nr_changed_cb(GBDATA *gb_awar, int *, GB_CB_TYPE) {
424    // this callback keeps 'current_user_code_nr' synced with AWAR_PROTEIN_TYPE
425    GBDATA         *gb_main = GB_get_root(gb_awar);
426    GB_transaction  ta(gb_main);
427    current_user_code_nr    = GB_read_int(gb_awar);
428}
429
430#define CACHED_TRANSLATORS 4
431
432#if defined(DEBUG)
433// #define DUMP_TRANSLATOR_ALLOC
434#endif // DEBUG
435
436
437AWT_translator *AWT_get_translator(int code_nr) {
438    static AWT_translator *cached[CACHED_TRANSLATORS];
439
440    if (!cached[0] || cached[0]->CodeNr() != code_nr) { // most recent != requested
441        AWT_translator *translator = 0;
442        int             i;
443
444        for (i = 1; i<CACHED_TRANSLATORS; i++) {
445            if (cached[i] && cached[i]->CodeNr() == code_nr) {
446                // found existing translator
447                translator = cached[i];
448                cached[i]  = 0;
449                break;
450            }
451        }
452
453        if (!translator) {
454            translator = new AWT_translator(code_nr);
455
456#if defined(DUMP_TRANSLATOR_ALLOC)
457            static int allocCount = 0;
458            allocCount++;
459            printf("Alloc translator for code_nr=%i (allocCount=%i)\n", translator->CodeNr(), allocCount);
460#endif // DUMP_TRANSLATOR_ALLOC
461
462        }
463
464        // insert new or found translator at front and shift existing to higher indices:
465        for (i = 0; i<CACHED_TRANSLATORS && translator; i++) {
466            AWT_translator *move = cached[i];
467            cached[i]  = translator;
468            translator = move;
469        }
470
471        // delete oldest translator,  if no empty array position was found:
472        if (translator) {
473#if defined(DUMP_TRANSLATOR_ALLOC)
474            static int freeCount = 0;
475            freeCount++;
476            printf("Free translator for code_nr=%i (freeCount=%i)\n", translator->CodeNr(), freeCount);
477#endif // DUMP_TRANSLATOR_ALLOC
478
479            delete translator;
480        }
481    }
482
483    awt_assert(cached[0]->CodeNr() == code_nr);
484    return cached[0];
485}
486
487int AWT_default_protein_type(GBDATA *gb_main) {
488    // returns protein code selected in AWAR_PROTEIN_TYPE
489
490    if (current_user_code_nr == -1) { // user protein code in AWAR_PROTEIN_TYPE not traced yet
491        awt_assert(gb_main != 0); // either pass gb_main here or call once with valid gb_main at program startup
492
493        {
494            GB_transaction ta(gb_main);
495            GBDATA *awar = GB_search(gb_main, AWAR_PROTEIN_TYPE, GB_INT);
496            GB_add_callback(awar, GB_CB_CHANGED, user_code_nr_changed_cb, 0); // bind a callback that traces AWAR_PROTEIN_TYPE
497            user_code_nr_changed_cb(awar, 0, GB_CB_CHANGED);
498        }
499
500        awt_assert(current_user_code_nr != -1);
501    }
502
503    return current_user_code_nr;
504}
505
506AWT_translator *AWT_get_user_translator(GBDATA *gb_main) {
507    return AWT_get_translator(AWT_default_protein_type(gb_main));
508}
Note: See TracBrowser for help on using the browser.