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

Last change on this file was 14803, checked in by westram, 3 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.1 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : AP_codon_table.cxx                                //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Coded by Ralf Westram (coder@reallysoft.de) in January 2010   //
7//   Institute of Microbiology (Technical University Munich)       //
8//   http://www.arb-home.de/                                       //
9//                                                                 //
10// =============================================================== //
11
12#include "AP_codon_table.hxx"
13#include "iupac.h"
14
15#include <arbdb.h>
16
17#include <cctype>
18
19#define pn_assert(cond) arb_assert(cond)
20
21#define EMBL_BACTERIAL_TABLE_INDEX 11
22
23// Info about translation codes was taken from
24// http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi
25
26static AWT_Codon_Code_Definition AWT_codon_def[AWT_CODON_TABLES+1] =
27    {
28        //   0000000001111111111222222222233333333334444444444555555555566666
29        //   1234567890123456789012345678901234567890123456789012345678901234
30
31        //  "TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG",  base1
32        //  "TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG",  base2
33        //  "TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG"   base3
34        {
35            " (1) Standard code",
36            "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG", // The first code in this table has to be 'Standard code'!
37            "---M---------------M---------------M----------------------------",
38            1
39        },
40        {
41            " (2) Vertebrate mitochondrial code",
42            "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSS**VVVVAAAADDEEGGGG",
43            "--------------------------------MMMM---------------M------------",
44            2
45        },
46        {
47            " (3) Yeast mitochondrial code",
48            "FFLLSSSSYY**CCWWTTTTPPPPHHQQRRRRIIMMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
49            "----------------------------------MM----------------------------",
50            3
51        },
52        {
53            " (4) Mold/Protozoan/Coelenterate mito. + Mycoplasma/Spiroplasma code",
54            "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
55            "--MM---------------M------------MMMM---------------M------------",
56            4
57        },
58        {
59            " (5) Invertebrate mitochondrial code",
60            "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSSSVVVVAAAADDEEGGGG",
61            "---M----------------------------MMMM---------------M------------",
62            5
63        },
64        {
65            " (6) Ciliate, Dasycladacean and Hexamita nuclear code",
66            "FFLLSSSSYYQQCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
67            "-----------------------------------M----------------------------",
68            6
69        },
70        {
71            " (9) Echinoderm and Flatworm mitochondrial code",
72            "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNNKSSSSVVVVAAAADDEEGGGG",
73            "-----------------------------------M---------------M------------",
74            9
75        },
76        {
77            "(10) Euplotid nuclear code",
78            "FFLLSSSSYY**CCCWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
79            "-----------------------------------M----------------------------",
80            10
81        },
82        //   0000000001111111111222222222233333333334444444444555555555566666
83        //   1234567890123456789012345678901234567890123456789012345678901234
84
85        //  "TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG",  base1
86        //  "TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG",  base2
87        //  "TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG"   base3
88        {
89            "(11) Bacterial and Plant Plastid code",
90            "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
91            "---M---------------M------------MMMM---------------M------------",
92            11
93        },
94        {
95            "(12) Alternative Yeast nuclear code",
96            "FFLLSSSSYY**CC*WLLLSPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
97            "-------------------M---------------M----------------------------",
98            12
99        },
100        {
101            "(13) Ascidian mitochondrial code",
102            "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSGGVVVVAAAADDEEGGGG",
103            "---M------------------------------MM---------------M------------",
104            13
105        },
106        {
107            "(14) Alternative Flatworm mitochondrial code",
108            "FFLLSSSSYYY*CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNNKSSSSVVVVAAAADDEEGGGG",
109            "-----------------------------------M----------------------------",
110            14
111        },
112        {
113            "(15) Blepharisma nuclear code",
114            "FFLLSSSSYY*QCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
115            "-----------------------------------M----------------------------",
116            15
117        },
118        {
119            "(16) Chlorophycean mitochondrial code",
120            "FFLLSSSSYY*LCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
121            "-----------------------------------M----------------------------",
122            16
123        },
124        {
125            "(21) Trematode mitochondrial code",
126            "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNNKSSSSVVVVAAAADDEEGGGG",
127            "-----------------------------------M---------------M------------",
128            21
129        },
130        {
131            "(22) Scenedesmus obliquus mitochondrial code",
132            "FFLLSS*SYY*LCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
133            "-----------------------------------M----------------------------",
134            22
135        },
136        {
137            "(23) Thraustochytrium mitochondrial code",
138            "FF*LSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
139            "--------------------------------M--M---------------M------------",
140            23
141        },
142
143        { 0, 0, 0, 0 } // end of table-marker
144    };
145
146#define MAX_EMBL_TRANSL_TABLE_VALUE 23 // maximum known EMBL transl_table value
147
148int AWT_embl_transl_table_2_arb_code_nr(int embl_code_nr) {
149    // returns -1 if embl_code_nr is not known by ARB
150
151    static bool initialized = false;
152    static int  arb_code_nr_table[MAX_EMBL_TRANSL_TABLE_VALUE+1];                 // key: embl_code_nr, value: arb_code_nr or -1
153
154    if (!initialized) {
155        for (int embl = 0; embl <= MAX_EMBL_TRANSL_TABLE_VALUE; ++embl) {
156            arb_code_nr_table[embl] = -1; // illegal table
157        }
158        for (int arb_code_nr = 0; arb_code_nr < AWT_CODON_TABLES; ++arb_code_nr) {
159            arb_code_nr_table[AWT_codon_def[arb_code_nr].embl_feature_transl_table] = arb_code_nr;
160        }
161        // should be index of 'Bacterial and Plant Plastid code'
162        // (otherwise maybe AWAR_PROTEIN_TYPE_bacterial_code_index  is wrong)
163        pn_assert(arb_code_nr_table[EMBL_BACTERIAL_TABLE_INDEX] == AWAR_PROTEIN_TYPE_bacterial_code_index);
164        pn_assert(arb_code_nr_table[1] == 0); // Standard code has to be on index zero!
165
166        initialized = true;
167    }
168
169    if (embl_code_nr<0 || embl_code_nr>MAX_EMBL_TRANSL_TABLE_VALUE) return -1;
170
171    int arb_code_nr = arb_code_nr_table[embl_code_nr];
172#ifdef DEBUG
173    if (arb_code_nr != -1) {
174        pn_assert(arb_code_nr >= 0 && arb_code_nr < AWT_CODON_TABLES);
175        pn_assert(AWT_arb_code_nr_2_embl_transl_table(arb_code_nr) == embl_code_nr);
176    }
177#endif
178    return arb_code_nr;
179}
180
181int AWT_arb_code_nr_2_embl_transl_table(int arb_code_nr) {
182    pn_assert(arb_code_nr >= 0 && arb_code_nr<AWT_CODON_TABLES);
183    return AWT_codon_def[arb_code_nr].embl_feature_transl_table;
184}
185
186
187static bool codon_tables_initialized = false;
188static char definite_translation[AWT_MAX_CODONS]; // contains 0 if ambiguous, otherwise it contains the definite translation
189static char *ambiguous_codons[AWT_MAX_CODONS]; // for each ambiguous codon: contains all translations (each only once)
190
191void AP_initialize_codon_tables() {
192    if (codon_tables_initialized) return;
193
194    int codon_nr;
195    int code_nr;
196
197    for (codon_nr=0; codon_nr<AWT_MAX_CODONS; codon_nr++) {
198        ambiguous_codons[codon_nr] = 0;
199    }
200
201    pn_assert(AWT_CODON_TABLES>=1);
202    memcpy(definite_translation, AWT_codon_def[0].aa, AWT_MAX_CODONS); // only one translation is really definite
203
204    pn_assert(AWT_codon_def[AWT_CODON_TABLES].aa==NULL); // Error in AWT_codon_def or AWT_CODON_CODES
205
206    for (code_nr=1; code_nr<AWT_CODON_TABLES; code_nr++) {
207        const char *translation = AWT_codon_def[code_nr].aa;
208
209        for (codon_nr=0; codon_nr<AWT_MAX_CODONS; codon_nr++) {
210            if (definite_translation[codon_nr]!='?') { // is definite till now
211                if (definite_translation[codon_nr]!=translation[codon_nr]) { // we found a different translation
212                    // create ambiguous_codons:
213                    char *amb = (char*)GB_calloc(AWT_MAX_CODONS+1, sizeof(char));
214                    amb[0] = definite_translation[codon_nr];
215                    amb[1] = translation[codon_nr];
216
217                    ambiguous_codons[codon_nr] = amb;
218                    definite_translation[codon_nr] = '?';
219#if defined(DEBUG) && 0
220                    printf("amb[%i]='%s'\n", codon_nr, amb);
221#endif
222                }
223            }
224            else { // is ambiguous
225                if (strchr(ambiguous_codons[codon_nr], translation[codon_nr])==0) { // not listed in ambiguous codons
226                    // append another ambiguous codon:
227                    char *amb = ambiguous_codons[codon_nr];
228                    amb[strlen(amb)] = translation[codon_nr];
229#if defined(DEBUG) && 0
230                    printf("amb[%i]='%s'\n", codon_nr, amb);
231#endif
232                }
233            }
234        }
235    }
236
237    codon_tables_initialized = true;
238}
239
240// return 0..3 (ok) or 4 (failure)
241inline int dna2idx(char c) {
242    switch (c) {
243        case 'T': case 't':
244        case 'U': case 'u': return 0;
245        case 'C': case 'c': return 1;
246        case 'A': case 'a': return 2;
247        case 'G': case 'g': return 3;
248    }
249    return 4;
250}
251
252inline char idx2dna(int idx) {
253    pn_assert(idx>=0 && idx<4);
254    return "TCAG"[idx];
255}
256
257inline int calc_codon_nr(const char *dna) {
258    int i1 = dna2idx(dna[0]); if (i1 == 4) return AWT_MAX_CODONS; // is not a codon
259    int i2 = dna2idx(dna[1]); if (i2 == 4) return AWT_MAX_CODONS;
260    int i3 = dna2idx(dna[2]); if (i3 == 4) return AWT_MAX_CODONS;
261
262    int codon_nr = i1*16 + i2*4 + i3;
263    pn_assert(codon_nr>=0 && codon_nr<=AWT_MAX_CODONS);
264    return codon_nr;
265}
266
267inline void build_codon(int codon_nr, char *to_buffer) {
268    pn_assert(codon_nr>=0 && codon_nr<AWT_MAX_CODONS);
269
270    to_buffer[0] = idx2dna((codon_nr>>4)&3);
271    to_buffer[1] = idx2dna((codon_nr>>2)&3);
272    to_buffer[2] = idx2dna(codon_nr&3);
273}
274
275const char* AWT_get_codon_code_name(int code) {
276    pn_assert(code>=0 && code<AWT_CODON_TABLES);
277    return AWT_codon_def[code].name;
278}
279
280static const char *protein_name[26+1] = {
281    "Ala", // A
282    "Asx", // B
283    "Cys", // C
284    "Asp", // D
285    "Glu", // E
286    "Phe", // F
287    "Gly", // G
288    "His", // H
289    "Ile", // I
290    0,     // J
291    "Lys", // K
292    "Leu", // L
293    "Met", // M
294    "Asn", // N
295    0,     // O
296    "Pro", // P
297    "Gln", // Q
298    "Arg", // R
299    "Ser", // S
300    "Thr", // T
301    0,     // U
302    "Val", // V
303    "Trp", // W
304    "Xxx", // X
305    "Tyr", // Y
306    "Glx", // Z
307    0
308};
309
310const char *AP_get_protein_name(char protein) {
311    if (protein=='*') return "End";
312    if (protein=='-') return "---";
313
314    pn_assert(protein>='A' && protein<='Z');
315    pn_assert(protein_name[protein-'A']!=0);
316    return protein_name[protein-'A'];
317}
318
319#ifdef DEBUG
320
321inline char nextBase(char c) {
322    switch (c) {
323        case 'T': return 'C';
324        case 'C': return 'A';
325        case 'A': return 'G';
326        case 'G': return 0;
327        default: pn_assert(0);
328    }
329    return 0;
330}
331
332void AWT_dump_codons() {
333    AWT_allowedCode allowed_code;
334
335    for (char c='*'; c<='Z'; c++) {
336        printf("Codes for '%c': ", c);
337        int first_line = 1;
338        int found = 0;
339        for (char b1='T'; b1; b1=nextBase(b1)) {
340            for (char b2='T'; b2; b2=nextBase(b2)) {
341                for (char b3='T'; b3; b3=nextBase(b3)) {
342                    char dna[4];
343                    dna[0]=b1;
344                    dna[1]=b2;
345                    dna[2]=b3;
346                    dna[3]=0;
347
348                    AWT_allowedCode allowed_code_left;
349                    if (AWT_is_codon(c, dna, allowed_code, allowed_code_left)) {
350                        if (!first_line) printf("\n               ");
351                        first_line = 0;
352                        printf("%s (", dna);
353
354                        int first=1;
355                        for (int code=0; code<AWT_CODON_TABLES; code++) {
356                            if (allowed_code_left.is_allowed(code)) {
357                                if (!first) printf(",");
358                                first=0;
359                                printf("%i", code);
360                            }
361                        }
362                        printf(") ");
363
364                        found = 1;
365                    }
366                }
367            }
368        }
369        if (!found) printf("none");
370        printf("\n");
371        if (c=='*') c='A'-1;
372    }
373}
374#endif
375
376char AWT_is_start_codon(const char *dna, int arb_code_nr) {
377    // if dna[0]..dna[2] is defined as start codon for 'arb_code_nr'
378    //                  return 'M' (or whatever is defined in tables)
379    // return 0 otherwise
380
381    char is_start_codon = 0;
382    int  codon_nr       = calc_codon_nr(dna);
383
384    pn_assert(arb_code_nr >= 0 && arb_code_nr<AWT_CODON_TABLES);
385
386    if (codon_nr != AWT_MAX_CODONS) { // dna is a clean codon (it contains no iupac-codes)
387        const char *starts = AWT_codon_def[arb_code_nr].starts;
388
389        is_start_codon = starts[codon_nr];
390        if (is_start_codon == '-') is_start_codon = 0;
391    }
392
393    return is_start_codon;
394}
395
396
397bool AWT_is_codon(char protein, const char *dna, const AWT_allowedCode& allowed_code, AWT_allowedCode& allowed_code_left, const char **fail_reason_ptr) {
398    // return TRUE if 'dna' contains a codon of 'protein' ('dna' must not contain any gaps)
399    // allowed_code contains 1 for each allowed code and 0 otherwise
400    // allowed_code_left contains a copy of allowed_codes with all impossible codes set to zero
401
402    pn_assert(codon_tables_initialized);
403
404    const char *fail_reason = 0;
405    bool        is_codon    = false;
406
407    if (fail_reason_ptr) *fail_reason_ptr = 0;
408
409    protein = toupper(protein);
410    if (protein=='B') {         // B is a shortcut for Asp(=D) or Asn(=N)
411        is_codon = AWT_is_codon('D', dna, allowed_code, allowed_code_left, &fail_reason);
412        if (!is_codon) {
413            pn_assert(fail_reason != 0); // if failed there should always be a failure-reason
414            char *fail1 = strdup(fail_reason);
415            is_codon    = AWT_is_codon('N', dna, allowed_code, allowed_code_left, &fail_reason);
416            if (!is_codon) {
417                char *fail2 = strdup(fail_reason);
418                fail_reason = GBS_global_string("%s and %s", fail1, fail2);
419                free(fail2);
420            }
421            free(fail1);
422        }
423    }
424    else if (protein=='Z') {    // Z is a shortcut for Glu(=E) or Gln(=Q)
425        is_codon = AWT_is_codon('E', dna, allowed_code, allowed_code_left, &fail_reason);
426        if (!is_codon) {
427            pn_assert(fail_reason != 0); // if failed there should always be a failure-reason
428            char *fail1 = strdup(fail_reason);
429            is_codon    = AWT_is_codon('Q', dna, allowed_code, allowed_code_left, &fail_reason);
430            if (!is_codon) {
431                char *fail2 = strdup(fail_reason);
432                fail_reason = GBS_global_string("%s and %s", fail1, fail2);
433                free(fail2);
434            }
435            free(fail1);
436        }
437    }
438    else {
439        int codon_nr = calc_codon_nr(dna);
440        if (codon_nr==AWT_MAX_CODONS) { // dna is not a clean codon (it contains iupac-codes)
441            int  error_positions = 0;
442            int  first_error_pos = -1;
443            bool too_short       = false;
444            {
445                int iupac_pos;
446                for (iupac_pos=0; iupac_pos<3 && !too_short; iupac_pos++) {
447                    if (!dna[iupac_pos]) {
448                        too_short = true;
449                    }
450                    else if (strchr("ACGTU", dna[iupac_pos]) == 0) {
451                        if (first_error_pos==-1) first_error_pos = iupac_pos;
452                        error_positions++;
453                    }
454                }
455            }
456
457            if (too_short) {
458                fail_reason = GBS_global_string("Not enough nucleotides (got '%s')", dna);
459            }
460            else {
461                pn_assert(error_positions);
462                if (error_positions==3) { // don't accept codons with 3 errors
463                    fail_reason = GBS_global_string("Three consecutive IUPAC codes '%c%c%c'", dna[0], dna[1], dna[2]);
464                }
465                else {
466                    const char *decoded_iupac = iupac::decode(dna[first_error_pos], GB_AT_DNA, 0);
467
468                    if (!decoded_iupac[0]) { // no valid IUPAC
469                        allowed_code_left.forbidAll();
470                        fail_reason = GBS_global_string("Not a valid IUPAC code:'%c'", dna[first_error_pos]);
471                    }
472                    else {
473                        char dna_copy[4];
474                        memcpy(dna_copy, dna, 3);
475                        dna_copy[3] = 0;
476
477#if defined(DEBUG) && 0
478                        printf("Check if '%s' is a codon for '%c'\n", dna_copy, protein);
479#endif
480
481                        int all_are_codons = 1;
482                        AWT_allowedCode allowed_code_copy;
483                        allowed_code_copy = allowed_code;
484
485                        for (int i=0; decoded_iupac[i]; i++) {
486                            dna_copy[first_error_pos] = decoded_iupac[i];
487                            if (!AWT_is_codon(protein, dna_copy, allowed_code_copy, allowed_code_left)) {
488                                all_are_codons = 0;
489                                break;
490                            }
491                            allowed_code_copy = allowed_code_left;
492                        }
493
494                        if (all_are_codons) {
495                            allowed_code_left = allowed_code_copy;
496                            is_codon          = true;
497                        }
498                        else {
499                            allowed_code_left.forbidAll();
500                            fail_reason = GBS_global_string("Not all IUPAC-combinations of '%s' translate", dna_copy);
501                        }
502#if defined(DEBUG) && 0
503                        printf("result      = %i\n", all_are_codons);
504#endif
505                    }
506                }
507            }
508        }
509        else if (definite_translation[codon_nr]!='?') {
510            int ok = definite_translation[codon_nr]==protein;
511
512            if (ok) {
513                allowed_code_left = allowed_code;
514                is_codon          = true;
515            }
516            else {
517                allowed_code_left.forbidAll();
518                fail_reason = GBS_global_string("'%c%c%c' does never translate to '%c' (1)", dna[0], dna[1], dna[2], protein);
519            }
520        }
521        else if (strchr(ambiguous_codons[codon_nr], protein)==0) {
522            allowed_code_left.forbidAll();
523            fail_reason = GBS_global_string("'%c%c%c' does never translate to '%c' (2)", dna[0], dna[1], dna[2], protein);
524        }
525        else {
526#if defined(ASSERTION_USED)
527            bool correct_disallowed_translation = false;
528#endif
529
530            // search for allowed correct translation possibility:
531            for (int code_nr=0; code_nr<AWT_CODON_TABLES; code_nr++) {
532                if (AWT_codon_def[code_nr].aa[codon_nr] == protein) { // does it translate correct?
533                    if (allowed_code.is_allowed(code_nr)) { // is this code allowed?
534                        allowed_code_left.allow(code_nr);
535                        is_codon = true;
536                    }
537                    else {
538                        allowed_code_left.forbid(code_nr); // otherwise forbid code in future
539#if defined(ASSERTION_USED)
540                        correct_disallowed_translation = true;
541#endif
542                    }
543                }
544                else {
545                    allowed_code_left.forbid(code_nr); // otherwise forbid code in future
546                }
547            }
548
549            if (!is_codon) {
550                pn_assert(correct_disallowed_translation); // should be true because otherwise we shouldn't run into this else-branch
551                char  left_tables[AWT_CODON_TABLES*3+1];
552                char *ltp   = left_tables;
553                bool  first = true;
554                for (int code_nr=0; code_nr<AWT_CODON_TABLES; code_nr++) {
555                    if (allowed_code.is_allowed(code_nr)) {
556                        if (!first) *ltp++ = ',';
557                        ltp   += sprintf(ltp, "%i", code_nr);
558                        first  = false;
559                    }
560                }
561                fail_reason = GBS_global_string("'%c%c%c' does not translate to '%c' for any of the leftover trans-tables (%s)",
562                                                dna[0], dna[1], dna[2], protein, left_tables);
563            }
564        }
565    }
566
567    if (!is_codon) {
568        pn_assert(fail_reason);
569        if (fail_reason_ptr) *fail_reason_ptr = fail_reason; // set failure-reason if requested
570    }
571    return is_codon;
572}
573
574// -------------------------------------------------------------------------------- Codon_Group
575
576class Codon_Group
577{
578    char codon[64]; // index is calculated with calc_codon_nr
579
580public:
581    Codon_Group(char protein, int code_nr);
582    ~Codon_Group() {}
583
584    Codon_Group& operator += (const Codon_Group& other);
585    int expand(char *to_buffer) const;
586};
587
588Codon_Group::Codon_Group(char protein, int code_nr) {
589    protein = toupper(protein);
590    pn_assert(protein=='*' || isalpha(protein));
591    pn_assert(code_nr>=0 && code_nr<AWT_CODON_TABLES);
592
593    const char *amino_table = AWT_codon_def[code_nr].aa;
594    for (int i=0; i<AWT_MAX_CODONS; i++) {
595        codon[i] = amino_table[i]==protein;
596    }
597}
598
599Codon_Group& Codon_Group::operator+=(const Codon_Group& other) {
600    for (int i=0; i<AWT_MAX_CODONS; i++) {
601        codon[i] = codon[i] || other.codon[i];
602    }
603    return *this;
604}
605
606inline int legal_dna_no(int i) { return i>=0 && i<4; }
607inline void my_memcpy(char *dest, const char *source, size_t length) { for (size_t l=0; l<length; l++) { dest[l] = source[l]; } }
608
609inline const char *buildMixedCodon(const char *con1, const char *con2) {
610    int mismatches = 0;
611    int mismatch_index = -1;
612    static char buf[4];
613
614    for (int i=0; i<3; i++) {
615        if (con1[i]!=con2[i]) {
616            mismatches++;
617            mismatch_index = i;
618        }
619        else {
620            buf[i] = con1[i];
621        }
622    }
623
624    if (mismatches==1) { // exactly one position differs between codons
625        pn_assert(mismatch_index!=-1);
626        buf[mismatch_index] = iupac::combine(con1[mismatch_index], con2[mismatch_index], GB_AT_DNA);
627        buf[3] = 0;
628        return buf;
629    }
630    return 0;
631}
632
633static int expandMore(const char *bufferStart, int no_of_condons, char*&to_buffer) {
634    int i, j;
635    const char *con1, *con2;
636    int added = 0;
637
638    for (i=0; i<no_of_condons; i++) {
639        con1 = bufferStart+3*i;
640
641        for (j=i+1; j<no_of_condons; j++) {
642            con2 = bufferStart+3*j;
643            const char *result = buildMixedCodon(con1, con2);
644            if (result) {
645                to_buffer[0] = 0;
646                // do we already have this codon?
647                const char *found;
648                const char *startSearch = bufferStart;
649                for (;;) {
650                    found = strstr(startSearch, result);
651                    if (!found) break;
652                    int pos = (found-bufferStart);
653                    if ((pos%3)==0) break; // yes already here!
654                    startSearch = found+1; // was misaligned -> try behind
655                }
656
657                if (!found) {
658                    my_memcpy(to_buffer, result, 3); to_buffer+=3;
659                    added++;
660                }
661            }
662        }
663    }
664    return no_of_condons+added;
665}
666
667int Codon_Group::expand(char *to_buffer) const {
668    int count = 0;
669    int i;
670    char *org_to_buffer = to_buffer;
671
672    for (i=0; i<AWT_MAX_CODONS; i++) {
673        if (codon[i]) {
674            build_codon(i, to_buffer);
675            to_buffer += 3;
676            count++;
677        }
678    }
679
680#if defined(DEBUG) && 0
681    to_buffer[0] = 0;
682    printf("codons = '%s'\n", org_to_buffer);
683#endif
684
685    for (;;) {
686        int new_count = expandMore(org_to_buffer, count, to_buffer);
687        if (new_count==count) break; // nothing expanded -> done
688        count = new_count;
689#if defined(DEBUG) && 0
690        to_buffer[0] = 0;
691        printf("codons (expandedMore) = '%s'\n", org_to_buffer);
692#endif
693    }
694
695    pn_assert(count==(int(to_buffer-org_to_buffer)/3));
696
697    return count;
698}
699
700// --------------------------------------------------------------------------------
701
702static Codon_Group *get_Codon_Group(char protein, int code_nr) {
703    pn_assert(code_nr>=0 && code_nr<AWT_CODON_TABLES);
704    protein = toupper(protein);
705    pn_assert(isalpha(protein) || protein=='*');
706    pn_assert(codon_tables_initialized);
707
708    Codon_Group *cgroup = 0;
709
710    if (protein=='B') {
711        cgroup = new Codon_Group('D', code_nr);
712        Codon_Group N('N', code_nr);
713        *cgroup += N;
714    }
715    else if (protein=='Z') {
716        cgroup = new Codon_Group('E', code_nr);
717        Codon_Group Q('Q', code_nr);
718        *cgroup += Q;
719    }
720    else {
721        cgroup = new Codon_Group(protein, code_nr);
722    }
723
724    pn_assert(cgroup);
725
726    return cgroup;
727}
728
729#define MAX_CODON_LIST_LENGTH (70*3)
730
731// get a list of all codons ("xyzxyzxyz...") encoding 'protein' in case we use Codon-Code 'code_nr'
732// (includes all completely contained IUPAC-encoded codons at the end of list)
733const char *AP_get_codons(char protein, int code_nr) {
734    Codon_Group *cgroup = get_Codon_Group(protein, code_nr);
735
736    static char buffer[MAX_CODON_LIST_LENGTH+1];
737    int offset = 3*cgroup->expand(buffer);
738    pn_assert(offset<MAX_CODON_LIST_LENGTH);
739    buffer[offset] = 0;
740
741    delete cgroup;
742
743    return buffer;
744}
745
Note: See TracBrowser for help on using the repository browser.