Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

TranslateRealign.cxx

Visit:

Last change on this file was 18125, checked in by westram, 5 years ago
full update from 'trunk' into 'nameserver' adds: log:branches/cleanup@17627:17650,18105:18112 log:branches/clone@17813:17835 log:branches/config@17733:17737 log:branches/dbui@17881:17889 log:branches/fix@17597:17617 log:branches/fts@17653:18088,18090:18094 log:branches/gcc@17850:17876 log:branches/gui@17754:18092 log:branches/mix@17631:17633,17637:17640 log:trunk@17594:18123
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 80.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: