source: tags/ms_r16q2/NTREE/ad_transpro.cxx

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