Context Navigation

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

Visit:

Last change on this file was 6143, checked in by westram, 16 years ago
backport [6141] (parts changing code, but only strings and comments)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 106.2 KB

Line
1	// =============================================================== //
2	// //
3	// File : AWTC_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	#include <island_hopping.h>
13
14	#include <awtc_next_neighbours.hxx>
15	#include <awtc_seq_search.hxx>
16	#include "ClustalV.hxx"
17	#include "fast_aligner.hxx"
18
19	#include <arbdbt.h>
20	#include <aw_window.hxx>
21	#include <aw_awars.hxx>
22	#include <awt.hxx>
23	#include <awt_sel_boxes.hxx>
24
25	#include <cstdio>
26	#include <cstdlib>
27	#include <cctype>
28	#include <cmath>
29	#include <cstring>
30	#include <climits>
31
32	// --------------------------------------------------------------------------------
33
34	enum FA_report {
35	FA_NO_REPORT, // no report
36	FA_TEMP_REPORT, // report to temporary entries
37	FA_REPORT, // report to resident entries
38	};
39
40	enum FA_range {
41	FA_WHOLE_SEQUENCE, // align whole sequence
42	FA_AROUND_CURSOR, // align xxx positions around current cursor position
43	FA_SELECTED_RANGE, // align selected range
44	};
45
46	enum FA_turn {
47	FA_TURN_NEVER, // never try to turn sequence
48	FA_TURN_INTERACTIVE, // try to turn, but query user
49	FA_TURN_ALWAYS, // turn if score is better
50	};
51
52	enum FA_reference {
53	FA_REF_EXPLICIT, // reference sequence explicitely specified
54	FA_REF_CONSENSUS, // use group consensus as reference
55	FA_REF_RELATIVES, // search next relatives by PT server
56	};
57
58	enum FA_alignTarget {
59	FA_CURRENT, // align current species
60	FA_MARKED, // align all marked species
61	FA_SELECTED, // align selected species (= range)
62	};
63
64	struct AlignParams {
65	// int temporary; // // ==1 -> create only temporary aligment report into alignment (2=resident,0=none)
66	FA_report report;
67	bool showGapsMessages; // display messages about missing gaps in master?
68	int firstColumn; // first column of range to be aligned (0..len-1)
69	int lastColumn; // last column of range to be aligned (0..len-1, -1 = (len-1))
70	};
71
72	struct SearchRelativeParams {
73	int pt_server_id; // pt_server to search for next relative
74	GB_CSTR pt_server_alignment; // alignment used in pt_server (may differ from 'alignment')
75	int maxRelatives; // max # of relatives to use
76	int fam_oligo_len; // oligo length
77	int fam_mismatches; // allowed mismatches
78	bool fam_fast_mode; // fast family find (serch only oligos starting with 'A')
79	bool fam_rel_matches; // sort results by relative matches
80	};
81
82	// --------------------------------------------------------------------------------
83
84	#define GAP_CHAR '-'
85	#define QUALITY_NAME "ASC_ALIGNER_CLIENT_SCORE"
86	#define INSERTS_NAME "AMI_ALIGNER_MASTER_INSERTS"
87
88	#define FA_AWAR_ROOT "faligner/"
89	#define FA_AWAR_TO_ALIGN FA_AWAR_ROOT "what"
90	#define FA_AWAR_REFERENCE FA_AWAR_ROOT "against"
91	#define FA_AWAR_REFERENCE_NAME FA_AWAR_ROOT "sagainst"
92	#define FA_AWAR_RANGE FA_AWAR_ROOT "range"
93	#define FA_AWAR_PROTECTION FA_AWAR_ROOT "protection"
94	#define FA_AWAR_AROUND FA_AWAR_ROOT "around"
95	#define FA_AWAR_MIRROR FA_AWAR_ROOT "mirror"
96	#define FA_AWAR_REPORT FA_AWAR_ROOT "report"
97	#define FA_AWAR_SHOW_GAPS_MESSAGES FA_AWAR_ROOT "show_gaps"
98	#define FA_AWAR_USE_SECONDARY FA_AWAR_ROOT "use_secondary"
99	#define FA_AWAR_NEXT_RELATIVES FA_AWAR_ROOT "next_relatives"
100	#define FA_AWAR_PT_SERVER_ALIGNMENT "tmp/" FA_AWAR_ROOT "relative_ali"
101
102	#define FA_AWAR_ISLAND_HOPPING_ROOT "island_hopping/"
103	#define FA_AWAR_USE_ISLAND_HOPPING FA_AWAR_ISLAND_HOPPING_ROOT "use"
104	#define FA_AWAR_ESTIMATE_BASE_FREQ FA_AWAR_ISLAND_HOPPING_ROOT "estimate_base_freq"
105	#define FA_AWAR_BASE_FREQ_A FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_a"
106	#define FA_AWAR_BASE_FREQ_C FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_c"
107	#define FA_AWAR_BASE_FREQ_G FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_g"
108	#define FA_AWAR_BASE_FREQ_T FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_t"
109	#define FA_AWAR_SUBST_PARA_AC FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ac"
110	#define FA_AWAR_SUBST_PARA_AG FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ag"
111	#define FA_AWAR_SUBST_PARA_AT FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_at"
112	#define FA_AWAR_SUBST_PARA_CG FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_cg"
113	#define FA_AWAR_SUBST_PARA_CT FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ct"
114	#define FA_AWAR_SUBST_PARA_GT FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_gt"
115	#define FA_AWAR_EXPECTED_DISTANCE FA_AWAR_ISLAND_HOPPING_ROOT "expected_dist"
116	#define FA_AWAR_STRUCTURE_SUPPLEMENT FA_AWAR_ISLAND_HOPPING_ROOT "struct_suppl"
117	#define FA_AWAR_THRESHOLD FA_AWAR_ISLAND_HOPPING_ROOT "threshold"
118	#define FA_AWAR_GAP_A FA_AWAR_ISLAND_HOPPING_ROOT "gapa"
119	#define FA_AWAR_GAP_B FA_AWAR_ISLAND_HOPPING_ROOT "gapb"
120	#define FA_AWAR_GAP_C FA_AWAR_ISLAND_HOPPING_ROOT "gapc"
121
122	// --------------------------------------------------------------------------------
123
124	extern GBDATA *GLOBAL_gb_main;
125
126	static IslandHopping *island_hopper = 0;
127	static GB_alignment_type global_alignmentType = GB_AT_UNKNOWN; // type of actually aligned sequence
128
129	static int currentSequenceNumber; // used for counter
130	static int overallSequenceNumber;
131
132	// --------------------------------------------------------------------------------
133
134	static inline GB_ERROR species_not_found(GB_CSTR species_name) {
135	return GB_export_errorf("No species '%s' found!", species_name);
136	}
137
138	static GB_ERROR reverseComplement(GBDATA *gb_species, GB_CSTR ali, int max_protection) {
139	GBDATA *gbd = GBT_read_sequence(gb_species, ali);
140	GB_ERROR error = 0;
141
142	if (!gbd) {
143	error = GB_export_errorf("No 'data' found for species '%s'", GBT_read_name(gb_species));
144	}
145	else {
146	int my_protection = GB_read_security_write(gbd);
147
148	if (my_protection<=max_protection) { // ok
149	char *seq = GB_read_string(gbd);
150	int length = GB_read_string_count(gbd);
151	GB_alignment_type ali_type = GBT_get_alignment_type(GLOBAL_gb_main, ali);
152
153	char T_or_U;
154	error = GBT_determine_T_or_U(ali_type, &T_or_U, "reverse-complement");
155	if (!error) {
156	GBT_reverseComplementNucSequence(seq, length, T_or_U);
157	error = GB_write_string(gbd, seq);
158	}
159	}
160	else { // protection error
161	error = GB_export_errorf("Cannot reverse-complement species '%s' because of protection level", GBT_read_name(gb_species));
162	}
163
164	}
165
166	return error;
167	}
168
169	void AWTC_build_reverse_complement(AW_window *aw, AW_CL cd2) {
170	GB_push_transaction(GLOBAL_gb_main);
171
172	AW_root *root = aw->get_root();
173	FA_alignTarget revComplWhat = static_cast<FA_alignTarget>(root->awar(FA_AWAR_TO_ALIGN)->read_int());
174	char *default_alignment = GBT_get_default_alignment(GLOBAL_gb_main);
175	GB_CSTR alignment = root->awar_string(AWAR_EDITOR_ALIGNMENT, default_alignment)->read_string();
176	GB_ERROR error = 0;
177	int max_protection = root->awar(FA_AWAR_PROTECTION)->read_int();
178
179	switch (revComplWhat) {
180	case FA_CURRENT: { // current species
181	GB_CSTR species_name = root->awar(AWAR_SPECIES_NAME)->read_string();
182	GBDATA *gb_species = GBT_find_species(GLOBAL_gb_main, species_name);
183	if (!gb_species) error = species_not_found(species_name);
184	if (!error) error = reverseComplement(gb_species, alignment, max_protection);
185	break;
186	}
187	case FA_MARKED: { // marked species
188	GBDATA *gb_species = GBT_first_marked_species(GLOBAL_gb_main);
189
190	if (!gb_species) {
191	error = GB_export_error("There is no marked species");
192	}
193	while (gb_species) {
194	error = reverseComplement(gb_species, alignment, max_protection);
195	if (error) break;
196	gb_species = GBT_next_marked_species(gb_species);
197	}
198	break;
199	}
200	case FA_SELECTED: { // selected species
201	static struct AWTC_faligner_cd cd = (struct AWTC_faligner_cd )cd2;
202	AWTC_get_first_selected_species get_first_selected_species = cd->get_first_selected_species;
203	AWTC_get_next_selected_species get_next_selected_species = cd->get_next_selected_species;
204	int count = 0;
205	GBDATA *gb_species = get_first_selected_species(&count);
206
207	if (!gb_species) {
208	error = GB_export_error("There is no selected species!");
209	}
210	while (gb_species) {
211	error = reverseComplement(gb_species, alignment, max_protection);
212	if (error) break;
213	gb_species = get_next_selected_species();
214	}
215	break;
216	}
217	default: {
218	awtc_assert(0);
219	break;
220	}
221	}
222
223	GB_end_transaction_show_error(GLOBAL_gb_main, error, aw_message);
224	}
225
226	// --------------------------------------------------------------------------------
227
228	class UnalignedBasesList;
229	class UnalignedBases
230	{
231	int start, // absolute positions
232	end;
233	UnalignedBases *next;
234
235	friend class UnalignedBasesList;
236
237	public:
238	UnalignedBases(int pos1, int posn) { start = pos1; end = posn; next = 0; }
239	~UnalignedBases() { delete next; }
240
241	int get_start() const { return start; }
242	int get_end() const { return end; }
243	};
244
245	class UnalignedBasesList
246	{
247	UnalignedBases *head;
248	public:
249	UnalignedBasesList() { head = 0; }
250	~UnalignedBasesList() { delete head; }
251
252	void clear() { delete head; head = 0; }
253	int is_empty() const { return head==0; }
254
255	void memorize(int start, int end) {
256	#ifdef DEBUG
257	printf("memorize %i-%i\n", start, end);
258	#endif
259	awtc_assert(start<=end);
260	UnalignedBases *ub = new UnalignedBases(start, end);
261	ub->next = head;
262	head = ub;
263	}
264	void recall(int start, int end) {
265	UnalignedBases *ub = head;
266
267	if (start) *start = ub->start;
268	if (end) *end = ub->end;
269
270	#ifdef DEBUG
271	printf("recall %i-%i\n", ub->start, ub->end);
272	#endif
273
274	head = head->next;
275	ub->next = 0;
276	delete ub;
277	}
278
279	void add(UnalignedBasesList *ubl);
280	int add_and_recalc_positions(UnalignedBasesList ubl, AWTC_CompactedSubSequence oldSequence, AWTC_CompactedSubSequence *newSequence);
281	};
282
283	inline void UnalignedBasesList::add(UnalignedBasesList *ubl)
284	{
285	if (!ubl->is_empty()) {
286	if (is_empty()) {
287	head = ubl->head;
288	}
289	else {
290	UnalignedBases *last = head;
291
292	while (last->next) {
293	last = last->next;
294	}
295	last->next = ubl->head;
296	}
297	ubl->head = 0;
298	}
299	}
300
301	inline int UnalignedBasesList::add_and_recalc_positions(UnalignedBasesList ubl, AWTC_CompactedSubSequence oldSequence, AWTC_CompactedSubSequence *newSequence) {
302	// returns the number of unaligned bases
303	int bases = 0;
304
305	if (!ubl->is_empty()) {
306	UnalignedBases *toCorrect = ubl->head;
307
308	awtc_assert(oldSequence->length()==newSequence->length());
309
310	add(ubl);
311	while (toCorrect) {
312	int cs = oldSequence->compPosition(toCorrect->start); // compressed positions in oldSequence = compressed positions in newSequence
313	int ce = oldSequence->compPosition(toCorrect->end);
314
315	bases += ce-cs+1;
316
317	#ifdef DEBUG
318	printf("add_and_recalc_positions %i/%i -> ", toCorrect->start, toCorrect->end);
319	#endif
320
321	toCorrect->start = newSequence->expdPosition(cs) - newSequence->no_of_gaps_before(cs);
322	toCorrect->end = newSequence->expdPosition(ce) + newSequence->no_of_gaps_after(ce);
323
324	// if (cs>0) {
325	// toCorrect->start = newSequence->expdPosition(cs-1)+1; // first position after left neighbour (to get all gaps left of part)
326	// }
327	// else {
328	// toCorrect->start = newSequence->expdStartPosition(); //newSequence->expdPosition(cs);
329	// }
330
331	// if (ce<newSequence->length()) {
332	// toCorrect->end = newSequence->expdPosition(ce+1)-1; // last position before right neighbour (to get all gaps right of part)
333	// }
334	// else {
335	// toCorrect->end = newSequence->expdPosition(ce);
336	// }
337
338	#ifdef DEBUG
339	printf("%i/%i\n", toCorrect->start, toCorrect->end);
340	#endif
341
342	toCorrect = toCorrect->next;
343	}
344	}
345
346	return bases;
347	}
348
349
350	static UnalignedBasesList unaligned_bases; // if fast_align cannot align (no master bases) -> stores positions here
351
352	ATTRIBUTED(__ATTR__USERESULT, static inline GB_ERROR AWTC_memerr()) {
353	return "out of memory";
354	}
355
356	static const char AWTC_read_name(GBDATA gbd) {
357	return gbd ? GBT_read_name(gbd) : "GROUP-CONSENSUS";
358	}
359
360	static inline int relatedBases(char base1, char base2)
361	{
362	return AWTC_baseMatch(base1, base2)==1;
363	}
364
365	static inline char alignQuality(char slave, char master)
366	{
367	awtc_assert(slave);
368	awtc_assert(master);
369	char result = '#';
370	if (slave==master) result = '-'; // equal
371	else if (slave==GAP_CHAR) result = '+'; // inserted gap
372	else if (master==GAP_CHAR) result = '+'; // no gap in master
373	else if (relatedBases(slave,master)) result = '~'; // mutation (related bases)
374	return result; // mutation (non-related bases)
375	}
376
377	// --------------------------------------------------------------------------------
378	// Debugging stuff
379	// --------------------------------------------------------------------------------
380
381	#ifdef DEBUG
382	static char lstr(const char s, int len)
383	{
384	static int alloc;
385	static char *buffer;
386
387	if (alloc<(len+1)) {
388	if (alloc) free(buffer);
389	buffer = (char*)malloc(alloc=len+100);
390	}
391
392	memcpy(buffer, s, len);
393	buffer[len] = 0;
394
395	return buffer;
396	}
397
398	#define BUFLEN 120
399
400	static inline char compareChar(char base1, char base2)
401	{
402	return base1==base2 ? '=' : (relatedBases(base1,base2) ? 'x' : 'X');
403	}
404	static void dump_n_compare_one(const char seq1, const char seq2, long len, long offset)
405	{
406	awtc_assert(len<=BUFLEN);
407	char compare[BUFLEN+1];
408
409	for (long l=0; l<len; l++)
410	compare[l] = AWTC_is_gap(seq1[l]) \|\| AWTC_is_gap(seq2[l]) ? ' ' : compareChar(seq1[l],seq2[l]);
411
412	compare[len] = 0;
413
414	printf(" %li '%s'\n", offset, lstr(seq1,len));
415	printf(" %li '%s'\n", offset, lstr(seq2,len));
416	printf(" %li '%s'\n", offset, compare);
417	}
418
419	static inline void dump_rest(const char *seq, long len, int idx, long offset)
420	{
421	printf(" Rest von Sequenz %i:\n", idx);
422	while (len>BUFLEN)
423	{
424	printf(" %li '%s'\n", offset, lstr(seq, BUFLEN));
425	seq += BUFLEN;
426	len -= BUFLEN;
427	offset += BUFLEN;
428	}
429
430	awtc_assert(len>0);
431	printf(" '%s'\n", lstr(seq, len));
432	}
433
434	static void dump_n_compare(const char text, const char seq1, long len1, const char *seq2, long len2)
435	{
436	long offset = 0;
437
438	printf(" Comparing %s:\n", text);
439
440	while (len1>0 && len2>0)
441	{
442	long done = 0;
443
444	if (len1>=BUFLEN && len2>=BUFLEN)
445	{
446	dump_n_compare_one(seq1, seq2, done=BUFLEN, offset);
447	}
448	else
449	{
450	long min = len1<len2 ? len1 : len2;
451	dump_n_compare_one(seq1, seq2, done=min, offset);
452	}
453
454	seq1 += done;
455	seq2 += done;
456	len1 -= done;
457	len2 -= done;
458	offset += done;
459	}
460
461	if (len1>0) dump_rest(seq1, len1, 1, offset);
462	if (len2>0) dump_rest(seq2, len2, 2, offset);
463	printf(" -------------------\n");
464	}
465
466	static void dump_n_compare(const char *text, const AWTC_CompactedSubSequence& seq1, const AWTC_CompactedSubSequence& seq2)
467	{
468	dump_n_compare(text, seq1.text(), seq1.length(), seq2.text(), seq2.length());
469	}
470
471	#undef BUFLEN
472
473	static inline void dumpSeq(const char *seq, long len, long pos)
474	{
475	printf("'%s' ", lstr(seq, len));
476	printf("(Pos=%li,Len=%li)", pos, len);
477	}
478
479	#define dump() \
480	do \
481	{ \
482	double sig = partSignificance(sequence().length(), slaveSequence.length(), bestLength); \
483	\
484	printf(" Score = %li (Significance=%f)\n" \
485	" Master = ", bestScore, sig); \
486	dumpSeq(bestMasterLeft.text(), bestLength, bestMasterLeft.leftOf()); \
487	printf("\n" \
488	" Slave = "); \
489	dumpSeq(bestSlaveLeft.text(), bestLength, bestSlaveLeft.leftOf()); \
490	printf("\n"); \
491	} \
492	while(0)
493
494	#endif /DEBUG/
495
496
497	// --------------------------------------------------------------------------------
498	// INLINE-functions used in fast_align():
499	// --------------------------------------------------------------------------------
500
501	static inline double log3(double d)
502	{
503	return log(d)/log(3.0);
504	}
505	static inline double partSignificance(long seq1len, long seq2len, long partlen)
506	// returns log3 of significance of the part
507	// usage: partSignificance(...) < log3(maxAllowedSignificance)
508	{
509	return log3((seq1len-partlen)*(seq2len-partlen)) - partlen;
510	}
511
512	#if 0
513	static inline void recalcUsedBuffer(char *bufPtr, long lenPtr, long *usedPtr, long used)
514	{
515	*bufPtr += used;
516	*lenPtr -= used;
517	*usedPtr += used;
518	}
519	#endif
520
521	static inline GB_ERROR bufferTooSmall()
522	{
523	return GB_export_error("Cannot align - reserved buffer is to small");
524	}
525
526	static inline long insertsToNextBase(AWTC_alignBuffer *alignBuffer, const AWTC_SequencePosition& master)
527	{
528	int inserts;
529	int nextBase;
530
531	if (master.rightOf()>0) {
532	nextBase = master.expdPosition();
533	}
534	else {
535	nextBase = master.sequence().expdPosition(master.sequence().length());
536	}
537	inserts = nextBase-alignBuffer->offset();
538
539	return inserts;
540	}
541
542	static inline void insertBase(AWTC_alignBuffer *alignBuffer,
543	AWTC_SequencePosition& master, AWTC_SequencePosition& slave,
544	AWTC_fast_align_report *report)
545	{
546	char slaveBase = *slave.text();
547	char masterBase = *master.text();
548
549	alignBuffer->set(slaveBase, alignQuality(slaveBase, masterBase));
550	report->count_aligned_base(slaveBase!=masterBase);
551	++slave;
552	++master;
553	}
554
555	static inline void insertSlaveBases(AWTC_alignBuffer *alignBuffer,
556	AWTC_SequencePosition& slave,
557	int length,
558	AWTC_fast_align_report *report)
559	{
560	alignBuffer->copy(slave.text(), alignQuality(*slave.text(),GAP_CHAR), length);
561	report->count_unaligned_base(length);
562	slave += length;
563	}
564
565	static inline void insertGap(AWTC_alignBuffer *alignBuffer,
566	AWTC_SequencePosition& master,
567	AWTC_fast_align_report *report)
568	{
569	char masterBase = *master.text();
570
571	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, masterBase));
572	report->count_aligned_base(masterBase!=GAP_CHAR);
573	++master;
574	}
575
576	static inline GB_ERROR insertClustalValigned(AWTC_alignBuffer *alignBuffer,
577	AWTC_SequencePosition& master,
578	AWTC_SequencePosition& slave,
579	const char masterAlignment, const char slaveAlignment, long alignmentLength,
580	AWTC_fast_align_report *report)
581	// inserts bases of 'slave' to 'alignBuffer' according to alignment in 'masterAlignment' and 'slaveAlignment'
582	{
583	#define ACID '*' // contents of 'masterAlignment' and 'slaveAlignment'
584	#define GAP '-'
585
586	int pos;
587	int baseAtLeft = 0; // TRUE -> last position in alignBuffer contains a base
588
589	for (pos=0; pos<alignmentLength; pos++) {
590	long insert = insertsToNextBase(alignBuffer, master); // in expanded seq
591
592	if (masterAlignment[pos]==ACID) {
593	if (insert>0) {
594	if (insert>alignBuffer->free()) return bufferTooSmall();
595	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
596	awtc_assert(insertsToNextBase(alignBuffer,master)==0);
597	insert = 0;
598	}
599
600	if (!alignBuffer->free()) return bufferTooSmall();
601	if (slaveAlignment[pos]==ACID) {
602	insertBase(alignBuffer, master, slave, report);
603	baseAtLeft = 1;
604	}
605	else {
606	insertGap(alignBuffer, master, report);
607	baseAtLeft = 0;
608	}
609	}
610	else {
611	int slave_bases;
612	//int memo = 0;
613
614	awtc_assert(masterAlignment[pos]==GAP);
615	for (slave_bases=1; pos+slave_bases<alignmentLength && masterAlignment[pos+slave_bases]==GAP; slave_bases++) {
616	; // count gaps in master (= # of slave bases to insert)
617	}
618	if (!baseAtLeft && insert>slave_bases) {
619	int ins_gaps = insert-slave_bases;
620
621	awtc_assert(alignBuffer->free()>=ins_gaps);
622	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), ins_gaps);
623	insert -= ins_gaps;
624	}
625
626	if (insert<slave_bases) { // master has less gaps than slave bases to insert
627	report->memorize_insertion(master.expdPosition(), slave_bases-insert);
628	}
629	else if (insert>slave_bases) { // master has more gaps than slave bases to insert
630	awtc_assert(baseAtLeft);
631	}
632
633	int start = slave.expdPosition(); // memorize base positions without counterpart in master
634	int end = slave.expdPosition(slave_bases-1);
635	unaligned_bases.memorize(start, end);
636
637	if (slave_bases>alignBuffer->free()) return bufferTooSmall();
638	insertSlaveBases(alignBuffer, slave, slave_bases, report);
639	pos += slave_bases-1; // -1 compensates for()-increment
640	baseAtLeft = 1;
641	}
642	}
643
644	return 0;
645
646	#undef GAP
647	#undef ACID
648	}
649
650	static inline GB_ERROR insertAligned(AWTC_alignBuffer *alignBuffer,
651	AWTC_SequencePosition& master, AWTC_SequencePosition& slave, long partLength,
652	AWTC_fast_align_report *report)
653	// insert bases of 'slave' to 'alignBuffer' according to 'master'
654	{
655	if (partLength) {
656	long insert = insertsToNextBase(alignBuffer, master);
657
658	awtc_assert(partLength>=0);
659
660	if (insert<0) { // insert gaps into master
661	long min_insert = insert;
662
663	report->memorize_insertion(master.expdPosition(), -insert);
664
665	while (insert<0 && partLength) {
666	if (insert<min_insert) min_insert = insert;
667	if (!alignBuffer->free()) {
668	return bufferTooSmall();
669	}
670	insertBase(alignBuffer, master, slave, report);
671	partLength--;
672	insert = insertsToNextBase(alignBuffer,master);
673	}
674
675	awtc_assert(partLength>=0);
676	if (partLength==0) { // all inserted
677	return NULL;
678	}
679	}
680
681	awtc_assert(insert>=0);
682
683	if (insert>0) { // insert gaps into slave
684	if (insert>alignBuffer->free()) return bufferTooSmall();
685	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
686	awtc_assert(insertsToNextBase(alignBuffer,master)==0);
687	}
688
689	awtc_assert(partLength>=0);
690
691	while (partLength--) {
692	insert = insertsToNextBase(alignBuffer,master);
693
694	awtc_assert(insert>=0);
695	if (insert>0) {
696	if (insert>=alignBuffer->free()) return bufferTooSmall();
697	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
698	}
699	else {
700	if (!alignBuffer->free()) {
701	return bufferTooSmall();
702	}
703	}
704
705	insertBase(alignBuffer, master, slave, report);
706	}
707	}
708
709	return NULL;
710	}
711
712	static inline GB_ERROR cannot_fast_align(const AWTC_CompactedSubSequence& master, long moffset, long mlength,
713	const AWTC_CompactedSubSequence& slaveSequence, long soffset, long slength,
714	int max_seq_length,
715	AWTC_alignBuffer *alignBuffer,
716	AWTC_fast_align_report *report)
717	{
718	const char *mtext = master.text(moffset);
719	const char *stext = slaveSequence.text(soffset);
720	char maligned, saligned;
721	int len;
722	GB_ERROR error = 0;
723
724	if (slength) {
725	if (mlength) { // if slave- and master-sequences contain bases, we call the slow aligner
726	int score;
727
728	#if defined(DEBUG) && 1
729	printf("ClustalV-Align:\n");
730	printf(" mseq = '%s'\n", lstr(mtext, mlength));
731	printf(" sseq = '%s'\n", lstr(stext, slength));
732	#endif
733
734	int is_dna = -1;
735
736	switch (global_alignmentType) {
737	case GB_AT_RNA:
738	case GB_AT_DNA: is_dna = 1; break;
739	case GB_AT_AA: is_dna = 0; break;
740	default: error = GB_export_error("Unknown alignment type - aligner aborted"); break;
741	}
742
743	if (!error) {
744	error = AWTC_ClustalV_align(is_dna,
745	1,
746	mtext, mlength, stext, slength,
747	master.gapsBefore(moffset),
748	max_seq_length,
749	&maligned, &saligned, &len, &score);
750	}
751
752	if (!error) {
753	#if defined(DEBUG) && 1
754	printf("ClustalV returns:\n");
755	printf(" maligned = '%s'\n", lstr(maligned, len));
756	printf(" saligned = '%s'\n", lstr(saligned, len));
757	#endif
758
759	AWTC_SequencePosition masterPos(master, moffset);
760	AWTC_SequencePosition slavePos(slaveSequence, soffset);
761
762	error = insertClustalValigned(alignBuffer, masterPos, slavePos, maligned, saligned, len, report);
763
764	#if (defined(DEBUG) && 0)
765
766	AWTC_SequencePosition master2(master->sequence(), moffset);
767	AWTC_SequencePosition slave2(slaveSequence, soffset);
768	char *cmp = new char[len];
769
770	for (int l=0; l<len; l++) {
771	int gaps = 0;
772
773	if (maligned[l]=='*') {
774	maligned[l] = *master2.text();
775	++master2;
776	}
777	else {
778	gaps++;
779	}
780
781	if (saligned[l]=='*') {
782	saligned[l] = *slave2.text();
783	++slave2;
784	}
785	else {
786	gaps++;
787	}
788
789	cmp[l] = gaps \|\| maligned[l]==saligned[l] ? '=' : 'X';
790	}
791
792	printf(" master = '%s'\n", lstr(maligned,len));
793	printf(" slave = '%s'\n", lstr(saligned,len));
794	printf(" '%s'\n", lstr(cmp,len));
795
796	delete [] cmp;
797	#endif
798	}
799	}
800	else { // master is empty here, so we just copy in the slave bases
801	if (slength<=alignBuffer->free()) {
802	int start = slaveSequence.expdPosition(soffset);
803	int end = slaveSequence.expdPosition(soffset+slength-1);
804
805	unaligned_bases.memorize(start, end);
806	// int fill = start-alignBuffer->offset();
807	// if (fill>0) {
808	// alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), fill);
809	// }
810	alignBuffer->copy(slaveSequence.text(soffset), '?', slength); // '?' means not aligned (just inserted)
811	// corrected by at alignBuffer->correctUnalignedPositions()
812	report->count_unaligned_base(slength);
813	}
814	else {
815	error = bufferTooSmall();
816	}
817	}
818	}
819
820	return error;
821	}
822
823	// --------------------------------------------------------------------------------
824	// #define's for fast_align()
825	// --------------------------------------------------------------------------------
826
827	#define TEST_BETTER_SCORE() \
828	do \
829	{ \
830	if (score>bestScore) \
831	{ \
832	bestScore = score; \
833	bestLength = masterRight.text() - masterLeft.text(); \
834	bestMasterLeft = masterLeft; \
835	bestSlaveLeft = slaveLeft; \
836	/dump();/ \
837	} \
838	} \
839	while (0)
840
841	#define CAN_SCORE_LEFT() (masterLeft.leftOf() && slaveLeft.leftOf())
842	#define CAN_SCORE_RIGHT() (masterRight.rightOf() && slaveRight.rightOf())
843
844	#define SCORE_LEFT() \
845	do \
846	{ \
847	score += (--masterLeft).text()==(--slaveLeft).text() ? match : mismatch; \
848	TEST_BETTER_SCORE(); \
849	} \
850	while (0)
851
852	#define SCORE_RIGHT() \
853	do \
854	{ \
855	score += (++masterRight).text()==(++slaveRight).text() ? match : mismatch; \
856	TEST_BETTER_SCORE(); \
857	} \
858	while (0)
859
860	// --------------------------------------------------------------------------------
861	// AWTC_FastSearchSequence::fast_align()
862	// --------------------------------------------------------------------------------
863
864	GB_ERROR AWTC_FastSearchSequence::fast_align(const AWTC_CompactedSubSequence& slaveSequence,
865	AWTC_alignBuffer *alignBuffer,
866	int max_seq_length,
867	int match, int mismatch,
868	AWTC_fast_align_report *report) const
869	//
870	// aligns 'slaveSequence' to 'this'
871	//
872	// returns ==NULL -> all ok -> 'alignBuffer' contains aligned sequence
873	// !=NULL -> failure -> no results
874	//
875	{
876	GB_ERROR error = NULL;
877	int aligned = 0;
878
879	// set the following #if to zero to test ClustalV-Aligner (not fast_aligner)
880	#if 1
881
882	static double lowSignificance;
883	static int lowSignificanceInitialized;
884
885	if (!lowSignificanceInitialized) {
886	lowSignificance = log3(0.01);
887	//printf("lowSignificance=%f\n", lowSignificance);
888	lowSignificanceInitialized = 1;
889	}
890
891	AWTC_SequencePosition slave(slaveSequence);
892	long bestScore=0;
893	AWTC_SequencePosition bestMasterLeft(sequence());
894	AWTC_SequencePosition bestSlaveLeft(slaveSequence);
895	long bestLength=0;
896
897	while (slave.rightOf()>=3) { // with all triples of slaveSequence
898	AWTC_FastSearchOccurrence occurrence(*this, slave.text()); // "search" first
899	AWTC_SequencePosition rightmostSlave = slave;
900
901	while (occurrence.found()) { // with all occurrences of the triple
902	long score = match*3;
903	AWTC_SequencePosition masterLeft(occurrence.sequence(), occurrence.offset());
904	AWTC_SequencePosition masterRight(occurrence.sequence(), occurrence.offset()+3);
905	AWTC_SequencePosition slaveLeft(slave);
906	AWTC_SequencePosition slaveRight(slave,3);
907
908	while (score>0) {
909	if (CAN_SCORE_LEFT()) {
910	SCORE_LEFT();
911	}
912	else {
913	while (score>0 && CAN_SCORE_RIGHT()) {
914	SCORE_RIGHT();
915	}
916	break;
917	}
918
919	if (CAN_SCORE_RIGHT()) {
920	SCORE_RIGHT();
921	}
922	else {
923	while (score>0 && CAN_SCORE_LEFT()) {
924	SCORE_LEFT();
925	}
926	break;
927	}
928	}
929
930	occurrence.gotoNext(); // "search" next
931
932	if (rightmostSlave<slaveRight) {
933	rightmostSlave = slaveRight;
934	rightmostSlave -= 3;
935	}
936	}
937
938	if (rightmostSlave>slave) slave = rightmostSlave;
939	else ++slave;
940	}
941
942	if (bestLength) {
943	double sig = partSignificance(sequence().length(), slaveSequence.length(), bestLength);
944
945	//dump();
946
947	if (sig<lowSignificance) {
948	long masterLeftOf = bestMasterLeft.leftOf();
949	long masterRightStart = masterLeftOf+bestLength;
950	long masterRightOf = bestMasterLeft.rightOf()-bestLength;
951	long slaveLeftOf = bestSlaveLeft.leftOf();
952	long slaveRightStart = slaveLeftOf+bestLength;
953	long slaveRightOf = bestSlaveLeft.rightOf()-bestLength;
954
955	#define MIN_ALIGNMENT_RANGE 4
956
957	if (!error) {
958	if (masterLeftOf >= MIN_ALIGNMENT_RANGE && slaveLeftOf >= MIN_ALIGNMENT_RANGE) {
959	AWTC_CompactedSubSequence leftCompactedMaster(sequence(), 0, masterLeftOf);
960	AWTC_FastSearchSequence leftMaster(leftCompactedMaster);
961
962	error = leftMaster.fast_align(AWTC_CompactedSubSequence(slaveSequence, 0, slaveLeftOf),
963	alignBuffer, max_seq_length, match, mismatch, report);
964	}
965	else if (slaveLeftOf>0) {
966	error = cannot_fast_align(sequence(), 0, masterLeftOf,
967	slaveSequence, 0, slaveLeftOf,
968	max_seq_length, alignBuffer, report);
969	}
970
971	aligned = 1;
972	}
973
974	// align part of slave sequence according to master sequence:
975
976	if (!error) {
977	#if (defined(DEBUG) && 0)
978	long offset = alignBuffer->offset();
979	long used;
980	#endif
981	error = insertAligned(alignBuffer, bestMasterLeft, bestSlaveLeft, bestLength, report);
982	#if (defined(DEBUG) && 0)
983	used = alignBuffer->offset()-offset;
984	printf("aligned '%s' (len=%li, address=%li)\n", lstr(alignBuffer->text()+offset, used), used, long(alignBuffer));
985	#endif
986	aligned = 1;
987	}
988
989	if (!error) {
990	if (masterRightOf >= MIN_ALIGNMENT_RANGE && slaveRightOf >= MIN_ALIGNMENT_RANGE) {
991	AWTC_CompactedSubSequence rightCompactedMaster(sequence(), masterRightStart, masterRightOf);
992	AWTC_FastSearchSequence rightMaster(rightCompactedMaster);
993
994	error = rightMaster.fast_align(AWTC_CompactedSubSequence(slaveSequence, slaveRightStart, slaveRightOf),
995	alignBuffer,
996	max_seq_length, match, mismatch, report);
997	}
998	else if (slaveRightOf>0) {
999	error = cannot_fast_align(sequence(), masterRightStart, masterRightOf,
1000	slaveSequence, slaveRightStart, slaveRightOf,
1001	max_seq_length, alignBuffer, report);
1002	}
1003
1004	aligned = 1;
1005	}
1006
1007	}
1008	else {
1009	//printf("Not significant!\n");
1010	}
1011	}
1012
1013	#endif
1014
1015	if (!aligned && !error) {
1016	error = cannot_fast_align(sequence(), 0, sequence().length(),
1017	slaveSequence, 0, slaveSequence.length(),
1018	max_seq_length, alignBuffer, report);
1019	}
1020
1021	return error;
1022	}
1023
1024	#undef dump
1025	#undef TEST_BETTER_SCORE
1026	#undef CAN_SCORE_LEFT
1027	#undef CAN_SCORE_RIGHT
1028	#undef SCORE_LEFT
1029	#undef SCORE_RIGHT
1030
1031	// --------------------------------------------------------------------------------
1032	// Tools for AWTC_aligner()
1033	// --------------------------------------------------------------------------------
1034
1035	static long calcSequenceChecksum(const char *data, int length)
1036	{
1037	static char *gaps;
1038
1039	if (!gaps)
1040	{
1041	gaps = (char*)malloc(257);
1042
1043	int c = 1;
1044	int cnt = 0;
1045
1046	while (c<256)
1047	{
1048	if (AWTC_is_gap(toupper(c)))
1049	gaps[cnt++] = c;
1050	c++;
1051	}
1052
1053	gaps[cnt] = 0;
1054	}
1055
1056	long sum = GB_checksum(data, length, 1, gaps);
1057	return sum;
1058	}
1059
1060	static AWTC_CompactedSubSequence readCompactedSequence(GBDATA gb_species,
1061	const char *ali,
1062	GB_ERROR *errorPtr,
1063	char **dataPtr, // if dataPtr != NULL it will be set to the WHOLE uncompacted sequence
1064	long *seqChksum, // may be NULL (of part of sequence)
1065	int firstColumn, // return only part of the sequence
1066	int lastColumn) // (-1 -> till end of sequence)
1067	{
1068	GB_ERROR error = 0;
1069	GBDATA *gbd;
1070	AWTC_CompactedSubSequence *seq = 0;
1071
1072	gbd = GBT_read_sequence(gb_species, ali); // get sequence
1073
1074	if (gbd)
1075	{
1076	long length = GB_read_string_count(gbd);
1077	char *data = GB_read_string(gbd);
1078	long partLength;
1079	char *partData;
1080
1081	if (dataPtr) { // make a copy of the whole uncompacted sequence (returned to caller)
1082	*dataPtr = data;
1083	}
1084
1085	if (firstColumn!=0 \|\| lastColumn!=-1) { // take only part of sequence
1086	awtc_assert(firstColumn>=0);
1087	awtc_assert(firstColumn<=lastColumn);
1088
1089	// include all surrounding gaps
1090	while (firstColumn>0 && AWTC_is_gap(data[firstColumn-1])) {
1091	firstColumn--;
1092	}
1093	if (lastColumn!=-1) {
1094	while (lastColumn<(length+1) && AWTC_is_gap(data[lastColumn+1])) {
1095	lastColumn++;
1096	}
1097	}
1098
1099	partData = data+firstColumn;
1100	int slen = length-firstColumn;
1101
1102	awtc_assert(slen>=0);
1103	awtc_assert((size_t)slen==strlen(partData));
1104
1105	if (lastColumn==-1) { // take all till end of sequence
1106	partLength = slen;
1107	}
1108	else {
1109	partLength = lastColumn-firstColumn+1;
1110	if (partLength>slen) partLength = slen; // cut rest, if we have any
1111	}
1112	}
1113	else {
1114	partLength = length;
1115	partData = data;
1116	}
1117
1118	if (!error) {
1119	if (seqChksum) {
1120	*seqChksum = calcSequenceChecksum(partData, partLength);
1121	}
1122
1123	seq = new AWTC_CompactedSubSequence(partData, partLength, GBT_read_name(gb_species), firstColumn);
1124	}
1125
1126	if (!dataPtr) { // free sequence only if user has not requested to get it
1127	free(data);
1128	}
1129	}
1130	else {
1131	error = GB_export_errorf("No 'data' found for species '%s'", GBT_read_name(gb_species));
1132	if (dataPtr) *dataPtr = NULL; // (user must not care to free data if we fail)
1133	}
1134
1135	awtc_assert(errorPtr);
1136	*errorPtr = error;
1137
1138	return seq;
1139	}
1140
1141	static GB_ERROR writeStringToAlignment(GBDATA *gb_species, GB_CSTR alignment, GB_CSTR data_name, GB_CSTR str, bool temporary) {
1142	GBDATA *gb_ali = GB_search(gb_species, alignment, GB_DB);
1143	GB_ERROR error = NULL;
1144	GBDATA *gb_name = GB_search(gb_ali, data_name, GB_STRING);
1145
1146	if (gb_name) {
1147	awtc_assert(GB_check_father(gb_name, gb_ali));
1148	error = GB_write_string(gb_name, str);
1149	if (temporary && !error) error = GB_set_temporary(gb_name);
1150	}
1151	else {
1152	error = GB_export_errorf("Cannot create entry '%s' for '%s'", data_name, GBT_read_name(gb_species));
1153	}
1154
1155	return error;
1156	}
1157
1158	// --------------------------------------------------------------------------------
1159
1160	static GB_ERROR alignCompactedTo(AWTC_CompactedSubSequence *toAlignSequence,
1161	const AWTC_FastSearchSequence *alignTo,
1162	int max_seq_length,
1163	GB_CSTR alignment,
1164	long toAlignChksum,
1165	GBDATA *gb_toAlign,
1166	GBDATA *gb_alignTo, // may be NULL
1167	const AlignParams& ali_params)
1168	// if only part of the sequence should be aligned, then this functions already gets only the part
1169	// (i.o.w.: toAlignSequence, alignTo and toAlignChksum refer to the partial sequence)
1170	{
1171	AWTC_alignBuffer alignBuffer(max_seq_length);
1172	const char *master_name = AWTC_read_name(gb_alignTo);
1173
1174	AWTC_fast_align_report report(master_name, ali_params.showGapsMessages);
1175
1176	{
1177	const char *to_align_name = AWTC_read_name(gb_toAlign);
1178	const char *align_to_name = AWTC_read_name(gb_alignTo);
1179
1180	const char *stat_buf = GBS_global_string("Aligning #%i:%i %s (to %s)",
1181	currentSequenceNumber, overallSequenceNumber,
1182	to_align_name, align_to_name);
1183
1184	static GBDATA *last_gb_toAlign = 0;
1185	if (gb_toAlign!=last_gb_toAlign) {
1186	last_gb_toAlign = gb_toAlign;
1187	currentSequenceNumber++;
1188	}
1189
1190	aw_status(stat_buf);
1191	}
1192
1193	#if (defined(DEBUG) && 1)
1194	printf("alignCompactedTo(): master='%s' ", master_name);
1195	printf("slave='%s'\n", toAlignSequence->name());
1196	#endif
1197
1198	GB_ERROR error = GB_pop_transaction(gb_toAlign);
1199	if (!error) {
1200	if (island_hopper) {
1201	error = island_hopper->do_align();
1202	if (!error) {
1203	awtc_assert(island_hopper->was_aligned());
1204
1205	#if defined(DEBUG)
1206	printf("Island-Hopper returns:\n");
1207	printf("maligned = '%s'\n", lstr(island_hopper->get_result_ref(), island_hopper->get_result_length()));
1208	printf("saligned = '%s'\n", lstr(island_hopper->get_result(), island_hopper->get_result_length()));
1209	#endif // DEBUG
1210
1211	AWTC_SequencePosition masterPos(alignTo->sequence(), 0);
1212	AWTC_SequencePosition slavePos(*toAlignSequence, 0);
1213
1214	error = insertClustalValigned(&alignBuffer, masterPos, slavePos, island_hopper->get_result_ref(), island_hopper->get_result(), island_hopper->get_result_length(), &report);
1215	}
1216	}
1217	else {
1218	error = alignTo->fast_align(*toAlignSequence, &alignBuffer, max_seq_length, 2, -10, &report); // <- here we align
1219	}
1220	}
1221
1222	if (!error) {
1223	alignBuffer.correctUnalignedPositions();
1224	if (alignBuffer.free()) {
1225	alignBuffer.set('-', alignQuality(GAP_CHAR,GAP_CHAR), alignBuffer.free()); // rest of alignBuffer is set to '-'
1226	}
1227	alignBuffer.expandPoints(*toAlignSequence);
1228	}
1229
1230	#if (defined(DEBUG) && 1)
1231	if (!error) {
1232	AWTC_CompactedSubSequence alignedSlave(alignBuffer.text(), alignBuffer.length(), toAlignSequence->name());
1233	dump_n_compare("reference vs. aligned:", alignTo->sequence(), alignedSlave);
1234	}
1235	#endif
1236
1237	{
1238	GB_ERROR err = GB_push_transaction(gb_toAlign);
1239	if (!error) error = err;
1240	}
1241
1242	if (!error) {
1243	if (calcSequenceChecksum(alignBuffer.text(), alignBuffer.length())!=toAlignChksum) { // sequence-chksum changed
1244	error = GB_export_error("Internal aligner error (sequence checksum changed) -- aborted");
1245
1246	# ifdef DEBUG
1247	AWTC_CompactedSubSequence alignedSlave(alignBuffer.text(), alignBuffer.length(), toAlignSequence->name());
1248	dump_n_compare("Old Slave vs. new Slave", *toAlignSequence, alignedSlave);
1249	# endif
1250	}
1251	else {
1252	GB_push_my_security(gb_toAlign);
1253	{
1254	GBDATA *gbd = GBT_add_data(gb_toAlign, alignment, "data", GB_STRING);
1255
1256	if (!gbd) {
1257	error = GB_export_error("Can't find/create sequence data");
1258	}
1259	else {
1260	if (ali_params.firstColumn!=0 \|\| ali_params.lastColumn!=-1) { // we aligned just a part of the sequence
1261	char *buffer = GB_read_string(gbd); // so we read old sequence data
1262	long len = GB_read_string_count(gbd);
1263	if (!buffer) error = GB_await_error();
1264	else {
1265	int lenToCopy = ali_params.lastColumn-ali_params.firstColumn+1;
1266	long wholeChksum = calcSequenceChecksum(buffer,len);
1267	// long partChksum = calcSequenceChecksum(buffer+firstColumn, lenToCopy);
1268
1269	memcpy(buffer+ali_params.firstColumn, alignBuffer.text()+ali_params.firstColumn, lenToCopy); // copy in the aligned part
1270
1271	if (calcSequenceChecksum(buffer, len)!=wholeChksum) {
1272	error = GB_export_error("Internal aligner error (sequence checksum changed) -- aborted");
1273	}
1274	else {
1275	GB_write_string(gbd,"");
1276	error = GB_write_string(gbd, buffer);
1277	}
1278	}
1279
1280	free(buffer);
1281	}
1282	else {
1283	alignBuffer.point_ends_of();
1284	error = GBT_write_sequence(gbd, alignment, max_seq_length, alignBuffer.text()); // aligned all -> write all
1285	}
1286	}
1287	}
1288	GB_pop_my_security(gb_toAlign);
1289
1290	// if (!error) error = err;
1291
1292	if (!error && ali_params.report != FA_NO_REPORT) {
1293	// create temp-entry for slave containing alignment quality:
1294
1295	error = writeStringToAlignment(gb_toAlign, alignment, QUALITY_NAME, alignBuffer.quality(), ali_params.report==FA_TEMP_REPORT);
1296	if (!error) {
1297	// create temp-entry for master containing insert points:
1298
1299	int buflen = max_seq_length*2;
1300	char buffer = (char)malloc(buflen+1);
1301	char *afterLast = buffer;
1302
1303	if (!buffer) {
1304	error = AWTC_memerr();
1305	}
1306	else {
1307	memset(buffer, '-', buflen);
1308	buffer[buflen] = 0;
1309
1310	const AWTC_fast_align_insertion *inserts = report.insertion();
1311	while (inserts)
1312	{
1313	memset(buffer+inserts->offset(), '>', inserts->gaps());
1314	afterLast = buffer+inserts->offset()+inserts->gaps();
1315	inserts = inserts->next();
1316	}
1317
1318	*afterLast = 0;
1319	if (gb_alignTo) {
1320	error = writeStringToAlignment(gb_alignTo, alignment, INSERTS_NAME, buffer, ali_params.report==FA_TEMP_REPORT);
1321	}
1322	}
1323	}
1324	}
1325	}
1326	}
1327	return error;
1328	}
1329
1330	GB_ERROR AWTC_delete_temp_entries(GBDATA *gb_species, GB_CSTR alignment)
1331	{
1332	awtc_assert(gb_species);
1333	awtc_assert(alignment);
1334
1335	GBDATA *gb_ali = GB_search(gb_species, alignment, GB_FIND);
1336	GB_ERROR error = NULL;
1337
1338	if (gb_ali)
1339	{
1340	GBDATA *gb_name = GB_search(gb_ali, QUALITY_NAME, GB_FIND);
1341	if (gb_name) {
1342	error = GB_delete(gb_name);
1343	#if defined(DEBUG)
1344	printf("deleted '%s' of '%s' (gb_name=%p)\n", QUALITY_NAME, GBT_read_name(gb_species), gb_name);
1345	#endif
1346	}
1347
1348	if (!error) {
1349	gb_name = GB_search(gb_ali, INSERTS_NAME, GB_FIND);
1350	if (gb_name) {
1351	error = GB_delete(gb_name);
1352	#if defined(DEBUG)
1353	printf("deleted '%s' of '%s' (gb_name=%p)\n", INSERTS_NAME, GBT_read_name(gb_species), gb_name);
1354	#endif
1355	}
1356	}
1357	}
1358
1359	return error;
1360	}
1361
1362	ATTRIBUTED(__ATTR__USERESULT, static GB_ERROR GB_align_error(GB_ERROR olderr, GBDATA gb_toAlign, GBDATA gb_alignTo))
1363	// used by alignTo() and alignToNextRelative() to transform errors coming from subroutines
1364	// can be used by higher functions
1365	{
1366	const char *name_toAlign = AWTC_read_name(gb_toAlign);
1367	const char *name_alignTo = AWTC_read_name(gb_alignTo);
1368
1369	return GBS_global_string("Error while aligning '%s' to '%s':\n%s", name_toAlign, name_alignTo, olderr);
1370	}
1371
1372	static GB_ERROR alignTo(GBDATA *gb_toAlign,
1373	GB_CSTR alignment,
1374	const AWTC_FastSearchSequence *alignTo,
1375	GBDATA *gb_alignTo, // may be NULL
1376	int max_seq_length,
1377	const AlignParams& ali_params)
1378	{
1379	GB_ERROR error = NULL;
1380	long chksum;
1381	AWTC_CompactedSubSequence *toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, NULL, &chksum, ali_params.firstColumn, ali_params.lastColumn);
1382
1383	if (island_hopper) {
1384	GBDATA *gb_seq = GBT_read_sequence(gb_toAlign, alignment); // get sequence
1385	if (gb_seq) {
1386	long length = GB_read_string_count(gb_seq);
1387	const char *data = GB_read_char_pntr(gb_seq);
1388
1389	island_hopper->set_toAlign_sequence(data);
1390	island_hopper->set_alignment_length(length);
1391	}
1392	}
1393
1394
1395
1396	if (!error)
1397	{
1398	error = alignCompactedTo(toAlignSequence, alignTo, max_seq_length, alignment, chksum, gb_toAlign, gb_alignTo, ali_params);
1399	if (error) error = GB_align_error(error, gb_toAlign, gb_alignTo);
1400	delete toAlignSequence;
1401	}
1402
1403	return error;
1404	}
1405
1406	static GB_ERROR alignToGroupConsensus(GBDATA *gb_toAlign,
1407	GB_CSTR alignment,
1408	AWTC_get_consensus_func get_consensus,
1409	int max_seq_length,
1410	const AlignParams& ali_params)
1411	{
1412	GB_ERROR error = 0;
1413	const char *species_name = AWTC_read_name(gb_toAlign);
1414	char *consensus = get_consensus(species_name, ali_params.firstColumn, ali_params.lastColumn);
1415	size_t cons_len = strlen(consensus);
1416
1417	for (size_t i = 0; i<cons_len; ++i) { // translate consensus to be accepted by aligner
1418	switch (consensus[i]) {
1419	case '=': consensus[i] = '-'; break;
1420	default : break;
1421	}
1422	}
1423
1424	AWTC_CompactedSubSequence compacted(consensus, cons_len, "group consensus");
1425
1426	{
1427	AWTC_FastSearchSequence fast(compacted);
1428	error = alignTo(gb_toAlign, alignment, &fast, NULL, max_seq_length, ali_params);
1429	}
1430
1431	free(consensus);
1432
1433	return error;
1434	}
1435
1436	static void appendNameAndUsedBasePositions(char *toString, GBDATA gb_species, int usedBasePositions) {
1437	// if usedBasePositions == -1 -> unknown;
1438
1439	char *currInfo;
1440	if (usedBasePositions<0) {
1441	currInfo = strdup(GBT_read_name(gb_species));
1442	}
1443	else {
1444	awtc_assert(usedBasePositions>0); // otherwise it should NOT be announced here!
1445	currInfo = GBS_global_string_copy("%s:%i", GBT_read_name(gb_species), usedBasePositions);
1446	}
1447
1448	char *newString = 0;
1449	if (*toString) {
1450	newString = GBS_global_string_copy("%s, %s", *toString, currInfo);
1451	}
1452	else {
1453	newString = currInfo;
1454	currInfo = 0; // dont free
1455	}
1456
1457	freeset(*toString, newString);
1458	free(currInfo);
1459	}
1460
1461	inline int min(int i, int j) { return i<j ? i : j; }
1462
1463	static GB_ERROR alignToNextRelative(const SearchRelativeParams& relSearch,
1464	int max_seq_length,
1465	FA_turn turnAllowed,
1466	GB_CSTR alignment,
1467	GBDATA *gb_toAlign,
1468	const AlignParams& ali_params)
1469	{
1470	AWTC_CompactedSubSequence *toAlignSequence = NULL;
1471	bool use_different_pt_server_alignment = 0 != strcmp(relSearch.pt_server_alignment, alignment);
1472
1473	GB_ERROR error;
1474	long chksum;
1475	int restart = 1;
1476	int relativesToTest = relSearch.maxRelatives*2; // get more relatives from pt-server (needed when use_different_pt_server_alignment == true)
1477	char *nearestRelative = new char[relativesToTest+1];
1478	int next_relatives;
1479	int i;
1480
1481	if (use_different_pt_server_alignment) {
1482	turnAllowed = FA_TURN_NEVER; // makes no sense if we're using a different alignment for the pt_server
1483	}
1484
1485	for (next_relatives=0; next_relatives<relativesToTest; next_relatives++) {
1486	nearestRelative[next_relatives] = 0;
1487	}
1488	next_relatives = 0;
1489
1490	while (restart) {
1491	char *toAlignExpSequence = 0;
1492
1493	restart = 0;
1494
1495	if (use_different_pt_server_alignment) {
1496	toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, 0, &chksum, ali_params.firstColumn, ali_params.lastColumn);
1497
1498	GBDATA *gbd = GBT_read_sequence(gb_toAlign, relSearch.pt_server_alignment); // use a different alignment for next relative search
1499	if (!gbd) {
1500	error = GB_export_errorf("Species '%s' has no data in alignment '%s'", GBT_read_name(gb_toAlign), relSearch.pt_server_alignment);
1501	}
1502	else {
1503	toAlignExpSequence = GB_read_string(gbd);
1504	}
1505	}
1506	else {
1507	toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, &toAlignExpSequence, &chksum, ali_params.firstColumn, ali_params.lastColumn);
1508	}
1509
1510	if (error) {
1511	return error;
1512	}
1513
1514	while (next_relatives) {
1515	next_relatives--;
1516	freeset(nearestRelative[next_relatives], 0);
1517	}
1518
1519	{
1520	// find relatives
1521	AWTC_FIND_FAMILY family(GLOBAL_gb_main);
1522	double bestScore = 0;
1523
1524	aw_status("Searching relatives");
1525	error = family.findFamily(relSearch.pt_server_id,
1526	toAlignExpSequence,
1527	relSearch.fam_oligo_len,
1528	relSearch.fam_mismatches,
1529	relSearch.fam_fast_mode,
1530	relSearch.fam_rel_matches,
1531	FF_FORWARD,
1532	relativesToTest+1);
1533
1534	if (!error) {
1535	#if defined(DEBUG)
1536	double lastScore = -1;
1537	#endif // DEBUG
1538	for (const AWTC_FIND_FAMILY_MEMBER *fl = family.getFamilyList(); fl; fl=fl->next) {
1539	if (strcmp(toAlignSequence->name(), fl->name)!=0) {
1540	if (GBT_find_species(GLOBAL_gb_main, fl->name)) {
1541	double thisScore = relSearch.fam_rel_matches ? fl->rel_matches : fl->matches;
1542	#if defined(DEBUG)
1543	// check whether family list is sorted correctly
1544	awtc_assert(lastScore < 0 \|\| lastScore >= thisScore);
1545	lastScore = thisScore;
1546	#endif // DEBUG
1547	if (thisScore>=bestScore) bestScore = thisScore;
1548	if (next_relatives<(relativesToTest+1)) {
1549	nearestRelative[next_relatives] = strdup(fl->name);
1550	next_relatives++;
1551	}
1552	}
1553	}
1554	}
1555	}
1556
1557	if (!error && turnAllowed != FA_TURN_NEVER) { // test if mirrored sequence has better relatives
1558	char *mirroredSequence = strdup(toAlignExpSequence);
1559	long length = strlen(mirroredSequence);
1560	double bestMirroredScore = 0;
1561
1562	char T_or_U;
1563	error = GBT_determine_T_or_U(global_alignmentType, &T_or_U, "reverse-complement");
1564	GBT_reverseComplementNucSequence(mirroredSequence, length, T_or_U);
1565
1566	error = family.findFamily(relSearch.pt_server_id,
1567	mirroredSequence,
1568	relSearch.fam_oligo_len,
1569	relSearch.fam_mismatches,
1570	relSearch.fam_fast_mode,
1571	relSearch.fam_rel_matches,
1572	FF_FORWARD,
1573	relativesToTest+1);
1574	if (!error) {
1575	#if defined(DEBUG)
1576	double lastScore = -1;
1577	#endif // DEBUG
1578	for (const AWTC_FIND_FAMILY_MEMBER *fl = family.getFamilyList(); fl; fl = fl->next) {
1579	double thisScore = relSearch.fam_rel_matches ? fl->rel_matches : fl->matches;
1580	#if defined(DEBUG)
1581	// check whether family list is sorted correctly
1582	awtc_assert(lastScore < 0 \|\| lastScore >= thisScore);
1583	lastScore = thisScore;
1584	#endif // DEBUG
1585	if (thisScore >= bestMirroredScore) {
1586	if (strcmp(toAlignSequence->name(), fl->name)!=0) {
1587	if (GBT_find_species(GLOBAL_gb_main, fl->name)) bestMirroredScore = thisScore;
1588	}
1589	}
1590	}
1591	}
1592
1593	int turnIt = 0;
1594	if (bestMirroredScore>bestScore) {
1595	if (turnAllowed==FA_TURN_INTERACTIVE) {
1596	const char *message;
1597	if (relSearch.fam_rel_matches) {
1598	message = GBS_global_string("'%s' seems to be the other way round (score: %.1f%%, score if turned: %.1f%%)",
1599	toAlignSequence->name(), bestScore100, bestMirroredScore100);
1600	}
1601	else {
1602	message = GBS_global_string("'%s' seems to be the other way round (score: %li, score if turned: %li)",
1603	toAlignSequence->name(), long(bestScore+.5), long(bestMirroredScore+.5));
1604	}
1605	turnIt = aw_question(message, "Turn sequence,Leave sequence alone")==0;
1606	}
1607	else {
1608	awtc_assert(turnAllowed == FA_TURN_ALWAYS);
1609	turnIt = 1;
1610	}
1611	}
1612
1613	if (turnIt) { // write mirrored sequence
1614	GBDATA *gbd = GBT_read_sequence(gb_toAlign, alignment);
1615	GB_push_my_security(gbd);
1616	error = GB_write_string(gbd, mirroredSequence);
1617	GB_pop_my_security(gbd);
1618	if (!error) {
1619	delete toAlignSequence;
1620	restart = 1; // continue while loop
1621	}
1622	}
1623
1624	free(mirroredSequence);
1625	}
1626	}
1627	free(toAlignExpSequence);
1628	}
1629
1630	if (!error) {
1631	if (!next_relatives) {
1632	char warning[200];
1633	sprintf(warning, "No relative found for '%s'", toAlignSequence->name());
1634	aw_message(warning);
1635	}
1636	else {
1637	// assuming relatives are sorted! (nearest to farer)
1638
1639	// get data pointers
1640	typedef GBDATA *GBDATAP;
1641	GBDATAP *gb_reference = new GBDATAP[relSearch.maxRelatives];
1642
1643	{
1644	for (i=0; i<relSearch.maxRelatives && i<next_relatives; i++) {
1645	GBDATA *gb_species = GBT_find_species(GLOBAL_gb_main, nearestRelative[i]);
1646	if (!gb_species) { // pt-server seems not up to date!
1647	error = species_not_found(nearestRelative[i]);
1648	break;
1649	}
1650
1651	GBDATA *gb_sequence = GBT_read_sequence(gb_species, alignment);
1652	if (gb_sequence) { // if relative has sequence data in the current alignment ..
1653	gb_reference[i] = gb_species;
1654	}
1655	else { // remove this relative
1656	free(nearestRelative[i]);
1657	for (int j = i+1; j<next_relatives; ++j) {
1658	nearestRelative[j-1] = nearestRelative[j];
1659	}
1660	next_relatives--;
1661	nearestRelative[next_relatives] = 0;
1662	i--; // redo same index
1663	}
1664	}
1665
1666	// delete superfluous relatives
1667	for (; i<next_relatives; ++i) freeset(nearestRelative[i], 0);
1668
1669	if (next_relatives>relSearch.maxRelatives) {
1670	next_relatives = relSearch.maxRelatives;
1671	}
1672	}
1673
1674	// align
1675
1676	if (!error) {
1677	AWTC_CompactedSubSequence *alignToSequence = readCompactedSequence(gb_reference[0], alignment, &error, NULL, NULL, ali_params.firstColumn, ali_params.lastColumn);
1678	awtc_assert(alignToSequence != 0);
1679
1680	if (island_hopper) {
1681	GBDATA *gb_ref = GBT_read_sequence(gb_reference[0], alignment); // get reference sequence
1682	GBDATA *gb_align = GBT_read_sequence(gb_toAlign, alignment); // get sequence to align
1683
1684	if (gb_ref && gb_align) {
1685	long length_ref = GB_read_string_count(gb_ref);
1686	// long length_align = GB_read_string_count(gb_align);
1687	const char *data;
1688
1689	data = GB_read_char_pntr(gb_ref);
1690	island_hopper->set_ref_sequence(data);
1691
1692	data = GB_read_char_pntr(gb_align);
1693	island_hopper->set_toAlign_sequence(data);
1694
1695	island_hopper->set_alignment_length(length_ref);
1696	}
1697	}
1698
1699	{
1700	AWTC_FastSearchSequence referenceFastSeq(*alignToSequence);
1701
1702	error = alignCompactedTo(toAlignSequence, &referenceFastSeq,
1703	max_seq_length, alignment, chksum,
1704	gb_toAlign, gb_reference[0], ali_params);
1705	}
1706
1707	if (error) {
1708	error = GB_align_error(error, gb_toAlign, gb_reference[0]);
1709	}
1710	else {
1711	char *used_relatives = 0;
1712
1713	if (!unaligned_bases.is_empty()) {
1714	if (island_hopper) {
1715	appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], -1);
1716	if (next_relatives>1) error = "Island hopping uses only one relative";
1717	}
1718	else {
1719	UnalignedBasesList ubl;
1720	UnalignedBasesList ubl_for_next_relative;
1721
1722	int unaligned_positions;
1723	{
1724	AWTC_CompactedSubSequence *alignedSequence = readCompactedSequence(gb_toAlign, alignment, &error, 0, 0, ali_params.firstColumn, ali_params.lastColumn);
1725
1726	unaligned_positions = ubl.add_and_recalc_positions(&unaligned_bases, toAlignSequence, alignedSequence);
1727	// now ubl holds the unaligned (and recalculated) parts from last relative
1728	delete alignedSequence;
1729	}
1730
1731	int toalign_positions = toAlignSequence->length();
1732
1733	if (unaligned_positions<toalign_positions) {
1734	appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], toalign_positions-unaligned_positions);
1735	}
1736
1737	for (i=1; i<next_relatives && !error; i++) {
1738	ubl.add(&ubl_for_next_relative);
1739	int unaligned_positions_for_next = 0;
1740
1741	while (!ubl.is_empty() && !error) {
1742	int start, end;
1743	ubl.recall(&start, &end);
1744
1745	awtc_assert(ali_params.firstColumn<=start && start<=end && (end<=ali_params.lastColumn \|\| ali_params.lastColumn==-1));
1746
1747	AWTC_CompactedSubSequence *alignToPart = readCompactedSequence(gb_reference[i], alignment, &error, 0, 0, start, end);
1748	if (error) break;
1749
1750	long part_chksum;
1751	AWTC_CompactedSubSequence *toAlignPart = readCompactedSequence(gb_toAlign, alignment, &error, 0, &part_chksum, start, end);
1752	if (error) break;
1753
1754	AlignParams ubl_ali_params = {
1755	ali_params.report,
1756	ali_params.showGapsMessages,
1757	start,
1758	end
1759	};
1760
1761	AWTC_FastSearchSequence referenceFastSeq(*alignToPart);
1762	error = alignCompactedTo(toAlignPart, &referenceFastSeq,
1763	max_seq_length, alignment, part_chksum,
1764	gb_toAlign, gb_reference[i], ubl_ali_params);
1765
1766	AWTC_CompactedSubSequence *alignedPart = readCompactedSequence(gb_toAlign, alignment, &error, 0, 0, start, end);
1767
1768	unaligned_positions_for_next += ubl_for_next_relative.add_and_recalc_positions(&unaligned_bases, toAlignPart, alignedPart);
1769
1770	delete alignedPart;
1771	delete alignToPart;
1772	delete toAlignPart;
1773	}
1774
1775	awtc_assert(unaligned_positions_for_next <= unaligned_positions); // means: number of unaligned positions has increased by use of relative
1776	if (unaligned_positions_for_next<unaligned_positions) {
1777	appendNameAndUsedBasePositions(&used_relatives, gb_reference[i], unaligned_positions-unaligned_positions_for_next);
1778	unaligned_positions = unaligned_positions_for_next;
1779	}
1780	}
1781	}
1782	}
1783
1784	if (!error) { // write used relatives to db-field
1785	error = GBT_write_string(gb_toAlign, "used_rels", used_relatives);
1786	}
1787	}
1788
1789	delete alignToSequence;
1790	}
1791
1792	delete [] gb_reference;
1793	}
1794	}
1795
1796	delete toAlignSequence;
1797
1798	for (i=0; i<next_relatives; i++) freeset(nearestRelative[i], 0);
1799	delete [] nearestRelative;
1800
1801	return error;
1802	}
1803
1804	// --------------------------------------------------------------------------------
1805	// AWTC_aligner()
1806	// --------------------------------------------------------------------------------
1807
1808	static GB_ERROR AWTC_aligner(
1809	// define alignment target(s):
1810	FA_alignTarget alignWhat,
1811	GB_CSTR alignment, // name of alignment to use (==NULL -> use default)
1812	GB_CSTR toalign, // name of species to align (used if alignWhat == FA_CURRENT)
1813	AWTC_get_first_selected_species get_first_selected_species, // used if alignWhat == FA_SELECTED
1814	AWTC_get_next_selected_species get_next_selected_species, // --- "" ---
1815
1816	// define reference sequence(s):
1817	GB_CSTR reference, // name of reference species
1818	AWTC_get_consensus_func get_consensus, // function to get consensus seq
1819	const SearchRelativeParams& relSearch, // params to search for relatives
1820
1821	// general params:
1822	FA_turn turnAllowed,
1823	const AlignParams& ali_params,
1824	int maxProtection) // protection level
1825	{
1826	// (reference==NULL && get_consensus==NULL -> use next relative for (each) sequence)
1827
1828	GB_ERROR error = GB_push_transaction(GLOBAL_gb_main);
1829	int search_by_pt_server = reference==NULL && get_consensus==NULL;
1830	int err_count = 0;
1831	int wasNotAllowedToAlign = 0; // incremented for every sequence which has higher protection level (and was not aligned)
1832
1833	awtc_assert(reference==NULL \|\| get_consensus==NULL); // can't do both modes
1834
1835	if (turnAllowed != FA_TURN_NEVER) {
1836	if ((ali_params.firstColumn!=0 \|\| ali_params.lastColumn!=-1) \|\| !search_by_pt_server) {
1837	// if not selected 'Range/Whole sequence' or not selected 'Reference/Auto search..'
1838	turnAllowed = FA_TURN_NEVER; // then disable mirroring for the actual call
1839	}
1840	}
1841
1842	if (!error && !alignment) {
1843	alignment = (GB_CSTR)GBT_get_default_alignment(GLOBAL_gb_main); // get default alignment
1844	if (!alignment) {
1845	error = GB_export_error("No default alignment");
1846	}
1847	}
1848
1849	if (!error && alignment) {
1850	global_alignmentType = GBT_get_alignment_type(GLOBAL_gb_main, alignment);
1851	if (search_by_pt_server) {
1852	GB_alignment_type pt_server_alignmentType = GBT_get_alignment_type(GLOBAL_gb_main, relSearch.pt_server_alignment);
1853
1854	if (pt_server_alignmentType != GB_AT_RNA &&
1855	pt_server_alignmentType != GB_AT_DNA) {
1856	error = GB_export_error("pt_servers only support RNA/DNA sequences.\n"
1857	"In the aligner window you may specify a RNA/DNA alignment \n"
1858	"and use a pt_server build on that alignment.");
1859	}
1860	}
1861	}
1862
1863	if (!error) {
1864	GBDATA *gb_species_data = GB_search(GLOBAL_gb_main, "species_data", GB_CREATE_CONTAINER);
1865	int max_seq_length = GBT_get_alignment_len(GLOBAL_gb_main, alignment);
1866
1867	if (reference) { // align to reference sequence
1868	GBDATA *gb_reference = GBT_find_species_rel_species_data(gb_species_data, reference);
1869
1870	if (!gb_reference) {
1871	error = species_not_found(reference);
1872	}
1873	else {
1874	long referenceChksum;
1875	AWTC_CompactedSubSequence *referenceSeq = readCompactedSequence(gb_reference, alignment, &error, NULL, &referenceChksum, ali_params.firstColumn, ali_params.lastColumn);
1876	if (island_hopper) {
1877	GBDATA *gb_seq = GBT_read_sequence(gb_reference, alignment); // get sequence
1878	if (gb_seq) {
1879	long length = GB_read_string_count(gb_seq);
1880	const char *data = GB_read_char_pntr(gb_seq);
1881
1882	island_hopper->set_ref_sequence(data);
1883	island_hopper->set_alignment_length(length);
1884	}
1885	}
1886
1887
1888	if (!error) {
1889	AWTC_FastSearchSequence referenceFastSeq(*referenceSeq);
1890
1891	switch (alignWhat) {
1892	case FA_CURRENT: { // align one sequence
1893	awtc_assert(toalign);
1894	GBDATA *gb_toalign = GBT_find_species_rel_species_data(gb_species_data, toalign);
1895	if (!gb_toalign) {
1896	error = species_not_found(toalign);
1897	}
1898	else {
1899	currentSequenceNumber = overallSequenceNumber = 1;
1900	int myProtection = GB_read_security_write(GBT_read_sequence(gb_toalign, alignment));
1901	error = 0;
1902	if (myProtection<=maxProtection) {
1903	error = alignTo(gb_toalign, alignment, &referenceFastSeq, gb_reference, max_seq_length, ali_params);
1904	}
1905	else {
1906	wasNotAllowedToAlign++;
1907	}
1908	}
1909	break;
1910	}
1911	case FA_MARKED: { // align all marked sequences
1912	int count = GBT_count_marked_species(GLOBAL_gb_main);
1913	int done = 0;
1914	GBDATA *gb_species = GBT_first_marked_species_rel_species_data(gb_species_data);
1915
1916	currentSequenceNumber = 1;
1917	overallSequenceNumber = count;
1918
1919	while (gb_species) {
1920	int myProtection = GB_read_security_write(GBT_read_sequence(gb_species, alignment));
1921	error = 0;
1922	if (myProtection<=maxProtection) {
1923	error = alignTo(gb_species, alignment, &referenceFastSeq, gb_reference, max_seq_length, ali_params);
1924	}
1925	else {
1926	wasNotAllowedToAlign++;
1927	}
1928
1929	if (error) {
1930	aw_message(error);
1931	error = 0;
1932	err_count++;
1933	GB_abort_transaction(GLOBAL_gb_main);
1934	}
1935	else {
1936	GB_pop_transaction(GLOBAL_gb_main);
1937	}
1938
1939	done++;
1940	GB_status(double(done)/double(count));
1941
1942	GB_push_transaction(GLOBAL_gb_main);
1943	gb_species = GBT_next_marked_species(gb_species);
1944	}
1945	break;
1946	}
1947	case FA_SELECTED: { // align all selected species
1948	int count;
1949	int done = 0;
1950	GBDATA *gb_species = get_first_selected_species(&count);
1951
1952	currentSequenceNumber = 1;
1953	overallSequenceNumber = count;
1954
1955	if (!gb_species) {
1956	aw_message("There is no selected species!");
1957	}
1958
1959	while (gb_species) {
1960	int myProtection = GB_read_security_write(GBT_read_sequence(gb_species, alignment));
1961
1962	error = 0;
1963	if (myProtection<=maxProtection) {
1964	error = alignTo(gb_species, alignment, &referenceFastSeq, gb_reference, max_seq_length, ali_params);
1965	}
1966	else {
1967	wasNotAllowedToAlign++;
1968	}
1969
1970	if (error) {
1971	aw_message(error);
1972	error = 0;
1973	err_count++;
1974	GB_abort_transaction(GLOBAL_gb_main);
1975	}
1976	else {
1977	GB_pop_transaction(GLOBAL_gb_main);
1978	}
1979
1980	done++;
1981	GB_status(double(done)/double(count));
1982
1983	GB_push_transaction(GLOBAL_gb_main);
1984	gb_species = get_next_selected_species();
1985	}
1986
1987	break;
1988	}
1989	}
1990
1991	delete referenceSeq;
1992	}
1993	}
1994	}
1995	else if (get_consensus) { // align to group consensi
1996	switch (alignWhat) {
1997	case FA_CURRENT: { // align one sequence
1998	awtc_assert(toalign);
1999	GBDATA *gb_toalign = GBT_find_species_rel_species_data(gb_species_data, toalign);
2000
2001	currentSequenceNumber = overallSequenceNumber = 1;
2002
2003	if (!gb_toalign) {
2004	error = species_not_found(toalign);
2005	}
2006	else {
2007	int myProtection = GB_read_security_write(GBT_read_sequence(gb_toalign, alignment));
2008
2009	error = 0;
2010	if (myProtection<=maxProtection) {
2011	error = alignToGroupConsensus(gb_toalign, alignment, get_consensus, max_seq_length, ali_params);
2012	}
2013	else {
2014	wasNotAllowedToAlign++;
2015	}
2016	}
2017	break;
2018	}
2019	case FA_MARKED: { // align all marked sequences
2020	int count = GBT_count_marked_species(GLOBAL_gb_main);
2021	int done = 0;
2022	GBDATA *gb_species = GBT_first_marked_species_rel_species_data(gb_species_data);
2023
2024	currentSequenceNumber = 1;
2025	overallSequenceNumber = count;
2026
2027	while (gb_species) {
2028	int myProtection = GB_read_security_write(GBT_read_sequence(gb_species, alignment));
2029
2030	error = 0;
2031	if (myProtection<=maxProtection) {
2032	error = alignToGroupConsensus(gb_species, alignment, get_consensus, max_seq_length, ali_params);
2033	}
2034	else {
2035	wasNotAllowedToAlign++;
2036	}
2037
2038	if (error) {
2039	aw_message(error);
2040	error = 0;
2041	err_count++;
2042	GB_abort_transaction(GLOBAL_gb_main);
2043	}
2044	else {
2045	GB_pop_transaction(GLOBAL_gb_main);
2046	}
2047
2048	done++;
2049	GB_status(double(done)/double(count));
2050
2051	GB_push_transaction(GLOBAL_gb_main);
2052	gb_species = GBT_next_marked_species(gb_species);
2053	}
2054	break;
2055	}
2056	case FA_SELECTED: { // align all selected species
2057	int count;
2058	int done = 0;
2059	GBDATA *gb_species = get_first_selected_species(&count);
2060
2061	currentSequenceNumber = 1;
2062	overallSequenceNumber = count;
2063
2064	if (!gb_species) {
2065	aw_message("There is no selected species!");
2066	}
2067
2068	while (gb_species) {
2069	int myProtection = GB_read_security_write(GBT_read_sequence(gb_species, alignment));
2070
2071	if (myProtection<=maxProtection) {
2072	error = alignToGroupConsensus(gb_species, alignment, get_consensus, max_seq_length, ali_params);
2073	}
2074	else {
2075	wasNotAllowedToAlign++;
2076	}
2077
2078	if (error) {
2079	aw_message(error);
2080	error = 0;
2081	err_count++;
2082	GB_abort_transaction(GLOBAL_gb_main);
2083	}
2084	else {
2085	GB_pop_transaction(GLOBAL_gb_main);
2086	}
2087
2088	done++;
2089	GB_status(double(done)/double(count));
2090
2091	GB_push_transaction(GLOBAL_gb_main);
2092	gb_species = get_next_selected_species();
2093	}
2094	break;
2095	}
2096	}
2097	}
2098	else { // align to next relative
2099	switch (alignWhat) {
2100	case FA_CURRENT: { // align one sequence
2101	awtc_assert(toalign);
2102	GBDATA *gb_toalign = GBT_find_species_rel_species_data(gb_species_data, toalign);
2103
2104	currentSequenceNumber = overallSequenceNumber = 1;
2105
2106	if (!gb_toalign) {
2107	error = species_not_found(toalign);
2108	}
2109	else {
2110	int myProtection = GB_read_security_write(GBT_read_sequence(gb_toalign, alignment));
2111
2112	error = 0;
2113	if (myProtection<=maxProtection) {
2114	error = alignToNextRelative(relSearch, max_seq_length, turnAllowed, alignment, gb_toalign, ali_params);
2115	}
2116	else {
2117	wasNotAllowedToAlign++;
2118	}
2119	}
2120	break;
2121	}
2122	case FA_MARKED: { // align all marked sequences
2123	int count = GBT_count_marked_species(GLOBAL_gb_main);
2124	int done = 0;
2125	GBDATA *gb_species = GBT_first_marked_species_rel_species_data(gb_species_data);
2126
2127	currentSequenceNumber = 1;
2128	overallSequenceNumber = count;
2129
2130	while (gb_species) {
2131	int myProtection = GB_read_security_write(GBT_read_sequence(gb_species, alignment));
2132
2133	error = 0;
2134	if (myProtection<=maxProtection) {
2135	error = alignToNextRelative(relSearch, max_seq_length, turnAllowed, alignment, gb_species, ali_params);
2136	}
2137	else {
2138	wasNotAllowedToAlign++;
2139	}
2140
2141	if (error) {
2142	aw_message(error);
2143	error = 0;
2144	err_count++;
2145	GB_abort_transaction(GLOBAL_gb_main);
2146	}
2147	else {
2148	GB_pop_transaction(GLOBAL_gb_main);
2149	}
2150
2151	done++;
2152	GB_status(double(done)/double(count));
2153
2154	GB_push_transaction(GLOBAL_gb_main);
2155	gb_species = GBT_next_marked_species(gb_species);
2156	}
2157	break;
2158	}
2159	case FA_SELECTED: { // align all selected species
2160	int count;
2161	int done = 0;
2162	GBDATA *gb_species = get_first_selected_species(&count);
2163
2164	currentSequenceNumber = 1;
2165	overallSequenceNumber = count;
2166
2167	if (!gb_species) {
2168	aw_message("There is no selected species!");
2169	}
2170
2171	while (gb_species) {
2172	int myProtection = GB_read_security_write(GBT_read_sequence(gb_species, alignment));
2173
2174	if (myProtection<=maxProtection) {
2175	error = alignToNextRelative(relSearch, max_seq_length, turnAllowed, alignment, gb_species, ali_params);
2176	}
2177	else {
2178	wasNotAllowedToAlign++;
2179	}
2180
2181	if (error) {
2182	aw_message(error);
2183	error = 0;
2184	err_count++;
2185	GB_abort_transaction(GLOBAL_gb_main);
2186	}
2187	else {
2188	GB_pop_transaction(GLOBAL_gb_main);
2189	}
2190
2191	done++;
2192	GB_status(double(done)/double(count));
2193
2194	GB_push_transaction(GLOBAL_gb_main);
2195	gb_species = get_next_selected_species();
2196	}
2197	break;
2198	}
2199	}
2200	}
2201	}
2202
2203	AWTC_ClustalV_align_exit();
2204	unaligned_bases.clear();
2205
2206	error = GB_end_transaction(GLOBAL_gb_main, error);
2207
2208	if (wasNotAllowedToAlign>0) {
2209	const char *mess = GB_export_errorf("%i species were not aligned (because of protection level)", wasNotAllowedToAlign);
2210	aw_popup_ok(mess);
2211	}
2212
2213	if (err_count) {
2214	if (error) aw_message(error);
2215	error = GB_export_errorf("Aligner produced %i error%c", err_count, err_count==1 ? '\0' : 's');
2216	}
2217
2218	return error;
2219	}
2220
2221	// --------------------------------------------------------------------------------
2222	// Parameter-Window for Fast-Aligner
2223	// --------------------------------------------------------------------------------
2224
2225	void AWTC_start_faligning(AW_window *aw, AW_CL cd2)
2226	{
2227	AW_root *root = aw->get_root();
2228	char *reference = NULL; // align against next relatives
2229	char *toalign = NULL; // align marked species
2230	GB_ERROR error = NULL;
2231	static struct AWTC_faligner_cd cd = (struct AWTC_faligner_cd )cd2;
2232	int get_consensus = 0;
2233	int pt_server_id = -1;
2234
2235	AWTC_get_first_selected_species get_first_selected_species = 0;
2236	AWTC_get_next_selected_species get_next_selected_species = 0;
2237
2238	awtc_assert(island_hopper == 0);
2239	if (root->awar(FA_AWAR_USE_ISLAND_HOPPING)->read_int()) {
2240	island_hopper = new IslandHopping();
2241	if (root->awar(FA_AWAR_USE_SECONDARY)->read_int()) {
2242	if (cd->helix_string) island_hopper->set_helix(cd->helix_string);
2243	else error = "Warning: No HELIX found. Can't use secondary structure";
2244	}
2245
2246	if (!error) {
2247	island_hopper->set_parameters(root->awar(FA_AWAR_ESTIMATE_BASE_FREQ)->read_int(),
2248	root->awar(FA_AWAR_BASE_FREQ_T)->read_float(),
2249	root->awar(FA_AWAR_BASE_FREQ_C)->read_float(),
2250	root->awar(FA_AWAR_BASE_FREQ_A)->read_float(),
2251	root->awar(FA_AWAR_BASE_FREQ_C)->read_float(),
2252	root->awar(FA_AWAR_SUBST_PARA_CT)->read_float(),
2253	root->awar(FA_AWAR_SUBST_PARA_AT)->read_float(),
2254	root->awar(FA_AWAR_SUBST_PARA_GT)->read_float(),
2255	root->awar(FA_AWAR_SUBST_PARA_AC)->read_float(),
2256	root->awar(FA_AWAR_SUBST_PARA_CG)->read_float(),
2257	root->awar(FA_AWAR_SUBST_PARA_AG)->read_float(),
2258	root->awar(FA_AWAR_EXPECTED_DISTANCE)->read_float(),
2259	root->awar(FA_AWAR_STRUCTURE_SUPPLEMENT)->read_float(),
2260	root->awar(FA_AWAR_GAP_A)->read_float(),
2261	root->awar(FA_AWAR_GAP_B)->read_float(),
2262	root->awar(FA_AWAR_GAP_C)->read_float(),
2263	root->awar(FA_AWAR_THRESHOLD)->read_float()
2264	);
2265	}
2266	}
2267
2268	FA_alignTarget alignWhat = static_cast<FA_alignTarget>(root->awar(FA_AWAR_TO_ALIGN)->read_int());
2269	if (!error) {
2270	switch (alignWhat) {
2271	case FA_CURRENT: { // align current species
2272	toalign = root->awar(AWAR_SPECIES_NAME)->read_string();
2273	break;
2274	}
2275	case FA_MARKED: { // align marked species
2276	break;
2277	}
2278	case FA_SELECTED: { // align selected species
2279	get_first_selected_species = cd->get_first_selected_species;
2280	get_next_selected_species = cd->get_next_selected_species;
2281	break;
2282	}
2283	default: {
2284	awtc_assert(0);
2285	break;
2286	}
2287	}
2288
2289	switch (static_cast<FA_reference>(root->awar(FA_AWAR_REFERENCE)->read_int())) {
2290	case FA_REF_EXPLICIT: // align against specified species
2291
2292	reference = root->awar(FA_AWAR_REFERENCE_NAME)->read_string();
2293	break;
2294
2295	case FA_REF_CONSENSUS: // align against group consensus
2296
2297	if (cd->get_group_consensus) {
2298	get_consensus = 1;
2299	}
2300	else {
2301	error = "Can't get group consensus here.";
2302	}
2303	break;
2304
2305	case FA_REF_RELATIVES: // align against species searched via pt_server
2306
2307	pt_server_id = root->awar(AWAR_PT_SERVER)->read_int();
2308	if (pt_server_id<0) {
2309	error = "No pt_server selected";
2310	}
2311	break;
2312
2313	default: awtc_assert(0);
2314	break;
2315	}
2316	}
2317
2318	int firstColumn = 0;
2319	int lastColumn = -1;
2320
2321	if (!error) {
2322	switch (static_cast<FA_range>(root->awar(FA_AWAR_RANGE)->read_int())) {
2323	case FA_WHOLE_SEQUENCE:
2324	break;
2325	case FA_AROUND_CURSOR: {
2326	int curpos = root->awar(AWAR_CURSOR_POSITION_LOCAL)->read_int();
2327	int size = root->awar(FA_AWAR_AROUND)->read_int();
2328
2329	if ((curpos-size)>0) firstColumn = curpos-size;
2330	lastColumn = curpos+size;
2331	break;
2332	}
2333	case FA_SELECTED_RANGE: {
2334	if (!cd->get_selected_range(&firstColumn, &lastColumn)) {
2335	error = "There is no selected species!";
2336	}
2337	#ifdef DEBUG
2338	else {
2339	printf("Selected range: %i .. %i\n", firstColumn, lastColumn);
2340	}
2341	#endif
2342	break;
2343	}
2344	default: { awtc_assert(0); break; }
2345	}
2346	}
2347
2348	if (!error) {
2349	char *editor_alignment = 0;
2350	{
2351	GB_transaction dummy(GLOBAL_gb_main);
2352	char *default_alignment = GBT_get_default_alignment(GLOBAL_gb_main);
2353
2354	editor_alignment = root->awar_string(AWAR_EDITOR_ALIGNMENT, default_alignment)->read_string();
2355	free(default_alignment);
2356	}
2357	char *pt_server_alignment = root->awar(FA_AWAR_PT_SERVER_ALIGNMENT)->read_string();
2358
2359	if (island_hopper) {
2360	island_hopper->set_range(firstColumn, lastColumn);
2361	firstColumn = 0;
2362	lastColumn = -1;
2363	}
2364
2365	struct SearchRelativeParams relSearch = {
2366	pt_server_id,
2367	pt_server_alignment,
2368	root->awar(FA_AWAR_NEXT_RELATIVES)->read_int(),
2369	root->awar(AWAR_NN_OLIGO_LEN)->read_int(),
2370	root->awar(AWAR_NN_MISMATCHES)->read_int(),
2371	root->awar(AWAR_NN_FAST_MODE)->read_int(),
2372	root->awar(AWAR_NN_REL_MATCHES)->read_int(),
2373	};
2374
2375	struct AlignParams ali_params = {
2376	static_cast<FA_report>(root->awar(FA_AWAR_REPORT)->read_int()),
2377	root->awar(FA_AWAR_SHOW_GAPS_MESSAGES)->read_int(),
2378	firstColumn,
2379	lastColumn,
2380	};
2381
2382	{
2383	GB_transaction ta(GLOBAL_gb_main);
2384	aw_openstatus("FastAligner");
2385	error = AWTC_aligner(
2386	alignWhat,
2387	editor_alignment,
2388	toalign,
2389	get_first_selected_species,
2390	get_next_selected_species,
2391
2392	reference,
2393	get_consensus ? cd->get_group_consensus : NULL,
2394	relSearch,
2395
2396	static_cast<FA_turn>(root->awar(FA_AWAR_MIRROR)->read_int()),
2397	ali_params,
2398	root->awar(FA_AWAR_PROTECTION)->read_int()
2399	);
2400	aw_closestatus();
2401	error = ta.close(error);
2402	}
2403
2404	free(pt_server_alignment);
2405	free(editor_alignment);
2406	}
2407
2408	if (island_hopper) {
2409	delete island_hopper;
2410	island_hopper = 0;
2411	}
2412
2413	if (toalign) free(toalign);
2414	if (error) aw_message(error);
2415	if (cd->do_refresh) cd->refresh_display();
2416	}
2417
2418
2419
2420	void AWTC_create_faligner_variables(AW_root *root,AW_default db1)
2421	{
2422	root->awar_string(FA_AWAR_REFERENCE_NAME, "<undef>", db1);
2423
2424	root->awar_int(FA_AWAR_TO_ALIGN, FA_CURRENT, db1);
2425	root->awar_int(FA_AWAR_REFERENCE, FA_REF_CONSENSUS, db1);
2426	root->awar_int(FA_AWAR_RANGE, FA_WHOLE_SEQUENCE, db1);
2427
2428	#if defined(DEVEL_RALF)
2429	root->awar_int(FA_AWAR_PROTECTION, 0, db1)->write_int(6);
2430	#else
2431	root->awar_int(FA_AWAR_PROTECTION, 0, db1)->write_int(0);
2432	#endif
2433
2434	root->awar_int(FA_AWAR_AROUND, 25, db1);
2435	root->awar_int(FA_AWAR_MIRROR, FA_TURN_INTERACTIVE, db1);
2436	root->awar_int(FA_AWAR_REPORT, FA_NO_REPORT, db1);
2437	root->awar_int(FA_AWAR_SHOW_GAPS_MESSAGES, 1, db1);
2438	root->awar_int(FA_AWAR_USE_SECONDARY, 0, db1);
2439	root->awar_int(AWAR_PT_SERVER, 0, db1);
2440	root->awar_int(FA_AWAR_NEXT_RELATIVES, 1, db1)->set_minmax(1,100);
2441
2442	root->awar_string(FA_AWAR_PT_SERVER_ALIGNMENT, root->awar(AWAR_DEFAULT_ALIGNMENT)->read_char_pntr(), 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 AWTC_awar_set_current_sequence(AW_root *root, AW_default db1)
2474	{
2475	root->awar_int(FA_AWAR_TO_ALIGN, FA_CURRENT, db1);
2476	}
2477
2478	// sets the aligner reference species to current species
2479	void AWTC_set_reference_species_name(AW_window /aww*/, AW_CL cl_AW_root)
2480	{
2481	AW_root root = (AW_root)cl_AW_root;
2482	char *specName = root->awar(AWAR_SPECIES_NAME)->read_string();
2483
2484	root->awar(FA_AWAR_REFERENCE_NAME)->write_string(specName);
2485	free(specName);
2486	}
2487
2488	AW_window AWTC_create_island_hopping_window(AW_root root, AW_CL ) {
2489	AW_window_simple *aws = new AW_window_simple;
2490
2491	aws->init( root, "ISLAND_HOPPING_PARA", "Parameters for Island Hopping");
2492	aws->load_xfig("awtc/islandhopping.fig");
2493
2494	aws->at( "close" );
2495	aws->callback ( (AW_CB0)AW_POPDOWN );
2496	aws->create_button( "CLOSE", "CLOSE", "O" );
2497
2498	aws->at( "help" );
2499	aws->callback ( AW_POPUP_HELP, (AW_CL) "islandhopping.hlp" );
2500	aws->create_button( "HELP", "HELP" );
2501
2502	aws->at("use_secondary");
2503	aws->label("Use secondary structure (only for re-align)");
2504	aws->create_toggle(FA_AWAR_USE_SECONDARY);
2505
2506	aws->at("freq");
2507	aws->create_toggle_field(FA_AWAR_ESTIMATE_BASE_FREQ,"Base freq.","B");
2508	aws->insert_default_toggle("Estimate","E",1);
2509	aws->insert_toggle("Define here: ","D",0);
2510	aws->update_toggle_field();
2511
2512	aws->at("freq_a"); aws->label("A:"); aws->create_input_field(FA_AWAR_BASE_FREQ_A, 4);
2513	aws->at("freq_c"); aws->label("C:"); aws->create_input_field(FA_AWAR_BASE_FREQ_C, 4);
2514	aws->at("freq_g"); aws->label("G:"); aws->create_input_field(FA_AWAR_BASE_FREQ_G, 4);
2515	aws->at("freq_t"); aws->label("T:"); aws->create_input_field(FA_AWAR_BASE_FREQ_T, 4);
2516
2517	int xpos[4], ypos[4];
2518	{
2519	aws->button_length(1);
2520
2521	int dummy;
2522	aws->at("h_a"); aws->get_at_position(&xpos[0], &dummy); aws->create_button("A","A");
2523	aws->at("h_c"); aws->get_at_position(&xpos[1], &dummy); aws->create_button("C","C");
2524	aws->at("h_g"); aws->get_at_position(&xpos[2], &dummy); aws->create_button("G","G");
2525	aws->at("h_t"); aws->get_at_position(&xpos[3], &dummy); aws->create_button("T","T");
2526
2527	aws->at("v_a"); aws->get_at_position(&dummy, &ypos[0] ); aws->create_button("A","A");
2528	aws->at("v_c"); aws->get_at_position(&dummy, &ypos[1] ); aws->create_button("C","C");
2529	aws->at("v_g"); aws->get_at_position(&dummy, &ypos[2] ); aws->create_button("G","G");
2530	aws->at("v_t"); aws->get_at_position(&dummy, &ypos[3] ); aws->create_button("T","T");
2531	}
2532
2533	aws->at("subst"); aws->create_button("subst_para", "Substitution rate parameters:");
2534
2535	#define XOFF -25
2536	#define YOFF 0
2537
2538	aws->at(xpos[1]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AC, 4);
2539	aws->at(xpos[2]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AG, 4);
2540	aws->at(xpos[3]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AT, 4);
2541	aws->at(xpos[2]+XOFF, ypos[1]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_CG, 4);
2542	aws->at(xpos[3]+XOFF, ypos[1]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_CT, 4);
2543	aws->at(xpos[3]+XOFF, ypos[2]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_GT, 4);
2544
2545	#undef XOFF
2546	#undef YOFF
2547
2548	aws->label_length(22);
2549
2550	aws->at("dist");
2551	aws->label("Expected distance");
2552	aws->create_input_field(FA_AWAR_EXPECTED_DISTANCE, 5);
2553
2554	aws->at("supp");
2555	aws->label("Structure supplement");
2556	aws->create_input_field(FA_AWAR_STRUCTURE_SUPPLEMENT, 5);
2557
2558	aws->at("thres");
2559	aws->label("Threshold");
2560	aws->create_input_field(FA_AWAR_THRESHOLD, 5);
2561
2562	aws->label_length(10);
2563
2564	aws->at("gapA");
2565	aws->label("Gap A");
2566	aws->create_input_field(FA_AWAR_GAP_A, 5);
2567
2568	aws->at("gapB");
2569	aws->label("Gap B");
2570	aws->create_input_field(FA_AWAR_GAP_B, 5);
2571
2572	aws->at("gapC");
2573	aws->label("Gap C");
2574	aws->create_input_field(FA_AWAR_GAP_C, 5);
2575
2576	return (AW_window *)aws;
2577	}
2578
2579	AW_window AWTC_create_family_settings_window(AW_root root) {
2580	static AW_window_simple *aws = 0;
2581
2582	if (!aws) {
2583	aws = new AW_window_simple;
2584
2585	aws->init(root, "FAMILY_PARAMS", "Family search paramaters");
2586	aws->load_xfig("awtc/family_settings.fig");
2587
2588	aws->at( "close" );
2589	aws->callback ( (AW_CB0)AW_POPDOWN );
2590	aws->create_button( "CLOSE", "CLOSE", "O" );
2591
2592	aws->at( "help" );
2593	aws->callback ( AW_POPUP_HELP, (AW_CL) "next_neighbours_common.hlp" );
2594	aws->create_button( "HELP", "HELP" );
2595
2596	AWTC_create_common_next_neighbour_fields(aws);
2597	}
2598
2599	return aws;
2600	}
2601
2602	AW_window AWTC_create_faligner_window(AW_root root, AW_CL cd2)
2603	{
2604	AW_window_simple *aws = new AW_window_simple;
2605
2606	aws->init( root, "INTEGRATED_ALIGNERS", INTEGRATED_ALIGNERS_TITLE);
2607	aws->load_xfig("awtc/faligner.fig");
2608
2609	aws->label_length( 10 );
2610	aws->button_length( 10 );
2611
2612	aws->at( "close" );
2613	aws->callback ( (AW_CB0)AW_POPDOWN );
2614	aws->create_button( "CLOSE", "CLOSE", "O" );
2615
2616	aws->at( "help" );
2617	aws->callback ( AW_POPUP_HELP, (AW_CL) "faligner.hlp" );
2618	aws->create_button( "HELP", "HELP" );
2619
2620	aws->at("aligner");
2621	aws->create_toggle_field(FA_AWAR_USE_ISLAND_HOPPING,"Aligner","A");
2622	aws->insert_default_toggle("Fast aligner", "F", 0);
2623	aws->sens_mask(AWM_EXP);
2624	aws->insert_toggle ("Island Hopping", "I", 1);
2625	aws->sens_mask(AWM_ALL);
2626	aws->update_toggle_field();
2627
2628	aws->button_length(12);
2629	aws->at("island_para");
2630	aws->callback(AW_POPUP, (AW_CL)AWTC_create_island_hopping_window, (AW_CL)0);
2631	aws->sens_mask(AWM_EXP);
2632	aws->create_button("island_para", "Parameters", "");
2633	aws->sens_mask(AWM_ALL);
2634
2635	aws->button_length(10);
2636
2637	aws->at( "rev_compl" );
2638	aws->callback(AWTC_build_reverse_complement, cd2);
2639	aws->create_button( "reverse_complement", "Turn now!", "");
2640
2641	aws->at("what");
2642	aws->create_toggle_field(FA_AWAR_TO_ALIGN,"Align what?","A");
2643	aws->insert_toggle ("Current Species:", "A", FA_CURRENT);
2644	aws->insert_default_toggle("Marked Species", "M", FA_MARKED);
2645	aws->insert_toggle ("Selected Species", "S", FA_SELECTED);
2646	aws->update_toggle_field();
2647
2648	aws->at( "swhat" );
2649	aws->create_input_field( AWAR_SPECIES_NAME, 2);
2650
2651	aws->at("against");
2652	aws->create_toggle_field(FA_AWAR_REFERENCE,"Reference","");
2653	aws->insert_toggle ("Species by name:", "S", FA_REF_EXPLICIT);
2654	aws->insert_toggle ("Group consensus", "K", FA_REF_CONSENSUS);
2655	aws->insert_default_toggle("Auto search by pt_server:", "A", FA_REF_RELATIVES);
2656	aws->update_toggle_field();
2657
2658	aws->at( "sagainst" );
2659	aws->create_input_field(FA_AWAR_REFERENCE_NAME, 2);
2660
2661	aws->at("copy");
2662	aws->callback(AWTC_set_reference_species_name, (AW_CL)root);
2663	aws->create_button("Copy", "Copy", "");
2664
2665	aws->label_length(0);
2666	aws->at( "pt_server" );
2667	awt_create_selection_list_on_pt_servers(aws, AWAR_PT_SERVER, true);
2668
2669	aws->at("use_ali");
2670	aws->create_input_field(FA_AWAR_PT_SERVER_ALIGNMENT, 15);
2671
2672	aws->at("relatives");
2673	aws->label("Number of relatives to use:");
2674	aws->create_input_field(FA_AWAR_NEXT_RELATIVES, 3);
2675
2676	aws->at("relSett");
2677	aws->callback(AW_POPUP, (AW_CL)AWTC_create_family_settings_window, (AW_CL)root);
2678	aws->create_button("Settings", "Settings", "");
2679
2680	// Range
2681
2682	aws->label_length( 10 );
2683	aws->at("range");
2684	aws->create_toggle_field(FA_AWAR_RANGE, "Range", "");
2685	aws->insert_default_toggle("Whole sequence", "", FA_WHOLE_SEQUENCE);
2686	aws->insert_toggle ("Positions around cursor: ", "", FA_AROUND_CURSOR);
2687	aws->insert_toggle ("Selected Range", "", FA_SELECTED_RANGE);
2688	aws->update_toggle_field();
2689
2690	aws->at("around");
2691	aws->create_input_field(FA_AWAR_AROUND, 2);
2692
2693	aws->at("protection");
2694	aws->create_option_menu(FA_AWAR_PROTECTION, "Protection", "");
2695	aws->insert_default_option("0", 0, 0);
2696	aws->insert_option("1", 0, 1);
2697	aws->insert_option("2", 0, 2);
2698	aws->insert_option("3", 0, 3);
2699	aws->insert_option("4", 0, 4);
2700	aws->insert_option("5", 0, 5);
2701	aws->insert_option("6", 0, 6);
2702	aws->update_option_menu();
2703
2704	// MirrorCheck
2705
2706	aws->at("mirror");
2707	aws->create_option_menu(FA_AWAR_MIRROR, "Turn check", "");
2708	aws->insert_option ("Never turn sequence", "", FA_TURN_NEVER);
2709	aws->insert_default_option("User acknowledgment", "", FA_TURN_INTERACTIVE);
2710	aws->insert_option ("Automatically turn sequence", "", FA_TURN_ALWAYS);
2711	aws->update_option_menu();
2712
2713	// Report
2714
2715	aws->at("insert");
2716	aws->create_option_menu(FA_AWAR_REPORT, "Report", "");
2717	aws->insert_option ("No report", "", FA_NO_REPORT);
2718	aws->sens_mask(AWM_EXP);
2719	aws->insert_default_option("Report to temporary entries", "", FA_TEMP_REPORT);
2720	aws->insert_option ("Report to resident entries", "", FA_REPORT);
2721	aws->sens_mask(AWM_ALL);
2722	aws->update_option_menu();
2723
2724	aws->at("gaps");
2725	aws->create_toggle(FA_AWAR_SHOW_GAPS_MESSAGES);
2726
2727	aws->at( "align" );
2728	aws->callback ( AWTC_start_faligning, cd2);
2729	aws->highlight();
2730	aws->create_button( "GO", "GO", "G");
2731
2732	return (AW_window *)aws;
2733	}
2734
2735
2736

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

Download in other formats: