Context Navigation

source: branches/profile/SL/FAST_ALIGNER/fast_aligner.cxx

Visit:

Last change on this file was 12831, checked in by westram, 10 years ago
merge [12806:12810] from 'fix' into 'trunk' eliminates GBT_write_sequence adds: log:branches/fix@12807:12810
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 134.6 KB

Line
1	// =============================================================== //
2	// //
3	// File : fast_aligner.cxx //
4	// Purpose : A fast aligner (not a multiple aligner!) //
5	// //
6	// Coded by Ralf Westram (coder@reallysoft.de) in 1998 //
7	// Institute of Microbiology (Technical University Munich) //
8	// http://www.arb-home.de/ //
9	// //
10	// =============================================================== //
11
12
13	#include "fast_aligner.hxx"
14	#include "ClustalV.hxx"
15	#include "seq_search.hxx"
16
17	#include <island_hopping.h>
18
19	#include <awtc_next_neighbours.hxx>
20	#include <awt_sel_boxes.hxx>
21
22	#include <aw_awars.hxx>
23	#include <aw_root.hxx>
24	#include <aw_question.hxx>
25
26	#include <arbdbt.h>
27	#include <ad_cb.h>
28
29	#include <arb_defs.h>
30	#include <arb_progress.h>
31	#include <RangeList.h>
32
33	#include <cctype>
34	#include <cmath>
35	#include <climits>
36
37	#include <list>
38
39	// --------------------------------------------------------------------------------
40
41	#if defined(DEBUG)
42	// #define TRACE_CLUSTAL_DATA
43	// #define TRACE_ISLANDHOPPER_DATA
44	// #define TRACE_COMPRESSED_ALIGNMENT
45	// #define TRACE_RELATIVES
46	#endif // DEBUG
47
48	// --------------------------------------------------------------------------------
49
50	enum FA_report {
51	FA_NO_REPORT, // no report
52	FA_TEMP_REPORT, // report to temporary entries
53	FA_REPORT, // report to resident entries
54	};
55
56	enum FA_range {
57	FA_WHOLE_SEQUENCE, // align whole sequence
58	FA_AROUND_CURSOR, // align xxx positions around current cursor position
59	FA_SELECTED_RANGE, // align selected range
60	FA_SAI_MULTI_RANGE, // align versus multi range define by SAI
61	};
62
63	enum FA_turn {
64	FA_TURN_NEVER, // never try to turn sequence
65	FA_TURN_INTERACTIVE, // try to turn, but query user
66	FA_TURN_ALWAYS, // turn if score is better
67	};
68
69	enum FA_reference {
70	FA_REF_EXPLICIT, // reference sequence explicitly specified
71	FA_REF_CONSENSUS, // use group consensus as reference
72	FA_REF_RELATIVES, // search next relatives by PT server
73	};
74
75	enum FA_alignTarget {
76	FA_CURRENT, // align current species
77	FA_MARKED, // align all marked species
78	FA_SELECTED, // align selected species (= range)
79	};
80
81	enum FA_errorAction {
82	FA_NO_ACTION, // do nothing
83	FA_MARK_FAILED, // mark failed species (unmark rest)
84	FA_MARK_ALIGNED, // mark aligned species (unmark rest)
85	};
86
87	struct AlignParams {
88	FA_report report;
89	bool showGapsMessages; // display messages about missing gaps in master?
90	PosRange range; // range to be aligned
91	};
92
93	struct SearchRelativeParams : virtual Noncopyable {
94	FamilyFinder *ff;
95	char *pt_server_alignment; // alignment used in pt_server (may differ from 'alignment')
96	int maxRelatives; // max # of relatives to use
97
98	SearchRelativeParams(FamilyFinder ff_, const char pt_server_alignment_, int maxRelatives_)
99	: ff(ff_),
100	pt_server_alignment(strdup(pt_server_alignment_)),
101	maxRelatives(maxRelatives_)
102	{}
103
104	~SearchRelativeParams() {
105	free(pt_server_alignment);
106	delete(ff);
107	}
108
109	FamilyFinder *getFamilyFinder() { return ff; }
110	};
111
112	// --------------------------------------------------------------------------------
113
114	#define GAP_CHAR '-'
115	#define QUALITY_NAME "ASC_ALIGNER_CLIENT_SCORE"
116	#define INSERTS_NAME "AMI_ALIGNER_MASTER_INSERTS"
117
118	#define FA_AWAR_ROOT "faligner/"
119	#define FA_AWAR_TO_ALIGN FA_AWAR_ROOT "what"
120	#define FA_AWAR_REFERENCE FA_AWAR_ROOT "against"
121	#define FA_AWAR_REFERENCE_NAME FA_AWAR_ROOT "sagainst"
122	#define FA_AWAR_RANGE FA_AWAR_ROOT "range"
123	#define FA_AWAR_PROTECTION FA_AWAR_ROOT "protection"
124	#define FA_AWAR_AROUND FA_AWAR_ROOT "around"
125	#define FA_AWAR_MIRROR FA_AWAR_ROOT "mirror"
126	#define FA_AWAR_REPORT FA_AWAR_ROOT "report"
127	#define FA_AWAR_SHOW_GAPS_MESSAGES FA_AWAR_ROOT "show_gaps"
128	#define FA_AWAR_CONTINUE_ON_ERROR FA_AWAR_ROOT "continue_on_error"
129	#define FA_AWAR_ACTION_ON_ERROR FA_AWAR_ROOT "action_on_error"
130	#define FA_AWAR_USE_SECONDARY FA_AWAR_ROOT "use_secondary"
131	#define FA_AWAR_NEXT_RELATIVES FA_AWAR_ROOT "next_relatives"
132	#define FA_AWAR_RELATIVE_RANGE FA_AWAR_ROOT "relrange"
133	#define FA_AWAR_PT_SERVER_ALIGNMENT "tmp/" FA_AWAR_ROOT "relative_ali"
134	#define FA_AWAR_SAI_RANGE_NAME FA_AWAR_ROOT "sai/sainame"
135	#define FA_AWAR_SAI_RANGE_CHARS FA_AWAR_ROOT "sai/chars"
136
137	#define FA_AWAR_ISLAND_HOPPING_ROOT "island_hopping/"
138	#define FA_AWAR_USE_ISLAND_HOPPING FA_AWAR_ISLAND_HOPPING_ROOT "use"
139	#define FA_AWAR_ESTIMATE_BASE_FREQ FA_AWAR_ISLAND_HOPPING_ROOT "estimate_base_freq"
140	#define FA_AWAR_BASE_FREQ_A FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_a"
141	#define FA_AWAR_BASE_FREQ_C FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_c"
142	#define FA_AWAR_BASE_FREQ_G FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_g"
143	#define FA_AWAR_BASE_FREQ_T FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_t"
144	#define FA_AWAR_SUBST_PARA_AC FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ac"
145	#define FA_AWAR_SUBST_PARA_AG FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ag"
146	#define FA_AWAR_SUBST_PARA_AT FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_at"
147	#define FA_AWAR_SUBST_PARA_CG FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_cg"
148	#define FA_AWAR_SUBST_PARA_CT FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ct"
149	#define FA_AWAR_SUBST_PARA_GT FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_gt"
150	#define FA_AWAR_EXPECTED_DISTANCE FA_AWAR_ISLAND_HOPPING_ROOT "expected_dist"
151	#define FA_AWAR_STRUCTURE_SUPPLEMENT FA_AWAR_ISLAND_HOPPING_ROOT "struct_suppl"
152	#define FA_AWAR_THRESHOLD FA_AWAR_ISLAND_HOPPING_ROOT "threshold"
153	#define FA_AWAR_GAP_A FA_AWAR_ISLAND_HOPPING_ROOT "gapa"
154	#define FA_AWAR_GAP_B FA_AWAR_ISLAND_HOPPING_ROOT "gapb"
155	#define FA_AWAR_GAP_C FA_AWAR_ISLAND_HOPPING_ROOT "gapc"
156
157	// --------------------------------------------------------------------------------
158
159	static IslandHopping *island_hopper = NULL; // NULL -> use fast aligner; else use island hopper
160
161	static GB_alignment_type global_alignmentType = GB_AT_UNKNOWN; // type of actually aligned sequence
162
163	static int currentSequenceNumber; // used for counter
164	static int overallSequenceNumber;
165
166	// --------------------------------------------------------------------------------
167
168	inline ARB_ERROR species_not_found(GB_CSTR species_name) {
169	return GBS_global_string("No species '%s' found!", species_name);
170	}
171
172	static ARB_ERROR reverseComplement(GBDATA *gb_species, GB_CSTR ali, int max_protection) {
173	GBDATA *gbd = GBT_find_sequence(gb_species, ali);
174	ARB_ERROR error = 0;
175
176	if (!gbd) {
177	error = GBS_global_string("No 'data' found for species '%s'", GBT_read_name(gb_species));
178	}
179	else {
180	int my_protection = GB_read_security_write(gbd);
181
182	if (my_protection<=max_protection) { // ok
183	char *seq = GB_read_string(gbd);
184	int length = GB_read_string_count(gbd);
185	GBDATA *gb_main = GB_get_root(gb_species);
186	GB_alignment_type ali_type = GBT_get_alignment_type(gb_main, ali);
187
188	char T_or_U;
189	error = GBT_determine_T_or_U(ali_type, &T_or_U, "reverse-complement");
190	if (!error) {
191	GBT_reverseComplementNucSequence(seq, length, T_or_U);
192	error = GB_write_string(gbd, seq);
193	}
194	}
195	else { // protection error
196	error = GBS_global_string("Cannot reverse-complement species '%s' because of protection level", GBT_read_name(gb_species));
197	}
198
199	}
200
201	return error;
202	}
203
204	static void build_reverse_complement(AW_window *aw, AW_CL cl_AlignDataAccess) {
205	const AlignDataAccess data_access = (const AlignDataAccess )cl_AlignDataAccess;
206	GBDATA *gb_main = data_access->gb_main;
207
208	GB_push_transaction(gb_main);
209
210	AW_root *root = aw->get_root();
211	FA_alignTarget revComplWhat = static_cast<FA_alignTarget>(root->awar(FA_AWAR_TO_ALIGN)->read_int());
212	char *default_alignment = GBT_get_default_alignment(gb_main);
213	GB_CSTR alignment = root->awar_string(AWAR_EDITOR_ALIGNMENT, default_alignment)->read_string();
214	ARB_ERROR error = 0;
215	int max_protection = root->awar(FA_AWAR_PROTECTION)->read_int();
216
217	switch (revComplWhat) {
218	case FA_CURRENT: { // current species
219	GB_CSTR species_name = root->awar(AWAR_SPECIES_NAME)->read_string();
220	GBDATA *gb_species = GBT_find_species(gb_main, species_name);
221	if (!gb_species) error = species_not_found(species_name);
222	if (!error) error = reverseComplement(gb_species, alignment, max_protection);
223	break;
224	}
225	case FA_MARKED: { // marked species
226	GBDATA *gb_species = GBT_first_marked_species(gb_main);
227
228	if (!gb_species) {
229	error = "There is no marked species";
230	}
231	while (gb_species) {
232	error = reverseComplement(gb_species, alignment, max_protection);
233	if (error) break;
234	gb_species = GBT_next_marked_species(gb_species);
235	}
236	break;
237	}
238	case FA_SELECTED: { // selected species (editor selection!)
239
240	Aligner_get_first_selected_species get_first_selected_species = data_access->get_first_selected_species;
241	Aligner_get_next_selected_species get_next_selected_species = data_access->get_next_selected_species;
242
243	int count = 0;
244	GBDATA *gb_species = get_first_selected_species(&count);
245
246	if (!gb_species) {
247	error = "There is no selected species!";
248	}
249	while (gb_species) {
250	error = reverseComplement(gb_species, alignment, max_protection);
251	if (error) break;
252	gb_species = get_next_selected_species();
253	}
254	break;
255	}
256	default: {
257	fa_assert(0);
258	break;
259	}
260	}
261
262	GB_end_transaction_show_error(gb_main, error, aw_message);
263	}
264
265	// --------------------------------------------------------------------------------
266
267	class AliChange { // describes a local alignment change
268	const CompactedSubSequence& Old;
269	const CompactedSubSequence& New;
270
271	public:
272	AliChange(const CompactedSubSequence& old_, const CompactedSubSequence& new_)
273	: Old(old_), New(new_)
274	{
275	fa_assert(Old.may_refer_to_same_part_as(New));
276	}
277
278	int follow(ExplicitRange& range) const {
279	ExplicitRange compressed(Old.compPosition(range.start()),
280	Old.compPosition(range.end()));
281
282	int exp_start = New.expdPosition(compressed.start());
283	int exp_end = New.expdPosition(compressed.end());
284
285	int gaps_before = New.no_of_gaps_before(compressed.start());
286	int gaps_after = New.no_of_gaps_after(compressed.end());
287
288	fa_assert(gaps_before >= 0);
289	fa_assert(gaps_after >= 0);
290
291	range = ExplicitRange(exp_start-gaps_before, exp_end+gaps_after);
292
293	return compressed.size(); // number of bases
294	}
295	};
296
297	class LooseBases {
298	typedef std::list<ExplicitRange> Ranges;
299
300	Ranges ranges;
301
302	public:
303
304	bool is_empty() const { return ranges.empty(); }
305	void clear() { ranges.clear(); }
306
307	void memorize(ExplicitRange range) {
308	ranges.push_front(range);
309	}
310	ExplicitRange recall() {
311	ExplicitRange range = ranges.front();
312	ranges.pop_front();
313	return range;
314	}
315	int follow_ali_change(const AliChange& change) {
316	// transform positions according to changed alignment (oldSequence -> newSequence)
317	// returns the number of bases contained in 'this'
318	int basecount = 0;
319	if (!is_empty()) {
320	for (Ranges::iterator r = ranges.begin(); r != ranges.end(); ++r) {
321	basecount += change.follow(*r);
322	}
323	}
324	return basecount;
325	}
326	void append(LooseBases& loose) { ranges.splice(ranges.end(), loose.ranges); }
327	int follow_ali_change_and_append(LooseBases& loose, const AliChange& change) {
328	// returns the number of loose bases (added from 'loose')
329	int basecount = loose.follow_ali_change(change);
330	append(loose);
331	return basecount;
332	}
333	};
334
335	static LooseBases unaligned_bases; // if fast_align cannot align (no master bases) -> stores positions here
336
337
338	static const char read_name(GBDATA gbd) {
339	return gbd ? GBT_read_name(gbd) : "GROUP-CONSENSUS";
340	}
341
342	inline int relatedBases(char base1, char base2) {
343	return baseMatch(base1, base2)==1;
344	}
345
346	inline char alignQuality(char slave, char master) {
347	fa_assert(slave);
348	fa_assert(master);
349	char result = '#';
350	if (slave==master) result = '-'; // equal
351	else if (slave==GAP_CHAR) result = '+'; // inserted gap
352	else if (master==GAP_CHAR) result = '+'; // no gap in master
353	else if (relatedBases(slave, master)) result = '~'; // mutation (related bases)
354	return result; // mutation (non-related bases)
355	}
356
357	// -------------------------
358	// Debugging stuff
359
360	#ifdef DEBUG
361	static char lstr(const char s, int len) {
362	static int alloc;
363	static char *buffer;
364
365	if (alloc<(len+1)) {
366	if (alloc) free(buffer);
367	buffer = (char*)malloc(alloc=len+100);
368	}
369
370	memcpy(buffer, s, len);
371	buffer[len] = 0;
372
373	return buffer;
374	}
375
376	#define BUFLEN 120
377
378	inline char compareChar(char base1, char base2) {
379	return base1==base2 ? '=' : (relatedBases(base1, base2) ? 'x' : 'X');
380	}
381
382	#if defined(TRACE_COMPRESSED_ALIGNMENT)
383
384	static void dump_n_compare_one(const char seq1, const char seq2, long len, long offset) {
385	fa_assert(len<=BUFLEN);
386	char compare[BUFLEN+1];
387
388	for (long l=0; l<len; l++) {
389	compare[l] = (is_ali_gap(seq1[l]) && is_ali_gap(seq2[l])) ? ' ' : compareChar(seq1[l], seq2[l]);
390	}
391
392	compare[len] = 0;
393
394	printf(" %li '%s'\n", offset, lstr(seq1, len));
395	printf(" %li '%s'\n", offset, lstr(seq2, len));
396	printf(" %li '%s'\n", offset, compare);
397	}
398
399	inline void dump_rest(const char *seq, long len, int idx, long offset) {
400	printf(" Rest von Sequenz %i:\n", idx);
401	while (len>BUFLEN) {
402	printf(" %li '%s'\n", offset, lstr(seq, BUFLEN));
403	seq += BUFLEN;
404	len -= BUFLEN;
405	offset += BUFLEN;
406	}
407
408	fa_assert(len>0);
409	printf(" '%s'\n", lstr(seq, len));
410	}
411
412	static void dump_n_compare(const char text, const char seq1, long len1, const char *seq2, long len2) {
413	long offset = 0;
414
415	printf(" Comparing %s:\n", text);
416
417	while (len1>0 && len2>0) {
418	long done = 0;
419
420	if (len1>=BUFLEN && len2>=BUFLEN) {
421	dump_n_compare_one(seq1, seq2, done=BUFLEN, offset);
422	}
423	else {
424	long min = len1<len2 ? len1 : len2;
425	dump_n_compare_one(seq1, seq2, done=min, offset);
426	}
427
428	seq1 += done;
429	seq2 += done;
430	len1 -= done;
431	len2 -= done;
432	offset += done;
433	}
434
435	if (len1>0) dump_rest(seq1, len1, 1, offset);
436	if (len2>0) dump_rest(seq2, len2, 2, offset);
437	printf(" -------------------\n");
438	}
439
440	static void dump_n_compare(const char *text, const CompactedSubSequence& seq1, const CompactedSubSequence& seq2) {
441	dump_n_compare(text, seq1.text(), seq1.length(), seq2.text(), seq2.length());
442	}
443	#endif // TRACE_COMPRESSED_ALIGNMENT
444
445	#undef BUFLEN
446
447	inline void dumpSeq(const char *seq, long len, long pos) {
448	printf("'%s' ", lstr(seq, len));
449	printf("(Pos=%li,Len=%li)", pos, len);
450	}
451
452	#define dump() \
453	do { \
454	double sig = partSignificance(sequence().length(), slaveSequence.length(), bestLength); \
455	\
456	printf(" Score = %li (Significance=%f)\n" \
457	" Master = ", bestScore, sig); \
458	dumpSeq(bestMasterLeft.text(), bestLength, bestMasterLeft.leftOf()); \
459	printf("\n" \
460	" Slave = "); \
461	dumpSeq(bestSlaveLeft.text(), bestLength, bestSlaveLeft.leftOf()); \
462	printf("\n"); \
463	} while (0)
464
465	#endif //DEBUG
466
467
468	// ------------------------------------------------
469	// INLINE-functions used in fast_align():
470
471	inline double log3(double d) {
472	return log(d)/log(3.0);
473	}
474	inline double partSignificance(long seq1len, long seq2len, long partlen) {
475	// returns log3 of significance of the part
476	// usage: partSignificance(...) < log3(maxAllowedSignificance)
477	return log3((seq1len-partlen)*(seq2len-partlen)) - partlen;
478	}
479
480	inline ARB_ERROR bufferTooSmall() {
481	return "Cannot align - reserved buffer is to small";
482	}
483
484	inline long insertsToNextBase(AlignBuffer *alignBuffer, const SequencePosition& master) {
485	int inserts;
486	int nextBase;
487
488	if (master.rightOf()>0) {
489	nextBase = master.expdPosition();
490	}
491	else {
492	nextBase = master.sequence().expdPosition(master.sequence().length());
493	}
494	inserts = nextBase-alignBuffer->offset();
495
496	return inserts;
497	}
498
499	inline void insertBase(AlignBuffer *alignBuffer,
500	SequencePosition& master, SequencePosition& slave,
501	FastAlignReport *report)
502	{
503	char slaveBase = *slave.text();
504	char masterBase = *master.text();
505
506	alignBuffer->set(slaveBase, alignQuality(slaveBase, masterBase));
507	report->count_aligned_base(slaveBase!=masterBase);
508	++slave;
509	++master;
510	}
511
512	inline void insertSlaveBases(AlignBuffer *alignBuffer,
513	SequencePosition& slave,
514	int length,
515	FastAlignReport *report)
516	{
517	alignBuffer->copy(slave.text(), alignQuality(*slave.text(), GAP_CHAR), length);
518	report->count_unaligned_base(length);
519	slave += length;
520	}
521
522	inline void insertGap(AlignBuffer *alignBuffer,
523	SequencePosition& master,
524	FastAlignReport *report)
525	{
526	char masterBase = *master.text();
527
528	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, masterBase));
529	report->count_aligned_base(masterBase!=GAP_CHAR);
530	++master;
531	}
532
533	static ARB_ERROR insertClustalValigned(AlignBuffer *alignBuffer,
534	SequencePosition& master,
535	SequencePosition& slave,
536	const char masterAlignment, const char slaveAlignment, long alignmentLength,
537	FastAlignReport *report)
538	{
539	// inserts bases of 'slave' to 'alignBuffer' according to alignment in 'masterAlignment' and 'slaveAlignment'
540	#define ACID '*' // contents of 'masterAlignment' and 'slaveAlignment'
541	#define GAP '-'
542
543	int pos;
544	int baseAtLeft = 0; // TRUE -> last position in alignBuffer contains a base
545
546	for (pos=0; pos<alignmentLength; pos++) {
547	long insert = insertsToNextBase(alignBuffer, master); // in expanded seq
548
549	if (masterAlignment[pos]==ACID) {
550	if (insert>0) {
551	if (insert>alignBuffer->free()) return bufferTooSmall();
552	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
553	fa_assert(insertsToNextBase(alignBuffer, master)==0);
554	insert = 0;
555	}
556
557	if (!alignBuffer->free()) return bufferTooSmall();
558	if (slaveAlignment[pos]==ACID) {
559	insertBase(alignBuffer, master, slave, report);
560	baseAtLeft = 1;
561	}
562	else {
563	insertGap(alignBuffer, master, report);
564	baseAtLeft = 0;
565	}
566	}
567	else {
568	int slave_bases;
569
570	fa_assert(masterAlignment[pos]==GAP);
571	for (slave_bases=1; pos+slave_bases<alignmentLength && masterAlignment[pos+slave_bases]==GAP; slave_bases++) {
572	; // count gaps in master (= # of slave bases to insert)
573	}
574	if (!baseAtLeft && insert>slave_bases) {
575	int ins_gaps = insert-slave_bases;
576
577	fa_assert(alignBuffer->free()>=ins_gaps);
578	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), ins_gaps);
579	insert -= ins_gaps;
580	}
581
582	if (insert<slave_bases) { // master has less gaps than slave bases to insert
583	report->memorize_insertion(master.expdPosition(), slave_bases-insert);
584	}
585	else if (insert>slave_bases) { // master has more gaps than slave bases to insert
586	fa_assert(baseAtLeft);
587	}
588
589	unaligned_bases.memorize(ExplicitRange(slave.expdPosition(), // memorize base positions without counterpart in master
590	slave.expdPosition(slave_bases-1)));
591
592	if (slave_bases>alignBuffer->free()) return bufferTooSmall();
593	insertSlaveBases(alignBuffer, slave, slave_bases, report);
594	pos += slave_bases-1; // -1 compensates for()-increment
595	baseAtLeft = 1;
596	}
597	}
598
599	return 0;
600
601	#undef GAP
602	#undef ACID
603	}
604
605	static ARB_ERROR insertAligned(AlignBuffer *alignBuffer,
606	SequencePosition& master, SequencePosition& slave, long partLength,
607	FastAlignReport *report)
608	{
609	// insert bases of 'slave' to 'alignBuffer' according to 'master'
610	if (partLength) {
611	long insert = insertsToNextBase(alignBuffer, master);
612
613	fa_assert(partLength>=0);
614
615	if (insert<0) { // insert gaps into master
616	long min_insert = insert;
617
618	report->memorize_insertion(master.expdPosition(), -insert);
619
620	while (insert<0 && partLength) {
621	if (insert<min_insert) min_insert = insert;
622	if (!alignBuffer->free()) {
623	return bufferTooSmall();
624	}
625	insertBase(alignBuffer, master, slave, report);
626	partLength--;
627	insert = insertsToNextBase(alignBuffer, master);
628	}
629
630	fa_assert(partLength>=0);
631	if (partLength==0) { // all inserted
632	return NULL;
633	}
634	}
635
636	fa_assert(insert>=0);
637
638	if (insert>0) { // insert gaps into slave
639	if (insert>alignBuffer->free()) return bufferTooSmall();
640	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
641	fa_assert(insertsToNextBase(alignBuffer, master)==0);
642	}
643
644	fa_assert(partLength>=0);
645
646	while (partLength--) {
647	insert = insertsToNextBase(alignBuffer, master);
648
649	fa_assert(insert>=0);
650	if (insert>0) {
651	if (insert>=alignBuffer->free()) return bufferTooSmall();
652	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
653	}
654	else {
655	if (!alignBuffer->free()) {
656	return bufferTooSmall();
657	}
658	}
659
660	insertBase(alignBuffer, master, slave, report);
661	}
662	}
663
664	return NULL;
665	}
666
667	static ARB_ERROR cannot_fast_align(const CompactedSubSequence& master, long moffset, long mlength,
668	const CompactedSubSequence& slaveSequence, long soffset, long slength,
669	int max_seq_length,
670	AlignBuffer *alignBuffer,
671	FastAlignReport *report)
672	{
673	const char *mtext = master.text(moffset);
674	const char *stext = slaveSequence.text(soffset);
675	ARB_ERROR error = 0;
676
677	if (slength) {
678	if (mlength) { // if slave- and master-sequences contain bases, we call the slow aligner
679	#ifdef TRACE_CLUSTAL_DATA
680	printf("ClustalV-Align:\n");
681	printf(" mseq = '%s'\n", lstr(mtext, mlength));
682	printf(" sseq = '%s'\n", lstr(stext, slength));
683	#endif // TRACE_CLUSTAL_DATA
684
685	int is_dna = -1;
686
687	switch (global_alignmentType) {
688	case GB_AT_RNA:
689	case GB_AT_DNA: is_dna = 1; break;
690	case GB_AT_AA: is_dna = 0; break;
691	default: error = "Unknown alignment type - aligner aborted"; break;
692	}
693
694	const char maligned, saligned;
695	int len;
696	if (!error) {
697	int score; // unused
698	error = ClustalV_align(is_dna,
699	1,
700	mtext, mlength, stext, slength,
701	master.gapsBefore(moffset),
702	max_seq_length,
703	maligned, saligned, len, score);
704	}
705
706	if (!error) {
707	#ifdef TRACE_CLUSTAL_DATA
708	printf("ClustalV returns:\n");
709	printf(" maligned = '%s'\n", lstr(maligned, len));
710	printf(" saligned = '%s'\n", lstr(saligned, len));
711	#endif // TRACE_CLUSTAL_DATA
712
713	SequencePosition masterPos(master, moffset);
714	SequencePosition slavePos(slaveSequence, soffset);
715
716	error = insertClustalValigned(alignBuffer, masterPos, slavePos, maligned, saligned, len, report);
717
718	#if (defined(DEBUG) && 0)
719
720	SequencePosition master2(master->sequence(), moffset);
721	SequencePosition slave2(slaveSequence, soffset);
722	char *cmp = new char[len];
723
724	for (int l=0; l<len; l++) {
725	int gaps = 0;
726
727	if (maligned[l]=='*') {
728	maligned[l] = *master2.text();
729	++master2;
730	}
731	else {
732	gaps++;
733	}
734
735	if (saligned[l]=='*') {
736	saligned[l] = *slave2.text();
737	++slave2;
738	}
739	else {
740	gaps++;
741	}
742
743	cmp[l] = gaps \|\| maligned[l]==saligned[l] ? '=' : 'X';
744	}
745
746	printf(" master = '%s'\n", lstr(maligned, len));
747	printf(" slave = '%s'\n", lstr(saligned, len));
748	printf(" '%s'\n", lstr(cmp, len));
749
750	delete [] cmp;
751	#endif
752	}
753	}
754	else { // master is empty here, so we just copy in the slave bases
755	if (slength<=alignBuffer->free()) {
756	unaligned_bases.memorize(ExplicitRange(slaveSequence.expdPosition(soffset),
757	slaveSequence.expdPosition(soffset+slength-1)));
758
759	alignBuffer->copy(slaveSequence.text(soffset), '?', slength); // '?' means not aligned (just inserted)
760	// corrected by at alignBuffer->correctUnalignedPositions()
761	report->count_unaligned_base(slength);
762	}
763	else {
764	error = bufferTooSmall();
765	}
766	}
767	}
768
769	return error;
770	}
771
772	// ------------------------------------
773	// #define's for fast_align()
774
775	#define TEST_BETTER_SCORE() \
776	do { \
777	if (score>bestScore) { \
778	bestScore = score; \
779	bestLength = masterRight.text() - masterLeft.text(); \
780	bestMasterLeft = masterLeft; \
781	bestSlaveLeft = slaveLeft; \
782	} \
783	} while (0)
784
785	#define CAN_SCORE_LEFT() (masterLeft.leftOf() && slaveLeft.leftOf())
786	#define CAN_SCORE_RIGHT() (masterRight.rightOf() && slaveRight.rightOf())
787
788	#define SCORE_LEFT() \
789	do { \
790	score += (--masterLeft).text()==(--slaveLeft).text() ? match : mismatch; \
791	TEST_BETTER_SCORE(); \
792	} while (0)
793
794	#define SCORE_RIGHT() \
795	do { \
796	score += (++masterRight).text()==(++slaveRight).text() ? match : mismatch; \
797	TEST_BETTER_SCORE(); \
798	} while (0)
799
800
801	ARB_ERROR FastSearchSequence::fast_align(const CompactedSubSequence& slaveSequence,
802	AlignBuffer *alignBuffer,
803	int max_seq_length,
804	int match, int mismatch,
805	FastAlignReport *report) const
806	{
807	// aligns 'slaveSequence' to 'this'
808	//
809	// returns
810	// == NULL -> all ok -> 'alignBuffer' contains aligned sequence
811	// != NULL -> failure -> no results
812
813	ARB_ERROR error = NULL;
814	int aligned = 0;
815
816	// set the following #if to zero to test ClustalV-Aligner (not fast_aligner)
817	#if 1
818
819	static double lowSignificance;
820	static int lowSignificanceInitialized;
821
822	if (!lowSignificanceInitialized) {
823	lowSignificance = log3(0.01);
824	lowSignificanceInitialized = 1;
825	}
826
827	SequencePosition slave(slaveSequence);
828	long bestScore=0;
829	SequencePosition bestMasterLeft(sequence());
830	SequencePosition bestSlaveLeft(slaveSequence);
831	long bestLength=0;
832
833	while (slave.rightOf()>=3) { // with all triples of slaveSequence
834	FastSearchOccurrence occurrence(*this, slave.text()); // "search" first
835	SequencePosition rightmostSlave = slave;
836
837	while (occurrence.found()) { // with all occurrences of the triple
838	long score = match*3;
839	SequencePosition masterLeft(occurrence.sequence(), occurrence.offset());
840	SequencePosition masterRight(occurrence.sequence(), occurrence.offset()+3);
841	SequencePosition slaveLeft(slave);
842	SequencePosition slaveRight(slave, 3);
843
844	while (score>0) {
845	if (CAN_SCORE_LEFT()) {
846	SCORE_LEFT();
847	}
848	else {
849	while (score>0 && CAN_SCORE_RIGHT()) {
850	SCORE_RIGHT();
851	}
852	break;
853	}
854
855	if (CAN_SCORE_RIGHT()) {
856	SCORE_RIGHT();
857	}
858	else {
859	while (score>0 && CAN_SCORE_LEFT()) {
860	SCORE_LEFT();
861	}
862	break;
863	}
864	}
865
866	occurrence.gotoNext(); // "search" next
867
868	if (rightmostSlave<slaveRight) {
869	rightmostSlave = slaveRight;
870	rightmostSlave -= 3;
871	}
872	}
873
874	if (rightmostSlave>slave) slave = rightmostSlave;
875	else ++slave;
876	}
877
878	if (bestLength) {
879	double sig = partSignificance(sequence().length(), slaveSequence.length(), bestLength);
880
881	if (sig<lowSignificance) {
882	long masterLeftOf = bestMasterLeft.leftOf();
883	long masterRightStart = masterLeftOf+bestLength;
884	long masterRightOf = bestMasterLeft.rightOf()-bestLength;
885	long slaveLeftOf = bestSlaveLeft.leftOf();
886	long slaveRightStart = slaveLeftOf+bestLength;
887	long slaveRightOf = bestSlaveLeft.rightOf()-bestLength;
888
889	#define MIN_ALIGNMENT_RANGE 4
890
891	if (!error) {
892	if (masterLeftOf >= MIN_ALIGNMENT_RANGE && slaveLeftOf >= MIN_ALIGNMENT_RANGE) {
893	CompactedSubSequence leftCompactedMaster(sequence(), 0, masterLeftOf);
894	FastSearchSequence leftMaster(leftCompactedMaster);
895
896	error = leftMaster.fast_align(CompactedSubSequence(slaveSequence, 0, slaveLeftOf),
897	alignBuffer, max_seq_length, match, mismatch, report);
898	}
899	else if (slaveLeftOf>0) {
900	error = cannot_fast_align(sequence(), 0, masterLeftOf,
901	slaveSequence, 0, slaveLeftOf,
902	max_seq_length, alignBuffer, report);
903	}
904
905	aligned = 1;
906	}
907
908	// align part of slave sequence according to master sequence:
909
910	if (!error) {
911	#if (defined(DEBUG) && 0)
912	long offset = alignBuffer->offset();
913	long used;
914	#endif
915	error = insertAligned(alignBuffer, bestMasterLeft, bestSlaveLeft, bestLength, report);
916	#if (defined(DEBUG) && 0)
917	used = alignBuffer->offset()-offset;
918	printf("aligned '%s' (len=%li, address=%li)\n", lstr(alignBuffer->text()+offset, used), used, long(alignBuffer));
919	#endif
920	aligned = 1;
921	}
922
923	if (!error) {
924	if (masterRightOf >= MIN_ALIGNMENT_RANGE && slaveRightOf >= MIN_ALIGNMENT_RANGE) {
925	CompactedSubSequence rightCompactedMaster(sequence(), masterRightStart, masterRightOf);
926	FastSearchSequence rightMaster(rightCompactedMaster);
927
928	error = rightMaster.fast_align(CompactedSubSequence(slaveSequence, slaveRightStart, slaveRightOf),
929	alignBuffer,
930	max_seq_length, match, mismatch, report);
931	}
932	else if (slaveRightOf>0) {
933	error = cannot_fast_align(sequence(), masterRightStart, masterRightOf,
934	slaveSequence, slaveRightStart, slaveRightOf,
935	max_seq_length, alignBuffer, report);
936	}
937
938	aligned = 1;
939	}
940
941	}
942	}
943
944	#endif
945
946	if (!aligned && !error) {
947	error = cannot_fast_align(sequence(), 0, sequence().length(),
948	slaveSequence, 0, slaveSequence.length(),
949	max_seq_length, alignBuffer, report);
950	}
951
952	return error;
953	}
954
955	#undef dump
956	#undef TEST_BETTER_SCORE
957	#undef CAN_SCORE_LEFT
958	#undef CAN_SCORE_RIGHT
959	#undef SCORE_LEFT
960	#undef SCORE_RIGHT
961
962	inline long calcSequenceChecksum(const char *data, long length) {
963	return GB_checksum(data, length, 1, GAP_CHARS);
964	}
965
966	#if defined(WARN_TODO)
967	#warning firstColumn + lastColumn -> PosRange
968	#endif
969
970	static CompactedSubSequence readCompactedSequence(GBDATA gb_species,
971	const char *ali,
972	ARB_ERROR *errorPtr,
973	char **dataPtr, // if dataPtr != NULL it will be set to the WHOLE uncompacted sequence
974	long *seqChksum, // may be NULL (of part of sequence)
975	PosRange range) // if !range.is_whole() -> return only part of the sequence
976	{
977	ARB_ERROR error = 0;
978	GBDATA *gbd;
979	CompactedSubSequence *seq = 0;
980
981
982
983	gbd = GBT_find_sequence(gb_species, ali); // get sequence
984
985	if (gbd) {
986	long length = GB_read_string_count(gbd);
987	char *data = GB_read_string(gbd);
988	long partLength;
989	char *partData;
990
991	if (dataPtr) { // make a copy of the whole uncompacted sequence (returned to caller)
992	*dataPtr = data;
993	}
994
995	int firstColumn = range.start();
996	if (range.is_part()) { // take only part of sequence
997	int lastColumn = range.end();
998
999	fa_assert(firstColumn>=0);
1000	fa_assert(firstColumn<=lastColumn);
1001
1002	// include all surrounding gaps (@@@ this might cause problems with partial alignment)
1003	while (firstColumn>0 && is_ali_gap(data[firstColumn-1])) {
1004	firstColumn--;
1005	}
1006	if (lastColumn!=-1) {
1007	while (lastColumn<(length-1) && is_ali_gap(data[lastColumn+1])) lastColumn++;
1008	fa_assert(lastColumn<length);
1009	}
1010
1011	partData = data+firstColumn;
1012	int slen = length-firstColumn;
1013
1014	fa_assert(slen>=0);
1015	fa_assert((size_t)slen==strlen(partData));
1016
1017	if (lastColumn==-1) { // take all till end of sequence
1018	partLength = slen;
1019	}
1020	else {
1021	partLength = lastColumn-firstColumn+1;
1022	if (partLength>slen) partLength = slen; // cut rest, if we have any
1023	}
1024	}
1025	else {
1026	partLength = length;
1027	partData = data;
1028	}
1029
1030	if (!error) {
1031	if (seqChksum) {
1032	*seqChksum = calcSequenceChecksum(partData, partLength);
1033	}
1034
1035	seq = new CompactedSubSequence(partData, partLength, GBT_read_name(gb_species), firstColumn);
1036	}
1037
1038	if (!dataPtr) { // free sequence only if user has not requested to get it
1039	free(data);
1040	}
1041	}
1042	else {
1043	error = GBS_global_string("No 'data' found for species '%s'", GBT_read_name(gb_species));
1044	if (dataPtr) *dataPtr = NULL; // (user must not care to free data if we fail)
1045	}
1046
1047	fa_assert(errorPtr);
1048	*errorPtr = error;
1049
1050	return seq;
1051	}
1052
1053	static ARB_ERROR writeStringToAlignment(GBDATA *gb_species, GB_CSTR alignment, GB_CSTR data_name, GB_CSTR str, bool temporary) {
1054	GBDATA *gb_ali = GB_search(gb_species, alignment, GB_DB);
1055	ARB_ERROR error = NULL;
1056	GBDATA *gb_name = GB_search(gb_ali, data_name, GB_STRING);
1057
1058	if (gb_name) {
1059	fa_assert(GB_check_father(gb_name, gb_ali));
1060	error = GB_write_string(gb_name, str);
1061	if (temporary && !error) error = GB_set_temporary(gb_name);
1062	}
1063	else {
1064	error = GBS_global_string("Cannot create entry '%s' for '%s'", data_name, GBT_read_name(gb_species));
1065	}
1066
1067	return error;
1068	}
1069
1070	// --------------------------------------------------------------------------------
1071
1072	static ARB_ERROR alignCompactedTo(CompactedSubSequence *toAlignSequence,
1073	const FastSearchSequence *alignTo,
1074	int max_seq_length,
1075	GB_CSTR alignment,
1076	long toAlignChksum,
1077	GBDATA *gb_toAlign,
1078	GBDATA *gb_alignTo, // may be NULL
1079	const AlignParams& ali_params)
1080	{
1081	// if only part of the sequence should be aligned, then this functions already gets only the part
1082	// (i.o.w.: toAlignSequence, alignTo and toAlignChksum refer to the partial sequence)
1083	AlignBuffer alignBuffer(max_seq_length);
1084	if (ali_params.range.start()>0) {
1085	alignBuffer.set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), ali_params.range.start());
1086	}
1087
1088	const char *master_name = read_name(gb_alignTo);
1089
1090	FastAlignReport report(master_name, ali_params.showGapsMessages);
1091
1092	{
1093	static GBDATA *last_gb_toAlign = 0;
1094	if (gb_toAlign!=last_gb_toAlign) {
1095	last_gb_toAlign = gb_toAlign;
1096	currentSequenceNumber++;
1097	}
1098	}
1099
1100	#ifdef TRACE_COMPRESSED_ALIGNMENT
1101	printf("alignCompactedTo(): master='%s' ", master_name);
1102	printf("slave='%s'\n", toAlignSequence->name());
1103	#endif // TRACE_COMPRESSED_ALIGNMENT
1104
1105	ARB_ERROR error = GB_pop_transaction(gb_toAlign);
1106	if (!error) {
1107	if (island_hopper) {
1108	error = island_hopper->do_align();
1109	if (!error) {
1110	fa_assert(island_hopper->was_aligned());
1111
1112	#ifdef TRACE_ISLANDHOPPER_DATA
1113	printf("Island-Hopper returns:\n");
1114	printf("maligned = '%s'\n", lstr(island_hopper->get_result_ref(), island_hopper->get_result_length()));
1115	printf("saligned = '%s'\n", lstr(island_hopper->get_result(), island_hopper->get_result_length()));
1116	#endif // TRACE_ISLANDHOPPER_DATA
1117
1118	SequencePosition masterPos(alignTo->sequence(), 0);
1119	SequencePosition slavePos(*toAlignSequence, 0);
1120
1121	error = insertClustalValigned(&alignBuffer, masterPos, slavePos, island_hopper->get_result_ref(), island_hopper->get_result(), island_hopper->get_result_length(), &report);
1122	}
1123	}
1124	else {
1125	error = alignTo->fast_align(*toAlignSequence, &alignBuffer, max_seq_length, 2, -10, &report); // <- here we align
1126	}
1127	}
1128
1129	if (!error) {
1130	alignBuffer.correctUnalignedPositions();
1131	if (alignBuffer.free()) {
1132	alignBuffer.set('-', alignQuality(GAP_CHAR, GAP_CHAR), alignBuffer.free()); // rest of alignBuffer is set to '-'
1133	}
1134	alignBuffer.restoreDots(*toAlignSequence);
1135	}
1136
1137	#ifdef TRACE_COMPRESSED_ALIGNMENT
1138	if (!error) {
1139	CompactedSubSequence alignedSlave(alignBuffer.text(), alignBuffer.length(), toAlignSequence->name());
1140	dump_n_compare("reference vs. aligned:", alignTo->sequence(), alignedSlave);
1141	}
1142	#endif // TRACE_COMPRESSED_ALIGNMENT
1143
1144	{
1145	GB_ERROR err = GB_push_transaction(gb_toAlign);
1146	if (!error) error = err;
1147	}
1148
1149	if (!error) {
1150	if (calcSequenceChecksum(alignBuffer.text(), alignBuffer.length())!=toAlignChksum) { // sequence-chksum changed
1151	error = "Internal aligner error (sequence checksum changed) -- aborted";
1152
1153	#ifdef TRACE_COMPRESSED_ALIGNMENT
1154	CompactedSubSequence alignedSlave(alignBuffer.text(), alignBuffer.length(), toAlignSequence->name());
1155	dump_n_compare("Old Slave vs. new Slave", *toAlignSequence, alignedSlave);
1156	#endif // TRACE_COMPRESSED_ALIGNMENT
1157	}
1158	else {
1159	GB_push_my_security(gb_toAlign);
1160	{
1161	GBDATA *gbd = GBT_add_data(gb_toAlign, alignment, "data", GB_STRING);
1162
1163	if (!gbd) {
1164	error = "Can't find/create sequence data";
1165	}
1166	else {
1167	if (ali_params.range.is_part()) { // we aligned just a part of the sequence
1168	char *buffer = GB_read_string(gbd); // so we read old sequence data
1169	long len = GB_read_string_count(gbd);
1170	if (!buffer) error = GB_await_error();
1171	else {
1172	int lenToCopy = ali_params.range.size();
1173	long wholeChksum = calcSequenceChecksum(buffer, len);
1174
1175	memcpy(buffer+ali_params.range.start(), alignBuffer.text()+ali_params.range.start(), lenToCopy); // copy in the aligned part
1176	// @@@ genau um 1 position zu niedrig
1177
1178	if (calcSequenceChecksum(buffer, len) != wholeChksum) {
1179	error = "Internal aligner error (sequence checksum changed) -- aborted";
1180	# ifdef TRACE_COMPRESSED_ALIGNMENT
1181	char *buffer_org = GB_read_string(gbd);
1182	dump_n_compare("Old seq vs. new seq (slave)", buffer_org, len, buffer, len);
1183	free(buffer_org);
1184	# endif // TRACE_COMPRESSED_ALIGNMENT
1185	}
1186	else {
1187	GB_write_string(gbd, "");
1188	error = GB_write_string(gbd, buffer);
1189	}
1190	}
1191
1192	free(buffer);
1193	}
1194	else {
1195	alignBuffer.setDotsAtEOSequence();
1196	error = GB_write_string(gbd, alignBuffer.text()); // aligned all -> write all
1197	}
1198	}
1199	}
1200	GB_pop_my_security(gb_toAlign);
1201
1202	if (!error && ali_params.report != FA_NO_REPORT) {
1203	// create temp-entry for slave containing alignment quality:
1204
1205	error = writeStringToAlignment(gb_toAlign, alignment, QUALITY_NAME, alignBuffer.quality(), ali_params.report==FA_TEMP_REPORT);
1206	if (!error) {
1207	// create temp-entry for master containing insert dots:
1208
1209	int buflen = max_seq_length*2;
1210	char buffer = (char)malloc(buflen+1);
1211	char *afterLast = buffer;
1212
1213	if (!buffer) {
1214	error = "out of memory";
1215	}
1216	else {
1217	memset(buffer, '-', buflen);
1218	buffer[buflen] = 0;
1219
1220	const FastAlignInsertion *inserts = report.insertion();
1221	while (inserts) {
1222	memset(buffer+inserts->offset(), '>', inserts->gaps());
1223	afterLast = buffer+inserts->offset()+inserts->gaps();
1224	inserts = inserts->next();
1225	}
1226
1227	*afterLast = 0;
1228	if (gb_alignTo) {
1229	error = writeStringToAlignment(gb_alignTo, alignment, INSERTS_NAME, buffer, ali_params.report==FA_TEMP_REPORT);
1230	}
1231	}
1232	}
1233	}
1234	}
1235	}
1236	return error;
1237	}
1238
1239	ARB_ERROR FastAligner_delete_temp_entries(GBDATA gb_species, const char alignment) {
1240	fa_assert(gb_species);
1241	fa_assert(alignment);
1242
1243	GBDATA *gb_ali = GB_search(gb_species, alignment, GB_FIND);
1244	ARB_ERROR error = NULL;
1245
1246	if (gb_ali) {
1247	GBDATA *gb_name = GB_search(gb_ali, QUALITY_NAME, GB_FIND);
1248	if (gb_name) {
1249	error = GB_delete(gb_name);
1250	#if defined(DEBUG)
1251	printf("deleted '%s' of '%s' (gb_name=%p)\n", QUALITY_NAME, GBT_read_name(gb_species), gb_name);
1252	#endif
1253	}
1254
1255	if (!error) {
1256	gb_name = GB_search(gb_ali, INSERTS_NAME, GB_FIND);
1257	if (gb_name) {
1258	error = GB_delete(gb_name);
1259	#if defined(DEBUG)
1260	printf("deleted '%s' of '%s' (gb_name=%p)\n", INSERTS_NAME, GBT_read_name(gb_species), gb_name);
1261	#endif
1262	}
1263	}
1264	}
1265
1266	return error;
1267	}
1268
1269	static ARB_ERROR align_error(ARB_ERROR olderr, GBDATA gb_toAlign, GBDATA gb_alignTo) {
1270	// used by alignTo() and alignToNextRelative() to transform errors coming from subroutines
1271	// can be used by higher functions
1272
1273	const char *name_toAlign = read_name(gb_toAlign);
1274	const char *name_alignTo = read_name(gb_alignTo);
1275
1276	fa_assert(olderr);
1277
1278	return GBS_global_string("Error while aligning '%s' to '%s':\n%s",
1279	name_toAlign, name_alignTo, olderr.deliver());
1280	}
1281
1282	static ARB_ERROR alignTo(GBDATA *gb_toAlign,
1283	GB_CSTR alignment,
1284	const FastSearchSequence *alignTo,
1285	GBDATA *gb_alignTo, // may be NULL
1286	int max_seq_length,
1287	const AlignParams& ali_params)
1288	{
1289	ARB_ERROR error = NULL;
1290	long chksum;
1291
1292	CompactedSubSequence *toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, NULL, &chksum, ali_params.range);
1293
1294	if (island_hopper) {
1295	GBDATA *gb_seq = GBT_find_sequence(gb_toAlign, alignment); // get sequence
1296	if (gb_seq) {
1297	long length = GB_read_string_count(gb_seq);
1298	const char *data = GB_read_char_pntr(gb_seq);
1299
1300	island_hopper->set_toAlign_sequence(data);
1301	island_hopper->set_alignment_length(length);
1302	}
1303	}
1304
1305
1306
1307	if (!error) {
1308	error = alignCompactedTo(toAlignSequence, alignTo, max_seq_length, alignment, chksum, gb_toAlign, gb_alignTo, ali_params);
1309	if (error) error = align_error(error, gb_toAlign, gb_alignTo);
1310	delete toAlignSequence;
1311	}
1312
1313	return error;
1314	}
1315
1316	static ARB_ERROR alignToGroupConsensus(GBDATA *gb_toAlign,
1317	GB_CSTR alignment,
1318	Aligner_get_consensus_func get_consensus,
1319	int max_seq_length,
1320	const AlignParams& ali_params)
1321	{
1322	ARB_ERROR error = 0;
1323	char *consensus = get_consensus(read_name(gb_toAlign), ali_params.range);
1324	size_t cons_len = strlen(consensus);
1325
1326	fa_assert(cons_len);
1327
1328	for (size_t i = 0; i<cons_len; ++i) { // translate consensus to be accepted by aligner
1329	switch (consensus[i]) {
1330	case '=': consensus[i] = '-'; break;
1331	default: break;
1332	}
1333	}
1334
1335	CompactedSubSequence compacted(consensus, cons_len, "group consensus", ali_params.range.start());
1336
1337	{
1338	FastSearchSequence fast(compacted);
1339	error = alignTo(gb_toAlign, alignment, &fast, NULL, max_seq_length, ali_params);
1340	}
1341
1342	free(consensus);
1343
1344	return error;
1345	}
1346
1347	static void appendNameAndUsedBasePositions(char *toString, GBDATA gb_species, int usedBasePositions) {
1348	// if usedBasePositions == -1 -> unknown;
1349
1350	char *currInfo;
1351	if (usedBasePositions<0) {
1352	currInfo = strdup(GBT_read_name(gb_species));
1353	}
1354	else {
1355	fa_assert(usedBasePositions>0); // otherwise it should NOT be announced here!
1356	currInfo = GBS_global_string_copy("%s:%i", GBT_read_name(gb_species), usedBasePositions);
1357	}
1358
1359	char *newString = 0;
1360	if (*toString) {
1361	newString = GBS_global_string_copy("%s, %s", *toString, currInfo);
1362	}
1363	else {
1364	newString = currInfo;
1365	currInfo = 0; // don't free
1366	}
1367
1368	freeset(*toString, newString);
1369	free(currInfo);
1370	}
1371
1372	inline int min(int i, int j) { return i<j ? i : j; }
1373
1374	static ARB_ERROR alignToNextRelative(SearchRelativeParams& relSearch,
1375	int max_seq_length,
1376	FA_turn turnAllowed,
1377	GB_CSTR alignment,
1378	GBDATA *gb_toAlign,
1379	const AlignParams& ali_params)
1380	{
1381	CompactedSubSequence *toAlignSequence = NULL;
1382	bool use_different_pt_server_alignment = 0 != strcmp(relSearch.pt_server_alignment, alignment);
1383
1384	ARB_ERROR error;
1385	long chksum;
1386	int relativesToTest = relSearch.maxRelatives*2; // get more relatives from pt-server (needed when use_different_pt_server_alignment == true)
1387	char *nearestRelative = new char[relativesToTest+1];
1388	int next_relatives;
1389	int i;
1390	GBDATA *gb_main = GB_get_root(gb_toAlign);
1391
1392	if (use_different_pt_server_alignment) {
1393	turnAllowed = FA_TURN_NEVER; // makes no sense if we're using a different alignment for the pt_server
1394	}
1395
1396	for (next_relatives=0; next_relatives<relativesToTest; next_relatives++) {
1397	nearestRelative[next_relatives] = 0;
1398	}
1399	next_relatives = 0;
1400
1401	bool restart = true;
1402	while (restart) {
1403	restart = false;
1404
1405	char *findRelsBySeq = 0;
1406	if (use_different_pt_server_alignment) {
1407	toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, 0, &chksum, ali_params.range);
1408
1409	GBDATA *gbd = GBT_find_sequence(gb_toAlign, relSearch.pt_server_alignment); // use a different alignment for next relative search
1410	if (!gbd) {
1411	error = GBS_global_string("Species '%s' has no data in alignment '%s'", GBT_read_name(gb_toAlign), relSearch.pt_server_alignment);
1412	}
1413	else {
1414	findRelsBySeq = GB_read_string(gbd);
1415	}
1416	}
1417	else {
1418	toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, &findRelsBySeq, &chksum, ali_params.range);
1419	}
1420
1421	if (error) {
1422	delete toAlignSequence;
1423	return error; // @@@ leaks ?
1424	}
1425
1426	while (next_relatives) {
1427	next_relatives--;
1428	freenull(nearestRelative[next_relatives]);
1429	}
1430
1431	{
1432	// find relatives
1433	FamilyFinder *familyFinder = relSearch.getFamilyFinder();
1434	const PosRange& range = familyFinder->get_TargetRange();
1435
1436	if (range.is_part()) {
1437	range.copy_corresponding_part(findRelsBySeq, findRelsBySeq, strlen(findRelsBySeq));
1438	turnAllowed = FA_TURN_NEVER; // makes no sense if we're using partial relative search
1439	}
1440
1441	error = familyFinder->searchFamily(findRelsBySeq, FF_FORWARD, relativesToTest+1, 0); // @@@ make min_score configurable
1442
1443	// @@@ case where no relative found (due to min score) handle how ? abort ? warn ?
1444
1445	double bestScore = 0;
1446	if (!error) {
1447	#if defined(DEBUG)
1448	double lastScore = -1;
1449	#if defined(TRACE_RELATIVES)
1450	fprintf(stderr, "List of relatives used for '%s':\n", GBT_read_name(gb_toAlign));
1451	#endif // TRACE_RELATIVES
1452	#endif // DEBUG
1453	for (const FamilyList *fl = familyFinder->getFamilyList(); fl; fl=fl->next) {
1454	if (strcmp(toAlignSequence->name(), fl->name)!=0) {
1455	if (GBT_find_species(gb_main, fl->name)) { // @@@
1456	double thisScore = familyFinder->uses_rel_matches() ? fl->rel_matches : fl->matches;
1457	#if defined(DEBUG)
1458	// check whether family list is sorted correctly
1459	fa_assert(lastScore < 0 \|\| lastScore >= thisScore);
1460	lastScore = thisScore;
1461	#if defined(TRACE_RELATIVES)
1462	fprintf(stderr, "- %s (%5.2f)\n", fl->name, thisScore);
1463	#endif // TRACE_RELATIVES
1464	#endif // DEBUG
1465	if (thisScore>=bestScore) bestScore = thisScore;
1466	if (next_relatives<(relativesToTest+1)) {
1467	nearestRelative[next_relatives] = strdup(fl->name);
1468	next_relatives++;
1469	}
1470	}
1471	}
1472	}
1473	}
1474
1475	if (!error && turnAllowed != FA_TURN_NEVER) { // test if mirrored sequence has better relatives
1476	char *mirroredSequence = strdup(findRelsBySeq);
1477	long length = strlen(mirroredSequence);
1478	double bestMirroredScore = 0;
1479
1480	char T_or_U;
1481	error = GBT_determine_T_or_U(global_alignmentType, &T_or_U, "reverse-complement");
1482	if (!error) {
1483	GBT_reverseComplementNucSequence(mirroredSequence, length, T_or_U);
1484	error = familyFinder->searchFamily(mirroredSequence, FF_FORWARD, relativesToTest+1, 0); // @@@ make min_score configurable
1485	}
1486	if (!error) {
1487	#if defined(DEBUG)
1488	double lastScore = -1;
1489	#if defined(TRACE_RELATIVES)
1490	fprintf(stderr, "List of relatives used for '%s' (turned seq):\n", GBT_read_name(gb_toAlign));
1491	#endif // TRACE_RELATIVES
1492	#endif // DEBUG
1493	for (const FamilyList *fl = familyFinder->getFamilyList(); fl; fl = fl->next) {
1494	double thisScore = familyFinder->uses_rel_matches() ? fl->rel_matches : fl->matches;
1495	#if defined(DEBUG)
1496	// check whether family list is sorted correctly
1497	fa_assert(lastScore < 0 \|\| lastScore >= thisScore);
1498	lastScore = thisScore;
1499	#if defined(TRACE_RELATIVES)
1500	fprintf(stderr, "- %s (%5.2f)\n", fl->name, thisScore);
1501	#endif // TRACE_RELATIVES
1502	#endif // DEBUG
1503	if (thisScore >= bestMirroredScore) {
1504	if (strcmp(toAlignSequence->name(), fl->name)!=0) {
1505	if (GBT_find_species(gb_main, fl->name)) bestMirroredScore = thisScore; // @@@
1506	}
1507	}
1508	}
1509	}
1510
1511	int turnIt = 0;
1512	if (bestMirroredScore>bestScore) {
1513	if (turnAllowed==FA_TURN_INTERACTIVE) {
1514	const char *message;
1515	if (familyFinder->uses_rel_matches()) {
1516	message = GBS_global_string("'%s' seems to be the other way round (score: %.1f%%, score if turned: %.1f%%)",
1517	toAlignSequence->name(), bestScore100, bestMirroredScore100);
1518	}
1519	else {
1520	message = GBS_global_string("'%s' seems to be the other way round (score: %li, score if turned: %li)",
1521	toAlignSequence->name(), long(bestScore+.5), long(bestMirroredScore+.5));
1522	}
1523	turnIt = aw_question("fastali_turn_sequence", message, "Turn sequence,Leave sequence alone")==0;
1524	}
1525	else {
1526	fa_assert(turnAllowed == FA_TURN_ALWAYS);
1527	turnIt = 1;
1528	#if defined(TRACE_RELATIVES)
1529	fprintf(stderr, "Using turned sequence!\n");
1530	#endif // TRACE_RELATIVES
1531	}
1532	}
1533
1534	if (turnIt) { // write mirrored sequence
1535	GBDATA *gbd = GBT_find_sequence(gb_toAlign, alignment);
1536	GB_push_my_security(gbd);
1537	error = GB_write_string(gbd, mirroredSequence);
1538	GB_pop_my_security(gbd);
1539	if (!error) {
1540	delete toAlignSequence;
1541	restart = true; // continue while loop
1542	}
1543	}
1544
1545	free(mirroredSequence);
1546	}
1547	}
1548	free(findRelsBySeq);
1549	}
1550
1551	if (!error) {
1552	if (!next_relatives) {
1553	char warning[200];
1554	sprintf(warning, "No relative found for '%s'", toAlignSequence->name());
1555	aw_message(warning);
1556	}
1557	else {
1558	// assuming relatives are sorted! (nearest to farthest)
1559
1560	// get data pointers
1561	typedef GBDATA *GBDATAP;
1562	GBDATAP *gb_reference = new GBDATAP[relSearch.maxRelatives];
1563
1564	{
1565	for (i=0; i<relSearch.maxRelatives && i<next_relatives; i++) {
1566	GBDATA *gb_species = GBT_find_species(gb_main, nearestRelative[i]);
1567	if (!gb_species) { // pt-server seems not up to date!
1568	error = species_not_found(nearestRelative[i]);
1569	break;
1570	}
1571
1572	GBDATA *gb_sequence = GBT_find_sequence(gb_species, alignment);
1573	if (gb_sequence) { // if relative has sequence data in the current alignment ..
1574	gb_reference[i] = gb_species;
1575	}
1576	else { // remove this relative
1577	free(nearestRelative[i]);
1578	for (int j = i+1; j<next_relatives; ++j) {
1579	nearestRelative[j-1] = nearestRelative[j];
1580	}
1581	next_relatives--;
1582	nearestRelative[next_relatives] = 0;
1583	i--; // redo same index
1584	}
1585	}
1586
1587	// delete superfluous relatives
1588	for (; i<next_relatives; ++i) freenull(nearestRelative[i]);
1589
1590	if (next_relatives>relSearch.maxRelatives) {
1591	next_relatives = relSearch.maxRelatives;
1592	}
1593	}
1594
1595	// align
1596
1597	if (!error) {
1598	CompactedSubSequence *alignToSequence = readCompactedSequence(gb_reference[0], alignment, &error, NULL, NULL, ali_params.range);
1599	fa_assert(alignToSequence != 0);
1600
1601	if (island_hopper) {
1602	GBDATA *gb_ref = GBT_find_sequence(gb_reference[0], alignment); // get reference sequence
1603	GBDATA *gb_align = GBT_find_sequence(gb_toAlign, alignment); // get sequence to align
1604
1605	if (gb_ref && gb_align) {
1606	long length_ref = GB_read_string_count(gb_ref);
1607	const char *data;
1608
1609	data = GB_read_char_pntr(gb_ref);
1610	island_hopper->set_ref_sequence(data);
1611
1612	data = GB_read_char_pntr(gb_align);
1613	island_hopper->set_toAlign_sequence(data);
1614
1615	island_hopper->set_alignment_length(length_ref);
1616	}
1617	}
1618
1619	{
1620	FastSearchSequence referenceFastSeq(*alignToSequence);
1621
1622	error = alignCompactedTo(toAlignSequence, &referenceFastSeq,
1623	max_seq_length, alignment, chksum,
1624	gb_toAlign, gb_reference[0], ali_params);
1625	}
1626
1627	if (error) {
1628	error = align_error(error, gb_toAlign, gb_reference[0]);
1629	}
1630	else {
1631	char *used_relatives = 0;
1632
1633	if (unaligned_bases.is_empty()) {
1634	appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], -1);
1635	}
1636	else {
1637	if (island_hopper) {
1638	appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], -1);
1639	if (next_relatives>1) error = "Island hopping uses only one relative";
1640	}
1641	else {
1642	LooseBases loose;
1643	LooseBases loose_for_next_relative;
1644
1645	int unaligned_positions;
1646	{
1647	CompactedSubSequence *alignedSequence = readCompactedSequence(gb_toAlign, alignment, &error, 0, 0, ali_params.range);
1648
1649	unaligned_positions = loose.follow_ali_change_and_append(unaligned_bases, AliChange(toAlignSequence, alignedSequence));
1650	// now loose holds the unaligned (and recalculated) parts from last relative
1651	delete alignedSequence;
1652	}
1653
1654	if (!error) {
1655	int toalign_positions = toAlignSequence->length();
1656	if (unaligned_positions<toalign_positions) {
1657	appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], toalign_positions-unaligned_positions);
1658	}
1659	}
1660
1661	for (i=1; i<next_relatives && !error; i++) {
1662	loose.append(loose_for_next_relative);
1663	int unaligned_positions_for_next = 0;
1664
1665	while (!loose.is_empty() && !error) {
1666	ExplicitRange partRange(intersection(loose.recall(), ali_params.range));
1667	CompactedSubSequence *alignToPart = readCompactedSequence(gb_reference[i], alignment, &error, 0, 0, partRange);
1668
1669	if (!error) {
1670	long part_chksum;
1671	CompactedSubSequence *toAlignPart = readCompactedSequence(gb_toAlign, alignment, &error, 0, &part_chksum, partRange);
1672
1673	fa_assert(contradicted(error, toAlignPart));
1674
1675	if (!error) {
1676	AlignParams loose_ali_params = { ali_params.report, ali_params.showGapsMessages, partRange };
1677
1678	FastSearchSequence referenceFastSeq(*alignToPart);
1679	error = alignCompactedTo(toAlignPart, &referenceFastSeq,
1680	max_seq_length, alignment, part_chksum,
1681	gb_toAlign, gb_reference[i], loose_ali_params);
1682	if (!error) {
1683	CompactedSubSequence *alignedPart = readCompactedSequence(gb_toAlign, alignment, &error, 0, 0, partRange);
1684	if (!error) {
1685	unaligned_positions_for_next += loose_for_next_relative.follow_ali_change_and_append(unaligned_bases, AliChange(toAlignPart, alignedPart));
1686	}
1687	delete alignedPart;
1688	}
1689	}
1690	delete toAlignPart;
1691	}
1692	delete alignToPart;
1693	}
1694
1695	if (!error) {
1696	fa_assert(unaligned_positions_for_next <= unaligned_positions); // means: number of unaligned positions has increased by use of relative
1697	if (unaligned_positions_for_next<unaligned_positions) {
1698	appendNameAndUsedBasePositions(&used_relatives, gb_reference[i], unaligned_positions-unaligned_positions_for_next);
1699	unaligned_positions = unaligned_positions_for_next;
1700	}
1701	}
1702	}
1703	}
1704	}
1705
1706	if (!error) { // write used relatives to db-field
1707	error = GBT_write_string(gb_toAlign, "used_rels", used_relatives);
1708	}
1709	free(used_relatives);
1710	}
1711
1712	delete alignToSequence;
1713	}
1714
1715	delete [] gb_reference;
1716	}
1717	}
1718
1719	delete toAlignSequence;
1720
1721	for (i=0; i<next_relatives; i++) freenull(nearestRelative[i]);
1722	delete [] nearestRelative;
1723
1724	return error;
1725	}
1726
1727	// ------------------------
1728	// AlignmentReference
1729
1730	class AlignmentReference : virtual Noncopyable {
1731	GB_CSTR alignment;
1732	int max_seq_length;
1733	const AlignParams& ali_params;
1734
1735	public:
1736	AlignmentReference(GB_CSTR alignment_,
1737	int max_seq_length_,
1738	const AlignParams& ali_params_)
1739	: alignment(alignment_),
1740	max_seq_length(max_seq_length_),
1741	ali_params(ali_params_)
1742	{}
1743	virtual ~AlignmentReference() {}
1744
1745	virtual ARB_ERROR align_to(GBDATA *gb_toalign) const = 0;
1746
1747	GB_CSTR get_alignment() const { return alignment; }
1748	int get_max_seq_length() const { return max_seq_length; }
1749	const AlignParams& get_ali_params() const { return ali_params; }
1750	};
1751
1752
1753	#if defined(WARN_TODO)
1754	#warning make alignTo a member of ExplicitReference (or of AlignmentReference)
1755	#warning let alignToGroupConsensus and alignToNextRelative use ExplicitReference
1756	#endif
1757
1758	class ExplicitReference: public AlignmentReference { // derived from a Noncopyable
1759	const FastSearchSequence *targetSequence;
1760	GBDATA *gb_alignTo;
1761
1762	public:
1763	ExplicitReference(GB_CSTR alignment_,
1764	const FastSearchSequence *targetSequence_,
1765	GBDATA *gb_alignTo_,
1766	int max_seq_length_,
1767	const AlignParams& ali_params_)
1768	: AlignmentReference(alignment_, max_seq_length_, ali_params_),
1769	targetSequence(targetSequence_),
1770	gb_alignTo(gb_alignTo_)
1771	{}
1772
1773	ARB_ERROR align_to(GBDATA *gb_toalign) const OVERRIDE {
1774	return alignTo(gb_toalign, get_alignment(), targetSequence, gb_alignTo, get_max_seq_length(), get_ali_params());
1775	}
1776	};
1777
1778	#if defined(WARN_TODO)
1779	#warning make alignToGroupConsensus a member of ConsensusReference
1780	#endif
1781
1782	class ConsensusReference: public AlignmentReference {
1783	Aligner_get_consensus_func get_consensus;
1784
1785	public:
1786	ConsensusReference(GB_CSTR alignment_,
1787	Aligner_get_consensus_func get_consensus_,
1788	int max_seq_length_,
1789	const AlignParams& ali_params_)
1790	: AlignmentReference(alignment_, max_seq_length_, ali_params_),
1791	get_consensus(get_consensus_)
1792	{}
1793
1794	ARB_ERROR align_to(GBDATA *gb_toalign) const OVERRIDE {
1795	return alignToGroupConsensus(gb_toalign, get_alignment(), get_consensus, get_max_seq_length(), get_ali_params());
1796	}
1797	};
1798
1799	#if defined(WARN_TODO)
1800	#warning make alignToNextRelative a member of SearchRelativesReference
1801	#endif
1802
1803	class SearchRelativesReference: public AlignmentReference {
1804	SearchRelativeParams& relSearch;
1805	FA_turn turnAllowed;
1806
1807	public:
1808	SearchRelativesReference(SearchRelativeParams& relSearch_,
1809	int max_seq_length_,
1810	FA_turn turnAllowed_,
1811	GB_CSTR alignment_,
1812	const AlignParams& ali_params_)
1813	: AlignmentReference(alignment_, max_seq_length_, ali_params_),
1814	relSearch(relSearch_),
1815	turnAllowed(turnAllowed_)
1816	{}
1817
1818	ARB_ERROR align_to(GBDATA *gb_toalign) const OVERRIDE {
1819	return alignToNextRelative(relSearch, get_max_seq_length(), turnAllowed, get_alignment(), gb_toalign, get_ali_params());
1820	}
1821	};
1822
1823
1824	// ----------------
1825	// Aligner
1826
1827	class Aligner : virtual Noncopyable {
1828	GBDATA *gb_main;
1829
1830	// define alignment target(s):
1831	FA_alignTarget alignWhat;
1832	GB_CSTR alignment; // name of alignment to use (==NULL -> use default)
1833	GB_CSTR toalign; // name of species to align (used if alignWhat == FA_CURRENT)
1834	Aligner_get_first_selected_species get_first_selected_species; // used if alignWhat == FA_SELECTED
1835	Aligner_get_next_selected_species get_next_selected_species; // --- "" ---
1836
1837	// define reference sequence(s):
1838	GB_CSTR reference; // name of reference species
1839	Aligner_get_consensus_func get_consensus; // function to get consensus seq
1840	SearchRelativeParams& relSearch; // params to search for relatives
1841
1842	// general params:
1843	FA_turn turnAllowed;
1844	const AlignParams& ali_params;
1845	int maxProtection; // protection level
1846
1847	// -------------------- new members
1848	int wasNotAllowedToAlign; // number of failures caused by wrong protection
1849	int err_count; // count errors
1850	bool continue_on_error; /* true -> run single alignments in separate transactions.
1851	* If one target fails, continue with rest.
1852	* false -> run all in one transaction
1853	* One fails -> all fail!
1854	*/
1855	FA_errorAction error_action;
1856
1857	typedef std::list<GBDATA*> GBDATAlist;
1858	GBDATAlist species_to_mark; // species that will be marked after aligning
1859
1860	ARB_ERROR alignToReference(GBDATA *gb_toalign, const AlignmentReference& ref);
1861	ARB_ERROR alignTargetsToReference(const AlignmentReference& ref, GBDATA *gb_species_data);
1862
1863	ARB_ERROR alignToExplicitReference(GBDATA *gb_species_data, int max_seq_length);
1864	ARB_ERROR alignToConsensus(GBDATA *gb_species_data, int max_seq_length);
1865	ARB_ERROR alignToRelatives(GBDATA *gb_species_data, int max_seq_length);
1866
1867	void triggerAction(GBDATA *gb_species, bool has_been_aligned) {
1868	bool mark = false;
1869	switch (error_action) {
1870	case FA_MARK_FAILED: mark = !has_been_aligned; break;
1871	case FA_MARK_ALIGNED: mark = has_been_aligned; break;
1872	case FA_NO_ACTION: mark = false; break;
1873	}
1874	if (mark) species_to_mark.push_back(gb_species);
1875	}
1876
1877	public:
1878
1879	#if defined(WARN_TODO)
1880	#warning pass AlignmentReference from caller (replacing reference parameters)
1881	#endif
1882
1883	Aligner(GBDATA *gb_main_,
1884
1885	// define alignment target(s):
1886	FA_alignTarget alignWhat_,
1887	GB_CSTR alignment_, // name of alignment to use (==NULL -> use default)
1888	GB_CSTR toalign_, // name of species to align (used if alignWhat == FA_CURRENT)
1889	Aligner_get_first_selected_species get_first_selected_species_, // used if alignWhat == FA_SELECTED
1890	Aligner_get_next_selected_species get_next_selected_species_, // --- "" ---
1891
1892	// define reference sequence(s):
1893	GB_CSTR reference_, // name of reference species
1894	Aligner_get_consensus_func get_consensus_, // function to get consensus seq
1895	SearchRelativeParams& relSearch_, // params to search for relatives
1896
1897	// general params:
1898	FA_turn turnAllowed_,
1899	const AlignParams& ali_params_,
1900	int maxProtection_, // protection level
1901	bool continue_on_error_,
1902	FA_errorAction error_action_)
1903	: gb_main(gb_main_),
1904	alignWhat(alignWhat_),
1905	alignment(alignment_),
1906	toalign(toalign_),
1907	get_first_selected_species(get_first_selected_species_),
1908	get_next_selected_species(get_next_selected_species_),
1909	reference(reference_),
1910	get_consensus(get_consensus_),
1911	relSearch(relSearch_),
1912	turnAllowed(turnAllowed_),
1913	ali_params(ali_params_),
1914	maxProtection(maxProtection_),
1915	wasNotAllowedToAlign(0),
1916	err_count(0),
1917	continue_on_error(continue_on_error_),
1918	error_action(continue_on_error ? error_action_ : FA_NO_ACTION)
1919	{}
1920
1921	ARB_ERROR run();
1922	};
1923
1924	ARB_ERROR Aligner::alignToReference(GBDATA *gb_toalign, const AlignmentReference& ref) {
1925	int myProtection = GB_read_security_write(GBT_find_sequence(gb_toalign, alignment));
1926	ARB_ERROR error;
1927
1928	if (myProtection<=maxProtection) {
1929	// Depending on 'continue_on_error' we either
1930	// * stop aligning if an error occurs or
1931	// * run the alignment of each species in its own transaction, ignore but report errors and
1932	// optionally mark aligned or failed species.
1933
1934	if (continue_on_error) {
1935	fa_assert(GB_get_transaction_level(gb_main) == 1);
1936	error = GB_end_transaction(gb_main, error); // end global transaction
1937	}
1938
1939	if (!error) {
1940	error = GB_push_transaction(gb_main);
1941	if (!error) error = ref.align_to(gb_toalign);
1942	error = GB_end_transaction(gb_main, error);
1943
1944	if (error) err_count++;
1945	triggerAction(gb_toalign, !error);
1946	}
1947
1948	if (continue_on_error) {
1949	if (error) {
1950	GB_warning(error.deliver());
1951	error = NULL;
1952	}
1953	error = GB_begin_transaction(gb_main); // re-open global transaction
1954	}
1955	}
1956	else {
1957	wasNotAllowedToAlign++;
1958	triggerAction(gb_toalign, false);
1959	}
1960
1961	return error;
1962	}
1963
1964	ARB_ERROR Aligner::alignTargetsToReference(const AlignmentReference& ref, GBDATA *gb_species_data) {
1965	ARB_ERROR error;
1966
1967	fa_assert(GB_get_transaction_level(gb_main) == 1);
1968
1969	switch (alignWhat) {
1970	case FA_CURRENT: { // align one sequence
1971	fa_assert(toalign);
1972
1973	GBDATA *gb_toalign = GBT_find_species_rel_species_data(gb_species_data, toalign);
1974	if (!gb_toalign) {
1975	error = species_not_found(toalign);
1976	}
1977	else {
1978	currentSequenceNumber = overallSequenceNumber = 1;
1979	error = alignToReference(gb_toalign, ref);
1980	}
1981	break;
1982	}
1983	case FA_MARKED: { // align all marked sequences
1984	int count = GBT_count_marked_species(gb_main);
1985	GBDATA *gb_species = GBT_first_marked_species_rel_species_data(gb_species_data);
1986
1987	arb_progress progress("Aligning marked species", count);
1988	progress.auto_subtitles("Species");
1989
1990	currentSequenceNumber = 1;
1991	overallSequenceNumber = count;
1992
1993	while (gb_species && !error) {
1994	error = alignToReference(gb_species, ref);
1995	progress.inc_and_check_user_abort(error);
1996	gb_species = GBT_next_marked_species(gb_species);
1997	}
1998	break;
1999	}
2000	case FA_SELECTED: { // align all selected species
2001	int count;
2002	GBDATA *gb_species = get_first_selected_species(&count);
2003
2004
2005	currentSequenceNumber = 1;
2006	overallSequenceNumber = count;
2007
2008	if (!gb_species) {
2009	aw_message("There is no selected species!");
2010	}
2011	else {
2012	arb_progress progress("Aligning selected species", count);
2013	progress.auto_subtitles("Species");
2014
2015	while (gb_species && !error) {
2016	error = alignToReference(gb_species, ref);
2017	progress.inc_and_check_user_abort(error);
2018	gb_species = get_next_selected_species();
2019	}
2020	}
2021	break;
2022	}
2023	}
2024
2025	fa_assert(GB_get_transaction_level(gb_main) == 1);
2026	return error;
2027	}
2028
2029	ARB_ERROR Aligner::alignToExplicitReference(GBDATA *gb_species_data, int max_seq_length) {
2030	ARB_ERROR error;
2031	GBDATA *gb_reference = GBT_find_species_rel_species_data(gb_species_data, reference);
2032
2033	if (!gb_reference) {
2034	error = species_not_found(reference);
2035	}
2036	else {
2037	long referenceChksum;
2038	CompactedSubSequence *referenceSeq = readCompactedSequence(gb_reference, alignment, &error, NULL, &referenceChksum, ali_params.range);
2039
2040	if (island_hopper) {
2041	#if defined(WARN_TODO)
2042	#warning setting island_hopper reference has to be done in called function (seems that it is NOT done for alignToConsensus and alignToRelatives). First get tests in place!
2043	#endif
2044	GBDATA *gb_seq = GBT_find_sequence(gb_reference, alignment); // get sequence
2045	if (gb_seq) {
2046	long length = GB_read_string_count(gb_seq);
2047	const char *data = GB_read_char_pntr(gb_seq);
2048
2049	island_hopper->set_ref_sequence(data);
2050	island_hopper->set_alignment_length(length);
2051	}
2052	}
2053
2054
2055	if (!error) {
2056	#if defined(WARN_TODO)
2057	#warning do not pass FastSearchSequence to ExplicitReference, instead pass sequence and length (ExplicitReference shall create it itself)
2058	#endif
2059
2060	FastSearchSequence referenceFastSeq(*referenceSeq);
2061	ExplicitReference target(alignment, &referenceFastSeq, gb_reference, max_seq_length, ali_params);
2062
2063	error = alignTargetsToReference(target, gb_species_data);
2064	}
2065	delete referenceSeq;
2066	}
2067	return error;
2068	}
2069
2070	ARB_ERROR Aligner::alignToConsensus(GBDATA *gb_species_data, int max_seq_length) {
2071	return alignTargetsToReference(ConsensusReference(alignment, get_consensus, max_seq_length, ali_params),
2072	gb_species_data);
2073	}
2074
2075	ARB_ERROR Aligner::alignToRelatives(GBDATA *gb_species_data, int max_seq_length) {
2076
2077	return alignTargetsToReference(SearchRelativesReference(relSearch, max_seq_length, turnAllowed, alignment, ali_params),
2078	gb_species_data);
2079	}
2080
2081	ARB_ERROR Aligner::run() {
2082	// (reference == NULL && get_consensus==NULL -> use next relative for (each) sequence)
2083
2084	fa_assert(GB_get_transaction_level(gb_main) == 0); // no open transaction allowed here!
2085	ARB_ERROR error = GB_push_transaction(gb_main);
2086	bool search_by_pt_server = reference==NULL && get_consensus==NULL;
2087
2088	err_count = 0;
2089	wasNotAllowedToAlign = 0; // incremented for every sequence which has higher protection level (and was not aligned)
2090	species_to_mark.clear();
2091
2092	fa_assert(reference==NULL \|\| get_consensus==NULL); // can't do both modes
2093
2094	if (turnAllowed != FA_TURN_NEVER) {
2095	if ((ali_params.range.is_part()) \|\| !search_by_pt_server) {
2096	// if not selected 'Range/Whole sequence' or not selected 'Reference/Auto search..'
2097	turnAllowed = FA_TURN_NEVER; // then disable mirroring for the current call
2098	}
2099	}
2100
2101	if (!error && !alignment) {
2102	alignment = (GB_CSTR)GBT_get_default_alignment(gb_main); // get default alignment
2103	if (!alignment) error = "No default alignment";
2104	}
2105
2106	if (!error && alignment) {
2107	global_alignmentType = GBT_get_alignment_type(gb_main, alignment);
2108	if (search_by_pt_server) {
2109	GB_alignment_type pt_server_alignmentType = GBT_get_alignment_type(gb_main, relSearch.pt_server_alignment);
2110
2111	if (pt_server_alignmentType != GB_AT_RNA &&
2112	pt_server_alignmentType != GB_AT_DNA) {
2113	error = "pt_servers only support RNA/DNA sequences.\n"
2114	"In the aligner window you may specify a RNA/DNA alignment \n"
2115	"and use a pt_server build on that alignment.";
2116	}
2117	}
2118	}
2119
2120	if (!error) {
2121	GBDATA *gb_species_data = GBT_get_species_data(gb_main);
2122	int max_seq_length = GBT_get_alignment_len(gb_main, alignment);
2123
2124	if (reference) error = alignToExplicitReference(gb_species_data, max_seq_length);
2125	else if (get_consensus) error = alignToConsensus(gb_species_data, max_seq_length);
2126	else error = alignToRelatives(gb_species_data, max_seq_length);
2127	}
2128
2129	ClustalV_exit();
2130	unaligned_bases.clear();
2131
2132	error = GB_end_transaction(gb_main, error); // close global transaction
2133
2134	if (wasNotAllowedToAlign>0) {
2135	const char *mess = GBS_global_string("%i species were not aligned (because of protection level)", wasNotAllowedToAlign);
2136	aw_message(mess);
2137	}
2138
2139	if (err_count) {
2140	aw_message_if(error);
2141	error = GBS_global_string("Aligner produced %i error%c", err_count, err_count==1 ? '\0' : 's');
2142	}
2143
2144	if (error_action != FA_NO_ACTION) {
2145	fa_assert(continue_on_error);
2146
2147	GB_transaction ta(gb_main);
2148	GBT_mark_all(gb_main, 0);
2149	for (GBDATAlist::iterator sp = species_to_mark.begin(); sp != species_to_mark.end(); ++sp) {
2150	GB_write_flag(*sp, 1);
2151	}
2152
2153	const char *whatsMarked = (error_action == FA_MARK_ALIGNED) ? "aligned" : "failed";
2154	size_t markCount = species_to_mark.size();
2155	if (markCount>0) {
2156	const char *msg = GBS_global_string("%zu %s species %s been marked",
2157	markCount,
2158	whatsMarked,
2159	(markCount == 1) ? "has" : "have");
2160	aw_message(msg);
2161	}
2162	}
2163
2164	return error;
2165	}
2166
2167	void FastAligner_start(AW_window *aw, AW_CL cl_AlignDataAccess) {
2168	AW_root *root = aw->get_root();
2169	char *reference = NULL; // align against next relatives
2170	char *toalign = NULL; // align marked species
2171	ARB_ERROR error = NULL;
2172	const AlignDataAccess data_access = (const AlignDataAccess )cl_AlignDataAccess;
2173	int get_consensus = 0;
2174	int pt_server_id = -1;
2175
2176	Aligner_get_first_selected_species get_first_selected_species = 0;
2177	Aligner_get_next_selected_species get_next_selected_species = 0;
2178
2179	fa_assert(island_hopper == 0);
2180	if (root->awar(FA_AWAR_USE_ISLAND_HOPPING)->read_int()) {
2181	island_hopper = new IslandHopping;
2182	if (root->awar(FA_AWAR_USE_SECONDARY)->read_int()) {
2183	if (data_access->helix_string) island_hopper->set_helix(data_access->helix_string);
2184	else error = "Warning: No HELIX found. Can't use secondary structure";
2185	}
2186
2187	if (!error) {
2188	island_hopper->set_parameters(root->awar(FA_AWAR_ESTIMATE_BASE_FREQ)->read_int(),
2189	root->awar(FA_AWAR_BASE_FREQ_T)->read_float(),
2190	root->awar(FA_AWAR_BASE_FREQ_C)->read_float(),
2191	root->awar(FA_AWAR_BASE_FREQ_A)->read_float(),
2192	root->awar(FA_AWAR_BASE_FREQ_C)->read_float(),
2193	root->awar(FA_AWAR_SUBST_PARA_CT)->read_float(),
2194	root->awar(FA_AWAR_SUBST_PARA_AT)->read_float(),
2195	root->awar(FA_AWAR_SUBST_PARA_GT)->read_float(),
2196	root->awar(FA_AWAR_SUBST_PARA_AC)->read_float(),
2197	root->awar(FA_AWAR_SUBST_PARA_CG)->read_float(),
2198	root->awar(FA_AWAR_SUBST_PARA_AG)->read_float(),
2199	root->awar(FA_AWAR_EXPECTED_DISTANCE)->read_float(),
2200	root->awar(FA_AWAR_STRUCTURE_SUPPLEMENT)->read_float(),
2201	root->awar(FA_AWAR_GAP_A)->read_float(),
2202	root->awar(FA_AWAR_GAP_B)->read_float(),
2203	root->awar(FA_AWAR_GAP_C)->read_float(),
2204	root->awar(FA_AWAR_THRESHOLD)->read_float()
2205	);
2206	}
2207	}
2208
2209	FA_alignTarget alignWhat = static_cast<FA_alignTarget>(root->awar(FA_AWAR_TO_ALIGN)->read_int());
2210	if (!error) {
2211	switch (alignWhat) {
2212	case FA_CURRENT: { // align current species
2213	toalign = root->awar(AWAR_SPECIES_NAME)->read_string();
2214	break;
2215	}
2216	case FA_MARKED: { // align marked species
2217	break;
2218	}
2219	case FA_SELECTED: { // align selected species
2220	get_first_selected_species = data_access->get_first_selected_species;
2221	get_next_selected_species = data_access->get_next_selected_species;
2222	break;
2223	}
2224	default: {
2225	fa_assert(0);
2226	break;
2227	}
2228	}
2229
2230	switch (static_cast<FA_reference>(root->awar(FA_AWAR_REFERENCE)->read_int())) {
2231	case FA_REF_EXPLICIT: // align against specified species
2232	reference = root->awar(FA_AWAR_REFERENCE_NAME)->read_string();
2233	break;
2234
2235	case FA_REF_CONSENSUS: // align against group consensus
2236	if (data_access->get_group_consensus) {
2237	get_consensus = 1;
2238	}
2239	else {
2240	error = "Can't get group consensus here.";
2241	}
2242	break;
2243
2244	case FA_REF_RELATIVES: // align against species searched via pt_server
2245	pt_server_id = root->awar(AWAR_PT_SERVER)->read_int();
2246	if (pt_server_id<0) {
2247	error = "No pt_server selected";
2248	}
2249	break;
2250
2251	default: fa_assert(0);
2252	break;
2253	}
2254	}
2255
2256	RangeList ranges;
2257	bool autoRestrictRange4nextRelSearch = true;
2258
2259	if (!error) {
2260	switch (static_cast<FA_range>(root->awar(FA_AWAR_RANGE)->read_int())) {
2261	case FA_WHOLE_SEQUENCE:
2262	autoRestrictRange4nextRelSearch = false;
2263	ranges.add(PosRange::whole());
2264	break;
2265
2266	case FA_AROUND_CURSOR: {
2267	int curpos = root->awar(AWAR_CURSOR_POSITION_LOCAL)->read_int();
2268	int size = root->awar(FA_AWAR_AROUND)->read_int();
2269
2270	ranges.add(PosRange(curpos-size, curpos+size));
2271	break;
2272	}
2273	case FA_SELECTED_RANGE: {
2274	PosRange range;
2275	if (!data_access->get_selected_range(range)) {
2276	error = "There is no selected species!";
2277	}
2278	else {
2279	ranges.add(range);
2280	}
2281	break;
2282	}
2283
2284	case FA_SAI_MULTI_RANGE: {
2285	GB_transaction ta(data_access->gb_main);
2286
2287	char *sai_name = root->awar(FA_AWAR_SAI_RANGE_NAME)->read_string();
2288	char *aliuse = GBT_get_default_alignment(data_access->gb_main);
2289
2290	GBDATA *gb_sai = GBT_expect_SAI(data_access->gb_main, sai_name);
2291	if (!gb_sai) error = GB_await_error();
2292	else {
2293	GBDATA *gb_data = GBT_find_sequence(gb_sai, aliuse);
2294	if (!gb_data) {
2295	error = GB_have_error()
2296	? GB_await_error()
2297	: GBS_global_string("SAI '%s' has no data in alignment '%s'", sai_name, aliuse);
2298	}
2299	else {
2300	char *sai_data = GB_read_string(gb_data);
2301	char *set_bits = root->awar(FA_AWAR_SAI_RANGE_CHARS)->read_string();
2302
2303	ranges = build_RangeList_from_string(sai_data, set_bits, false);
2304
2305	free(set_bits);
2306	free(sai_data);
2307	}
2308	}
2309	free(aliuse);
2310	free(sai_name);
2311	break;
2312	}
2313	}
2314	}
2315
2316	if (!error) {
2317	for (RangeList::iterator r = ranges.begin(); r != ranges.end() && !error; ++r) {
2318	PosRange range = *r;
2319
2320	GBDATA *gb_main = data_access->gb_main;
2321	char *editor_alignment = 0;
2322	{
2323	GB_transaction ta(gb_main);
2324	char *default_alignment = GBT_get_default_alignment(gb_main);
2325
2326	editor_alignment = root->awar_string(AWAR_EDITOR_ALIGNMENT, default_alignment)->read_string();
2327	free(default_alignment);
2328	}
2329
2330	char *pt_server_alignment = root->awar(FA_AWAR_PT_SERVER_ALIGNMENT)->read_string();
2331	PosRange relRange = PosRange::whole(); // unrestricted!
2332
2333	if (autoRestrictRange4nextRelSearch) {
2334	AW_awar *awar_relrange = root->awar(FA_AWAR_RELATIVE_RANGE);
2335	const char *relrange = awar_relrange->read_char_pntr();
2336	if (relrange[0]) {
2337	int region_plus = atoi(relrange);
2338
2339	fa_assert(range.is_part());
2340
2341	relRange = PosRange(range.start()-region_plus, range.end()+region_plus); // restricted
2342	awar_relrange->write_as_string(GBS_global_string("%i", region_plus)); // set awar to detected value (avoid misbehavior when it contains ' ' or similar)
2343	}
2344	}
2345
2346	if (island_hopper) {
2347	island_hopper->set_range(range);
2348	range = PosRange::whole();
2349	}
2350
2351	SearchRelativeParams relSearch(new PT_FamilyFinder(gb_main,
2352	pt_server_id,
2353	root->awar(AWAR_NN_OLIGO_LEN)->read_int(),
2354	root->awar(AWAR_NN_MISMATCHES)->read_int(),
2355	root->awar(AWAR_NN_FAST_MODE)->read_int(),
2356	root->awar(AWAR_NN_REL_MATCHES)->read_int(),
2357	RSS_BOTH_MIN), // old scaling as b4 [8520] @@@ make configurable
2358	pt_server_alignment,
2359	root->awar(FA_AWAR_NEXT_RELATIVES)->read_int());
2360
2361	relSearch.getFamilyFinder()->restrict_2_region(relRange);
2362
2363	struct AlignParams ali_params = {
2364	static_cast<FA_report>(root->awar(FA_AWAR_REPORT)->read_int()),
2365	bool(root->awar(FA_AWAR_SHOW_GAPS_MESSAGES)->read_int()),
2366	range
2367	};
2368
2369	{
2370	arb_progress progress("FastAligner");
2371	progress.allow_title_reuse();
2372
2373	int cont_on_error = root->awar(FA_AWAR_CONTINUE_ON_ERROR)->read_int();
2374
2375	Aligner aligner(gb_main,
2376	alignWhat,
2377	editor_alignment,
2378	toalign,
2379	get_first_selected_species,
2380	get_next_selected_species,
2381
2382	reference,
2383	get_consensus ? data_access->get_group_consensus : NULL,
2384	relSearch,
2385
2386	static_cast<FA_turn>(root->awar(FA_AWAR_MIRROR)->read_int()),
2387	ali_params,
2388	root->awar(FA_AWAR_PROTECTION)->read_int(),
2389	cont_on_error,
2390	(FA_errorAction)root->awar(FA_AWAR_ACTION_ON_ERROR)->read_int());
2391	error = aligner.run();
2392
2393	if (error && cont_on_error) {
2394	aw_message_if(error);
2395	error = NULL;
2396	}
2397	}
2398
2399	free(pt_server_alignment);
2400	free(editor_alignment);
2401	}
2402	}
2403
2404	if (island_hopper) {
2405	delete island_hopper;
2406	island_hopper = 0;
2407	}
2408
2409	if (toalign) free(toalign);
2410	aw_message_if(error);
2411	if (data_access->do_refresh) data_access->refresh_display();
2412	}
2413
2414
2415
2416	void FastAligner_create_variables(AW_root *root, AW_default db1) {
2417	root->awar_string(FA_AWAR_REFERENCE_NAME, "", db1);
2418
2419	root->awar_int(FA_AWAR_TO_ALIGN, FA_CURRENT, db1);
2420	root->awar_int(FA_AWAR_REFERENCE, FA_REF_EXPLICIT, db1);
2421	root->awar_int(FA_AWAR_RANGE, FA_WHOLE_SEQUENCE, db1);
2422
2423	AW_awar *ali_protect = root->awar_int(FA_AWAR_PROTECTION, 0, db1);
2424	if (ARB_in_novice_mode(root)) {
2425	ali_protect->write_int(0); // reset protection for noobs
2426	}
2427
2428	root->awar_int(FA_AWAR_AROUND, 25, db1);
2429	root->awar_int(FA_AWAR_MIRROR, FA_TURN_INTERACTIVE, db1);
2430	root->awar_int(FA_AWAR_REPORT, FA_NO_REPORT, db1);
2431	root->awar_int(FA_AWAR_SHOW_GAPS_MESSAGES, 1, db1);
2432	root->awar_int(FA_AWAR_CONTINUE_ON_ERROR, 1, db1);
2433	root->awar_int(FA_AWAR_ACTION_ON_ERROR, FA_NO_ACTION, db1);
2434	root->awar_int(FA_AWAR_USE_SECONDARY, 0, db1);
2435	root->awar_int(AWAR_PT_SERVER, 0, db1);
2436	root->awar_int(FA_AWAR_NEXT_RELATIVES, 1, db1)->set_minmax(1, 100);
2437
2438	root->awar_string(FA_AWAR_RELATIVE_RANGE, "", db1);
2439	root->awar_string(FA_AWAR_PT_SERVER_ALIGNMENT, root->awar(AWAR_DEFAULT_ALIGNMENT)->read_char_pntr(), db1);
2440
2441	root->awar_string(FA_AWAR_SAI_RANGE_NAME, "", db1);
2442	root->awar_string(FA_AWAR_SAI_RANGE_CHARS, "x1", db1);
2443
2444	// island hopping:
2445
2446	root->awar_int(FA_AWAR_USE_ISLAND_HOPPING, 0, db1);
2447
2448	root->awar_int(FA_AWAR_ESTIMATE_BASE_FREQ, 1, db1);
2449
2450	root->awar_float(FA_AWAR_BASE_FREQ_A, 0.25, db1);
2451	root->awar_float(FA_AWAR_BASE_FREQ_C, 0.25, db1);
2452	root->awar_float(FA_AWAR_BASE_FREQ_G, 0.25, db1);
2453	root->awar_float(FA_AWAR_BASE_FREQ_T, 0.25, db1);
2454
2455	root->awar_float(FA_AWAR_SUBST_PARA_AC, 1.0, db1);
2456	root->awar_float(FA_AWAR_SUBST_PARA_AG, 4.0, db1);
2457	root->awar_float(FA_AWAR_SUBST_PARA_AT, 1.0, db1);
2458	root->awar_float(FA_AWAR_SUBST_PARA_CG, 1.0, db1);
2459	root->awar_float(FA_AWAR_SUBST_PARA_CT, 4.0, db1);
2460	root->awar_float(FA_AWAR_SUBST_PARA_GT, 1.0, db1);
2461
2462	root->awar_float(FA_AWAR_EXPECTED_DISTANCE, 0.3, db1);
2463	root->awar_float(FA_AWAR_STRUCTURE_SUPPLEMENT, 0.5, db1);
2464	root->awar_float(FA_AWAR_THRESHOLD, 0.005, db1);
2465
2466	root->awar_float(FA_AWAR_GAP_A, 8.0, db1);
2467	root->awar_float(FA_AWAR_GAP_B, 4.0, db1);
2468	root->awar_float(FA_AWAR_GAP_C, 7.0, db1);
2469
2470	AWTC_create_common_next_neighbour_vars(root);
2471	}
2472
2473	void FastAligner_set_align_current(AW_root *root, AW_default db1) {
2474	root->awar_int(FA_AWAR_TO_ALIGN, FA_CURRENT, db1);
2475	}
2476
2477	void FastAligner_set_reference_species(AW_window * /* aww */, AW_CL cl_AW_root) {
2478	// sets the aligner reference species to current species
2479	AW_root root = (AW_root)cl_AW_root;
2480	char *specName = root->awar(AWAR_SPECIES_NAME)->read_string();
2481
2482	root->awar(FA_AWAR_REFERENCE_NAME)->write_string(specName);
2483	free(specName);
2484	}
2485
2486	static AW_window create_island_hopping_window(AW_root root, AW_CL) {
2487	AW_window_simple *aws = new AW_window_simple;
2488
2489	aws->init(root, "ISLAND_HOPPING_PARA", "Parameters for Island Hopping");
2490	aws->load_xfig("faligner/islandhopping.fig");
2491
2492	aws->at("close");
2493	aws->callback ((AW_CB0)AW_POPDOWN);
2494	aws->create_button("CLOSE", "CLOSE", "O");
2495
2496	aws->at("help");
2497	aws->callback(makeHelpCallback("islandhopping.hlp"));
2498	aws->create_button("HELP", "HELP");
2499
2500	aws->at("use_secondary");
2501	aws->label("Use secondary structure (only for re-align)");
2502	aws->create_toggle(FA_AWAR_USE_SECONDARY);
2503
2504	aws->at("freq");
2505	aws->create_toggle_field(FA_AWAR_ESTIMATE_BASE_FREQ, "Base freq.", "B");
2506	aws->insert_default_toggle("Estimate", "E", 1);
2507	aws->insert_toggle("Define here: ", "D", 0);
2508	aws->update_toggle_field();
2509
2510	aws->at("freq_a"); aws->label("A:"); aws->create_input_field(FA_AWAR_BASE_FREQ_A, 4);
2511	aws->at("freq_c"); aws->label("C:"); aws->create_input_field(FA_AWAR_BASE_FREQ_C, 4);
2512	aws->at("freq_g"); aws->label("G:"); aws->create_input_field(FA_AWAR_BASE_FREQ_G, 4);
2513	aws->at("freq_t"); aws->label("T:"); aws->create_input_field(FA_AWAR_BASE_FREQ_T, 4);
2514
2515	int xpos[4], ypos[4];
2516	{
2517	aws->button_length(1);
2518
2519	int dummy;
2520	aws->at("h_a"); aws->get_at_position(&xpos[0], &dummy); aws->create_button(NULL, "A");
2521	aws->at("h_c"); aws->get_at_position(&xpos[1], &dummy); aws->create_button(NULL, "C");
2522	aws->at("h_g"); aws->get_at_position(&xpos[2], &dummy); aws->create_button(NULL, "G");
2523	aws->at("h_t"); aws->get_at_position(&xpos[3], &dummy); aws->create_button(NULL, "T");
2524
2525	aws->at("v_a"); aws->get_at_position(&dummy, &ypos[0]); aws->create_button(NULL, "A");
2526	aws->at("v_c"); aws->get_at_position(&dummy, &ypos[1]); aws->create_button(NULL, "C");
2527	aws->at("v_g"); aws->get_at_position(&dummy, &ypos[2]); aws->create_button(NULL, "G");
2528	aws->at("v_t"); aws->get_at_position(&dummy, &ypos[3]); aws->create_button(NULL, "T");
2529	}
2530
2531	aws->at("subst"); aws->create_button(NULL, "Substitution rate parameters:");
2532
2533	#define XOFF -25
2534	#define YOFF 0
2535
2536	aws->at(xpos[1]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AC, 4);
2537	aws->at(xpos[2]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AG, 4);
2538	aws->at(xpos[3]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AT, 4);
2539	aws->at(xpos[2]+XOFF, ypos[1]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_CG, 4);
2540	aws->at(xpos[3]+XOFF, ypos[1]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_CT, 4);
2541	aws->at(xpos[3]+XOFF, ypos[2]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_GT, 4);
2542
2543	#undef XOFF
2544	#undef YOFF
2545
2546	aws->label_length(22);
2547
2548	aws->at("dist");
2549	aws->label("Expected distance");
2550	aws->create_input_field(FA_AWAR_EXPECTED_DISTANCE, 5);
2551
2552	aws->at("supp");
2553	aws->label("Structure supplement");
2554	aws->create_input_field(FA_AWAR_STRUCTURE_SUPPLEMENT, 5);
2555
2556	aws->at("thres");
2557	aws->label("Threshold");
2558	aws->create_input_field(FA_AWAR_THRESHOLD, 5);
2559
2560	aws->label_length(10);
2561
2562	aws->at("gapA");
2563	aws->label("Gap A");
2564	aws->create_input_field(FA_AWAR_GAP_A, 5);
2565
2566	aws->at("gapB");
2567	aws->label("Gap B");
2568	aws->create_input_field(FA_AWAR_GAP_B, 5);
2569
2570	aws->at("gapC");
2571	aws->label("Gap C");
2572	aws->create_input_field(FA_AWAR_GAP_C, 5);
2573
2574	return (AW_window *)aws;
2575	}
2576
2577	static AW_window create_family_settings_window(AW_root root) {
2578	static AW_window_simple *aws = 0;
2579
2580	if (!aws) {
2581	aws = new AW_window_simple;
2582
2583	aws->init(root, "FAMILY_PARAMS", "Family search parameters");
2584	aws->load_xfig("faligner/family_settings.fig");
2585
2586	aws->at("close");
2587	aws->callback ((AW_CB0)AW_POPDOWN);
2588	aws->create_button("CLOSE", "CLOSE", "O");
2589
2590	aws->at("help");
2591	aws->callback(makeHelpCallback("next_neighbours_common.hlp"));
2592	aws->create_button("HELP", "HELP");
2593
2594	AWTC_create_common_next_neighbour_fields(aws);
2595	}
2596
2597	return aws;
2598	}
2599
2600	AW_window FastAligner_create_window(AW_root root, const AlignDataAccess *data_access) {
2601	AW_window_simple *aws = new AW_window_simple;
2602
2603	aws->init(root, "INTEGRATED_ALIGNERS", INTEGRATED_ALIGNERS_TITLE);
2604	aws->load_xfig("faligner/faligner.fig");
2605
2606	aws->label_length(10);
2607	aws->button_length(10);
2608
2609	aws->at("close");
2610	aws->callback ((AW_CB0)AW_POPDOWN);
2611	aws->create_button("CLOSE", "CLOSE", "O");
2612
2613	aws->at("help");
2614	aws->callback(makeHelpCallback("faligner.hlp"));
2615	aws->create_button("HELP", "HELP");
2616
2617	aws->at("aligner");
2618	aws->create_toggle_field(FA_AWAR_USE_ISLAND_HOPPING, "Aligner", "A");
2619	aws->insert_default_toggle("Fast aligner", "F", 0);
2620	aws->sens_mask(AWM_EXP);
2621	aws->insert_toggle ("Island Hopping", "I", 1);
2622	aws->sens_mask(AWM_ALL);
2623	aws->update_toggle_field();
2624
2625	aws->button_length(12);
2626	aws->at("island_para");
2627	aws->callback(AW_POPUP, (AW_CL)create_island_hopping_window, 0);
2628	aws->sens_mask(AWM_EXP);
2629	aws->create_button("island_para", "Parameters", "");
2630	aws->sens_mask(AWM_ALL);
2631
2632	aws->button_length(10);
2633
2634	aws->at("rev_compl");
2635	aws->callback(build_reverse_complement, (AW_CL)data_access);
2636	aws->create_button("reverse_complement", "Turn now!", "");
2637
2638	aws->at("what");
2639	aws->create_toggle_field(FA_AWAR_TO_ALIGN, "Align what?", "A");
2640	aws->insert_toggle ("Current Species:", "A", FA_CURRENT);
2641	aws->insert_default_toggle("Marked Species", "M", FA_MARKED);
2642	aws->insert_toggle ("Selected Species", "S", FA_SELECTED);
2643	aws->update_toggle_field();
2644
2645	aws->at("swhat");
2646	aws->create_input_field(AWAR_SPECIES_NAME, 2);
2647
2648	aws->at("against");
2649	aws->create_toggle_field(FA_AWAR_REFERENCE, "Reference", "");
2650	aws->insert_toggle ("Species by name:", "S", FA_REF_EXPLICIT);
2651	aws->insert_toggle ("Group consensus", "K", FA_REF_CONSENSUS);
2652	aws->insert_default_toggle("Auto search by pt_server:", "A", FA_REF_RELATIVES);
2653	aws->update_toggle_field();
2654
2655	aws->at("sagainst");
2656	aws->create_input_field(FA_AWAR_REFERENCE_NAME, 2);
2657
2658	aws->at("copy");
2659	aws->callback(FastAligner_set_reference_species, (AW_CL)root);
2660	aws->create_button("Copy", "Copy", "");
2661
2662	aws->label_length(0);
2663	aws->at("pt_server");
2664	awt_create_PTSERVER_selection_button(aws, AWAR_PT_SERVER);
2665
2666	aws->at("relrange");
2667	aws->label("Data from range only, plus");
2668	aws->create_input_field(FA_AWAR_RELATIVE_RANGE, 3);
2669
2670	aws->at("use_ali");
2671	aws->label("Alignment");
2672	aws->create_input_field(FA_AWAR_PT_SERVER_ALIGNMENT, 12);
2673
2674	aws->at("relatives");
2675	aws->label("Number of relatives to use");
2676	aws->create_input_field(FA_AWAR_NEXT_RELATIVES, 3);
2677
2678	aws->at("relSett");
2679	aws->callback(AW_POPUP, (AW_CL)create_family_settings_window, (AW_CL)root);
2680	aws->create_autosize_button("Settings", "More settings", "");
2681
2682	// Range
2683
2684	aws->label_length(10);
2685	aws->at("range");
2686	aws->create_toggle_field(FA_AWAR_RANGE, "Range", "");
2687	aws->insert_default_toggle("Whole sequence", "", FA_WHOLE_SEQUENCE);
2688	aws->insert_toggle ("Positions around cursor: ", "", FA_AROUND_CURSOR);
2689	aws->insert_toggle ("Selected range", "", FA_SELECTED_RANGE);
2690	aws->insert_toggle ("Multi-Range by SAI", "", FA_SAI_MULTI_RANGE);
2691	aws->update_toggle_field();
2692
2693	aws->at("around");
2694	aws->create_input_field(FA_AWAR_AROUND, 2);
2695
2696	aws->at("sai");
2697	awt_create_SAI_selection_button(data_access->gb_main, aws, FA_AWAR_SAI_RANGE_NAME);
2698
2699	aws->at("rchars");
2700	aws->create_input_field(FA_AWAR_SAI_RANGE_CHARS, 2);
2701
2702	// Protection
2703
2704	aws->at("protection");
2705	aws->label("Protection");
2706	aws->create_option_menu(FA_AWAR_PROTECTION, true);
2707	aws->insert_default_option("0", 0, 0);
2708	aws->insert_option("1", 0, 1);
2709	aws->insert_option("2", 0, 2);
2710	aws->insert_option("3", 0, 3);
2711	aws->insert_option("4", 0, 4);
2712	aws->insert_option("5", 0, 5);
2713	aws->insert_option("6", 0, 6);
2714	aws->update_option_menu();
2715
2716	// MirrorCheck
2717
2718	aws->at("mirror");
2719	aws->label("Turn check");
2720	aws->create_option_menu(FA_AWAR_MIRROR, true);
2721	aws->insert_option ("Never turn sequence", "", FA_TURN_NEVER);
2722	aws->insert_default_option("User acknowledgment ", "", FA_TURN_INTERACTIVE);
2723	aws->insert_option ("Automatically turn sequence", "", FA_TURN_ALWAYS);
2724	aws->update_option_menu();
2725
2726	// Report
2727
2728	aws->at("insert");
2729	aws->label("Report");
2730	aws->create_option_menu(FA_AWAR_REPORT, true);
2731	aws->insert_option ("No report", "", FA_NO_REPORT);
2732	aws->sens_mask(AWM_EXP);
2733	aws->insert_default_option("Report to temporary entries", "", FA_TEMP_REPORT);
2734	aws->insert_option ("Report to resident entries", "", FA_REPORT);
2735	aws->sens_mask(AWM_ALL);
2736	aws->update_option_menu();
2737
2738	aws->at("gaps");
2739	aws->create_toggle(FA_AWAR_SHOW_GAPS_MESSAGES);
2740
2741	aws->at("continue");
2742	aws->create_toggle(FA_AWAR_CONTINUE_ON_ERROR);
2743
2744	aws->at("on_failure");
2745	aws->label("On failure");
2746	aws->create_option_menu(FA_AWAR_ACTION_ON_ERROR, true);
2747	aws->insert_default_option("do nothing", "", FA_NO_ACTION);
2748	aws->insert_option ("mark failed", "", FA_MARK_FAILED);
2749	aws->insert_option ("mark aligned", "", FA_MARK_ALIGNED);
2750	aws->update_option_menu();
2751
2752	aws->at("align");
2753	aws->callback(FastAligner_start, (AW_CL)data_access);
2754	aws->highlight();
2755	aws->create_button("GO", "GO", "G");
2756
2757	return aws;
2758	}
2759
2760	// --------------------------------------------------------------------------------
2761
2762	#ifdef UNIT_TESTS
2763
2764	#include <test_unit.h>
2765
2766	// ---------------------
2767	// OligoCounter
2768
2769	#include <map>
2770	#include <string>
2771
2772	using std::map;
2773	using std::string;
2774
2775	typedef map<string, size_t> OligoCount;
2776
2777	class OligoCounter {
2778	size_t oligo_len;
2779	size_t datasize;
2780
2781	mutable OligoCount occurance;
2782
2783	static string removeGaps(const char *seq) {
2784	size_t len = strlen(seq);
2785	string nogaps;
2786	nogaps.reserve(len);
2787
2788	for (size_t p = 0; p<len; ++p) {
2789	char c = seq[p];
2790	if (c != '-' && c != '.') nogaps.append(1, c);
2791	}
2792	return nogaps;
2793	}
2794
2795	void count_oligos(const string& seq) {
2796	occurance.clear();
2797	size_t max_pos = seq.length()-oligo_len;
2798	for (size_t p = 0; p <= max_pos; ++p) {
2799	string oligo(seq, p, oligo_len);
2800	occurance[oligo]++;
2801	}
2802	}
2803
2804	public:
2805	OligoCounter()
2806	: oligo_len(0),
2807	datasize(0)
2808	{}
2809	OligoCounter(const char *seq, size_t oligo_len_)
2810	: oligo_len(oligo_len_)
2811	{
2812	string seq_nogaps = removeGaps(seq);
2813	datasize = seq_nogaps.length();
2814	count_oligos(seq_nogaps);
2815	}
2816
2817	OligoCounter(const OligoCounter& other)
2818	: oligo_len(other.oligo_len),
2819	datasize(other.datasize),
2820	occurance(other.occurance)
2821	{}
2822
2823	size_t oligo_count(const char *oligo) {
2824	fa_assert(strlen(oligo) == oligo_len);
2825	return occurance[oligo];
2826	}
2827
2828	size_t similarity_score(const OligoCounter& other) const {
2829	size_t score = 0;
2830	if (oligo_len == other.oligo_len) {
2831	for (OligoCount::const_iterator o = occurance.begin(); o != occurance.end(); ++o) {
2832	const string& oligo = o->first;
2833	size_t count = o->second;
2834
2835	score += min(count, other.occurance[oligo]);
2836	}
2837	}
2838	return score;
2839	}
2840
2841	size_t getDataSize() const { return datasize; }
2842	};
2843
2844	void TEST_OligoCounter() {
2845	OligoCounter oc1("CCAGGT", 3);
2846	OligoCounter oc2("GGTCCA", 3);
2847	OligoCounter oc2_gaps("..GGT--CCA..", 3);
2848	OligoCounter oc3("AGGTCC", 3);
2849	OligoCounter oc4("AGGTCCAGG", 3);
2850
2851	TEST_EXPECT_EQUAL(oc1.oligo_count("CCA"), 1);
2852	TEST_EXPECT_ZERO(oc1.oligo_count("CCG"));
2853	TEST_EXPECT_EQUAL(oc4.oligo_count("AGG"), 2);
2854
2855	int sc1_2 = oc1.similarity_score(oc2);
2856	int sc2_1 = oc2.similarity_score(oc1);
2857	TEST_EXPECT_EQUAL(sc1_2, sc2_1);
2858
2859	int sc1_2gaps = oc1.similarity_score(oc2_gaps);
2860	TEST_EXPECT_EQUAL(sc1_2, sc1_2gaps);
2861
2862	int sc1_3 = oc1.similarity_score(oc3);
2863	int sc2_3 = oc2.similarity_score(oc3);
2864	int sc3_4 = oc3.similarity_score(oc4);
2865
2866	TEST_EXPECT_EQUAL(sc1_2, 2); // common oligos (CCA GGT)
2867	TEST_EXPECT_EQUAL(sc1_3, 2); // common oligos (AGG GGT)
2868	TEST_EXPECT_EQUAL(sc2_3, 3); // common oligos (GGT GTC TCC)
2869
2870	TEST_EXPECT_EQUAL(sc3_4, 4);
2871	}
2872
2873	// -------------------------
2874	// FakeFamilyFinder
2875
2876	class FakeFamilyFinder: public FamilyFinder { // derived from a Noncopyable
2877	// used by unit tests to detect next relatives instead of asking the pt-server
2878
2879	GBDATA *gb_main;
2880	string ali_name;
2881	map<string, OligoCounter> oligos_counted; // key = species name
2882	PosRange counted_for_range;
2883	size_t oligo_len;
2884
2885	public:
2886	FakeFamilyFinder(GBDATA *gb_main_, string ali_name_, bool rel_matches_, size_t oligo_len_)
2887	: FamilyFinder(rel_matches_, RSS_BOTH_MIN),
2888	gb_main(gb_main_),
2889	ali_name(ali_name_),
2890	counted_for_range(PosRange::whole()),
2891	oligo_len(oligo_len_)
2892	{}
2893
2894	GB_ERROR searchFamily(const char *sequence, FF_complement compl_mode, int max_results, double min_score) OVERRIDE { // @@@ use min_score
2895	// 'sequence' has to contain full sequence or part corresponding to 'range'
2896
2897	TEST_EXPECT_EQUAL(compl_mode, FF_FORWARD); // not fit for other modes
2898
2899	delete_family_list();
2900
2901	OligoCounter seq_oligo_count(sequence, oligo_len);
2902
2903	if (range != counted_for_range) {
2904	oligos_counted.clear(); // forget results for different range
2905	counted_for_range = range;
2906	}
2907
2908	char *buffer = 0;
2909	int buffersize = 0;
2910
2911	bool partial_match = range.is_part();
2912
2913	GB_transaction ta(gb_main);
2914	int results = 0;
2915
2916	for (GBDATA *gb_species = GBT_first_species(gb_main);
2917	gb_species && results<max_results;
2918	gb_species = GBT_next_species(gb_species))
2919	{
2920	string name = GBT_get_name(gb_species);
2921	if (oligos_counted.find(name) == oligos_counted.end()) {
2922	GBDATA *gb_data = GBT_find_sequence(gb_species, ali_name.c_str());
2923	const char *spec_seq = GB_read_char_pntr(gb_data);
2924
2925	if (partial_match) {
2926	int spec_seq_len = GB_read_count(gb_data);
2927	int range_len = ExplicitRange(range, spec_seq_len).size();
2928
2929	if (buffersize<range_len) {
2930	delete [] buffer;
2931	buffersize = range_len;
2932	buffer = new char[buffersize+1];
2933	}
2934
2935	range.copy_corresponding_part(buffer, spec_seq, spec_seq_len);
2936	oligos_counted[name] = OligoCounter(buffer, oligo_len);
2937	}
2938	else {
2939	oligos_counted[name] = OligoCounter(spec_seq, oligo_len);
2940	}
2941	}
2942
2943	const OligoCounter& spec_oligo_count = oligos_counted[name];
2944	size_t score = seq_oligo_count.similarity_score(spec_oligo_count);
2945
2946	FamilyList *newMember = new FamilyList;
2947
2948	newMember->name = strdup(name.c_str());
2949	newMember->matches = score;
2950	newMember->rel_matches = score/spec_oligo_count.getDataSize();
2951	newMember->next = NULL;
2952
2953	family_list = newMember->insertSortedBy_matches(family_list);
2954	results++;
2955	}
2956
2957	delete [] buffer;
2958
2959	return NULL;
2960	}
2961	};
2962
2963	// ----------------------------
2964	// test_alignment_data
2965
2966	#include <arb_unit_test.h>
2967
2968	static const char *test_aliname = "ali_test";
2969
2970	static const char get_aligned_data_of(GBDATA gb_main, const char *species_name) {
2971	GB_transaction ta(gb_main);
2972	ARB_ERROR error;
2973	const char *data = NULL;
2974
2975	GBDATA *gb_species = GBT_find_species(gb_main, species_name);
2976	if (!gb_species) error = GB_await_error();
2977	else {
2978	GBDATA *gb_data = GBT_find_sequence(gb_species, test_aliname);
2979	if (!gb_data) error = GB_await_error();
2980	else {
2981	data = GB_read_char_pntr(gb_data);
2982	if (!data) error = GB_await_error();
2983	}
2984	}
2985
2986	TEST_EXPECT_NULL(error.deliver());
2987
2988	return data;
2989	}
2990
2991	static const char get_used_rels_for(GBDATA gb_main, const char *species_name) {
2992	GB_transaction ta(gb_main);
2993	const char *result = NULL;
2994	GBDATA *gb_species = GBT_find_species(gb_main, species_name);
2995	if (!gb_species) result = GBS_global_string("<No such species '%s'>", species_name);
2996	else {
2997	GBDATA *gb_used_rels = GB_search(gb_species, "used_rels", GB_FIND);
2998	if (!gb_used_rels) result = "<No such field 'used_rels'>";
2999	else result = GB_read_char_pntr(gb_used_rels);
3000	}
3001	return result;
3002	}
3003
3004	static GB_ERROR forget_used_relatives(GBDATA *gb_main) {
3005	GB_ERROR error = NULL;
3006	GB_transaction ta(gb_main);
3007	for (GBDATA *gb_species = GBT_first_species(gb_main);
3008	gb_species && !error;
3009	gb_species = GBT_next_species(gb_species))
3010	{
3011	GBDATA *gb_used_rels = GB_search(gb_species, "used_rels", GB_FIND);
3012	if (gb_used_rels) {
3013	error = GB_delete(gb_used_rels);
3014	}
3015	}
3016	return error;
3017	}
3018
3019
3020	#define ALIGNED_DATA_OF(name) get_aligned_data_of(gb_main, name)
3021	#define USED_RELS_FOR(name) get_used_rels_for(gb_main, name)
3022
3023	// ----------------------------------------
3024
3025	static GBDATA *selection_fake_gb_main = NULL;
3026	static GBDATA *selection_fake_gb_last = NULL;
3027
3028	static GBDATA fake_first_selected(int count) {
3029	selection_fake_gb_last = NULL;
3030	*count = 2; // we fake two species as selected ("c1" and "c2")
3031	return GBT_find_species(selection_fake_gb_main, "c1");
3032	}
3033	static GBDATA *fake_next_selected() {
3034	if (!selection_fake_gb_last) {
3035	selection_fake_gb_last = GBT_find_species(selection_fake_gb_main, "c2");
3036	}
3037	else {
3038	selection_fake_gb_last = NULL;
3039	}
3040	return selection_fake_gb_last;
3041	}
3042
3043	static char fake_get_consensus(const char, PosRange range) {
3044	const char *data = get_aligned_data_of(selection_fake_gb_main, "s1");
3045	if (range.is_whole()) return strdup(data);
3046	return GB_strpartdup(data+range.start(), data+range.end());
3047	}
3048
3049	static void test_install_fakes(GBDATA *gb_main) {
3050	selection_fake_gb_main = gb_main;
3051	}
3052
3053
3054	// ----------------------------------------
3055
3056	static AlignParams test_ali_params = {
3057	FA_NO_REPORT,
3058	false, // showGapsMessages
3059	PosRange::whole()
3060	};
3061
3062	static AlignParams test_ali_params_partial = {
3063	FA_NO_REPORT,
3064	false, // showGapsMessages
3065	PosRange(25, 60)
3066	};
3067
3068	// ----------------------------------------
3069
3070	static struct arb_unit_test::test_alignment_data TestAlignmentData_TargetAndReferenceHandling[] = {
3071	{ 0, "s1", ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........." }, // reference
3072	{ 0, "s2", "AUCUCCUAAACCCAACCGUAGUUCGAAUUGAGGACUGUAACUC......................................................" }, // align single sequence (same data as reference)
3073	{ 1, "m1", "UAGAGGAUUUGGGUUGGCAUCAAGCUUAACUCCUGACAUUGAG......................................................" }, // align marked sequences.. (complement of reference)
3074	{ 1, "m2", "...UCCUAAACCAACCCGUAGUUCGAAUUGAGGACUGUAA........................................................." },
3075	{ 1, "m3", "AUC---UAAACCAACCCGUAGUUCGAAUUGAGGACUG---CUC......................................................" },
3076	{ 0, "c1", "AUCUCCUAAACCCAACC--------AAUUGAGGACUGUAACUC......................................................" }, // align selected sequences..
3077	{ 0, "c2", "AUCUCCU------AACCGUAGUUCCCCGAA------ACUGUAACUC..................................................." },
3078	{ 0, "r1", "GAGUUACAGUCCUCAAUUCGGGGAACUACGGUUGGGUUUAGGAGAU..................................................." }, // align by faked pt_server
3079	};
3080
3081	void TEST_Aligner_TargetAndReferenceHandling() {
3082	// performs some alignments to check logic of target and reference handling
3083
3084	GB_shell shell;
3085	ARB_ERROR error;
3086	GBDATA *gb_main = TEST_CREATE_DB(error, test_aliname, TestAlignmentData_TargetAndReferenceHandling, false);
3087
3088	TEST_EXPECT_NULL(error.deliver());
3089
3090	SearchRelativeParams search_relative_params(new FakeFamilyFinder(gb_main, test_aliname, false, 8),
3091	test_aliname,
3092	2);
3093
3094	test_install_fakes(gb_main);
3095	arb_suppress_progress silence;
3096
3097	// bool cont_on_err = true;
3098	bool cont_on_err = false;
3099
3100	TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 3); // we got 3 marked species
3101	{
3102	Aligner aligner(gb_main,
3103	FA_CURRENT,
3104	test_aliname,
3105	"s2", // toalign
3106	NULL, // get_first_selected_species
3107	NULL, // get_next_selected_species
3108	"s1", // reference species
3109	NULL, // get_consensus
3110	search_relative_params, // relative search
3111	FA_TURN_ALWAYS,
3112	test_ali_params,
3113	0,
3114	cont_on_err,
3115	FA_NO_ACTION);
3116	error = aligner.run();
3117	TEST_EXPECT_NULL(error.deliver());
3118	}
3119	TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 3); // we still got 3 marked species
3120	{
3121	Aligner aligner(gb_main,
3122	FA_MARKED,
3123	test_aliname,
3124	NULL, // toalign
3125	NULL, // get_first_selected_species
3126	NULL, // get_next_selected_species
3127	"s1", // reference species
3128	NULL, // get_consensus
3129	search_relative_params, // relative search
3130	FA_TURN_ALWAYS,
3131	test_ali_params,
3132	0,
3133	cont_on_err,
3134	FA_MARK_FAILED);
3135	error = aligner.run();
3136	TEST_EXPECT_NULL(error.deliver());
3137
3138	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 0); // FA_MARK_FAILED (none failed -> none marked)
3139	}
3140	{
3141	Aligner aligner(gb_main,
3142	FA_SELECTED,
3143	test_aliname,
3144	NULL, // toalign
3145	fake_first_selected, // get_first_selected_species
3146	fake_next_selected, // get_next_selected_species
3147	NULL, // reference species
3148	fake_get_consensus, // get_consensus
3149	search_relative_params, // relative search
3150	FA_TURN_ALWAYS,
3151	test_ali_params,
3152	0,
3153	cont_on_err,
3154	FA_MARK_ALIGNED);
3155	error = aligner.run();
3156	TEST_EXPECT_NULL(error.deliver());
3157
3158	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 2); // FA_MARK_ALIGNED (2 selected were aligned)
3159	}
3160	{
3161	Aligner aligner(gb_main,
3162	FA_CURRENT,
3163	test_aliname,
3164	"r1", // toalign
3165	NULL, // get_first_selected_species
3166	NULL, // get_next_selected_species
3167	NULL, // reference species
3168	NULL, // get_consensus
3169	search_relative_params, // relative search
3170	FA_TURN_ALWAYS,
3171	test_ali_params,
3172	0,
3173	cont_on_err,
3174	FA_MARK_ALIGNED);
3175
3176	error = aligner.run();
3177	TEST_EXPECT_NULL(error.deliver());
3178
3179	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 1);
3180	}
3181
3182	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s2"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C...........");
3183
3184	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m1"), ".......UAG--AGG-A------U-U-UGGGU-UG-G-C-A-U-CAA-GCU--------UAA-C-UCCUG-AC--A-UUGAG...............");
3185	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m2"), "..............U-C------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA................");
3186	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m3"), ".........A--U----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----CU-C...........");
3187
3188	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-------------------AA-U-UGAGG-AC--U-GUAA-CU-C...........");
3189	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c2"), ".........A--UCU-C------C-U-AA---------C-C-G-UAG-UUC------------C-CCGAA-AC--U-GUAA-CU-C...........");
3190
3191	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("r1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUCCCC-----GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // here sequence shall be turned!
3192
3193
3194	TEST_EXPECT_EQUAL(USED_RELS_FOR("r1"), "s2:43, s1:1");
3195
3196	// ----------------------------------------------
3197	// now align all others vs next relative
3198
3199	search_relative_params.maxRelatives = 5;
3200	TEST_EXPECT_NO_ERROR(forget_used_relatives(gb_main));
3201
3202	int species_count = ARRAY_ELEMS(TestAlignmentData_TargetAndReferenceHandling);
3203	for (int sp = 0; sp<species_count; ++sp) {
3204	const char *name = TestAlignmentData_TargetAndReferenceHandling[sp].name;
3205	if (strcmp(name, "r1") != 0) { // skip "r1" (already done above)
3206	Aligner aligner(gb_main,
3207	FA_CURRENT,
3208	test_aliname,
3209	name, // toalign
3210	NULL, // get_first_selected_species
3211	NULL, // get_next_selected_species
3212	NULL, // reference species
3213	NULL, // get_consensus
3214	search_relative_params, // relative search
3215	FA_TURN_ALWAYS,
3216	test_ali_params,
3217	0,
3218	cont_on_err,
3219	FA_MARK_ALIGNED);
3220
3221	error = aligner.run();
3222	TEST_EXPECT_NULL(error.deliver());
3223
3224	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 1);
3225	}
3226	}
3227
3228	TEST_EXPECT_EQUAL(USED_RELS_FOR("s1"), "s2");
3229	TEST_EXPECT_EQUAL(USED_RELS_FOR("s2"), "s1"); // same as done manually
3230	TEST_EXPECT_EQUAL(USED_RELS_FOR("m1"), "r1:42");
3231	TEST_EXPECT_EQUAL(USED_RELS_FOR("m2"), "m3");
3232	TEST_EXPECT_EQUAL(USED_RELS_FOR("m3"), "m2");
3233	TEST_EXPECT_EQUAL(USED_RELS_FOR("c1"), "r1");
3234	TEST_EXPECT_EQUAL(USED_RELS_FOR("c2"), "r1");
3235
3236	// range aligned below (see test_ali_params_partial)
3237	// "-------------------------XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX------------------------------------"
3238	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // 1st aligning of 's1'
3239	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s2"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above (again aligned vs 's1')
3240
3241	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m1"), ".........U--AGA-G------G---AUUUG-GG-U-U-G-G-CAU-CAAGCU-----UAA-C-UCCUG-AC--A-UUGAG---------------"); // changed; @@@ bug: no dots at end
3242
3243	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m2"), ".........U--C----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----UA-A..........."); // changed (1st align vs 's1', this align vs 'm3')
3244	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m3"), ".........A--U----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----CU-C..........."); // same_as_above (1st align vs 's1', this align vs 'm2')
3245	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-------------------AA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3246	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c2"), ".........A--UCU-C------C-U--------A-A-C-C-G-UAG-UUCCCC-----GA--------A-AC--U-GUAA-CU-C..........."); // changed
3247
3248	// --------------------------------------
3249	// test partial relative search
3250
3251	search_relative_params.getFamilyFinder()->restrict_2_region(test_ali_params_partial.range);
3252	TEST_EXPECT_NO_ERROR(forget_used_relatives(gb_main));
3253
3254	for (int sp = 0; sp<species_count; ++sp) {
3255	const char *name = TestAlignmentData_TargetAndReferenceHandling[sp].name;
3256	Aligner aligner(gb_main,
3257	FA_CURRENT,
3258	test_aliname,
3259	name, // toalign
3260	NULL, // get_first_selected_species
3261	NULL, // get_next_selected_species
3262	NULL, // reference species
3263	NULL, // get_consensus
3264	search_relative_params, // relative search
3265	FA_TURN_NEVER,
3266	test_ali_params_partial,
3267	0,
3268	cont_on_err,
3269	FA_MARK_ALIGNED);
3270
3271	error = aligner.run();
3272	TEST_EXPECT_NULL(error.deliver());
3273
3274	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 1);
3275	}
3276
3277	TEST_EXPECT_EQUAL(USED_RELS_FOR("s1"), "s2");
3278	TEST_EXPECT_EQUAL(USED_RELS_FOR("s2"), "s1");
3279	TEST_EXPECT_EQUAL(USED_RELS_FOR("m1"), "r1"); // (not really differs)
3280	TEST_EXPECT_EQUAL(USED_RELS_FOR("m2"), "m3");
3281	TEST_EXPECT_EQUAL(USED_RELS_FOR("m3"), "m2");
3282	TEST_EXPECT_EQUAL(USED_RELS_FOR("c1"), "r1");
3283	TEST_EXPECT_EQUAL(USED_RELS_FOR("c2"), "r1");
3284	TEST_EXPECT_EQUAL(USED_RELS_FOR("r1"), "s2:20, c2:3");
3285
3286	// aligned range (see test_ali_params_partial)
3287	// "-------------------------XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX------------------------------------"
3288	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3289	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s2"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3290
3291	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m1"), ".........U--AGA-G------G-A-UU-UG-GG-U-U-G-G-CAU-CAAGCU-----UAA-C-UCCUG-AC--A-UUGAG---------------"); // changed
3292	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m2"), ".........U--C----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----UA-A..........."); // same_as_above
3293	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m3"), ".........A--U----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----CU-C..........."); // same_as_above
3294
3295	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-------------------AA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3296	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c2"), ".........A--UCU-C------C---------UA-A-C-C-G-UAG-UUCCCC-----GA--------A-AC--U-GUAA-CU-C..........."); // changed
3297
3298	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("r1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUCCCC-----GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3299
3300	GB_close(gb_main);
3301	}
3302
3303	// ----------------------------------------
3304
3305	static struct arb_unit_test::test_alignment_data TestAlignmentData_checksumError[] = {
3306	{ 0, "MtnK1722", "...G-GGC-C-G............CCC-GG--------CAAUGGGGGCGGCCCGGCGGAC----GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCC....." }, // gets aligned
3307	{ 0, "MhnFormi", "...A-CGA-U-C------------CUUCGG--------GGUCG-U-GG-C-GU-A--C------GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCCU...." }, // 1st relative
3308	{ 0, "MhnT1916", "...A-CGA-A-C------------CUU-GU--------GUUCG-U-GG-C-GA-A--C------GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCCU...." }, // next relative
3309	{ 0, "MthVanni", "...U-GGU-U-U------------C-------------GGCCA-U-GG-C-GG-A--C------GG--C-UCAUU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCC....." }, // next relative
3310	{ 0, "ThcCeler", "...G-GGG-C-G...CC-U---U--------GC--G--CGCAC-C-GG-C-GG-A--C------GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUA-CGGC-UCGAUCACCUCCU...." }, // next relative
3311	};
3312
3313	void TEST_SLOW_Aligner_checksumError() {
3314	// @@@ SLOW cause this often gets terminated in nightly builds
3315	// no idea why. normally it runs a few ms
3316
3317	// this produces an internal aligner error
3318
3319	GB_shell shell;
3320	ARB_ERROR error;
3321	GBDATA *gb_main = TEST_CREATE_DB(error, test_aliname, TestAlignmentData_checksumError, false);
3322
3323	SearchRelativeParams search_relative_params(new FakeFamilyFinder(gb_main, test_aliname, false, 8),
3324	test_aliname,
3325	10); // use up to 10 relatives
3326
3327	test_install_fakes(gb_main);
3328	arb_suppress_progress silence;
3329
3330	bool cont_on_err = true;
3331	if (!error) {
3332	Aligner aligner(gb_main,
3333	FA_CURRENT,
3334	test_aliname,
3335	"MtnK1722", // toalign
3336	NULL, // get_first_selected_species
3337	NULL, // get_next_selected_species
3338	NULL, // reference species
3339	NULL, // get_consensus
3340	search_relative_params, // relative search
3341	FA_TURN_ALWAYS,
3342	test_ali_params,
3343	0,
3344	cont_on_err,
3345	FA_MARK_ALIGNED);
3346
3347	error = aligner.run();
3348	}
3349	{
3350	GB_ERROR err = error.deliver();
3351	TEST_EXPECT_NULL__BROKEN(err, "Aligner produced 1 error");
3352	}
3353	TEST_EXPECT_EQUAL__BROKEN(USED_RELS_FOR("MtnK1722"), "???", "<No such field 'used_rels'>"); // subsequent failure
3354
3355	GB_close(gb_main);
3356	}
3357
3358	static const char *asstr(LooseBases& ub) {
3359	LooseBases tmp;
3360	while (!ub.is_empty()) tmp.memorize(ub.recall());
3361
3362	const char *result = "";
3363	while (!tmp.is_empty()) {
3364	ExplicitRange r = tmp.recall();
3365	result = GBS_global_string("%s %i/%i", result, r.start(), r.end());
3366	ub.memorize(r);
3367	}
3368	return result;
3369	}
3370
3371	void TEST_BASIC_UnalignedBases() {
3372	LooseBases ub;
3373	TEST_EXPECT(ub.is_empty());
3374	TEST_EXPECT_EQUAL(asstr(ub), "");
3375
3376	// test add+remove
3377	ub.memorize(ExplicitRange(5, 7));
3378	TEST_REJECT(ub.is_empty());
3379	TEST_EXPECT_EQUAL(asstr(ub), " 5/7");
3380
3381	TEST_EXPECT(ub.recall() == ExplicitRange(5, 7));
3382	TEST_EXPECT(ub.is_empty());
3383
3384	ub.memorize(ExplicitRange(2, 4));
3385	TEST_EXPECT_EQUAL(asstr(ub), " 2/4");
3386
3387	ub.memorize(ExplicitRange(4, 9));
3388	TEST_EXPECT_EQUAL(asstr(ub), " 2/4 4/9");
3389
3390	ub.memorize(ExplicitRange(8, 10));
3391	ub.memorize(ExplicitRange(11, 14));
3392	ub.memorize(ExplicitRange(12, 17));
3393	TEST_EXPECT_EQUAL(asstr(ub), " 2/4 4/9 8/10 11/14 12/17");
3394	TEST_EXPECT_EQUAL(asstr(ub), " 2/4 4/9 8/10 11/14 12/17"); // check asstr has no side-effect
3395
3396	{
3397	LooseBases toaddNrecalc;
3398
3399	CompactedSubSequence Old("ACGTACGT", 8, "name1");
3400	CompactedSubSequence New("--A-C--G-T--A-C-G-T", 19, "name2");
3401	// ---------------------- 0123456789012345678
3402
3403	toaddNrecalc.memorize(ExplicitRange(1, 7));
3404	toaddNrecalc.memorize(ExplicitRange(3, 5));
3405	TEST_EXPECT_EQUAL(asstr(toaddNrecalc), " 1/7 3/5");
3406
3407	ub.follow_ali_change_and_append(toaddNrecalc, AliChange(Old, New));
3408
3409	TEST_EXPECT_EQUAL(asstr(ub), " 3/18 8/15 2/4 4/9 8/10 11/14 12/17");
3410	TEST_EXPECT(toaddNrecalc.is_empty());
3411
3412	LooseBases selfRecalc;
3413	selfRecalc.follow_ali_change_and_append(ub, AliChange(New, New));
3414	TEST_EXPECT_EQUAL__BROKEN(asstr(selfRecalc),
3415	" 3/18 8/15 0/6 3/11 8/11 10/15 10/17", // wanted behavior?
3416	" 3/18 8/17 0/6 3/11 8/13 10/15 10/18"); // doc wrong behavior @@@ "8/17", "8/13", "10/18" are wrong
3417
3418	ub.follow_ali_change_and_append(selfRecalc, AliChange(New, Old));
3419	TEST_EXPECT_EQUAL__BROKEN(asstr(ub),
3420	" 1/7 3/5 0/1 1/3 3/3 4/5 4/6", // wanted behavior? (from wanted behavior above)
3421	" 1/7 3/7 0/2 1/4 3/5 4/6 4/7"); // document wrong behavior
3422	TEST_EXPECT_EQUAL__BROKEN(asstr(ub),
3423	" 1/7 3/6 0/1 1/3 3/4 4/5 4/7", // wanted behavior? (from wrong result above)
3424	" 1/7 3/7 0/2 1/4 3/5 4/6 4/7"); // document wrong behavior
3425	}
3426	}
3427
3428
3429	#endif // UNIT_TESTS
3430
3431

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

Download in other formats: