source: branches/stable/SL/PRONUC/AP_pro_a_nucs.cxx

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