source: branches/help/SL/ALILINK/TranslateRealign.cxx

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