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