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