source: tags/ms_r18q1/SL/ALILINK/TranslateRealign.cxx

Last change on this file was 16986, checked in by westram, 6 years ago

Update: continued by [17178]
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 76.9 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : TranslateRealign.cxx                              //
4//   Purpose   : Translate and realign                             //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include <TranslateRealign.h>
12#include <Translate.hxx>
13#include <AP_codon_table.hxx>
14#include <AP_pro_a_nucs.hxx>
15#include <aw_question.hxx> // @@@ remove (this module should not ask questions!)
16#include <arb_progress.h>
17#include <arbdbt.h>
18#include <arb_defs.h>
19#include <string>
20
21#define ali_assert(cond) arb_assert(cond)
22
23template<typename T>
24class BufferPtr {
25    T *const  bstart;
26    T        *curr;
27public:
28    explicit BufferPtr(T *b) : bstart(b), curr(b) {}
29
30    const T* start() const { return bstart; }
31    size_t offset() const { return curr-bstart; }
32
33    T get() { return *curr++; }
34
35    void put(T c) { *curr++ = c; }
36    void put(T c1, T c2, T c3) { put(c1); put(c2); put(c3); }
37    void put(T c, size_t count) {
38        memset(curr, c, count*sizeof(T));
39        inc(count);
40    }
41    void copy(BufferPtr<const T>& source, size_t count) {
42        memcpy(curr, source, count*sizeof(T));
43        inc(count);
44        source.inc(count);
45    }
46
47    T operator[](int i) const {
48        ali_assert(i>=0 || size_t(-i)<=offset());
49        return curr[i];
50    }
51
52    operator const T*() const { return curr; }
53    operator T*() { return curr; }
54
55    void inc(int o) { curr += o; ali_assert(curr>=bstart); }
56
57    BufferPtr<T>& operator++() { curr++; return *this; }
58    BufferPtr<T>& operator--() { inc(-1); return *this; }
59};
60
61template<typename T>
62class SizedBufferPtr : public BufferPtr<T> {
63    size_t len;
64public:
65    SizedBufferPtr(T *b, size_t len_) : BufferPtr<T>(b), len(len_) {}
66    ~SizedBufferPtr() { ali_assert(valid()); }
67    bool valid() const { return this->offset()<=len; }
68    size_t restLength() const { ali_assert(valid()); return len-this->offset(); }
69    size_t length() const { return len; }
70};
71
72typedef SizedBufferPtr<const char> SizedReadBuffer;
73typedef SizedBufferPtr<char>       SizedWriteBuffer;
74
75// ----------------------------------
76//      Translate protein -> dna
77
78inline bool legal_ORF_pos(int p) { return p >= 0 && p<=2; }
79
80GB_ERROR ALI_translate_marked(GBDATA *gb_main, bool use_entries, bool save_entries, int selected_startpos, bool translate_all, const char *ali_source, const char *ali_dest) {
81    // if use_entries   == true -> use fields 'codon_start' and 'transl_table' for translation
82    //                             (selected_startpos and AWAR_PROTEIN_TYPE are only used if both fields are missing,
83    //                             if only one is missing, now an error occurs)
84    // if use_entries   == false -> always use selected_startpos and AWAR_PROTEIN_TYPE
85    // if translate_all == true -> a selected_startpos > 1 produces a leading 'X' in protein data
86    //                             (otherwise nucleotides in front of the starting pos are simply ignored)
87    // if selected_startpos == AUTODETECT_STARTPOS -> the start pos is chosen to minimise number of stop codons
88
89    ali_assert(legal_ORF_pos(selected_startpos) || selected_startpos == AUTODETECT_STARTPOS);
90
91    GB_ERROR  error   = NULp;
92    char     *to_free = NULp;
93
94    // check/create alignments
95    {
96        GBDATA *gb_source = GBT_get_alignment(gb_main, ali_source);
97        if (!gb_source) {
98            error = GBS_global_string("No valid source alignment (%s)", GB_await_error());
99        }
100        else {
101            GBDATA *gb_dest = GBT_get_alignment(gb_main, ali_dest);
102            if (!gb_dest) {
103                GB_clear_error();
104                const char *msg = GBS_global_string("You have not selected a destination alignment\n"
105                                                    "Shall I create one ('%s_pro') for you?", ali_source);
106                if (!aw_ask_sure("create_protein_ali", msg)) { // @@@ remove (pass answer as parameter and fail if needed)
107                    error = "Cancelled by user";
108                }
109                else {
110                    long slen = GBT_get_alignment_len(gb_main, ali_source);
111                    to_free   = GBS_global_string_copy("%s_pro", ali_source);
112                    ali_dest  = to_free;
113                    gb_dest   = GBT_create_alignment(gb_main, ali_dest, slen/3+1, 0, 1, "ami");
114
115                    if (!gb_dest) error = GB_await_error();
116                    else {
117                        char *fname = GBS_global_string_copy("%s/data", ali_dest);
118                        error       = GBT_add_new_changekey(gb_main, fname, GB_STRING);
119                        free(fname);
120                    }
121                }
122            }
123        }
124    }
125
126    int no_data             = 0;  // count species w/o data
127    int spec_no_transl_info = 0;  // counts species w/o or with illegal transl_table and/or codon_start
128    int count               = 0;  // count translated species
129    int stops               = 0;  // count overall stop codons
130    int selected_ttable     = -1;
131
132    if (!error) {
133        arb_progress progress("Translating", GBT_count_marked_species(gb_main));
134
135        bool table_used[AWT_CODON_TABLES];
136        memset(table_used, 0, sizeof(table_used));
137        selected_ttable = *GBT_read_int(gb_main, AWAR_PROTEIN_TYPE); // read selected table
138
139        if (use_entries) {
140            for (GBDATA *gb_species = GBT_first_marked_species(gb_main);
141                 gb_species && !error;
142                 gb_species = GBT_next_marked_species(gb_species))
143            {
144                int arb_table, codon_start;
145                error = AWT_getTranslationInfo(gb_species, arb_table, codon_start);
146
147                if (!error) {
148                    if (arb_table == -1) arb_table = selected_ttable; // no transl_table entry -> default to selected standard code
149                    table_used[arb_table] = true;
150                }
151            }
152        }
153        else {
154            table_used[selected_ttable] = true; // and mark it used
155        }
156
157        for (int table = 0; table<AWT_CODON_TABLES && !error; ++table) {
158            if (!table_used[table]) continue;
159
160            for (GBDATA *gb_species = GBT_first_marked_species(gb_main);
161                 gb_species && !error;
162                 gb_species = GBT_next_marked_species(gb_species))
163            {
164                bool found_transl_info = false;
165                int  startpos          = selected_startpos;
166
167                if (use_entries) {  // if entries are used, test if field 'transl_table' matches current table
168                    int sp_arb_table, sp_codon_start;
169
170                    error = AWT_getTranslationInfo(gb_species, sp_arb_table, sp_codon_start);
171
172                    ali_assert(!error); // should already have been handled after first call to AWT_getTranslationInfo above
173
174                    if (sp_arb_table == -1) { // no table in DB
175                        ali_assert(sp_codon_start == -1);    // either both should be defined or none
176                        sp_arb_table   = selected_ttable;   // use selected translation table as default (if 'transl_table' field is missing)
177                        sp_codon_start = selected_startpos; // use selected codon startpos (if 'codon_start' field is missing)
178                    }
179                    else {
180                        ali_assert(sp_codon_start != -1); // either both should be defined or none
181                        found_transl_info = true;
182                        ali_assert(legal_ORF_pos(sp_codon_start));
183                    }
184
185                    if (sp_arb_table != table) continue; // species has not current transl_table
186
187                    startpos = sp_codon_start;
188                }
189
190                GBDATA *gb_source = GB_entry(gb_species, ali_source);
191                if (!gb_source) { ++no_data; }
192                else {
193                    GBDATA *gb_source_data = GB_entry(gb_source, "data");
194                    if (!gb_source_data) { ++no_data; }
195                    else {
196                        char *data = GB_read_string(gb_source_data);
197                        size_t  data_size = GB_read_string_count(gb_source_data);
198                        if (!data) {
199                            GB_print_error(); // cannot read data (ignore species)
200                            ++no_data;
201                        }
202                        else {
203                            if (!found_transl_info) ++spec_no_transl_info; // count species with missing info
204
205                            if (startpos == AUTODETECT_STARTPOS) {
206                                int   cn;
207                                int   stop_codons;
208                                int   least_stop_codons = -1;
209                                char* trial_data[3]     = {data, ARB_strdup(data), ARB_strdup(data)};
210
211                                for (cn = 0 ; cn < 3 ; cn++) {
212                                    stop_codons = AWT_pro_a_nucs_convert(table, trial_data[cn], data_size, cn, translate_all, false, false, NULp); // do the translation
213
214                                    if ((stop_codons < least_stop_codons) ||
215                                        (least_stop_codons == -1))
216                                    {
217                                        least_stop_codons = stop_codons;
218                                        startpos          = cn;
219                                    }
220                                }
221
222                                for (cn = 0 ; cn < 3 ; cn++) {
223                                    if (cn != startpos) {
224                                        free(trial_data[cn]);
225                                    }
226                                }
227
228                                data   = trial_data[startpos];
229                                stops += least_stop_codons;
230
231                            }
232                            else {
233                                stops += AWT_pro_a_nucs_convert(table, data, data_size, startpos, translate_all, false, false, NULp); // do the translation
234                            }
235
236                            ali_assert(legal_ORF_pos(startpos));
237                            ++count;
238
239                            GBDATA *gb_dest_data     = GBT_add_data(gb_species, ali_dest, "data", GB_STRING);
240                            if (!gb_dest_data) error = GB_await_error();
241                            else    error            = GB_write_string(gb_dest_data, data);
242
243
244                            if (!error && save_entries && !found_transl_info) {
245                                error = AWT_saveTranslationInfo(gb_species, selected_ttable, startpos);
246                            }
247
248                            free(data);
249                        }
250                    }
251                }
252                progress.inc_and_check_user_abort(error);
253            }
254        }
255    }
256
257    if (!error) {
258        if (use_entries) { // use 'transl_table' and 'codon_start' fields ?
259            if (spec_no_transl_info) {
260                int embl_transl_table = AWT_arb_code_nr_2_embl_transl_table(selected_ttable);
261                GB_warning(GBS_global_string("%i taxa had no valid translation info (fields 'transl_table' and 'codon_start')\n"
262                                             "Defaults (%i and %i) have been used%s.",
263                                             spec_no_transl_info,
264                                             embl_transl_table, selected_startpos+1,
265                                             save_entries ? " and written to DB entries" : ""));
266            }
267            else { // all entries were present
268                GB_warning("codon_start and transl_table entries were found for all translated taxa");
269            }
270        }
271
272        if (no_data>0) {
273            GB_warning(GBS_global_string("%i taxa had no data in '%s'", no_data, ali_source));
274        }
275        if ((count+no_data) == 0) {
276            GB_warning("Please mark species to translate");
277        }
278        else {
279            GB_warning(GBS_global_string("%i taxa converted\n  %f stops per sequence found",
280                                         count, (double)stops/(double)count));
281        }
282    }
283
284    free(to_free);
285
286    return error;
287}
288
289// -----------------------------------------------------------
290//      Realign a dna alignment to a given protein source
291
292class Distributor {
293    int xcount;
294    int *dist;
295    int *left;
296
297    GB_ERROR error;
298
299    void fillFrom(int off) {
300        ali_assert(!error);
301        ali_assert(off<xcount);
302
303        do {
304            int leftX    = xcount-off;
305            int leftDNA  = left[off];
306            int minLeave = leftX-1;
307            int maxLeave = minLeave*3;
308            int minTake  = std::max(1, leftDNA-maxLeave);
309
310#if defined(ASSERTION_USED)
311            int maxTake  = std::min(3, leftDNA-minLeave);
312            ali_assert(minTake<=maxTake);
313#endif
314
315            dist[off]   = minTake;
316            left[off+1] = left[off]-dist[off];
317
318            off++;
319        } while (off<xcount);
320
321        ali_assert(left[xcount] == 0); // expect correct amount of dna has been used
322    }
323    bool incAt(int off) {
324        ali_assert(!error);
325        ali_assert(off<xcount);
326
327        if (dist[off] == 3) {
328            return false;
329        }
330
331        int leftX    = xcount-off;
332        int leftDNA  = left[off];
333        int minLeave = leftX-1;
334        int maxTake  = std::min(3, leftDNA-minLeave);
335
336        if (dist[off] == maxTake) {
337            return false;
338        }
339
340        dist[off]++;
341        left[off+1]--;
342        fillFrom(off+1);
343        return true;
344    }
345
346public:
347    Distributor(int xcount_, int dnacount) :
348        xcount(xcount_),
349        dist(new int[xcount]),
350        left(new int[xcount+1]),
351        error(NULp)
352    {
353        if (dnacount<xcount) {
354            error = "not enough nucleotides";
355        }
356        else if (dnacount>3*xcount) {
357            error = "too much nucleotides";
358        }
359        else {
360            left[0] = dnacount;
361            fillFrom(0);
362        }
363    }
364    Distributor(const Distributor& other)
365        : xcount(other.xcount),
366          dist(new int[xcount]),
367          left(new int[xcount+1]),
368          error(other.error)
369    {
370        memcpy(dist, other.dist, sizeof(*dist)*xcount);
371        memcpy(left, other.left, sizeof(*left)*(xcount+1));
372    }
373    DECLARE_ASSIGNMENT_OPERATOR(Distributor);
374    ~Distributor() {
375        delete [] dist;
376        delete [] left;
377    }
378
379    void reset() { *this = Distributor(xcount, left[0]); }
380
381    int operator[](int off) const {
382        ali_assert(!error);
383        ali_assert(off>=0 && off<xcount);
384        return dist[off];
385    }
386
387    int size() const { return xcount; }
388
389    GB_ERROR get_error() const { return error; }
390
391    bool next() {
392        for (int incPos = xcount-2; incPos>=0; --incPos) {
393            if (incAt(incPos)) return true;
394        }
395        return false;
396    }
397
398    bool mayFailTranslation() const {
399        for (int i = 0; i<xcount; ++i) {
400            if (dist[i] == 3) return true;
401        }
402        return false;
403    }
404    int get_score() const {
405        // rates balanced distributions high
406        int score = 1;
407        for (int i = 0; i<xcount; ++i) { // LOOP_VECTORIZED=4
408            score *= dist[i];
409        }
410        return score + 6 - dist[0] - dist[xcount-1]; // prefer border positions with less nucs
411    }
412
413    bool translates_to_Xs(const char *dna, TransTables allowed, TransTables& remaining) const {
414        /*! checks whether distribution of 'dna' translates to X's
415         * @param dna compressed dna
416         * @param allowed allowed translation tables
417         * @param remaining remaining translation tables
418         * @return true if 'dna' translates to X's
419         */
420        bool translates = true;
421        int  off        = 0;
422        for (int p = 0; translates && p<xcount; off += dist[p++]) {
423            if (dist[p] == 3) {
424                TransTables this_remaining;
425                translates = AWT_is_codon('X', dna+off, allowed, this_remaining);
426                if (translates) {
427                    ali_assert(this_remaining.is_subset_of(allowed));
428                    allowed = this_remaining;
429                }
430            }
431        }
432        if (translates) remaining = allowed;
433        return translates;
434    }
435};
436
437inline bool isGap(char c) { return c == '-' || c == '.'; }
438
439using std::string;
440
441class FailedAt {
442    string             reason;
443    RefPtr<const char> at_prot; // points into aligned protein seq
444    RefPtr<const char> at_dna;  // points into compressed seq
445
446    int cmp(const FailedAt& other) const {
447        ptrdiff_t d = at_prot - other.at_prot;
448        if (!d)   d = at_dna - other.at_dna;
449        return d<0 ? -1 : d>0 ? 1 : 0;
450    }
451
452public:
453    FailedAt() :
454        at_prot(NULp),
455        at_dna(NULp)
456    {}
457    FailedAt(GB_ERROR reason_, const char *at_prot_, const char *at_dna_)
458        : reason(reason_),
459          at_prot(at_prot_),
460          at_dna(at_dna_)
461    {
462        ali_assert(reason_);
463    }
464
465    GB_ERROR why() const { return reason.empty() ? NULp : reason.c_str(); }
466    const char *protein_at() const { return at_prot; }
467    const char *dna_at() const { return at_dna; }
468
469    operator bool() const { return !reason.empty(); }
470
471    void add_prefix(const char *prefix) {
472        ali_assert(!reason.empty());
473        reason = string(prefix)+reason;
474    }
475
476    bool operator>(const FailedAt& other) const { return cmp(other)>0; }
477};
478
479class RealignAttempt : virtual Noncopyable {
480    TransTables           allowed;
481    SizedReadBuffer       compressed_dna;
482    BufferPtr<const char> aligned_protein;
483    SizedWriteBuffer      target_dna;
484    FailedAt              fail;
485    bool                  cutoff_dna;
486
487    void perform();
488
489    bool sync_behind_X_and_distribute(const int x_count, char *const x_start, const char *const x_start_prot);
490
491public:
492    RealignAttempt(const TransTables& allowed_, const char *compressed_dna_, size_t compressed_len_, const char *aligned_protein_, char *target_dna_, size_t target_len_, bool cutoff_dna_)
493        : allowed(allowed_),
494          compressed_dna(compressed_dna_, compressed_len_),
495          aligned_protein(aligned_protein_),
496          target_dna(target_dna_, target_len_),
497          cutoff_dna(cutoff_dna_)
498    {
499        ali_assert(aligned_protein[0]);
500        perform();
501    }
502
503    const TransTables& get_remaining_tables() const { return allowed; }
504    const FailedAt& failed() const { return fail; }
505};
506
507static GB_ERROR distribute_xdata(SizedReadBuffer& dna, size_t xcount, char *xtarget_, bool gap_before, bool gap_after, const TransTables& allowed, TransTables& remaining) {
508    /*! distributes 'dna' to marked X-positions
509     * @param xtarget destination buffer (target positions are marked with '!')
510     * @param xcount number of X's encountered
511     * @param gap_before true if resulting realignment has a gap or the start of alignment before the X-positions
512     * @param gap_after analog to 'gap_before'
513     * @param allowed allowed translation tables
514     * @param remaining remaining allowed translation tables (with those tables disabled for which no distribution possible)
515     * @return error if dna distribution wasn't possible
516     */
517
518    BufferPtr<char> xtarget(xtarget_);
519    Distributor     dist(xcount, dna.length());
520    GB_ERROR        error = dist.get_error();
521    if (!error) {
522        Distributor best(dist);
523        TransTables best_remaining = allowed;
524
525        while (dist.next()) {
526            if (dist.get_score() > best.get_score()) {
527                if (!dist.mayFailTranslation() || best.mayFailTranslation()) {
528                    best           = dist;
529                    best_remaining = allowed;
530                    ali_assert(best_remaining.is_subset_of(allowed));
531                }
532            }
533        }
534
535        if (best.mayFailTranslation()) {
536            TransTables curr_remaining;
537            if (best.translates_to_Xs(dna, allowed, curr_remaining)) {
538                best_remaining = curr_remaining;
539                ali_assert(best_remaining.is_subset_of(allowed));
540            }
541            else {
542                ali_assert(!error);
543                error = "no translating X-distribution found";
544                dist.reset();
545                do {
546                    if (dist.translates_to_Xs(dna, allowed, curr_remaining)) {
547                        best           = dist;
548                        best_remaining = curr_remaining;
549                        error          = NULp;
550                        ali_assert(best_remaining.is_subset_of(allowed));
551                        break;
552                    }
553                } while (dist.next());
554
555                while (dist.next()) {
556                    if (dist.get_score() > best.get_score()) {
557                        if (dist.translates_to_Xs(dna, allowed, curr_remaining)) {
558                            best           = dist;
559                            best_remaining = curr_remaining;
560                            ali_assert(best_remaining.is_subset_of(allowed));
561                        }
562                    }
563                }
564            }
565        }
566
567        if (!error) { // now really distribute nucs
568            for (int x = 0; x<best.size(); ++x) {
569                while (xtarget[0] != '!') {
570                    ali_assert(xtarget[1] && xtarget[2]); // buffer overflow
571                    xtarget.inc(3);
572                }
573
574                switch (best[x]) {
575                    case 2: {
576                        enum { UNDECIDED, SPREAD, LEFT, RIGHT } mode = UNDECIDED;
577
578                        bool is_1st_X  = xtarget.offset() == 0;
579                        bool gaps_left = is_1st_X ? gap_before : isGap(xtarget[-1]);
580
581                        if (gaps_left) mode = LEFT;
582                        else { // definitely has no gap left!
583                            bool is_last_X  = x == best.size()-1;
584                            int  next_nucs  = is_last_X ? 0 : best[x+1];
585                            bool gaps_right = isGap(xtarget[3]) || next_nucs == 1 || (is_last_X && gap_after);
586
587                            if (gaps_right) mode = RIGHT;
588                            else {
589                                bool nogaps_right = next_nucs == 3 || (is_last_X && !gap_after);
590                                if (nogaps_right) { // we know, we have NO adjacent gaps
591                                    mode = is_last_X ? LEFT : (is_1st_X ? RIGHT : SPREAD);
592                                }
593                                else {
594                                    ali_assert(!is_last_X);
595                                    mode = RIGHT; // forward problem to next X
596                                }
597                            }
598                        }
599
600                        char d1 = dna.get();
601                        char d2 = dna.get();
602
603                        switch (mode) {
604                            case UNDECIDED: ali_assert(0); FALLTHROUGH; // in NDEBUG
605                            case SPREAD: xtarget.put(d1,  '-', d2);  break;
606                            case LEFT:   xtarget.put(d1,  d2,  '-'); break;
607                            case RIGHT:  xtarget.put('-', d1,  d2);  break;
608                        }
609
610                        break;
611                    }
612                    case 1: xtarget.put('-', dna.get(), '-'); break;
613                    case 3: xtarget.copy(dna, 3); break;
614                    default: ali_assert(0); break;
615                }
616                ali_assert(dna.valid());
617            }
618
619            ali_assert(!error);
620            remaining = best_remaining;
621            ali_assert(remaining.is_subset_of(allowed));
622        }
623    }
624
625    return error;
626}
627
628bool RealignAttempt::sync_behind_X_and_distribute(const int x_count, char *const x_start, const char *const x_start_prot) {
629    /*! brute-force search for sync behind 'X' and distribute dna onto X positions
630     * @param x_count number of X encountered
631     * @param x_start dna read position
632     * @param x_start_prot protein read position
633     * @return true if sync and distribution succeed
634     */
635
636    bool complete = false;
637
638    ali_assert(!failed());
639    ali_assert(aligned_protein.offset()>0);
640    const char p = aligned_protein[-1];
641
642    size_t compressed_rest_len = compressed_dna.restLength();
643    ali_assert(strlen(compressed_dna) == compressed_rest_len);
644
645    size_t min_dna = x_count;
646    size_t max_dna = std::min(size_t(x_count)*3, compressed_rest_len);
647
648    if (min_dna>max_dna) {
649        fail = FailedAt("not enough nucs for X's at sequence end", x_start_prot, compressed_dna);
650    }
651    else if (p) {
652        FailedAt foremost;
653        size_t   target_rest_len = target_dna.restLength();
654
655        for (size_t x_dna = min_dna; x_dna<=max_dna; ++x_dna) { // prefer low amounts of used dna
656            const char *dna_rest     = compressed_dna + x_dna;
657            size_t      dna_rest_len = compressed_rest_len     - x_dna;
658
659            ali_assert(strlen(dna_rest) == dna_rest_len);
660            ali_assert(compressed_rest_len>=x_dna);
661
662            RealignAttempt attemptRest(allowed, dna_rest, dna_rest_len, aligned_protein-1, target_dna, target_rest_len, cutoff_dna);
663            FailedAt       restFailed = attemptRest.failed();
664
665            if (!restFailed) {
666                SizedReadBuffer distrib_dna(compressed_dna, x_dna);
667
668                bool has_gap_before = x_start == target_dna.start() ? true : isGap(x_start[-1]);
669                bool has_gap_after  = isGap(dna_rest[0]);
670
671                TransTables remaining;
672                GB_ERROR    disterr = distribute_xdata(distrib_dna, x_count, x_start, has_gap_before, has_gap_after, attemptRest.get_remaining_tables(), remaining);
673                if (disterr) {
674                    restFailed = FailedAt(disterr, x_start_prot, dna_rest); // prot=start of Xs; dna=start of sync (behind Xs)
675                }
676                else {
677                    ali_assert(remaining.is_subset_of(allowed));
678                    ali_assert(remaining.is_subset_of(attemptRest.get_remaining_tables()));
679                    allowed = remaining;
680                }
681            }
682
683            if (restFailed) {
684                if (restFailed > foremost) foremost = restFailed; // track "best" failure (highest fail position)
685            }
686            else { // success
687                foremost = FailedAt();
688                complete = true;
689                break; // use first success and return
690            }
691        }
692
693        if (foremost) {
694            ali_assert(!complete);
695            fail = foremost;
696            if (!strstr(fail.why(), "Sync behind 'X'")) { // do not spam repetitive sync-failures
697                fail.add_prefix("Sync behind 'X' failed foremost with: ");
698            }
699        }
700        else {
701            ali_assert(complete);
702        }
703    }
704    else {
705        GB_ERROR fail_reason = "internal error: no distribution attempted";
706        ali_assert(min_dna>0);
707        size_t x_dna;
708        for (x_dna = max_dna; x_dna>=min_dna; --x_dna) {     // prefer high amounts of dna
709            SizedReadBuffer append_dna(compressed_dna, x_dna);
710            TransTables     remaining;
711            fail_reason = distribute_xdata(append_dna, x_count, x_start, false, true, allowed, remaining);
712            if (!fail_reason) { // found distribution -> done
713                ali_assert(remaining.is_subset_of(allowed));
714                allowed = remaining;
715                break;
716            }
717        }
718
719        if (fail_reason) {
720            fail = FailedAt(fail_reason, x_start_prot+1, compressed_dna); // report error at start of X's
721        }
722        else {
723            fail = FailedAt(); // clear
724            compressed_dna.inc(x_dna);
725        }
726    }
727
728    ali_assert(implicated(complete, allowed.any()));
729
730    return complete;
731}
732
733void RealignAttempt::perform() {
734    bool complete = false; // set to true, if recursive attempt succeeds
735
736    while (char p = toupper(aligned_protein.get())) {
737        if (p=='X') { // one X represents 1 to 3 DNAs (normally 1 or 2, but 'NNN' translates to 'X')
738            char       *x_start      = target_dna;
739            const char *x_start_prot = aligned_protein-1;
740            int         x_count      = 0;
741
742            for (;;) {
743                if      (p=='X')   { x_count++; target_dna.put('!', 3); } // fill X space with marker
744                else if (isGap(p)) target_dna.put(p, 3);
745                else break;
746
747                p = toupper(aligned_protein.get());
748            }
749
750            ali_assert(x_count);
751            ali_assert(!complete);
752            complete = sync_behind_X_and_distribute(x_count, x_start, x_start_prot);
753            if (!complete && !failed()) {
754                if (p) { // not all proteins were processed
755                    fail = FailedAt("internal error", aligned_protein-1, compressed_dna);
756                    ali_assert(0);
757                }
758            }
759            break; // done
760        }
761
762        if (isGap(p)) target_dna.put(p, 3);
763        else {
764            TransTables remaining;
765            size_t      compressed_rest_len = compressed_dna.restLength();
766
767            if (compressed_rest_len<3) {
768                fail = FailedAt(GBS_global_string("not enough nucs left for codon of '%c'", p), aligned_protein-1, compressed_dna);
769            }
770            else {
771                ali_assert(strlen(compressed_dna) == compressed_rest_len);
772                ali_assert(compressed_rest_len >= 3);
773                const char *why_fail;
774                if (!AWT_is_codon(p, compressed_dna, allowed, remaining, &why_fail)) {
775                    fail = FailedAt(why_fail, aligned_protein-1, compressed_dna);
776                }
777            }
778
779            if (failed()) break;
780
781            ali_assert(remaining.is_subset_of(allowed));
782            allowed = remaining;
783            target_dna.copy(compressed_dna, 3);
784        }
785    }
786
787    ali_assert(compressed_dna.valid());
788
789    if (!failed() && !complete) {
790        while (target_dna.offset()>0 && isGap(target_dna[-1])) --target_dna; // remove terminal gaps
791
792        if (!cutoff_dna) { // append leftover dna-data (data w/o corresponding aa)
793            size_t compressed_rest_len = compressed_dna.restLength();
794            size_t target_rest_len = target_dna.restLength();
795            if (compressed_rest_len<=target_rest_len) {
796                target_dna.copy(compressed_dna, compressed_rest_len);
797            }
798            else {
799                fail = FailedAt(GBS_global_string("too much trailing DNA (%zu nucs, but only %zu columns left)",
800                                                  compressed_rest_len, target_rest_len),
801                                aligned_protein-1, compressed_dna);
802            }
803        }
804
805        if (!failed()) target_dna.put('.', target_dna.restLength()); // fill rest of sequence with dots
806        *target_dna = 0;
807    }
808
809#if defined(ASSERTION_USED)
810    if (!failed()) {
811        ali_assert(strlen(target_dna.start()) == target_dna.length());
812    }
813#endif
814}
815
816inline char *unalign(const char *data, size_t len, size_t& compressed_len) {
817    // removes gaps from sequence
818    char *compressed = ARB_alloc<char>(len+1);
819    compressed_len = 0;
820    for (size_t p = 0; p<len && data[p]; ++p) {
821        if (!isGap(data[p])) {
822            compressed[compressed_len++] = data[p];
823        }
824    }
825    compressed[compressed_len] = 0;
826    return compressed;
827}
828
829class Realigner {
830    const char *ali_source;
831    const char *ali_dest;
832
833    size_t ali_len;        // of ali_dest
834    size_t needed_ali_len; // >ali_len if ali_dest is too short; 0 otherwise
835
836    const char *fail_reason;
837
838    GB_ERROR annotate_fail_position(const FailedAt& failed, const char *source, const char *dest, const char *compressed_dest) {
839        int source_fail_pos = failed.protein_at() - source;
840        int dest_fail_pos   = 0;
841        {
842            int fail_d_base_count = failed.dna_at() - compressed_dest;
843
844            const char *dp = dest;
845
846            for (;;) {
847                char c = *dp++;
848
849                if (!c) { // failure at end of sequence
850                    dest_fail_pos++; // report position behind last non-gap
851                    break;
852                }
853                if (!isGap(c)) {
854                    dest_fail_pos = (dp-1)-dest;
855                    if (!fail_d_base_count) break;
856                    fail_d_base_count--;
857                }
858            }
859        }
860        return GBS_global_string("%s at %s:%i / %s:%i",
861                                 failed.why(),
862                                 ali_source, info2bio(source_fail_pos),
863                                 ali_dest, info2bio(dest_fail_pos));
864    }
865
866
867    static void calc_needed_dna(const char *prot, size_t len, size_t& minDNA, size_t& maxDNA) {
868        minDNA = maxDNA = 0;
869        for (size_t o = 0; o<len; ++o) {
870            char p = toupper(prot[o]);
871            if (p == 'X') {
872                minDNA += 1;
873                maxDNA += 3;
874            }
875            else if (!isGap(p)) {
876                minDNA += 3;
877                maxDNA += 3;
878            }
879        }
880    }
881    static size_t countLeadingGaps(const char *buffer) {
882        size_t gaps = 0;
883        for (int o = 0; isGap(buffer[o]); ++o) ++gaps;
884        return gaps;
885    }
886
887public:
888    Realigner(const char *ali_source_, const char *ali_dest_, size_t ali_len_)
889        : ali_source(ali_source_),
890          ali_dest(ali_dest_),
891          ali_len(ali_len_),
892          needed_ali_len(0)
893    {
894        clear_failure();
895    }
896
897    size_t get_needed_dest_alilen() const { return needed_ali_len; }
898
899    void set_failure(const char *reason) { fail_reason = reason; }
900    void clear_failure() { fail_reason = NULp; }
901
902    const char *failure() const { return fail_reason; }
903
904    char *realign_seq(TransTables& allowed, const char *const source, size_t source_len, const char *const dest, size_t dest_len, bool cutoff_dna) {
905        ali_assert(!failure());
906
907        size_t  wanted_ali_len = source_len*3;
908        char   *buffer         = NULp;
909
910        if (ali_len<wanted_ali_len) {
911            fail_reason = GBS_global_string("Alignment '%s' is too short (increase its length to %zu)", ali_dest, wanted_ali_len);
912            if (wanted_ali_len>needed_ali_len) needed_ali_len = wanted_ali_len;
913        }
914        else {
915            // compress destination DNA (=remove align-characters):
916            size_t  compressed_len;
917            char   *compressed_dest = unalign(dest, dest_len, compressed_len);
918
919            ARB_alloc(buffer, ali_len+1);
920
921            RealignAttempt attempt(allowed, compressed_dest, compressed_len, source, buffer, ali_len, cutoff_dna);
922            FailedAt       failed = attempt.failed();
923
924            if (failed) {
925                // test for superfluous DNA at sequence start
926                size_t min_dna, max_dna;
927                calc_needed_dna(source, source_len, min_dna, max_dna);
928
929                if (min_dna<compressed_len) { // we have more DNA than we need
930                    size_t extra_dna = compressed_len-min_dna;
931                    for (size_t skip = 1; skip<=extra_dna; ++skip) {
932                        RealignAttempt attemptSkipped(allowed, compressed_dest+skip, compressed_len-skip, source, buffer, ali_len, cutoff_dna);
933                        if (!attemptSkipped.failed()) {
934                            failed = FailedAt(); // clear
935                            if (!cutoff_dna) {
936                                size_t start_gaps = countLeadingGaps(buffer);
937                                if (start_gaps<skip) {
938                                    failed = FailedAt(GBS_global_string("Not enough gaps to place %zu extra nucs at start of sequence",
939                                                                        skip), source, compressed_dest);
940                                }
941                                else { // success
942                                    memcpy(buffer+(start_gaps-skip), compressed_dest, skip); // copy-in skipped dna
943                                }
944                            }
945                            if (!failed) {
946                                ali_assert(attempt.get_remaining_tables().is_subset_of(allowed));
947                                allowed = attemptSkipped.get_remaining_tables();
948                            }
949                            break; // no need to skip more dna, when we already have too few leading gaps
950                        }
951                    }
952                }
953            }
954            else {
955                ali_assert(attempt.get_remaining_tables().is_subset_of(allowed));
956                allowed = attempt.get_remaining_tables();
957            }
958
959            if (failed) {
960                fail_reason = annotate_fail_position(failed, source, dest, compressed_dest);
961                freenull(buffer);
962            }
963            free(compressed_dest);
964        }
965        ali_assert(contradicted(buffer, fail_reason));
966        return buffer;
967    }
968};
969
970struct Data : virtual Noncopyable {
971    GBDATA *gb_data;
972    char   *data;
973    size_t  len;
974    char   *error;
975
976    Data(GBDATA *gb_species, const char *aliName) :
977        gb_data(NULp),
978        data(NULp),
979        len(0),
980        error(NULp)
981    {
982        GBDATA *gb_ali = GB_entry(gb_species, aliName);
983        if (gb_ali) {
984            gb_data = GB_entry(gb_ali, "data");
985            if (gb_data) {
986                data          = GB_read_string(gb_data);
987                if (data) len = GB_read_string_count(gb_data);
988                else error    = ARB_strdup(GB_await_error());
989                return;
990            }
991        }
992        error = GBS_global_string_copy("No data in alignment '%s'", aliName);
993    }
994    ~Data() {
995        free(data);
996        free(error);
997    }
998};
999
1000GB_ERROR ALI_realign_marked(GBDATA *gb_main, const char *ali_source, const char *ali_dest, size_t& neededLength, bool unmark_succeeded, bool cutoff_dna) {
1001    /*! realigns DNA alignment of marked sequences according to their protein alignment
1002     * @param ali_source protein source alignment
1003     * @param ali_dest modified DNA alignment
1004     * @param neededLength result: minimum alignment length needed in ali_dest (if too short) or 0 if long enough
1005     * @param unmark_succeeded unmark all species that were successfully realigned
1006     */
1007    AP_initialize_codon_tables();
1008
1009    ali_assert(GB_get_transaction_level(gb_main) == 0);
1010    GB_transaction ta(gb_main); // do not abort (otherwise sth goes wrong with species marks)
1011
1012    {
1013        GBDATA *gb_source = GBT_get_alignment(gb_main, ali_source); if (!gb_source) return "Please select a valid source alignment";
1014        GBDATA *gb_dest   = GBT_get_alignment(gb_main, ali_dest);   if (!gb_dest)   return "Please select a valid destination alignment";
1015    }
1016
1017    if (GBT_get_alignment_type(gb_main, ali_source) != GB_AT_AA)  return "Invalid source alignment type";
1018    if (GBT_get_alignment_type(gb_main, ali_dest)   != GB_AT_DNA) return "Invalid destination alignment type";
1019
1020    long ali_len = GBT_get_alignment_len(gb_main, ali_dest);
1021    ali_assert(ali_len>0);
1022
1023    GB_ERROR error = NULp;
1024
1025    long no_of_marked_species    = GBT_count_marked_species(gb_main);
1026    long no_of_realigned_species = 0; // count successfully realigned species
1027
1028    arb_progress progress("Re-aligner", no_of_marked_species);
1029    progress.auto_subtitles("Re-aligning species");
1030
1031    Realigner realigner(ali_source, ali_dest, ali_len);
1032
1033    for (GBDATA *gb_species = GBT_first_marked_species(gb_main);
1034         !error && gb_species;
1035         gb_species = GBT_next_marked_species(gb_species))
1036    {
1037        realigner.clear_failure();
1038
1039        Data source(gb_species, ali_source);
1040        Data dest(gb_species, ali_dest);
1041
1042        if      (source.error) realigner.set_failure(source.error);
1043        else if (dest.error)   realigner.set_failure(dest.error);
1044
1045        if (!realigner.failure()) {
1046            TransTables allowed; // default: all translation tables allowed
1047#if defined(ASSERTION_USED)
1048            bool has_valid_translation_info = false;
1049#endif
1050            {
1051                int arb_transl_table, codon_start;
1052                GB_ERROR local_error = AWT_getTranslationInfo(gb_species, arb_transl_table, codon_start);
1053                if (local_error) {
1054                    realigner.set_failure(GBS_global_string("Error while reading 'transl_table' (%s)", local_error));
1055                }
1056                else if (arb_transl_table >= 0) {
1057                    // we found a 'transl_table' entry -> restrict used code to the code stored there
1058                    allowed.forbidAllBut(arb_transl_table);
1059#if defined(ASSERTION_USED)
1060                    has_valid_translation_info = true;
1061#endif
1062                }
1063            }
1064
1065            if (!realigner.failure()) {
1066                char *buffer = realigner.realign_seq(allowed, source.data, source.len, dest.data, dest.len, cutoff_dna);
1067                if (buffer) { // re-alignment successful
1068                    error = GB_write_string(dest.gb_data, buffer);
1069
1070                    if (!error) {
1071                        int explicit_table_known = allowed.explicit_table();
1072
1073                        if (explicit_table_known >= 0) { // we know the exact code -> write codon_start and transl_table
1074                            const int codon_start  = 0; // by definition (after realignment)
1075                            error = AWT_saveTranslationInfo(gb_species, explicit_table_known, codon_start);
1076                        }
1077#if defined(ASSERTION_USED)
1078                        else { // we dont know the exact code -> can only happen if species has no translation info
1079                            ali_assert(allowed.any()); // bug in realigner
1080                            ali_assert(!has_valid_translation_info);
1081                        }
1082#endif
1083                    }
1084                    free(buffer);
1085                    if (!error && unmark_succeeded) GB_write_flag(gb_species, 0);
1086                }
1087            }
1088        }
1089
1090        if (realigner.failure()) {
1091            ali_assert(!error);
1092            GB_warningf("Automatic re-align failed for '%s'\nReason: %s", GBT_read_name(gb_species), realigner.failure());
1093        }
1094        else if (!error) {
1095            no_of_realigned_species++;
1096        }
1097
1098        progress.inc_and_check_user_abort(error);
1099    }
1100
1101    neededLength = realigner.get_needed_dest_alilen();
1102
1103    if (no_of_marked_species == 0) {
1104        GB_warning("Please mark some species to realign them");
1105    }
1106    else if (no_of_realigned_species != no_of_marked_species) {
1107        long failed = no_of_marked_species-no_of_realigned_species;
1108        ali_assert(failed>0);
1109        if (no_of_realigned_species) {
1110            GB_warningf("%li marked species failed to realign (%li succeeded)", failed, no_of_realigned_species);
1111        }
1112        else {
1113            GB_warning("All marked species failed to realign");
1114        }
1115    }
1116
1117    if (error) progress.done();
1118    else error = GBT_check_data(gb_main,ali_dest);
1119
1120    return error;
1121}
1122
1123
1124// --------------------------------------------------------------------------------
1125
1126#ifdef UNIT_TESTS
1127#ifndef TEST_UNIT_H
1128#include <test_unit.h>
1129#endif
1130
1131#include <arb_handlers.h>
1132
1133static std::string msgs;
1134
1135static void msg_to_string(const char *msg) {
1136    msgs += msg;
1137    msgs += '\n';
1138}
1139
1140static const char *translation_info(GBDATA *gb_species) {
1141    int      arb_transl_table;
1142    int      codon_start;
1143    GB_ERROR error = AWT_getTranslationInfo(gb_species, arb_transl_table, codon_start);
1144
1145    static SmartCharPtr result;
1146
1147    if (error) result = GBS_global_string_copy("Error: %s", error);
1148    else       result = GBS_global_string_copy("t=%i,cs=%i", arb_transl_table, codon_start);
1149
1150    return &*result;
1151}
1152
1153static arb_handlers test_handlers = {
1154    msg_to_string,
1155    msg_to_string,
1156    msg_to_string,
1157    active_arb_handlers->status,
1158};
1159
1160#define DNASEQ(name) GB_read_char_pntr(GBT_find_sequence(GBT_find_species(gb_main, name), "ali_dna"))
1161#define PROSEQ(name) GB_read_char_pntr(GBT_find_sequence(GBT_find_species(gb_main, name), "ali_pro"))
1162
1163#define TRANSLATION_INFO(name) translation_info(GBT_find_species(gb_main, name))
1164
1165void TEST_realign() {
1166    arb_handlers *old_handlers = active_arb_handlers;
1167    ARB_install_handlers(test_handlers);
1168
1169    GB_shell  shell;
1170    GBDATA   *gb_main = GB_open("TEST_realign.arb", "rw");
1171
1172    arb_suppress_progress here;
1173    enum TransResult { SAME, CHANGED };
1174
1175    {
1176        GB_ERROR error;
1177        size_t   neededLength = 0;
1178
1179        {
1180            struct transinfo_check {
1181                const char  *species_name;
1182                const char  *old_info;
1183                TransResult  changed;
1184                const char  *new_info;
1185            };
1186
1187            transinfo_check info[] = {
1188                { "BctFra12", "t=0,cs=1",  SAME,    NULp        }, // fails -> unchanged
1189                { "CytLyti6", "t=9,cs=1",  CHANGED, "t=9,cs=0"  },
1190                { "TaxOcell", "t=14,cs=1", CHANGED, "t=14,cs=0" },
1191                { "StrRamo3", "t=0,cs=1",  SAME,    NULp        }, // fails -> unchanged
1192                { "StrCoel9", "t=0,cs=0",  SAME,    NULp        }, // already correct
1193                { "MucRacem", "t=0,cs=1",  CHANGED, "t=0,cs=0"  },
1194                { "MucRace2", "t=0,cs=1",  CHANGED, "t=0,cs=0"  },
1195                { "MucRace3", "t=0,cs=0",  SAME,    NULp        }, // fails -> unchanged
1196                { "AbdGlauc", "t=0,cs=0",  SAME,    NULp        }, // already correct
1197                { "CddAlbic", "t=0,cs=0",  SAME,    NULp        }, // already correct
1198
1199                { NULp, NULp, SAME, NULp }
1200            };
1201
1202            {
1203                GB_transaction ta(gb_main);
1204
1205                for (int i = 0; info[i].species_name; ++i) {
1206                    const transinfo_check& I = info[i];
1207                    TEST_ANNOTATE(I.species_name);
1208                    TEST_EXPECT_EQUAL(TRANSLATION_INFO(I.species_name), I.old_info);
1209                }
1210            }
1211            TEST_ANNOTATE(NULp);
1212
1213            msgs  = "";
1214            error = ALI_realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1215            TEST_EXPECT_NO_ERROR(error);
1216            TEST_EXPECT_EQUAL(msgs,
1217                              "Automatic re-align failed for 'BctFra12'\nReason: not enough nucs for X's at sequence end at ali_pro:40 / ali_dna:109\n" // new correct report (got no nucs for 1 X)
1218                              "Automatic re-align failed for 'StrRamo3'\nReason: not enough nucs for X's at sequence end at ali_pro:36 / ali_dna:106\n" // new correct report (got 3 nucs for 4 Xs)
1219                              "Automatic re-align failed for 'MucRace3'\nReason: Sync behind 'X' failed foremost with: Not all IUPAC-combinations of 'NCC' translate to 'T' (for trans-table 0) at ali_pro:28 / ali_dna:78\n" // correct report
1220                              "3 marked species failed to realign (7 succeeded)\n"
1221                );
1222
1223            {
1224                GB_transaction ta(gb_main);
1225
1226                TEST_EXPECT_EQUAL(DNASEQ("BctFra12"),    "ATGGCTAAAGAGAAATTTGAACGTACCAAACCGCACGTAAACATTGGTACAATCGGTCACGTTGACCACGGTAAAACCACTTTGACTGCTGCTATCACTACTGTGTTG------------------"); // failed = > seq unchanged
1227                TEST_EXPECT_EQUAL(DNASEQ("CytLyti6"),    "-A-TGGCAAAGGAAACTTTTGATCGTTCCAAACCGCACTTAA---ATATAG---GTACTATTGGACACGTAGATCACGGTAAAACTACTTTAACTGCTGCTATTACAASAGTAT-T-----G....");
1228                TEST_EXPECT_EQUAL(DNASEQ("TaxOcell"),    "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G..........");
1229                TEST_EXPECT_EQUAL(DNASEQ("StrRamo3"),    "ATGTCCAAGACGGCATACGTGCGCACCAAACCGCATCTGAACATCGGCACGATGGGTCATGTCGACCACGGCAAGACCACGTTGACCGCCGCCATCACCAAGGTCCTC------------------"); // failed = > seq unchanged
1230                TEST_EXPECT_EQUAL(DNASEQ("StrCoel9"),    "ATGTCCAAGACGGCGTACGTCCGC-C--C--A-CC-TG--A----GGCACGATG-G-CC--C-GACCACGGCAAGACCACCCTGACCGCCGCCATCACCAAGGTC-C--T--------C.......");
1231                TEST_EXPECT_EQUAL(DNASEQ("MucRacem"),    "......ATGGGTAAAGAG---------AAGACTCACGTTAACGTCGTCGTCATTGGTCACGTCGATTCCGGTAAATCTACTACTACTGGTCACTTGATTTACAAGTGTGGTGGTATA-AA......");
1232                TEST_EXPECT_EQUAL(DNASEQ("MucRace2"),    "ATGGGTAAGGAG---------AAGACTCACGTTAACGTCGTCGTCATTGGTCACGTCGATTCCGGTAAATCTACTACTACTGGTCACTTGATTTACAAGTGTGGTGGT-ATNNNAT-AAA......");
1233                TEST_EXPECT_EQUAL(DNASEQ("MucRace3"),    "-----------ATGGGTAAAGAGAAGACTCACGTTRAYGTTGTCGTTATTGGTCACGTCRATTCCGGTAAGTCCACCNCCRCTGGTCACTTGATTTACAAGTGTGGTGGTATAA-A----------"); // failed = > seq unchanged
1234                TEST_EXPECT_EQUAL(DNASEQ("AbdGlauc"),    "ATGGGTAAA-G--A--A--A--A--G-AC--T-CACGTTAACGTCGTTGTCATTGGTCACGTCGATTCTGGTAAATCCACCACCACTGGTCATTTGATCTACAAGTGCGGTGGTATA-AA......");
1235                TEST_EXPECT_EQUAL(DNASEQ("CddAlbic"),    "ATG-GG-TAAA-GAA------------AAAACTCACGTTAACGTTGTTGTTATTGGTCACGTCGATTCCGGTAAATCTACTACCACCGGTCACTTAATTTACAAGTGTGGTGGTATA-AA......");
1236                // ------------------------------------- "123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123"
1237
1238                for (int i = 0; info[i].species_name; ++i) {
1239                    const transinfo_check& I = info[i];
1240                    TEST_ANNOTATE(I.species_name);
1241                    switch (I.changed) {
1242                        case SAME:
1243                            TEST_EXPECT_EQUAL(TRANSLATION_INFO(I.species_name), I.old_info);
1244                            TEST_EXPECT_NULL(static_cast<const char*>(I.new_info));
1245                            break;
1246                        case CHANGED:
1247                            TEST_EXPECT_EQUAL(TRANSLATION_INFO(I.species_name), I.new_info);
1248                            TEST_EXPECT_DIFFERENT(I.new_info, I.old_info);
1249                            break;
1250                    }
1251                }
1252                TEST_ANNOTATE(NULp);
1253            }
1254        }
1255
1256        // test translation of successful realignments (see previous section)
1257        {
1258            GB_transaction ta(gb_main);
1259
1260            struct translate_check {
1261                const char *species_name;
1262                const char *original_prot;
1263                TransResult retranslation;
1264                const char *changed_prot; // if changed by translation (NULp for SAME)
1265            };
1266
1267            translate_check trans[] = {
1268                { "CytLyti6", "XWQRKLLIVPNRT*-I*-VLLDT*ITVKLL*SSLLZZYX-X.",
1269                  CHANGED,    "XWQRKLLIVPNRT*-I*-VLLDT*ITVKLL*SSLLQZYX-X." }, // ok: one of the Zs near end translates to Q
1270                { "TaxOcell", "XG*SNFWPVQAARNHRHD--RSRGPRQBDSDRCYHHGAX-..",
1271                  CHANGED,    "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - only changes gaptype at EOS
1272                { "MucRacem", "..MGKE---KTHVNVVVIGHVDSGKSTTTGHLIYKCGGIX..", SAME, NULp },
1273                { "MucRace2", "MGKE---KTHVNVVVIGHVDSGKSTTTGHLIYKCGGXXXK--",
1274                  CHANGED,    "MGKE---KTHVNVVVIGHVDSGKSTTTGHLIYKCGGXXXK.." }, // ok - only changes gaptype at EOS
1275                { "AbdGlauc", "MGKXXXXXXXXHVNVVVIGHVDSGKSTTTGHLIYKCGGIX..", SAME, NULp },
1276                { "StrCoel9", "MSKTAYVRXXXXXX-GTMXXXDHGKTTLTAAITKVXX--X..", SAME, NULp },
1277                { "CddAlbic", "MXXXE----KTHVNVVVIGHVDSGKSTTTGHLIYKCGGIX..", SAME, NULp },
1278
1279                { NULp, NULp, SAME, NULp }
1280            };
1281
1282            // check original protein sequences
1283            for (int t = 0; trans[t].species_name; ++t) {
1284                const translate_check& T = trans[t];
1285                TEST_ANNOTATE(T.species_name);
1286                TEST_EXPECT_EQUAL(PROSEQ(T.species_name), T.original_prot);
1287            }
1288            TEST_ANNOTATE(NULp);
1289
1290            msgs  = "";
1291            error = ALI_translate_marked(gb_main, true, false, 0, true, "ali_dna", "ali_pro");
1292            TEST_EXPECT_NO_ERROR(error);
1293            TEST_EXPECT_EQUAL(msgs, "codon_start and transl_table entries were found for all translated taxa\n10 taxa converted\n  1.100000 stops per sequence found\n");
1294
1295            // check re-translated protein sequences
1296            for (int t = 0; trans[t].species_name; ++t) {
1297                const translate_check& T = trans[t];
1298                TEST_ANNOTATE(T.species_name);
1299                switch (T.retranslation) {
1300                    case SAME:
1301                        TEST_EXPECT_NULL(static_cast<const char*>(T.changed_prot));
1302                        TEST_EXPECT_EQUAL(PROSEQ(T.species_name), T.original_prot);
1303                        break;
1304                    case CHANGED:
1305                        TEST_REJECT_NULL(static_cast<const char*>(T.changed_prot));
1306                        TEST_EXPECT_DIFFERENT(T.original_prot, T.changed_prot);
1307                        TEST_EXPECT_EQUAL(PROSEQ(T.species_name), T.changed_prot);
1308                        break;
1309                }
1310            }
1311            TEST_ANNOTATE(NULp);
1312
1313            ta.close("dont commit");
1314        }
1315
1316        // -----------------------------
1317        //      provoke some errors
1318
1319        GBDATA *gb_TaxOcell;
1320        // unmark all but gb_TaxOcell
1321        {
1322            GB_transaction ta(gb_main);
1323
1324            gb_TaxOcell = GBT_find_species(gb_main, "TaxOcell");
1325            TEST_REJECT_NULL(gb_TaxOcell);
1326
1327            GBT_mark_all(gb_main, 0);
1328            GB_write_flag(gb_TaxOcell, 1);
1329        }
1330
1331        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1332
1333        // wrong alignment type
1334        {
1335            msgs  = "";
1336            error = ALI_realign_marked(gb_main, "ali_dna", "ali_pro", neededLength, false, false);
1337            TEST_EXPECT_ERROR_CONTAINS(error, "Invalid source alignment type");
1338            TEST_EXPECT_EQUAL(msgs, "");
1339        }
1340
1341        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1342
1343        GBDATA *gb_TaxOcell_amino;
1344        GBDATA *gb_TaxOcell_dna;
1345        {
1346            GB_transaction ta(gb_main);
1347            gb_TaxOcell_amino = GBT_find_sequence(gb_TaxOcell, "ali_pro");
1348            gb_TaxOcell_dna   = GBT_find_sequence(gb_TaxOcell, "ali_dna");
1349        }
1350        TEST_REJECT_NULL(gb_TaxOcell_amino);
1351        TEST_REJECT_NULL(gb_TaxOcell_dna);
1352
1353        // -----------------------------------------
1354        //      document some existing behavior
1355        {
1356            struct realign_check {
1357                const char  *seq;
1358                const char  *result;
1359                bool         cutoff;
1360                TransResult  retranslation;
1361                const char  *changed_prot; // if changed by translation (NULp for SAME)
1362            };
1363
1364            realign_check seq[] = {
1365                //"XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-.." // original aa sequence
1366                // { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..", "sdfjlksdjf" }, // templ
1367                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G..........", false, CHANGED, // original
1368                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - only changes gaptype at EOS
1369
1370                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHG.....", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCTG...........", false, CHANGED, // missing some AA at right end (extra DNA gets no longer truncated!)
1371                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - adds translation of extra DNA (DNA should never be modified by realigner!)
1372                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHG.....", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGT...............", true,  SAME, NULp }, // missing some AA at right end -> cutoff DNA
1373
1374                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYH-----..", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCTG...........", false, CHANGED,
1375                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - adds translation of extra DNA
1376                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCY---H....", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTAT---------CACCACGGTGCTG..", false, CHANGED, // rightmost possible position of 'H' (see failing test below)
1377                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCY---HHGAX" }, // ok - adds translation of extra DNA
1378
1379                { "---SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..", "-ATGGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G..........", false, CHANGED, // missing some AA at left end (extra DNA gets detected now)
1380                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - adds translation of extra DNA (start of alignment)
1381                { "...SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX...", ".........AGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G..........", true,  SAME, NULp }, // missing some AA at left end -> cutoff DNA
1382
1383
1384                { "XG*SNFXXXXXXAXXXNHRHDXXXXXXPRQNDSDRCYHHGAX", "AT-GGCTAAAGAAACTTT-TG-AC-CG-GT-CCAA-GCC-GC-ACGT-AAACATCGGCACGAT-CG-GT-CA-CG-TGGA-CCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G.", false, SAME, NULp },
1385                { "XG*SNFWPVQAARNHRHD-XXXXXX-PRQNDSDRCYHHGAX-", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT---CG-GT-CA-CG-TG-GA----CCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G....", false, CHANGED,
1386                  "XG*SNFWPVQAARNHRHD-XXXXXX-PRQNDSDRCYHHGAX." }, // ok - only changes gaptype at EOS
1387                { "XG*SNXLXRXQA-ARNHRHD-RXXVX-PRQNDSDRCYHHGAX", "AT-GGCTAAAGAAACTT-TTGAC-CGGTC-CAAGCC---GCACGTAAACATCGGCACGAT---CGG-TCAC-GTG-GA---CCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G.", false, SAME, NULp },
1388                { "XG*SXXFXDXVQAXT*TSARXRSXVX-PRQNDSDRCYHHGAX", "AT-GGCTAAAGA-A-AC-TTT-T-GACCG-GTCCAAGCCGC-ACGTAAACATCGGCACGA-T-CGGTCA-C-GTG-GA---CCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G.", false, SAME, NULp },
1389                // -------------------------------------------- "123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123"
1390
1391                { NULp, NULp, false, SAME, NULp }
1392            };
1393
1394            int   arb_transl_table, codon_start;
1395            char *org_dna;
1396            {
1397                GB_transaction ta(gb_main);
1398                TEST_EXPECT_NO_ERROR(AWT_getTranslationInfo(gb_TaxOcell, arb_transl_table, codon_start));
1399                TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), "t=14,cs=0");
1400                org_dna = GB_read_string(gb_TaxOcell_dna);
1401            }
1402
1403            for (int s = 0; seq[s].seq; ++s) {
1404                TEST_ANNOTATE(GBS_global_string("s=%i", s));
1405                realign_check& S = seq[s];
1406
1407                {
1408                    GB_transaction ta(gb_main);
1409                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_amino, S.seq));
1410                }
1411                msgs  = "";
1412                error = ALI_realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, S.cutoff);
1413                TEST_EXPECT_NO_ERROR(error);
1414                TEST_EXPECT_EQUAL(msgs, "");
1415                {
1416                    GB_transaction ta(gb_main);
1417                    TEST_EXPECT_EQUAL(GB_read_char_pntr(gb_TaxOcell_dna), S.result);
1418
1419                    // test retranslation:
1420                    msgs  = "";
1421                    error = ALI_translate_marked(gb_main, true, false, 0, true, "ali_dna", "ali_pro");
1422                    TEST_EXPECT_NO_ERROR(error);
1423                    if (s == 10) {
1424                        TEST_EXPECT_EQUAL(msgs, "codon_start and transl_table entries were found for all translated taxa\n1 taxa converted\n  2.000000 stops per sequence found\n");
1425                    }
1426                    else if (s == 6) {
1427                        TEST_EXPECT_EQUAL(msgs, "codon_start and transl_table entries were found for all translated taxa\n1 taxa converted\n  0.000000 stops per sequence found\n");
1428                    }
1429                    else {
1430                        TEST_EXPECT_EQUAL(msgs, "codon_start and transl_table entries were found for all translated taxa\n1 taxa converted\n  1.000000 stops per sequence found\n");
1431                    }
1432
1433                    switch (S.retranslation) {
1434                        case SAME:
1435                            TEST_EXPECT_NULL(S.changed_prot);
1436                            TEST_EXPECT_EQUAL(GB_read_char_pntr(gb_TaxOcell_amino), S.seq);
1437                            break;
1438                        case CHANGED:
1439                            TEST_REJECT_NULL(S.changed_prot);
1440                            TEST_EXPECT_EQUAL(GB_read_char_pntr(gb_TaxOcell_amino), S.changed_prot);
1441                            break;
1442                    }
1443
1444                    TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), "t=14,cs=0");
1445                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_dna, org_dna)); // restore changed DB entry
1446                }
1447            }
1448            TEST_ANNOTATE(NULp);
1449
1450            free(org_dna);
1451        }
1452
1453        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1454
1455        // ----------------------------------------------------
1456        //      write some aa sequences provoking failures
1457        {
1458            struct realign_fail {
1459                const char *seq;
1460                const char *failure;
1461            };
1462
1463#define ERRPREFIX     "Automatic re-align failed for 'TaxOcell'\nReason: "
1464#define ERRPREFIX_LEN 49
1465
1466#define FAILONE "All marked species failed to realign\n"
1467
1468            // dna of TaxOcell:
1469            // "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G----......"
1470
1471            realign_fail seq[] = {
1472                //"XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-.." // original aa sequence
1473                // { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..", "sdfjlksdjf" }, // templ
1474
1475                // wanted realign failures:
1476                { "XG*SNFXXXXXAXNHRHD--XXX-PRQNDSDRCYHHGAX-..", "Sync behind 'X' failed foremost with: 'GGA' translates to 'G', not to 'P' at ali_pro:25 / ali_dna:70\n" FAILONE },    // ok to fail: 5 Xs impossible
1477                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..", "Alignment 'ali_dna' is too short (increase its length to 252)\n" FAILONE }, // ok to fail: wrong alignment length
1478                { "XG*SNFWPVQAARNHRHD--XXX-PRQNDSDRCYHHGAX-..", "Sync behind 'X' failed foremost with: 'GGA' translates to 'G', not to 'P' at ali_pro:25 / ali_dna:70\n" FAILONE },    // ok to fail
1479                { "XG*SNX-A-X-ARNHRHD--XXX-PRQNDSDRCYHHGAX-..", "Sync behind 'X' failed foremost with: 'TGA' never translates to 'A' at ali_pro:8 / ali_dna:19\n" FAILONE },           // ok to fail
1480                { "XG*SXFXPXQAXRNHRHD--RSRGPRQNDSDRCYHHGAX-..", "Sync behind 'X' failed foremost with: 'ACG' translates to 'T', not to 'R' at ali_pro:13 / ali_dna:36\n" FAILONE },    // ok to fail
1481                { "XG*SNFWPVQAARNHRHD-----GPRQNDSDRCYHHGAX-..", "Sync behind 'X' failed foremost with: 'CGG' translates to 'R', not to 'G' at ali_pro:24 / ali_dna:61\n" FAILONE },    // ok to fail: some AA missing in the middle
1482                { "XG*SNFWPVQAARNHRHDRSRGPRQNDSDRCYHHGAXHHGA.", "Sync behind 'X' failed foremost with: not enough nucs left for codon of 'H' at ali_pro:38 / ali_dna:117\n" FAILONE }, // ok to fail: too many AA
1483                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCY----H...", "Sync behind 'X' failed foremost with: too much trailing DNA (10 nucs, but only 9 columns left) at ali_pro:43 / ali_dna:106\n" FAILONE }, // ok to fail: not enough space to place extra nucs behind 'H'
1484                { "--SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX--..", "Not enough gaps to place 8 extra nucs at start of sequence at ali_pro:1 / ali_dna:1\n" FAILONE }, // also see related, succeeding test above (which has same AA seq; just one more leading gap)
1485
1486                // failing realignments that should work:
1487
1488                { NULp, NULp }
1489            };
1490
1491            {
1492                GB_transaction ta(gb_main);
1493                TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), "t=14,cs=0");
1494            }
1495
1496            for (int s = 0; seq[s].seq; ++s) {
1497                TEST_ANNOTATE(GBS_global_string("s=%i", s));
1498                {
1499                    GB_transaction ta(gb_main);
1500                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_amino, seq[s].seq));
1501                }
1502                msgs  = "";
1503                error = ALI_realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1504                TEST_EXPECT_NO_ERROR(error);
1505                TEST_EXPECT_CONTAINS(msgs, ERRPREFIX);
1506                TEST_EXPECT_EQUAL(msgs.c_str()+ERRPREFIX_LEN, seq[s].failure);
1507
1508                {
1509                    GB_transaction ta(gb_main);
1510                    TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), "t=14,cs=0"); // should not change if error
1511                }
1512            }
1513            TEST_ANNOTATE(NULp);
1514        }
1515
1516        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1517
1518        // ----------------------------------------------
1519        //      some examples for given DNA/AA pairs
1520
1521        {
1522            struct explicit_realign {
1523                const char *acids;
1524                const char *dna;
1525                int         table;
1526                const char *info;
1527                const char *msgs;
1528            };
1529
1530            // YTR (=X(2,9,16), =L(else))
1531            //     CTA (=T(2),        =L(else))
1532            //     CTG (=T(2), =S(9), =L(else))
1533            //     TTA (=*(16),       =L(else))
1534            //     TTG (=L(always))
1535            //
1536            // AAR (=X(6,11,14), =K(else))
1537            //     AAA (=N(6,11,14), =K(else))
1538            //     AAG (=K(always))
1539            //
1540            // ATH (=X(1,2,4,10,14), =I(else))
1541            //     ATA (=M(1,2,4,10,14), =I(else))
1542            //     ATC (=I(always))
1543            //     ATT (=I(always))
1544
1545            const char*const NO_TI = "t=-1,cs=-1";
1546
1547            explicit_realign example[] = {
1548                { "LK", "TTGAAG", -1, NO_TI,        NULp }, // fine (for any table)
1549
1550                { "G",  "RGG",    -1, "t=10,cs=0",  NULp }, // correctly detects TI(10)
1551
1552                { "LK",  "YTRAAR",    2,  "t=2,cs=0",  "Not all IUPAC-combinations of 'YTR' translate to 'L' (for trans-table 2) at ali_pro:1 / ali_dna:1\n" }, // expected failure (CTA->T for table=2)
1553                { "LX",  "YTRAAR",    -1, NO_TI,       NULp }, // fine (AAR->X for table=6,11,14)
1554                { "LXX", "YTRAARATH", -1, "t=14,cs=0", NULp }, // correctly detects TI(14)
1555                { "LXI", "YTRAARATH", -1, NO_TI,       NULp }, // fine (for table=6,11)
1556
1557                { "LX", "YTRAAR", 2,  "t=2,cs=0",   "Not all IUPAC-combinations of 'YTR' translate to 'L' (for trans-table 2) at ali_pro:1 / ali_dna:1\n" }, // expected failure (AAR->K for table=2)
1558                { "LK", "YTRAAR", -1, NO_TI,        NULp }, // fine           (AAR->K for table!=6,11,14)
1559                { "LK", "YTRAAR", 6,  "t=6,cs=0",   "Not all IUPAC-combinations of 'AAR' translate to 'K' (for trans-table 6) at ali_pro:2 / ali_dna:4\n" }, // expected failure (AAA->N for table=6)
1560                { "XK", "YTRAAR", -1, NO_TI,        NULp }, // fine           (YTR->X for table=2,9,16)
1561
1562                { "XX",   "-YTRAAR",      0,  "t=0,cs=0", NULp },                                                                                             // does not fail because it realigns such that it translates back to 'XXX'
1563                { "XXL",  "YTRAARTTG",    0,  "t=0,cs=0", "Not enough gaps to place 2 extra nucs at start of sequence at ali_pro:1 / ali_dna:1\n" },          // expected failure (none can translate to X with table= 0, so it tries )
1564                { "-XXL", "-YTRA-AR-TTG", 0,  "t=0,cs=0", NULp },                                                                                             // does not fail because it realigns such that it translates back to 'XXXL'
1565                { "IXXL", "ATTYTRAARTTG", 0,  "t=0,cs=0", "Sync behind 'X' failed foremost with: 'RTT' never translates to 'L' (for trans-table 0) at ali_pro:4 / ali_dna:9\n" }, // expected failure (none of the 2 middle codons can translate to X with table= 0)
1566                { "XX",   "-YTRAAR",      -1, NO_TI,      NULp },                                                                                             // does not fail because it realigns such that it translates back to 'XXX'
1567                { "IXXL", "ATTYTRAARTTG", -1, NO_TI,      "Sync behind 'X' failed foremost with: 'RTT' never translates to 'L' at ali_pro:4 / ali_dna:9\n" }, // expected failure (not both 2 middle codons can translate to X with same table)
1568
1569                { "LX", "YTRATH", -1, NO_TI,        NULp }, // fine                (ATH->X for table=1,2,4,10,14)
1570                { "LX", "YTRATH", 2,  "t=2,cs=0",   "Not all IUPAC-combinations of 'YTR' translate to 'L' (for trans-table 2) at ali_pro:1 / ali_dna:1\n" }, // expected failure (YTR->X for table=2)
1571                { "XX", "YTRATH", 2,  "t=2,cs=0",   NULp }, // fine                (both->X for table=2)
1572                { "XX", "YTRATH", -1, "t=2,cs=0",   NULp }, // correctly detects TI(2)
1573
1574                { "XX", "AARATH", 14, "t=14,cs=0",  NULp }, // fine (both->X for table=14)
1575                { "XX", "AARATH", -1, "t=14,cs=0",  NULp }, // correctly detects TI(14)
1576                { "KI", "AARATH", -1, NO_TI,        NULp }, // fine (for table!=1,2,4,6,10,11,14)
1577                { "KI", "AARATH", 4,  "t=4,cs=0",   "Not all IUPAC-combinations of 'ATH' translate to 'I' (for trans-table 4) at ali_pro:2 / ali_dna:4\n" }, // expected failure (ATH->X for table=4)
1578                { "KX", "AARATH", 14, "t=14,cs=0",  "Not all IUPAC-combinations of 'AAR' translate to 'K' (for trans-table 14) at ali_pro:1 / ali_dna:1\n" }, // expected failure (AAR->X for table=14)
1579                { "KX", "AARATH", -1, NO_TI,        NULp }, // fine for table=1,2,4,10
1580                { "KX", "AARATH", 4,  "t=4,cs=0",   NULp }, // test table=4
1581                { "XI", "AARATH", 14, "t=14,cs=0",  "Sync behind 'X' failed foremost with: Not all IUPAC-combinations of 'ATH' translate to 'I' (for trans-table 14) at ali_pro:2 / ali_dna:4\n" }, // expected failure (ATH->X for table=14)
1582                { "KI", "AARATH", 14, "t=14,cs=0",  "Not all IUPAC-combinations of 'AAR' translate to 'K' (for trans-table 14) at ali_pro:1 / ali_dna:1\n" }, // expected failure (AAR->X for table=14)
1583
1584                { NULp, NULp, 0, NULp, NULp }
1585            };
1586
1587            for (int e = 0; example[e].acids; ++e) {
1588                const explicit_realign& E = example[e];
1589                TEST_ANNOTATE(GBS_global_string("%s <- %s (#%i)", E.acids, E.dna, E.table));
1590
1591                {
1592                    GB_transaction ta(gb_main);
1593                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_dna, E.dna));
1594                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_amino, E.acids));
1595                    if (E.table == -1) {
1596                        TEST_EXPECT_NO_ERROR(AWT_removeTranslationInfo(gb_TaxOcell));
1597                    }
1598                    else {
1599                        TEST_EXPECT_NO_ERROR(AWT_saveTranslationInfo(gb_TaxOcell, E.table, 0));
1600                    }
1601                }
1602
1603                msgs  = "";
1604                error = ALI_realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1605                TEST_EXPECT_NULL(error);
1606                if (E.msgs) {
1607                    TEST_EXPECT_CONTAINS(msgs, ERRPREFIX);
1608                    string wanted_msgs = string(E.msgs)+FAILONE;
1609                    TEST_EXPECT_EQUAL(msgs.c_str()+ERRPREFIX_LEN, wanted_msgs);
1610                }
1611                else {
1612                    TEST_EXPECT_EQUAL(msgs, "");
1613                }
1614
1615                GB_transaction ta(gb_main);
1616                if (!error) {
1617                    const char *dnaseq      = GB_read_char_pntr(gb_TaxOcell_dna);
1618                    size_t      expextedLen = strlen(E.dna);
1619                    size_t      seqlen      = strlen(dnaseq);
1620                    char       *firstPart   = ARB_strndup(dnaseq, expextedLen);
1621                    size_t      dna_behind;
1622                    char       *nothing     = unalign(dnaseq+expextedLen, seqlen-expextedLen, dna_behind);
1623
1624                    TEST_EXPECT_EQUAL(firstPart, E.dna);
1625                    TEST_EXPECT_EQUAL(dna_behind, 0);
1626                    TEST_EXPECT_EQUAL(nothing, "");
1627
1628                    free(nothing);
1629                    free(firstPart);
1630                }
1631                TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), E.info);
1632            }
1633        }
1634
1635        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1636
1637        // ----------------------------------
1638        //      invalid translation info
1639        {
1640            GB_transaction ta(gb_main);
1641
1642            TEST_EXPECT_NO_ERROR(AWT_saveTranslationInfo(gb_TaxOcell, 14, 0));
1643            GBDATA *gb_trans_table = GB_entry(gb_TaxOcell, "transl_table");
1644            TEST_EXPECT_NO_ERROR(GB_write_string(gb_trans_table, "666")); // evil translation table
1645        }
1646
1647        msgs  = "";
1648        error = ALI_realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1649        TEST_EXPECT_NO_ERROR(error);
1650        TEST_EXPECT_EQUAL(msgs, ERRPREFIX "Error while reading 'transl_table' (Illegal (or unsupported) value (666) in 'transl_table' (item='TaxOcell'))\n" FAILONE);
1651        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1652
1653        // ---------------------------------------
1654        //      source/dest alignment missing
1655        for (int i = 0; i<2; ++i) {
1656            TEST_ANNOTATE(GBS_global_string("i=%i", i));
1657
1658            {
1659                GB_transaction ta(gb_main);
1660
1661                GBDATA *gb_ali = GB_get_father(GBT_find_sequence(gb_TaxOcell, i ? "ali_pro" : "ali_dna"));
1662                GB_push_my_security(gb_main);
1663                TEST_EXPECT_NO_ERROR(GB_delete(gb_ali));
1664                GB_pop_my_security(gb_main);
1665            }
1666
1667            msgs  = "";
1668            error = ALI_realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1669            TEST_EXPECT_NO_ERROR(error);
1670            if (i) {
1671                TEST_EXPECT_EQUAL(msgs, ERRPREFIX "No data in alignment 'ali_pro'\n" FAILONE);
1672            }
1673            else {
1674                TEST_EXPECT_EQUAL(msgs, ERRPREFIX "No data in alignment 'ali_dna'\n" FAILONE);
1675            }
1676        }
1677        TEST_ANNOTATE(NULp);
1678
1679        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1680    }
1681
1682#undef ERRPREFIX
1683#undef ERRPREFIX_LEN
1684
1685    GB_close(gb_main);
1686    ARB_install_handlers(*old_handlers);
1687}
1688
1689static const char *permOf(const Distributor& dist) {
1690    const int   MAXDIST = 10;
1691    static char buffer[MAXDIST+1];
1692
1693    ali_assert(dist.size() <= MAXDIST);
1694    for (int p = 0; p<dist.size(); ++p) {
1695        buffer[p] = '0'+dist[p];
1696    }
1697    buffer[dist.size()] = 0;
1698
1699    return buffer;
1700}
1701
1702static arb_test::match_expectation stateOf(Distributor& dist, const char *expected_perm, bool hasNext) {
1703    using namespace arb_test;
1704
1705    expectation_group expected;
1706    expected.add(that(permOf(dist)).is_equal_to(expected_perm));
1707    expected.add(that(dist.next()).is_equal_to(hasNext));
1708    return all().ofgroup(expected);
1709}
1710
1711void TEST_distributor() {
1712    TEST_EXPECT_EQUAL(Distributor(3, 2).get_error(), "not enough nucleotides");
1713    TEST_EXPECT_EQUAL(Distributor(3, 10).get_error(), "too much nucleotides");
1714
1715    Distributor minDist(3, 3);
1716    TEST_EXPECTATION(stateOf(minDist, "111", false));
1717
1718    Distributor maxDist(3, 9);
1719    TEST_EXPECTATION(stateOf(maxDist, "333", false));
1720
1721    Distributor meanDist(3, 6);
1722    TEST_EXPECTATION(stateOf(meanDist, "123", true));
1723    TEST_EXPECTATION(stateOf(meanDist, "132", true));
1724    TEST_EXPECTATION(stateOf(meanDist, "213", true));
1725    TEST_EXPECTATION(stateOf(meanDist, "222", true));
1726    TEST_EXPECTATION(stateOf(meanDist, "231", true));
1727    TEST_EXPECTATION(stateOf(meanDist, "312", true));
1728    TEST_EXPECTATION(stateOf(meanDist, "321", false));
1729
1730    Distributor belowMax(4, 11);
1731    TEST_EXPECTATION(stateOf(belowMax, "2333", true));
1732    TEST_EXPECTATION(stateOf(belowMax, "3233", true));
1733    TEST_EXPECTATION(stateOf(belowMax, "3323", true));
1734    TEST_EXPECTATION(stateOf(belowMax, "3332", false));
1735
1736    Distributor aboveMin(4, 6);
1737    TEST_EXPECTATION(stateOf(aboveMin, "1113", true));
1738    TEST_EXPECTATION(stateOf(aboveMin, "1122", true));
1739    TEST_EXPECTATION(stateOf(aboveMin, "1131", true));
1740    TEST_EXPECTATION(stateOf(aboveMin, "1212", true));
1741    TEST_EXPECTATION(stateOf(aboveMin, "1221", true));
1742    TEST_EXPECTATION(stateOf(aboveMin, "1311", true));
1743    TEST_EXPECTATION(stateOf(aboveMin, "2112", true));
1744    TEST_EXPECTATION(stateOf(aboveMin, "2121", true));
1745    TEST_EXPECTATION(stateOf(aboveMin, "2211", true));
1746    TEST_EXPECTATION(stateOf(aboveMin, "3111", false));
1747
1748    Distributor check(6, 8);
1749    TEST_EXPECTATION(stateOf(check, "111113", true));
1750    TEST_EXPECTATION(stateOf(check, "111122", true));
1751    TEST_EXPECTATION(stateOf(check, "111131", true));
1752    TEST_EXPECTATION(stateOf(check, "111212", true));
1753    TEST_EXPECTATION(stateOf(check, "111221", true));
1754    TEST_EXPECTATION(stateOf(check, "111311", true));
1755    TEST_EXPECTATION(stateOf(check, "112112", true));
1756    TEST_EXPECTATION(stateOf(check, "112121", true));
1757    TEST_EXPECTATION(stateOf(check, "112211", true));
1758    TEST_EXPECTATION(stateOf(check, "113111", true));
1759    TEST_EXPECTATION(stateOf(check, "121112", true));
1760    TEST_EXPECTATION(stateOf(check, "121121", true));
1761    TEST_EXPECTATION(stateOf(check, "121211", true));
1762    TEST_EXPECTATION(stateOf(check, "122111", true));
1763    TEST_EXPECTATION(stateOf(check, "131111", true));
1764    TEST_EXPECTATION(stateOf(check, "211112", true));
1765    TEST_EXPECTATION(stateOf(check, "211121", true));
1766    TEST_EXPECTATION(stateOf(check, "211211", true));
1767    TEST_EXPECTATION(stateOf(check, "212111", true));
1768    TEST_EXPECTATION(stateOf(check, "221111", true));
1769    TEST_EXPECTATION(stateOf(check, "311111", false));
1770}
1771
1772#endif // UNIT_TESTS
1773
1774// --------------------------------------------------------------------------------
1775
Note: See TracBrowser for help on using the repository browser.