source: tags/ms_r16q3/NTREE/ad_transpro.cxx

Last change on this file was 15176, checked in by westram, 8 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 84.8 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, ARB_strdup(data), ARB_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    RefPtr<const char> at_prot; // points into aligned protein seq
568    RefPtr<const char> at_dna;  // points into 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
589    GB_ERROR why() const { return reason.empty() ? NULL : reason.c_str(); }
590    const char *protein_at() const { return at_prot; }
591    const char *dna_at() const { return at_dna; }
592
593    operator bool() const { return !reason.empty(); }
594
595    void add_prefix(const char *prefix) {
596        nt_assert(!reason.empty());
597        reason = string(prefix)+reason;
598    }
599
600    bool operator>(const FailedAt& other) const { return cmp(other)>0; }
601};
602
603class RealignAttempt : virtual Noncopyable {
604    TransTables           allowed;
605    SizedReadBuffer       compressed_dna;
606    BufferPtr<const char> aligned_protein;
607    SizedWriteBuffer      target_dna;
608    FailedAt              fail;
609    bool                  cutoff_dna;
610
611    void perform();
612
613    bool sync_behind_X_and_distribute(const int x_count, char *const x_start, const char *const x_start_prot);
614
615public:
616    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_)
617        : allowed(allowed_),
618          compressed_dna(compressed_dna_, compressed_len_),
619          aligned_protein(aligned_protein_),
620          target_dna(target_dna_, target_len_),
621          cutoff_dna(cutoff_dna_)
622    {
623        nt_assert(aligned_protein[0]);
624        perform();
625    }
626
627    const TransTables& get_remaining_tables() const { return allowed; }
628    const FailedAt& failed() const { return fail; }
629};
630
631static GB_ERROR distribute_xdata(SizedReadBuffer& dna, size_t xcount, char *xtarget_, bool gap_before, bool gap_after, const TransTables& allowed, TransTables& remaining) {
632    /*! distributes 'dna' to marked X-positions
633     * @param xtarget destination buffer (target positions are marked with '!')
634     * @param xcount number of X's encountered
635     * @param gap_before true if resulting realignment has a gap or the start of alignment before the X-positions
636     * @param gap_after analog to 'gap_before'
637     * @param allowed allowed translation tables
638     * @param remaining remaining allowed translation tables (with those tables disabled for which no distribution possible)
639     * @return error if dna distribution wasn't possible
640     */
641
642    BufferPtr<char> xtarget(xtarget_);
643    Distributor     dist(xcount, dna.length());
644    GB_ERROR        error = dist.get_error();
645    if (!error) {
646        Distributor best(dist);
647        TransTables best_remaining = allowed;
648
649        while (dist.next()) {
650            if (dist.get_score() > best.get_score()) {
651                if (!dist.mayFailTranslation() || best.mayFailTranslation()) {
652                    best           = dist;
653                    best_remaining = allowed;
654                    nt_assert(best_remaining.is_subset_of(allowed));
655                }
656            }
657        }
658
659        if (best.mayFailTranslation()) {
660            TransTables curr_remaining;
661            if (best.translates_to_Xs(dna, allowed, curr_remaining)) {
662                best_remaining = curr_remaining;
663                nt_assert(best_remaining.is_subset_of(allowed));
664            }
665            else {
666                nt_assert(!error);
667                error = "no translating X-distribution found";
668                dist.reset();
669                do {
670                    if (dist.translates_to_Xs(dna, allowed, curr_remaining)) {
671                        best           = dist;
672                        best_remaining = curr_remaining;
673                        error          = NULL;
674                        nt_assert(best_remaining.is_subset_of(allowed));
675                        break;
676                    }
677                } while (dist.next());
678
679                while (dist.next()) {
680                    if (dist.get_score() > best.get_score()) {
681                        if (dist.translates_to_Xs(dna, allowed, curr_remaining)) {
682                            best           = dist;
683                            best_remaining = curr_remaining;
684                            nt_assert(best_remaining.is_subset_of(allowed));
685                        }
686                    }
687                }
688            }
689        }
690
691        if (!error) { // now really distribute nucs
692            for (int x = 0; x<best.size(); ++x) {
693                while (xtarget[0] != '!') {
694                    nt_assert(xtarget[1] && xtarget[2]); // buffer overflow
695                    xtarget.inc(3);
696                }
697
698                switch (best[x]) {
699                    case 2: {
700                        enum { UNDECIDED, SPREAD, LEFT, RIGHT } mode = UNDECIDED;
701
702                        bool is_1st_X  = xtarget.offset() == 0;
703                        bool gaps_left = is_1st_X ? gap_before : isGap(xtarget[-1]);
704
705                        if (gaps_left) mode = LEFT;
706                        else { // definitely has no gap left!
707                            bool is_last_X  = x == best.size()-1;
708                            int  next_nucs  = is_last_X ? 0 : best[x+1];
709                            bool gaps_right = isGap(xtarget[3]) || next_nucs == 1 || (is_last_X && gap_after);
710
711                            if (gaps_right) mode = RIGHT;
712                            else {
713                                bool nogaps_right = next_nucs == 3 || (is_last_X && !gap_after);
714                                if (nogaps_right) { // we know, we have NO adjacent gaps
715                                    mode = is_last_X ? LEFT : (is_1st_X ? RIGHT : SPREAD);
716                                }
717                                else {
718                                    nt_assert(!is_last_X);
719                                    mode = RIGHT; // forward problem to next X
720                                }
721                            }
722                        }
723
724                        char d1 = dna.get();
725                        char d2 = dna.get();
726
727                        switch (mode) {
728                            case UNDECIDED: nt_assert(0); // fall-through in NDEBUG
729                            case SPREAD: xtarget.put(d1,  '-', d2);  break;
730                            case LEFT:   xtarget.put(d1,  d2,  '-'); break;
731                            case RIGHT:  xtarget.put('-', d1,  d2);  break;
732                        }
733
734                        break;
735                    }
736                    case 1: xtarget.put('-', dna.get(), '-'); break;
737                    case 3: xtarget.copy(dna, 3); break;
738                    default: nt_assert(0); break;
739                }
740                nt_assert(dna.valid());
741            }
742
743            nt_assert(!error);
744            remaining = best_remaining;
745            nt_assert(remaining.is_subset_of(allowed));
746        }
747    }
748
749    return error;
750}
751
752bool RealignAttempt::sync_behind_X_and_distribute(const int x_count, char *const x_start, const char *const x_start_prot) {
753    /*! brute-force search for sync behind 'X' and distribute dna onto X positions
754     * @param x_count number of X encountered
755     * @param x_start dna read position
756     * @param x_start_prot protein read position
757     * @return true if sync and distribution succeed
758     */
759
760    bool complete = false;
761
762    nt_assert(!failed());
763    nt_assert(aligned_protein.offset()>0);
764    const char p = aligned_protein[-1];
765
766    size_t compressed_rest_len = compressed_dna.restLength();
767    nt_assert(strlen(compressed_dna) == compressed_rest_len);
768
769    size_t min_dna = x_count;
770    size_t max_dna = std::min(size_t(x_count)*3, compressed_rest_len);
771
772    if (min_dna>max_dna) {
773        fail = FailedAt("not enough nucs for X's at sequence end", x_start_prot, compressed_dna);
774    }
775    else if (p) {
776        FailedAt foremost;
777        size_t   target_rest_len = target_dna.restLength();
778
779        for (size_t x_dna = min_dna; x_dna<=max_dna; ++x_dna) { // prefer low amounts of used dna
780            const char *dna_rest     = compressed_dna + x_dna;
781            size_t      dna_rest_len = compressed_rest_len     - x_dna;
782
783            nt_assert(strlen(dna_rest) == dna_rest_len);
784            nt_assert(compressed_rest_len>=x_dna);
785
786            RealignAttempt attemptRest(allowed, dna_rest, dna_rest_len, aligned_protein-1, target_dna, target_rest_len, cutoff_dna);
787            FailedAt       restFailed = attemptRest.failed();
788
789            if (!restFailed) {
790                SizedReadBuffer distrib_dna(compressed_dna, x_dna);
791
792                bool has_gap_before = x_start == target_dna.start() ? true : isGap(x_start[-1]);
793                bool has_gap_after  = isGap(dna_rest[0]);
794
795                TransTables remaining;
796                GB_ERROR    disterr = distribute_xdata(distrib_dna, x_count, x_start, has_gap_before, has_gap_after, attemptRest.get_remaining_tables(), remaining);
797                if (disterr) {
798                    restFailed = FailedAt(disterr, x_start_prot, dna_rest); // prot=start of Xs; dna=start of sync (behind Xs)
799                }
800                else {
801                    nt_assert(remaining.is_subset_of(allowed));
802                    nt_assert(remaining.is_subset_of(attemptRest.get_remaining_tables()));
803                    allowed = remaining;
804                }
805            }
806
807            if (restFailed) {
808                if (restFailed > foremost) foremost = restFailed; // track "best" failure (highest fail position)
809            }
810            else { // success
811                foremost = FailedAt();
812                complete = true;
813                break; // use first success and return
814            }
815        }
816
817        if (foremost) {
818            nt_assert(!complete);
819            fail = foremost;
820            if (!strstr(fail.why(), "Sync behind 'X'")) { // do not spam repetitive sync-failures
821                fail.add_prefix("Sync behind 'X' failed foremost with: ");
822            }
823        }
824        else {
825            nt_assert(complete);
826        }
827    }
828    else {
829        GB_ERROR fail_reason = "internal error: no distribution attempted";
830        nt_assert(min_dna>0);
831        size_t x_dna;
832        for (x_dna = max_dna; x_dna>=min_dna; --x_dna) {     // prefer high amounts of dna
833            SizedReadBuffer append_dna(compressed_dna, x_dna);
834            TransTables     remaining;
835            fail_reason = distribute_xdata(append_dna, x_count, x_start, false, true, allowed, remaining);
836            if (!fail_reason) { // found distribution -> done
837                nt_assert(remaining.is_subset_of(allowed));
838                allowed = remaining;
839                break;
840            }
841        }
842
843        if (fail_reason) {
844            fail = FailedAt(fail_reason, x_start_prot+1, compressed_dna); // report error at start of X's
845        }
846        else {
847            fail = FailedAt(); // clear
848            compressed_dna.inc(x_dna);
849        }
850    }
851
852    nt_assert(implicated(complete, allowed.any()));
853
854    return complete;
855}
856
857void RealignAttempt::perform() {
858    bool complete = false; // set to true, if recursive attempt succeeds
859
860    while (char p = toupper(aligned_protein.get())) {
861        if (p=='X') { // one X represents 1 to 3 DNAs (normally 1 or 2, but 'NNN' translates to 'X')
862            char       *x_start      = target_dna;
863            const char *x_start_prot = aligned_protein-1;
864            int         x_count      = 0;
865
866            for (;;) {
867                if      (p=='X')   { x_count++; target_dna.put('!', 3); } // fill X space with marker
868                else if (isGap(p)) target_dna.put(p, 3);
869                else break;
870
871                p = toupper(aligned_protein.get());
872            }
873
874            nt_assert(x_count);
875            nt_assert(!complete);
876            complete = sync_behind_X_and_distribute(x_count, x_start, x_start_prot);
877            if (!complete && !failed()) {
878                if (p) { // not all proteins were processed
879                    fail = FailedAt("internal error", aligned_protein-1, compressed_dna);
880                    nt_assert(0);
881                }
882            }
883            break; // done
884        }
885
886        if (isGap(p)) target_dna.put(p, 3);
887        else {
888            TransTables remaining;
889            size_t      compressed_rest_len = compressed_dna.restLength();
890
891            if (compressed_rest_len<3) {
892                fail = FailedAt(GBS_global_string("not enough nucs left for codon of '%c'", p), aligned_protein-1, compressed_dna);
893            }
894            else {
895                nt_assert(strlen(compressed_dna) == compressed_rest_len);
896                nt_assert(compressed_rest_len >= 3);
897                const char *why_fail;
898                if (!AWT_is_codon(p, compressed_dna, allowed, remaining, &why_fail)) {
899                    fail = FailedAt(why_fail, aligned_protein-1, compressed_dna);
900                }
901            }
902
903            if (failed()) break;
904
905            nt_assert(remaining.is_subset_of(allowed));
906            allowed = remaining;
907            target_dna.copy(compressed_dna, 3);
908        }
909    }
910
911    nt_assert(compressed_dna.valid());
912
913    if (!failed() && !complete) {
914        while (target_dna.offset()>0 && isGap(target_dna[-1])) --target_dna; // remove terminal gaps
915
916        if (!cutoff_dna) { // append leftover dna-data (data w/o corresponding aa)
917            size_t compressed_rest_len = compressed_dna.restLength();
918            size_t target_rest_len = target_dna.restLength();
919            if (compressed_rest_len<=target_rest_len) {
920                target_dna.copy(compressed_dna, compressed_rest_len);
921            }
922            else {
923                fail = FailedAt(GBS_global_string("too much trailing DNA (%zu nucs, but only %zu columns left)",
924                                                  compressed_rest_len, target_rest_len),
925                                aligned_protein-1, compressed_dna);
926            }
927        }
928
929        if (!failed()) target_dna.put('.', target_dna.restLength()); // fill rest of sequence with dots
930        *target_dna = 0;
931    }
932
933#if defined(ASSERTION_USED)
934    if (!failed()) {
935        nt_assert(strlen(target_dna.start()) == target_dna.length());
936    }
937#endif
938}
939
940inline char *unalign(const char *data, size_t len, size_t& compressed_len) {
941    // removes gaps from sequence
942    char *compressed = ARB_alloc<char>(len+1);
943    compressed_len = 0;
944    for (size_t p = 0; p<len && data[p]; ++p) {
945        if (!isGap(data[p])) {
946            compressed[compressed_len++] = data[p];
947        }
948    }
949    compressed[compressed_len] = 0;
950    return compressed;
951}
952
953class Realigner {
954    const char *ali_source;
955    const char *ali_dest;
956
957    size_t ali_len;        // of ali_dest
958    size_t needed_ali_len; // >ali_len if ali_dest is too short; 0 otherwise
959
960    const char *fail_reason;
961
962    GB_ERROR annotate_fail_position(const FailedAt& failed, const char *source, const char *dest, const char *compressed_dest) {
963        int source_fail_pos = failed.protein_at() - source;
964        int dest_fail_pos   = 0;
965        {
966            int fail_d_base_count = failed.dna_at() - compressed_dest;
967
968            const char *dp = dest;
969
970            for (;;) {
971                char c = *dp++;
972
973                if (!c) { // failure at end of sequence
974                    dest_fail_pos++; // report position behind last non-gap
975                    break;
976                }
977                if (!isGap(c)) {
978                    dest_fail_pos = (dp-1)-dest;
979                    if (!fail_d_base_count) break;
980                    fail_d_base_count--;
981                }
982            }
983        }
984        return GBS_global_string("%s at %s:%i / %s:%i",
985                                 failed.why(),
986                                 ali_source, info2bio(source_fail_pos),
987                                 ali_dest, info2bio(dest_fail_pos));
988    }
989
990
991    static void calc_needed_dna(const char *prot, size_t len, size_t& minDNA, size_t& maxDNA) {
992        minDNA = maxDNA = 0;
993        for (size_t o = 0; o<len; ++o) {
994            char p = toupper(prot[o]);
995            if (p == 'X') {
996                minDNA += 1;
997                maxDNA += 3;
998            }
999            else if (!isGap(p)) {
1000                minDNA += 3;
1001                maxDNA += 3;
1002            }
1003        }
1004    }
1005    static size_t countLeadingGaps(const char *buffer) {
1006        size_t gaps = 0;
1007        for (int o = 0; isGap(buffer[o]); ++o) ++gaps;
1008        return gaps;
1009    }
1010
1011public:
1012    Realigner(const char *ali_source_, const char *ali_dest_, size_t ali_len_)
1013        : ali_source(ali_source_),
1014          ali_dest(ali_dest_),
1015          ali_len(ali_len_),
1016          needed_ali_len(0)
1017    {
1018        clear_failure();
1019    }
1020
1021    size_t get_needed_dest_alilen() const { return needed_ali_len; }
1022
1023    void set_failure(const char *reason) { fail_reason = reason; }
1024    void clear_failure() { fail_reason = NULL; }
1025
1026    const char *failure() const { return fail_reason; }
1027
1028    char *realign_seq(TransTables& allowed, const char *const source, size_t source_len, const char *const dest, size_t dest_len, bool cutoff_dna) {
1029        nt_assert(!failure());
1030
1031        size_t  wanted_ali_len = source_len*3;
1032        char   *buffer         = NULL;
1033
1034        if (ali_len<wanted_ali_len) {
1035            fail_reason = GBS_global_string("Alignment '%s' is too short (increase its length to %zu)", ali_dest, wanted_ali_len);
1036            if (wanted_ali_len>needed_ali_len) needed_ali_len = wanted_ali_len;
1037        }
1038        else {
1039            // compress destination DNA (=remove align-characters):
1040            size_t  compressed_len;
1041            char   *compressed_dest = unalign(dest, dest_len, compressed_len);
1042
1043            ARB_alloc(buffer, ali_len+1);
1044
1045            RealignAttempt attempt(allowed, compressed_dest, compressed_len, source, buffer, ali_len, cutoff_dna);
1046            FailedAt       failed = attempt.failed();
1047
1048            if (failed) {
1049                // test for superfluous DNA at sequence start
1050                size_t min_dna, max_dna;
1051                calc_needed_dna(source, source_len, min_dna, max_dna);
1052
1053                if (min_dna<compressed_len) { // we have more DNA than we need
1054                    size_t extra_dna = compressed_len-min_dna;
1055                    for (size_t skip = 1; skip<=extra_dna; ++skip) {
1056                        RealignAttempt attemptSkipped(allowed, compressed_dest+skip, compressed_len-skip, source, buffer, ali_len, cutoff_dna);
1057                        if (!attemptSkipped.failed()) {
1058                            failed = FailedAt(); // clear
1059                            if (!cutoff_dna) {
1060                                size_t start_gaps = countLeadingGaps(buffer);
1061                                if (start_gaps<skip) {
1062                                    failed = FailedAt(GBS_global_string("Not enough gaps to place %zu extra nucs at start of sequence",
1063                                                                        skip), source, compressed_dest);
1064                                }
1065                                else { // success
1066                                    memcpy(buffer+(start_gaps-skip), compressed_dest, skip); // copy-in skipped dna
1067                                }
1068                            }
1069                            if (!failed) {
1070                                nt_assert(attempt.get_remaining_tables().is_subset_of(allowed));
1071                                allowed = attemptSkipped.get_remaining_tables();
1072                            }
1073                            break; // no need to skip more dna, when we already have too few leading gaps
1074                        }
1075                    }
1076                }
1077            }
1078            else {
1079                nt_assert(attempt.get_remaining_tables().is_subset_of(allowed));
1080                allowed = attempt.get_remaining_tables();
1081            }
1082
1083            if (failed) {
1084                fail_reason = annotate_fail_position(failed, source, dest, compressed_dest);
1085                freenull(buffer);
1086            }
1087            free(compressed_dest);
1088        }
1089        nt_assert(contradicted(buffer, fail_reason));
1090        return buffer;
1091    }
1092};
1093
1094struct Data : virtual Noncopyable {
1095    GBDATA *gb_data;
1096    char   *data;
1097    size_t  len;
1098    char   *error;
1099
1100    Data(GBDATA *gb_species, const char *aliName)
1101        : gb_data(NULL),
1102          data(NULL),
1103          len(0),
1104          error(NULL)
1105    {
1106        GBDATA *gb_ali = GB_entry(gb_species, aliName);
1107        if (gb_ali) {
1108            gb_data = GB_entry(gb_ali, "data");
1109            if (gb_data) {
1110                data          = GB_read_string(gb_data);
1111                if (data) len = GB_read_string_count(gb_data);
1112                else error    = ARB_strdup(GB_await_error());
1113                return;
1114            }
1115        }
1116        error = GBS_global_string_copy("No data in alignment '%s'", aliName);
1117    }
1118    ~Data() {
1119        free(data);
1120        free(error);
1121    }
1122};
1123
1124
1125static GB_ERROR realign_marked(GBDATA *gb_main, const char *ali_source, const char *ali_dest, size_t& neededLength, bool unmark_succeeded, bool cutoff_dna) {
1126    /*! realigns DNA alignment of marked sequences according to their protein alignment
1127     * @param ali_source protein source alignment
1128     * @param ali_dest modified DNA alignment
1129     * @param neededLength result: minimum alignment length needed in ali_dest (if too short) or 0 if long enough
1130     * @param unmark_succeeded unmark all species that were successfully realigned
1131     */
1132    AP_initialize_codon_tables();
1133
1134    nt_assert(GB_get_transaction_level(gb_main) == 0);
1135    GB_transaction ta(gb_main); // do not abort (otherwise sth goes wrong with species marks)
1136
1137    {
1138        GBDATA *gb_source = GBT_get_alignment(gb_main, ali_source); if (!gb_source) return "Please select a valid source alignment";
1139        GBDATA *gb_dest   = GBT_get_alignment(gb_main, ali_dest);   if (!gb_dest)   return "Please select a valid destination alignment";
1140    }
1141
1142    if (GBT_get_alignment_type(gb_main, ali_source) != GB_AT_AA)  return "Invalid source alignment type";
1143    if (GBT_get_alignment_type(gb_main, ali_dest)   != GB_AT_DNA) return "Invalid destination alignment type";
1144
1145    long ali_len = GBT_get_alignment_len(gb_main, ali_dest);
1146    nt_assert(ali_len>0);
1147
1148    GB_ERROR error = 0;
1149
1150    long no_of_marked_species    = GBT_count_marked_species(gb_main);
1151    long no_of_realigned_species = 0; // count successfully realigned species
1152
1153    arb_progress progress("Re-aligner", no_of_marked_species);
1154    progress.auto_subtitles("Re-aligning species");
1155
1156    Realigner realigner(ali_source, ali_dest, ali_len);
1157
1158    for (GBDATA *gb_species = GBT_first_marked_species(gb_main);
1159         !error && gb_species;
1160         gb_species = GBT_next_marked_species(gb_species))
1161    {
1162        realigner.clear_failure();
1163
1164        Data source(gb_species, ali_source);
1165        Data dest(gb_species, ali_dest);
1166
1167        if      (source.error) realigner.set_failure(source.error);
1168        else if (dest.error)   realigner.set_failure(dest.error);
1169
1170        if (!realigner.failure()) {
1171            TransTables allowed; // default: all translation tables allowed
1172#if defined(ASSERTION_USED)
1173            bool has_valid_translation_info = false;
1174#endif
1175            {
1176                int arb_transl_table, codon_start;
1177                GB_ERROR local_error = AWT_getTranslationInfo(gb_species, arb_transl_table, codon_start);
1178                if (local_error) {
1179                    realigner.set_failure(GBS_global_string("Error while reading 'transl_table' (%s)", local_error));
1180                }
1181                else if (arb_transl_table >= 0) {
1182                    // we found a 'transl_table' entry -> restrict used code to the code stored there
1183                    allowed.forbidAllBut(arb_transl_table);
1184#if defined(ASSERTION_USED)
1185                    has_valid_translation_info = true;
1186#endif
1187                }
1188            }
1189
1190            if (!realigner.failure()) {
1191                char *buffer = realigner.realign_seq(allowed, source.data, source.len, dest.data, dest.len, cutoff_dna);
1192                if (buffer) { // re-alignment successful
1193                    error = GB_write_string(dest.gb_data, buffer);
1194
1195                    if (!error) {
1196                        int explicit_table_known = allowed.explicit_table();
1197
1198                        if (explicit_table_known >= 0) { // we know the exact code -> write codon_start and transl_table
1199                            const int codon_start  = 0; // by definition (after realignment)
1200                            error = AWT_saveTranslationInfo(gb_species, explicit_table_known, codon_start);
1201                        }
1202#if defined(ASSERTION_USED)
1203                        else { // we dont know the exact code -> can only happen if species has no translation info
1204                            nt_assert(allowed.any()); // bug in realigner
1205                            nt_assert(!has_valid_translation_info);
1206                        }
1207#endif
1208                    }
1209                    free(buffer);
1210                    if (!error && unmark_succeeded) GB_write_flag(gb_species, 0);
1211                }
1212            }
1213        }
1214
1215        if (realigner.failure()) {
1216            nt_assert(!error);
1217            GB_warningf("Automatic re-align failed for '%s'\nReason: %s", GBT_read_name(gb_species), realigner.failure());
1218        }
1219        else if (!error) {
1220            no_of_realigned_species++;
1221        }
1222
1223        progress.inc_and_check_user_abort(error);
1224    }
1225
1226    neededLength = realigner.get_needed_dest_alilen();
1227
1228    if (no_of_marked_species == 0) {
1229        GB_warning("Please mark some species to realign them");
1230    }
1231    else if (no_of_realigned_species != no_of_marked_species) {
1232        long failed = no_of_marked_species-no_of_realigned_species;
1233        nt_assert(failed>0);
1234        if (no_of_realigned_species) {
1235            GB_warningf("%li marked species failed to realign (%li succeeded)", failed, no_of_realigned_species);
1236        }
1237        else {
1238            GB_warning("All marked species failed to realign");
1239        }
1240    }
1241
1242    if (error) progress.done();
1243    else error = GBT_check_data(gb_main,ali_dest);
1244
1245    return error;
1246}
1247
1248static void realign_event(AW_window *aww) {
1249    AW_root  *aw_root          = aww->get_root();
1250    char     *ali_source       = aw_root->awar(AWAR_TRANSPRO_DEST)->read_string();
1251    char     *ali_dest         = aw_root->awar(AWAR_TRANSPRO_SOURCE)->read_string();
1252    bool      unmark_succeeded = aw_root->awar(AWAR_REALIGN_UNMARK)->read_int();
1253    bool      cutoff_dna       = aw_root->awar(AWAR_REALIGN_CUTOFF)->read_int();
1254    size_t    neededLength     = 0;
1255    GBDATA   *gb_main          = GLOBAL.gb_main;
1256    GB_ERROR  error            = realign_marked(gb_main, ali_source, ali_dest, neededLength, unmark_succeeded, cutoff_dna);
1257
1258    if (!error && neededLength) {
1259        bool auto_inc_alisize = aw_root->awar(AWAR_REALIGN_INCALI)->read_int();
1260        if (auto_inc_alisize) {
1261            {
1262                GB_transaction ta(gb_main);
1263                error = ta.close(GBT_set_alignment_len(gb_main, ali_dest, neededLength));
1264            }
1265            if (!error) {
1266                aw_message(GBS_global_string("Alignment length of '%s' has been set to %zu\n"
1267                                             "running re-aligner again!",
1268                                             ali_dest, neededLength));
1269
1270                error = realign_marked(gb_main, ali_source, ali_dest, neededLength, unmark_succeeded, cutoff_dna);
1271                if (neededLength) {
1272                    error = GBS_global_string("internal error: neededLength=%zu (after autoinc)", neededLength);
1273                }
1274            }
1275        }
1276        else {
1277            GB_transaction ta(gb_main);
1278            long           destLen = GBT_get_alignment_len(gb_main, ali_dest);
1279            nt_assert(destLen>0 && size_t(destLen)<neededLength);
1280            error = GBS_global_string("Missing %zu columns in alignment '%s' (got=%li, need=%zu)\n"
1281                                      "(check toggle to permit auto-increment)",
1282                                      size_t(neededLength-destLen), ali_dest, destLen, neededLength);
1283        }
1284    }
1285
1286    if (error) aw_message(error);
1287    free(ali_dest);
1288    free(ali_source);
1289}
1290
1291AW_window *NT_create_realign_dna_window(AW_root *root) {
1292    AW_window_simple *aws = new AW_window_simple;
1293    aws->init(root, "REALIGN_DNA", "REALIGN DNA");
1294
1295    aws->load_xfig("transdna.fig");
1296
1297    aws->at("close");
1298    aws->callback(AW_POPDOWN);
1299    aws->create_button("CLOSE", "CLOSE", "C");
1300
1301    aws->callback(makeHelpCallback("realign_dna.hlp"));
1302    aws->at("help");
1303    aws->create_button("HELP", "HELP", "H");
1304
1305    aws->at("source");
1306#if defined(DEVEL_RALF)
1307    awt_create_ALI_selection_button(GLOBAL.gb_main, aws, AWAR_TRANSPRO_SOURCE, "dna=:rna="); // @@@ nonsense here - just testing awt_create_ALI_selection_button somewhere
1308#else // !defined(DEVEL_RALF)
1309    awt_create_ALI_selection_list(GLOBAL.gb_main, aws, AWAR_TRANSPRO_SOURCE, "dna=:rna=");
1310#endif
1311    aws->at("dest");
1312    awt_create_ALI_selection_list(GLOBAL.gb_main, aws, AWAR_TRANSPRO_DEST, "pro=:ami=");
1313
1314    aws->at("autolen"); aws->create_toggle(AWAR_REALIGN_INCALI);
1315    aws->at("unmark");  aws->create_toggle(AWAR_REALIGN_UNMARK);
1316    aws->at("cutoff");  aws->create_toggle(AWAR_REALIGN_CUTOFF);
1317
1318    aws->at("realign");
1319    aws->callback(realign_event);
1320    aws->create_autosize_button("REALIGN", "Realign marked species", "R");
1321
1322    return aws;
1323}
1324
1325
1326void NT_create_transpro_variables(AW_root *root, AW_default props) {
1327    root->awar_string(AWAR_TRANSPRO_SOURCE, "",         props);
1328    root->awar_string(AWAR_TRANSPRO_DEST,   "",         props);
1329    root->awar_string(AWAR_TRANSPRO_MODE,   "settings", props);
1330
1331    root->awar_int(AWAR_TRANSPRO_POS,    0, props);
1332    root->awar_int(AWAR_TRANSPRO_XSTART, 1, props);
1333    root->awar_int(AWAR_TRANSPRO_WRITE,  0, props);
1334    root->awar_int(AWAR_REALIGN_INCALI,  0, props);
1335    root->awar_int(AWAR_REALIGN_UNMARK,  0, props);
1336    root->awar_int(AWAR_REALIGN_CUTOFF,  0, props);
1337}
1338
1339// --------------------------------------------------------------------------------
1340
1341#ifdef UNIT_TESTS
1342#ifndef TEST_UNIT_H
1343#include <test_unit.h>
1344#endif
1345
1346#include <arb_handlers.h>
1347
1348static std::string msgs;
1349
1350static void msg_to_string(const char *msg) {
1351    msgs += msg;
1352    msgs += '\n';
1353}
1354
1355static const char *translation_info(GBDATA *gb_species) {
1356    int      arb_transl_table;
1357    int      codon_start;
1358    GB_ERROR error = AWT_getTranslationInfo(gb_species, arb_transl_table, codon_start);
1359
1360    static SmartCharPtr result;
1361
1362    if (error) result = GBS_global_string_copy("Error: %s", error);
1363    else       result = GBS_global_string_copy("t=%i,cs=%i", arb_transl_table, codon_start);
1364
1365    return &*result;
1366}
1367
1368static arb_handlers test_handlers = {
1369    msg_to_string,
1370    msg_to_string,
1371    msg_to_string,
1372    active_arb_handlers->status,
1373};
1374
1375#define DNASEQ(name) GB_read_char_pntr(GBT_find_sequence(GBT_find_species(gb_main, name), "ali_dna"))
1376#define PROSEQ(name) GB_read_char_pntr(GBT_find_sequence(GBT_find_species(gb_main, name), "ali_pro"))
1377
1378#define TRANSLATION_INFO(name) translation_info(GBT_find_species(gb_main, name))
1379
1380void TEST_realign() {
1381    arb_handlers *old_handlers = active_arb_handlers;
1382    ARB_install_handlers(test_handlers);
1383
1384    GB_shell  shell;
1385    GBDATA   *gb_main = GB_open("TEST_realign.arb", "rw");
1386
1387    arb_suppress_progress here;
1388    enum TransResult { SAME, CHANGED };
1389
1390    {
1391        GB_ERROR error;
1392        size_t   neededLength = 0;
1393
1394        {
1395            struct transinfo_check {
1396                const char  *species_name;
1397                const char  *old_info;
1398                TransResult  changed;
1399                const char  *new_info;
1400            };
1401
1402            transinfo_check info[] = {
1403                { "BctFra12", "t=0,cs=1",  SAME,    NULL },        // fails -> unchanged
1404                { "CytLyti6", "t=9,cs=1",  CHANGED, "t=9,cs=0" },
1405                { "TaxOcell", "t=14,cs=1", CHANGED, "t=14,cs=0" },
1406                { "StrRamo3", "t=0,cs=1",  SAME,    NULL },        // fails -> unchanged
1407                { "StrCoel9", "t=0,cs=0",  SAME,    NULL },        // already correct
1408                { "MucRacem", "t=0,cs=1",  CHANGED, "t=0,cs=0" },
1409                { "MucRace2", "t=0,cs=1",  CHANGED, "t=0,cs=0" },
1410                { "MucRace3", "t=0,cs=0",  SAME,    NULL },        // fails -> unchanged
1411                { "AbdGlauc", "t=0,cs=0",  SAME,    NULL },        // already correct
1412                { "CddAlbic", "t=0,cs=0",  SAME,    NULL },        // already correct
1413
1414                { NULL, NULL, SAME, NULL }
1415            };
1416
1417            {
1418                GB_transaction ta(gb_main);
1419
1420                for (int i = 0; info[i].species_name; ++i) {
1421                    const transinfo_check& I = info[i];
1422                    TEST_ANNOTATE(I.species_name);
1423                    TEST_EXPECT_EQUAL(TRANSLATION_INFO(I.species_name), I.old_info);
1424                }
1425            }
1426            TEST_ANNOTATE(NULL);
1427
1428            msgs  = "";
1429            error = realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1430            TEST_EXPECT_NO_ERROR(error);
1431            TEST_EXPECT_EQUAL(msgs,
1432                              "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)
1433                              "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)
1434                              "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
1435                              "3 marked species failed to realign (7 succeeded)\n"
1436                );
1437
1438            {
1439                GB_transaction ta(gb_main);
1440
1441                TEST_EXPECT_EQUAL(DNASEQ("BctFra12"),    "ATGGCTAAAGAGAAATTTGAACGTACCAAACCGCACGTAAACATTGGTACAATCGGTCACGTTGACCACGGTAAAACCACTTTGACTGCTGCTATCACTACTGTGTTG------------------"); // failed = > seq unchanged
1442                TEST_EXPECT_EQUAL(DNASEQ("CytLyti6"),    "-A-TGGCAAAGGAAACTTTTGATCGTTCCAAACCGCACTTAA---ATATAG---GTACTATTGGACACGTAGATCACGGTAAAACTACTTTAACTGCTGCTATTACAASAGTAT-T-----G....");
1443                TEST_EXPECT_EQUAL(DNASEQ("TaxOcell"),    "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G..........");
1444                TEST_EXPECT_EQUAL(DNASEQ("StrRamo3"),    "ATGTCCAAGACGGCATACGTGCGCACCAAACCGCATCTGAACATCGGCACGATGGGTCATGTCGACCACGGCAAGACCACGTTGACCGCCGCCATCACCAAGGTCCTC------------------"); // failed = > seq unchanged
1445                TEST_EXPECT_EQUAL(DNASEQ("StrCoel9"),    "ATGTCCAAGACGGCGTACGTCCGC-C--C--A-CC-TG--A----GGCACGATG-G-CC--C-GACCACGGCAAGACCACCCTGACCGCCGCCATCACCAAGGTC-C--T--------C.......");
1446                TEST_EXPECT_EQUAL(DNASEQ("MucRacem"),    "......ATGGGTAAAGAG---------AAGACTCACGTTAACGTCGTCGTCATTGGTCACGTCGATTCCGGTAAATCTACTACTACTGGTCACTTGATTTACAAGTGTGGTGGTATA-AA......");
1447                TEST_EXPECT_EQUAL(DNASEQ("MucRace2"),    "ATGGGTAAGGAG---------AAGACTCACGTTAACGTCGTCGTCATTGGTCACGTCGATTCCGGTAAATCTACTACTACTGGTCACTTGATTTACAAGTGTGGTGGT-ATNNNAT-AAA......");
1448                TEST_EXPECT_EQUAL(DNASEQ("MucRace3"),    "-----------ATGGGTAAAGAGAAGACTCACGTTRAYGTTGTCGTTATTGGTCACGTCRATTCCGGTAAGTCCACCNCCRCTGGTCACTTGATTTACAAGTGTGGTGGTATAA-A----------"); // failed = > seq unchanged
1449                TEST_EXPECT_EQUAL(DNASEQ("AbdGlauc"),    "ATGGGTAAA-G--A--A--A--A--G-AC--T-CACGTTAACGTCGTTGTCATTGGTCACGTCGATTCTGGTAAATCCACCACCACTGGTCATTTGATCTACAAGTGCGGTGGTATA-AA......");
1450                TEST_EXPECT_EQUAL(DNASEQ("CddAlbic"),    "ATG-GG-TAAA-GAA------------AAAACTCACGTTAACGTTGTTGTTATTGGTCACGTCGATTCCGGTAAATCTACTACCACCGGTCACTTAATTTACAAGTGTGGTGGTATA-AA......");
1451                // ------------------------------------- "123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123"
1452
1453                for (int i = 0; info[i].species_name; ++i) {
1454                    const transinfo_check& I = info[i];
1455                    TEST_ANNOTATE(I.species_name);
1456                    switch (I.changed) {
1457                        case SAME:
1458                            TEST_EXPECT_EQUAL(TRANSLATION_INFO(I.species_name), I.old_info);
1459                            TEST_EXPECT_NULL(I.new_info);
1460                            break;
1461                        case CHANGED:
1462                            TEST_EXPECT_EQUAL(TRANSLATION_INFO(I.species_name), I.new_info);
1463                            TEST_EXPECT_DIFFERENT(I.new_info, I.old_info);
1464                            break;
1465                    }
1466                }
1467                TEST_ANNOTATE(NULL);
1468            }
1469        }
1470
1471        // test translation of sucessful realignments (see previous section)
1472        {
1473            GB_transaction ta(gb_main);
1474
1475            struct translate_check {
1476                const char *species_name;
1477                const char *original_prot;
1478                TransResult retranslation;
1479                const char *changed_prot; // if changed by translation (NULL for SAME)
1480            };
1481
1482            translate_check trans[] = {
1483                { "CytLyti6", "XWQRKLLIVPNRT*-I*-VLLDT*ITVKLL*SSLLZZYX-X.",
1484                  CHANGED,    "XWQRKLLIVPNRT*-I*-VLLDT*ITVKLL*SSLLQZYX-X." }, // ok: one of the Zs near end translates to Q
1485                { "TaxOcell", "XG*SNFWPVQAARNHRHD--RSRGPRQBDSDRCYHHGAX-..",
1486                  CHANGED,    "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - only changes gaptype at EOS
1487                { "MucRacem", "..MGKE---KTHVNVVVIGHVDSGKSTTTGHLIYKCGGIX..", SAME, NULL },
1488                { "MucRace2", "MGKE---KTHVNVVVIGHVDSGKSTTTGHLIYKCGGXXXK--",
1489                  CHANGED,    "MGKE---KTHVNVVVIGHVDSGKSTTTGHLIYKCGGXXXK.." }, // ok - only changes gaptype at EOS
1490                { "AbdGlauc", "MGKXXXXXXXXHVNVVVIGHVDSGKSTTTGHLIYKCGGIX..", SAME, NULL },
1491                { "StrCoel9", "MSKTAYVRXXXXXX-GTMXXXDHGKTTLTAAITKVXX--X..", SAME, NULL },
1492                { "CddAlbic", "MXXXE----KTHVNVVVIGHVDSGKSTTTGHLIYKCGGIX..", SAME, NULL },
1493
1494                { NULL, NULL, SAME, NULL }
1495            };
1496
1497            // check original protein sequences
1498            for (int t = 0; trans[t].species_name; ++t) {
1499                const translate_check& T = trans[t];
1500                TEST_ANNOTATE(T.species_name);
1501                TEST_EXPECT_EQUAL(PROSEQ(T.species_name), T.original_prot);
1502            }
1503            TEST_ANNOTATE(NULL);
1504
1505            msgs  = "";
1506            error = arb_r2a(gb_main, true, false, 0, true, "ali_dna", "ali_pro");
1507            TEST_EXPECT_NO_ERROR(error);
1508            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");
1509
1510            // check re-translated protein sequences
1511            for (int t = 0; trans[t].species_name; ++t) {
1512                const translate_check& T = trans[t];
1513                TEST_ANNOTATE(T.species_name);
1514                switch (T.retranslation) {
1515                    case SAME:
1516                        TEST_EXPECT_NULL(T.changed_prot);
1517                        TEST_EXPECT_EQUAL(PROSEQ(T.species_name), T.original_prot);
1518                        break;
1519                    case CHANGED:
1520                        TEST_REJECT_NULL(T.changed_prot);
1521                        TEST_EXPECT_DIFFERENT(T.original_prot, T.changed_prot);
1522                        TEST_EXPECT_EQUAL(PROSEQ(T.species_name), T.changed_prot);
1523                        break;
1524                }
1525            }
1526            TEST_ANNOTATE(NULL);
1527
1528            ta.close("dont commit");
1529        }
1530
1531        // -----------------------------
1532        //      provoke some errors
1533
1534        GBDATA *gb_TaxOcell;
1535        // unmark all but gb_TaxOcell
1536        {
1537            GB_transaction ta(gb_main);
1538
1539            gb_TaxOcell = GBT_find_species(gb_main, "TaxOcell");
1540            TEST_REJECT_NULL(gb_TaxOcell);
1541
1542            GBT_mark_all(gb_main, 0);
1543            GB_write_flag(gb_TaxOcell, 1);
1544        }
1545
1546        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1547
1548        // wrong alignment type
1549        {
1550            msgs  = "";
1551            error = realign_marked(gb_main, "ali_dna", "ali_pro", neededLength, false, false);
1552            TEST_EXPECT_ERROR_CONTAINS(error, "Invalid source alignment type");
1553            TEST_EXPECT_EQUAL(msgs, "");
1554        }
1555
1556        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1557
1558        GBDATA *gb_TaxOcell_amino;
1559        GBDATA *gb_TaxOcell_dna;
1560        {
1561            GB_transaction ta(gb_main);
1562            gb_TaxOcell_amino = GBT_find_sequence(gb_TaxOcell, "ali_pro");
1563            gb_TaxOcell_dna   = GBT_find_sequence(gb_TaxOcell, "ali_dna");
1564        }
1565        TEST_REJECT_NULL(gb_TaxOcell_amino);
1566        TEST_REJECT_NULL(gb_TaxOcell_dna);
1567
1568        // -----------------------------------------
1569        //      document some existing behavior
1570        {
1571            struct realign_check {
1572                const char  *seq;
1573                const char  *result;
1574                bool         cutoff;
1575                TransResult  retranslation;
1576                const char  *changed_prot; // if changed by translation (NULL for SAME)
1577            };
1578
1579            realign_check seq[] = {
1580                //"XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-.." // original aa sequence
1581                // { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..", "sdfjlksdjf" }, // templ
1582                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G..........", false, CHANGED, // original
1583                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - only changes gaptype at EOS
1584
1585                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHG.....", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCTG...........", false, CHANGED, // missing some AA at right end (extra DNA gets no longer truncated!)
1586                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - adds translation of extra DNA (DNA should never be modified by realigner!)
1587                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHG.....", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGT...............", true,  SAME, NULL }, // missing some AA at right end -> cutoff DNA
1588
1589                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYH-----..", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCTG...........", false, CHANGED,
1590                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - adds translation of extra DNA
1591                { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCY---H....", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTAT---------CACCACGGTGCTG..", false, CHANGED, // rightmost possible position of 'H' (see failing test below)
1592                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCY---HHGAX" }, // ok - adds translation of extra DNA
1593
1594                { "---SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..", "-ATGGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G..........", false, CHANGED, // missing some AA at left end (extra DNA gets detected now)
1595                  "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX..." }, // ok - adds translation of extra DNA (start of alignment)
1596                { "...SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX...", ".........AGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G..........", true,  SAME, NULL }, // missing some AA at left end -> cutoff DNA
1597
1598
1599                { "XG*SNFXXXXXXAXXXNHRHDXXXXXXPRQNDSDRCYHHGAX", "AT-GGCTAAAGAAACTTT-TG-AC-CG-GT-CCAA-GCC-GC-ACGT-AAACATCGGCACGAT-CG-GT-CA-CG-TGGA-CCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G.", false, SAME, NULL },
1600                { "XG*SNFWPVQAARNHRHD-XXXXXX-PRQNDSDRCYHHGAX-", "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT---CG-GT-CA-CG-TG-GA----CCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G....", false, CHANGED,
1601                  "XG*SNFWPVQAARNHRHD-XXXXXX-PRQNDSDRCYHHGAX." }, // ok - only changes gaptype at EOS
1602                { "XG*SNXLXRXQA-ARNHRHD-RXXVX-PRQNDSDRCYHHGAX", "AT-GGCTAAAGAAACTT-TTGAC-CGGTC-CAAGCC---GCACGTAAACATCGGCACGAT---CGG-TCAC-GTG-GA---CCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G.", false, SAME, NULL },
1603                { "XG*SXXFXDXVQAXT*TSARXRSXVX-PRQNDSDRCYHHGAX", "AT-GGCTAAAGA-A-AC-TTT-T-GACCG-GTCCAAGCCGC-ACGTAAACATCGGCACGA-T-CGGTCA-C-GTG-GA---CCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G.", false, SAME, NULL },
1604                // -------------------------------------------- "123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123123"
1605
1606                { 0, 0, false, SAME, NULL }
1607            };
1608
1609            int   arb_transl_table, codon_start;
1610            char *org_dna;
1611            {
1612                GB_transaction ta(gb_main);
1613                TEST_EXPECT_NO_ERROR(AWT_getTranslationInfo(gb_TaxOcell, arb_transl_table, codon_start));
1614                TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), "t=14,cs=0");
1615                org_dna = GB_read_string(gb_TaxOcell_dna);
1616            }
1617
1618            for (int s = 0; seq[s].seq; ++s) {
1619                TEST_ANNOTATE(GBS_global_string("s=%i", s));
1620                realign_check& S = seq[s];
1621
1622                {
1623                    GB_transaction ta(gb_main);
1624                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_amino, S.seq));
1625                }
1626                msgs  = "";
1627                error = realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, S.cutoff);
1628                TEST_EXPECT_NO_ERROR(error);
1629                TEST_EXPECT_EQUAL(msgs, "");
1630                {
1631                    GB_transaction ta(gb_main);
1632                    TEST_EXPECT_EQUAL(GB_read_char_pntr(gb_TaxOcell_dna), S.result);
1633
1634                    // test retranslation:
1635                    msgs  = "";
1636                    error = arb_r2a(gb_main, true, false, 0, true, "ali_dna", "ali_pro");
1637                    TEST_EXPECT_NO_ERROR(error);
1638                    if (s == 10) {
1639                        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");
1640                    }
1641                    else if (s == 6) {
1642                        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");
1643                    }
1644                    else {
1645                        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");
1646                    }
1647
1648                    switch (S.retranslation) {
1649                        case SAME:
1650                            TEST_EXPECT_NULL(S.changed_prot);
1651                            TEST_EXPECT_EQUAL(GB_read_char_pntr(gb_TaxOcell_amino), S.seq);
1652                            break;
1653                        case CHANGED:
1654                            TEST_REJECT_NULL(S.changed_prot);
1655                            TEST_EXPECT_EQUAL(GB_read_char_pntr(gb_TaxOcell_amino), S.changed_prot);
1656                            break;
1657                    }
1658
1659                    TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), "t=14,cs=0");
1660                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_dna, org_dna)); // restore changed DB entry
1661                }
1662            }
1663            TEST_ANNOTATE(NULL);
1664
1665            free(org_dna);
1666        }
1667
1668        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1669
1670        // ----------------------------------------------------
1671        //      write some aa sequences provoking failures
1672        {
1673            struct realign_fail {
1674                const char *seq;
1675                const char *failure;
1676            };
1677
1678#define ERRPREFIX     "Automatic re-align failed for 'TaxOcell'\nReason: "
1679#define ERRPREFIX_LEN 49
1680
1681#define FAILONE "All marked species failed to realign\n"
1682
1683            // dna of TaxOcell:
1684            // "AT-GGCTAAAGAAACTTTTGACCGGTCCAAGCCGCACGTAAACATCGGCACGAT------CGGTCACGTGGACCACGGCAAAACGACTCTGACCGCTGCTATCACCACGGTGCT-G----......"
1685
1686            realign_fail seq[] = {
1687                //"XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-.." // original aa sequence
1688                // { "XG*SNFWPVQAARNHRHD--RSRGPRQNDSDRCYHHGAX-..", "sdfjlksdjf" }, // templ
1689
1690                // wanted realign failures:
1691                { "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
1692                { "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
1693                { "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
1694                { "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
1695                { "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
1696                { "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
1697                { "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
1698                { "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'
1699                { "--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)
1700
1701                // failing realignments that should work:
1702
1703                { 0, 0 }
1704            };
1705
1706            {
1707                GB_transaction ta(gb_main);
1708                TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), "t=14,cs=0");
1709            }
1710
1711            for (int s = 0; seq[s].seq; ++s) {
1712                TEST_ANNOTATE(GBS_global_string("s=%i", s));
1713                {
1714                    GB_transaction ta(gb_main);
1715                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_amino, seq[s].seq));
1716                }
1717                msgs  = "";
1718                error = realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1719                TEST_EXPECT_NO_ERROR(error);
1720                TEST_EXPECT_CONTAINS(msgs, ERRPREFIX);
1721                TEST_EXPECT_EQUAL(msgs.c_str()+ERRPREFIX_LEN, seq[s].failure);
1722
1723                {
1724                    GB_transaction ta(gb_main);
1725                    TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), "t=14,cs=0"); // should not change if error
1726                }
1727            }
1728            TEST_ANNOTATE(NULL);
1729        }
1730
1731        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1732
1733        // ----------------------------------------------
1734        //      some examples for given DNA/AA pairs
1735
1736        {
1737            struct explicit_realign {
1738                const char *acids;
1739                const char *dna;
1740                int         table;
1741                const char *info;
1742                const char *msgs;
1743            };
1744
1745            // YTR (=X(2,9,16), =L(else))
1746            //     CTA (=T(2),        =L(else))
1747            //     CTG (=T(2), =S(9), =L(else))
1748            //     TTA (=*(16),       =L(else))
1749            //     TTG (=L(always))
1750            //
1751            // AAR (=X(6,11,14), =K(else))
1752            //     AAA (=N(6,11,14), =K(else))
1753            //     AAG (=K(always))
1754            //
1755            // ATH (=X(1,2,4,10,14), =I(else))
1756            //     ATA (=M(1,2,4,10,14), =I(else))
1757            //     ATC (=I(always))
1758            //     ATT (=I(always))
1759
1760            const char*const NO_TI = "t=-1,cs=-1";
1761
1762            explicit_realign example[] = {
1763                { "LK", "TTGAAG", -1, NO_TI,        NULL }, // fine (for any table)
1764
1765                { "G",  "RGG",    -1, "t=10,cs=0",  NULL }, // correctly detects TI(10)
1766
1767                { "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)
1768                { "LX",  "YTRAAR",    -1, NO_TI,       NULL }, // fine (AAR->X for table=6,11,14)
1769                { "LXX", "YTRAARATH", -1, "t=14,cs=0", NULL }, // correctly detects TI(14)
1770                { "LXI", "YTRAARATH", -1, NO_TI,       NULL }, // fine (for table=6,11)
1771
1772                { "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)
1773                { "LK", "YTRAAR", -1, NO_TI,        NULL }, // fine           (AAR->K for table!=6,11,14)
1774                { "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)
1775                { "XK", "YTRAAR", -1, NO_TI,        NULL }, // fine           (YTR->X for table=2,9,16)
1776
1777                { "XX",   "-YTRAAR",      0,  "t=0,cs=0", NULL },                                                                                             // does not fail because it realigns such that it translates back to 'XXX'
1778                { "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 )
1779                { "-XXL", "-YTRA-AR-TTG", 0,  "t=0,cs=0", NULL },                                                                                             // does not fail because it realigns such that it translates back to 'XXXL'
1780                { "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)
1781                { "XX",   "-YTRAAR",      -1, NO_TI,      NULL },                                                                                             // does not fail because it realigns such that it translates back to 'XXX'
1782                { "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)
1783
1784                { "LX", "YTRATH", -1, NO_TI,        NULL }, // fine                (ATH->X for table=1,2,4,10,14)
1785                { "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)
1786                { "XX", "YTRATH", 2,  "t=2,cs=0",   NULL }, // fine                (both->X for table=2)
1787                { "XX", "YTRATH", -1, "t=2,cs=0",   NULL }, // correctly detects TI(2)
1788
1789                { "XX", "AARATH", 14, "t=14,cs=0",  NULL }, // fine (both->X for table=14)
1790                { "XX", "AARATH", -1, "t=14,cs=0",  NULL }, // correctly detects TI(14)
1791                { "KI", "AARATH", -1, NO_TI,        NULL }, // fine (for table!=1,2,4,6,10,11,14)
1792                { "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)
1793                { "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)
1794                { "KX", "AARATH", -1, NO_TI,        NULL }, // fine for table=1,2,4,10
1795                { "KX", "AARATH", 4,  "t=4,cs=0",   NULL }, // test table=4
1796                { "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)
1797                { "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)
1798
1799                { NULL, NULL, 0, NULL, NULL }
1800            };
1801
1802            for (int e = 0; example[e].acids; ++e) {
1803                const explicit_realign& E = example[e];
1804                TEST_ANNOTATE(GBS_global_string("%s <- %s (#%i)", E.acids, E.dna, E.table));
1805
1806                {
1807                    GB_transaction ta(gb_main);
1808                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_dna, E.dna));
1809                    TEST_EXPECT_NO_ERROR(GB_write_string(gb_TaxOcell_amino, E.acids));
1810                    if (E.table == -1) {
1811                        TEST_EXPECT_NO_ERROR(AWT_removeTranslationInfo(gb_TaxOcell));
1812                    }
1813                    else {
1814                        TEST_EXPECT_NO_ERROR(AWT_saveTranslationInfo(gb_TaxOcell, E.table, 0));
1815                    }
1816                }
1817
1818                msgs  = "";
1819                error = realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1820                TEST_EXPECT_NULL(error);
1821                if (E.msgs) {
1822                    TEST_EXPECT_CONTAINS(msgs, ERRPREFIX);
1823                    string wanted_msgs = string(E.msgs)+FAILONE;
1824                    TEST_EXPECT_EQUAL(msgs.c_str()+ERRPREFIX_LEN, wanted_msgs);
1825                }
1826                else {
1827                    TEST_EXPECT_EQUAL(msgs, "");
1828                }
1829
1830                GB_transaction ta(gb_main);
1831                if (!error) {
1832                    const char *dnaseq      = GB_read_char_pntr(gb_TaxOcell_dna);
1833                    size_t      expextedLen = strlen(E.dna);
1834                    size_t      seqlen      = strlen(dnaseq);
1835                    char       *firstPart   = ARB_strndup(dnaseq, expextedLen);
1836                    size_t      dna_behind;
1837                    char       *nothing     = unalign(dnaseq+expextedLen, seqlen-expextedLen, dna_behind);
1838
1839                    TEST_EXPECT_EQUAL(firstPart, E.dna);
1840                    TEST_EXPECT_EQUAL(dna_behind, 0);
1841                    TEST_EXPECT_EQUAL(nothing, "");
1842
1843                    free(nothing);
1844                    free(firstPart);
1845                }
1846                TEST_EXPECT_EQUAL(translation_info(gb_TaxOcell), E.info);
1847            }
1848        }
1849
1850        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1851
1852        // ----------------------------------
1853        //      invalid translation info
1854        {
1855            GB_transaction ta(gb_main);
1856
1857            TEST_EXPECT_NO_ERROR(AWT_saveTranslationInfo(gb_TaxOcell, 14, 0));
1858            GBDATA *gb_trans_table = GB_entry(gb_TaxOcell, "transl_table");
1859            TEST_EXPECT_NO_ERROR(GB_write_string(gb_trans_table, "666")); // evil translation table
1860        }
1861
1862        msgs  = "";
1863        error = realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1864        TEST_EXPECT_NO_ERROR(error);
1865        TEST_EXPECT_EQUAL(msgs, ERRPREFIX "Error while reading 'transl_table' (Illegal (or unsupported) value (666) in 'transl_table' (item='TaxOcell'))\n" FAILONE);
1866        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1867
1868        // ---------------------------------------
1869        //      source/dest alignment missing
1870        for (int i = 0; i<2; ++i) {
1871            TEST_ANNOTATE(GBS_global_string("i=%i", i));
1872
1873            {
1874                GB_transaction ta(gb_main);
1875
1876                GBDATA *gb_ali = GB_get_father(GBT_find_sequence(gb_TaxOcell, i ? "ali_pro" : "ali_dna"));
1877                GB_push_my_security(gb_main);
1878                TEST_EXPECT_NO_ERROR(GB_delete(gb_ali));
1879                GB_pop_my_security(gb_main);
1880            }
1881
1882            msgs  = "";
1883            error = realign_marked(gb_main, "ali_pro", "ali_dna", neededLength, false, false);
1884            TEST_EXPECT_NO_ERROR(error);
1885            if (i) {
1886                TEST_EXPECT_EQUAL(msgs, ERRPREFIX "No data in alignment 'ali_pro'\n" FAILONE);
1887            }
1888            else {
1889                TEST_EXPECT_EQUAL(msgs, ERRPREFIX "No data in alignment 'ali_dna'\n" FAILONE);
1890            }
1891        }
1892        TEST_ANNOTATE(NULL);
1893
1894        TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 1);
1895    }
1896
1897#undef ERRPREFIX
1898#undef ERRPREFIX_LEN
1899
1900    GB_close(gb_main);
1901    ARB_install_handlers(*old_handlers);
1902}
1903
1904static const char *permOf(const Distributor& dist) {
1905    const int   MAXDIST = 10;
1906    static char buffer[MAXDIST+1];
1907
1908    nt_assert(dist.size() <= MAXDIST);
1909    for (int p = 0; p<dist.size(); ++p) {
1910        buffer[p] = '0'+dist[p];
1911    }
1912    buffer[dist.size()] = 0;
1913
1914    return buffer;
1915}
1916
1917static arb_test::match_expectation stateOf(Distributor& dist, const char *expected_perm, bool hasNext) {
1918    using namespace arb_test;
1919
1920    expectation_group expected;
1921    expected.add(that(permOf(dist)).is_equal_to(expected_perm));
1922    expected.add(that(dist.next()).is_equal_to(hasNext));
1923    return all().ofgroup(expected);
1924}
1925
1926void TEST_distributor() {
1927    TEST_EXPECT_EQUAL(Distributor(3, 2).get_error(), "not enough nucleotides");
1928    TEST_EXPECT_EQUAL(Distributor(3, 10).get_error(), "too much nucleotides");
1929
1930    Distributor minDist(3, 3);
1931    TEST_EXPECTATION(stateOf(minDist, "111", false));
1932
1933    Distributor maxDist(3, 9);
1934    TEST_EXPECTATION(stateOf(maxDist, "333", false));
1935
1936    Distributor meanDist(3, 6);
1937    TEST_EXPECTATION(stateOf(meanDist, "123", true));
1938    TEST_EXPECTATION(stateOf(meanDist, "132", true));
1939    TEST_EXPECTATION(stateOf(meanDist, "213", true));
1940    TEST_EXPECTATION(stateOf(meanDist, "222", true));
1941    TEST_EXPECTATION(stateOf(meanDist, "231", true));
1942    TEST_EXPECTATION(stateOf(meanDist, "312", true));
1943    TEST_EXPECTATION(stateOf(meanDist, "321", false));
1944
1945    Distributor belowMax(4, 11);
1946    TEST_EXPECTATION(stateOf(belowMax, "2333", true));
1947    TEST_EXPECTATION(stateOf(belowMax, "3233", true));
1948    TEST_EXPECTATION(stateOf(belowMax, "3323", true));
1949    TEST_EXPECTATION(stateOf(belowMax, "3332", false));
1950
1951    Distributor aboveMin(4, 6);
1952    TEST_EXPECTATION(stateOf(aboveMin, "1113", true));
1953    TEST_EXPECTATION(stateOf(aboveMin, "1122", true));
1954    TEST_EXPECTATION(stateOf(aboveMin, "1131", true));
1955    TEST_EXPECTATION(stateOf(aboveMin, "1212", true));
1956    TEST_EXPECTATION(stateOf(aboveMin, "1221", true));
1957    TEST_EXPECTATION(stateOf(aboveMin, "1311", true));
1958    TEST_EXPECTATION(stateOf(aboveMin, "2112", true));
1959    TEST_EXPECTATION(stateOf(aboveMin, "2121", true));
1960    TEST_EXPECTATION(stateOf(aboveMin, "2211", true));
1961    TEST_EXPECTATION(stateOf(aboveMin, "3111", false));
1962
1963    Distributor check(6, 8);
1964    TEST_EXPECTATION(stateOf(check, "111113", true));
1965    TEST_EXPECTATION(stateOf(check, "111122", true));
1966    TEST_EXPECTATION(stateOf(check, "111131", true));
1967    TEST_EXPECTATION(stateOf(check, "111212", true));
1968    TEST_EXPECTATION(stateOf(check, "111221", true));
1969    TEST_EXPECTATION(stateOf(check, "111311", true));
1970    TEST_EXPECTATION(stateOf(check, "112112", true));
1971    TEST_EXPECTATION(stateOf(check, "112121", true));
1972    TEST_EXPECTATION(stateOf(check, "112211", true));
1973    TEST_EXPECTATION(stateOf(check, "113111", true));
1974    TEST_EXPECTATION(stateOf(check, "121112", true));
1975    TEST_EXPECTATION(stateOf(check, "121121", true));
1976    TEST_EXPECTATION(stateOf(check, "121211", true));
1977    TEST_EXPECTATION(stateOf(check, "122111", true));
1978    TEST_EXPECTATION(stateOf(check, "131111", true));
1979    TEST_EXPECTATION(stateOf(check, "211112", true));
1980    TEST_EXPECTATION(stateOf(check, "211121", true));
1981    TEST_EXPECTATION(stateOf(check, "211211", true));
1982    TEST_EXPECTATION(stateOf(check, "212111", true));
1983    TEST_EXPECTATION(stateOf(check, "221111", true));
1984    TEST_EXPECTATION(stateOf(check, "311111", false));
1985}
1986
1987#endif // UNIT_TESTS
1988
1989// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.