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

Last change on this file was 11889, checked in by westram, 5 years ago
  • add missing comment to closing #endif of test-sections
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.4 KB
Line 
1// ================================================================ //
2//                                                                  //
3//   File      : iupac.cxx                                          //
4//   Purpose   :                                                    //
5//                                                                  //
6//   Institute of Microbiology (Technical University Munich)        //
7//   http://www.arb-home.de/                                        //
8//                                                                  //
9// ================================================================ //
10
11#include "iupac.h"
12#include <cctype>
13
14#define awt_assert(bed) arb_assert(bed)
15
16#define IUPAC_EMPTY "" // changed from " " to "" (2009-03-19 -- ralf)
17#define ILL_CODE    char(26)
18
19using namespace iupac;
20
21namespace iupac {
22#define _____________ {0,0}
23
24    const Nuc_Group nuc_group[26][2] = {
25        { { "A",    1 }, { "A",    1 } }, // A
26        { { "CGT",  3 }, { "CGU",  3 } }, // B
27        { { "C",    1 }, { "C",    1 } }, // C
28        { { "AGT",  3 }, { "AGU",  3 } }, // D
29        { _____________, _____________ }, // E
30        { _____________, _____________ }, // F
31        { { "G",    1 }, { "G",    1 } }, // G
32        { { "ACT",  3 }, { "ACU",  3 } }, // H
33        { _____________, _____________ }, // I
34        { _____________, _____________ }, // J
35        { { "GT",   2 }, { "GU",   2 } }, // K
36        { _____________, _____________ }, // L
37        { { "AC",   2 }, { "AC",   2 } }, // M
38        { { "ACGT", 1 }, { "ACGU", 1 } }, // N
39        { _____________, _____________ }, // O
40        { _____________, _____________ }, // P
41        { _____________, _____________ }, // Q
42        { { "AG",   2 }, { "AG",   2 } }, // R
43        { { "CG",   2 }, { "CG",   2 } }, // S
44        { { "T",    1 }, { "U",    1 } }, // T
45        { { "T",    1 }, { "U",    1 } }, // U
46        { { "ACG",  3 }, { "ACG",  3 } }, // V
47        { { "AT",   2 }, { "AU",   2 } }, // W
48        { _____________, _____________ }, // X
49        { { "CT",   2 }, { "CU",   2 } }, // Y
50        { _____________, _____________ }, // Z
51
52        // (In each single string the nucs have to be sorted alphabetically)
53    };
54
55#undef _____________
56
57    static const char *aminoGroupMembers[AA_GROUP_COUNT] = {
58        "JOUX",   // AA_GROUP_NONE
59        "AGPST",  // AA_GROUP_ALPHA
60        "BDENQZ", // AA_GROUP_BETA
61        "HKR",    // AA_GROUP_GAMMA
62        "ILMV",   // AA_GROUP_DELTA
63        "FWY",    // AA_GROUP_EPSILON
64        "C",      // AA_GROUP_ZETA
65    };
66   
67    static Amino_Group amino_group[26];
68
69    class Setup { // setup static data in this module
70        void setup_amino_group() {
71            for (int c = 0; c<26; ++c) amino_group[c] = AA_GROUP_ILLEGAL;
72            for (Amino_Group g = AA_GROUP_NONE; g<AA_GROUP_COUNT; g = Amino_Group(g+1)) {
73                const char *members = aminoGroupMembers[g];
74                for (int p = 0; members[p]; ++p) {
75                    int c = members[p];
76                    awt_assert(c >= 'A' && c <= 'Z');
77                    amino_group[c-'A'] = g;
78                }
79            }
80        }
81    public:
82        Setup() { setup_amino_group(); }
83    };
84    static Setup perform;
85
86};
87
88Amino_Group iupac::get_amino_group_for(char aa) {
89    aa = toupper(aa);
90    if (aa<'A' || aa>'Z') { return AA_GROUP_ILLEGAL; }
91    return amino_group[aa-'A'];
92}
93
94char iupac::get_amino_consensus_char(Amino_Group ag) {
95    if (ag>AA_GROUP_NONE && ag<AA_GROUP_ILLEGAL) {
96        return aminoGroupMembers[ag][0];
97    }
98    return '?';
99}
100
101static char IUPAC_add[26][26]; // uses T
102static int IUPAC_add_initialized = 0;
103
104static void initialize_IUPAC_add()
105{
106    int c1, c2;
107
108    for (c1=0; c1<26; c1++) {
109        const char *decoded1 = nuc_group[c1][0].members;
110
111        if (!decoded1) {
112            for (c2=0; c2<=c1; c2++) {
113                IUPAC_add[c1][c2] = ILL_CODE;
114                IUPAC_add[c2][c1] = ILL_CODE;
115            }
116        }
117        else {
118            IUPAC_add[c1][c1] = char(c1);       // add char to same char
119
120            for (c2=0; c2<c1; c2++) {
121                if (strchr(decoded1, 'A'+c2)!=0) {      // char is already contained in this IUPAC
122                    IUPAC_add[c1][c2] = char(c1);
123                }
124                else {
125                    const char *decoded2 = nuc_group[c2][0].members;
126
127                    if (!decoded2) { // char is illegal
128                        IUPAC_add[c1][c2] = ILL_CODE;
129                    }
130                    else {
131#define MAX_MIXED 5
132                        char mixed[MAX_MIXED];
133                        char *mp = mixed;
134                        const char *d1 = decoded1;
135                        const char *d2 = decoded2;
136
137                        while (1) {
138                            char z1 = *d1;
139                            char z2 = *d2;
140
141                            if (!z1 && !z2) break;
142
143                            if (z1==z2) {
144                                *mp++ = z1;
145                                d1++;
146                                d2++;
147                            }
148                            else if (!z2 || (z1 && z1<z2)) {
149                                *mp++ = z1;
150                                d1++;
151                            }
152                            else {
153                                awt_assert(!z1 || (z2 && z2<z1));
154                                *mp++ = z2;
155                                d2++;
156                            }
157                        }
158
159                        awt_assert((mp-mixed)<MAX_MIXED);
160                        *mp++ = 0;
161
162#if !defined(NDEBUG) && 0
163                        printf("Mix '%s' + '%s' = '%s'\n", decoded1, decoded2, mixed);
164#endif
165
166                        int c3;
167
168                        for (c3=0; c3<26; c3++) {
169                            if (nuc_group[c3][0].members && strcmp(mixed, nuc_group[c3][0].members)==0) {
170                                IUPAC_add[c1][c2] = char(c3);
171                                break;
172                            }
173                        }
174
175                        if (c3>=26) {
176                            IUPAC_add[c1][c2] = 0;
177                        }
178                    }
179                }
180
181                if (IUPAC_add[c1][c2]==('U'-'A')) {
182                    IUPAC_add[c1][c2] = 'T'-'A';
183                }
184
185                IUPAC_add[c2][c1] = IUPAC_add[c1][c2];
186            }
187        }
188
189    }
190}
191
192char iupac::encode(const char bases[], GB_alignment_type ali)
193{
194    if (!IUPAC_add_initialized) {
195        initialize_IUPAC_add();
196        IUPAC_add_initialized = 1;
197    }
198
199    if (ali==GB_AT_AA) {
200        return '?';
201    }
202
203    int i = 0;
204    unsigned char c1;
205
206    while (1) {
207        c1 = bases[i++];
208        if (!c1) return '-';
209        if (isalpha(c1)) break;
210    }
211
212    c1 = toupper(c1)-'A';
213    unsigned char c = IUPAC_add[c1][c1];
214
215    while (c!=ILL_CODE && bases[i]) {
216        if (isalpha(bases[i])) {
217            int c2 = toupper(bases[i])-'A';
218            c = IUPAC_add[c][c2];
219        }
220        i++;
221    }
222
223    if (c==ILL_CODE) {
224        return '-';
225    }
226
227    c += 'A';
228    if (c=='T' && ali==GB_AT_RNA) c = 'U';
229    return c;
230}
231
232char iupac::combine(char c1, char c2, GB_alignment_type ali) {
233    char buffer[3];
234    buffer[0] = c1;
235    buffer[1] = c2;
236    buffer[2] = 0;
237    return iupac::encode(buffer, ali);
238}
239
240const char* iupac::decode(char iupac, GB_alignment_type ali, int decode_amino_iupac_groups) {
241    if (!isalpha(iupac)) {
242        return IUPAC_EMPTY;
243    }
244
245    if (ali==GB_AT_AA) {
246        if (decode_amino_iupac_groups) {
247            Amino_Group group = get_amino_group_for(iupac);
248            if (group == AA_GROUP_NONE) {
249                return "";
250            }
251            else {
252                awt_assert(group >= AA_GROUP_NONE && group<AA_GROUP_COUNT);
253                return aminoGroupMembers[group];
254            }
255        }
256        return "?";
257    }
258
259    const char *decoded = nuc_group[toupper(iupac)-'A'][ali==GB_AT_RNA ? 1 : 0].members;
260
261    return decoded ? decoded : IUPAC_EMPTY;
262}
263
264// --------------------------------------------------------------------------------
265
266#ifdef UNIT_TESTS
267
268#include <test_unit.h>
269
270void TEST_amino_groups() {
271    TEST_EXPECT_EQUAL(decode('x', GB_AT_AA, true), "");
272    TEST_EXPECT_EQUAL(decode('A', GB_AT_AA, true), "AGPST");
273
274    // each character (A-Z) shall be member of a group:
275    for (char c = 'A'; c <= 'Z'; ++c) {
276        Amino_Group group = get_amino_group_for(c);
277        TEST_EXPECT_LESS(group, AA_GROUP_COUNT);
278        TEST_EXPECT_DIFFERENT(group, AA_GROUP_ILLEGAL);
279        TEST_EXPECT_IN_RANGE(group, AA_GROUP_NONE, AA_GROUP_ZETA);
280
281        {
282            const char * const members = aminoGroupMembers[group];
283            const char * const found   = strchr(members, c);
284            TEST_EXPECT_EQUAL(found ? found[0] : 0, c); // char has to be a member of its group
285        }
286
287        if (group != AA_GROUP_NONE) {
288            const char * const decoded = decode(c, GB_AT_AA, true);
289            const char * const found   = strchr(decoded, c);
290
291            TEST_EXPECT_EQUAL(found ? found[0] : 0, c); // check char is
292        }
293    }
294
295    // no character may be member of two groups
296    for (Amino_Group group = AA_GROUP_NONE; group < AA_GROUP_COUNT; group = Amino_Group(group+1)) {
297        const char * const member = aminoGroupMembers[group];
298        for (int pos = 0; member[pos]; ++pos) {
299            Amino_Group groupOfChar = get_amino_group_for(member[pos]);
300            TEST_EXPECT_EQUAL(groupOfChar, group);
301        }
302    }
303}
304
305void TEST_nuc_groups() {
306    for (int base = 0; base<26; base++) {
307        TEST_EXPECT_EQUAL(nuc_group[base][0].count, nuc_group[base][1].count);
308        for (int alitype = 0; alitype<2; alitype++) {
309            const Nuc_Group& group = nuc_group[base][alitype];
310            if (group.members) {
311                if ((base+'A') == 'N') {
312                    TEST_EXPECT_EQUAL__BROKEN(group.count, strlen(group.members), 1);
313                    // @@@ fails because count is 1 for "ACGT" [N]
314                    // maybe be expected by resolve_IUPAC_target_string (fix this first)
315                    TEST_EXPECT_EQUAL(group.count, 1U); // fixture for behavior
316                }
317                else {
318                    TEST_EXPECT_EQUAL(group.count, strlen(group.members));
319                }
320               
321                for (size_t pos = 1; pos<group.count; ++pos) {
322                    TEST_EXPECT_LESS(group.members[pos-1], group.members[pos]);
323                }
324            }
325            else {
326                TEST_REJECT(group.count);
327            }
328        }
329    }
330}
331
332#endif // UNIT_TESTS
Note: See TracBrowser for help on using the repository browser.