Context Navigation

source: tags/old_import_filter/MERGE/MG_adapt_ali.cxx

Visit:

Last change on this file was 10113, checked in by westram, 11 years ago
changed dynamic value-initializations into default-initializations i believe/know they were used with ignorance of ctor-behavior.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 37.1 KB

Line
1	// =============================================================== //
2	// //
3	// File : MG_adapt_ali.cxx //
4	// Purpose : //
5	// //
6	// Institute of Microbiology (Technical University Munich) //
7	// http://www.arb-home.de/ //
8	// //
9	// =============================================================== //
10
11	#include "merge.hxx"
12	#include "MG_adapt_ali.hxx"
13
14	#include <aw_msg.hxx>
15	#include <arbdbt.h>
16	#include <arb_strbuf.h>
17
18	#include <algorithm>
19	#include <list>
20
21	#if defined(DEBUG)
22	// #define DUMP_MAPPING
23	// #define DUMP_SOFTMAPPING
24	#endif // DEBUG
25
26	// -----------------
27	// MG_remap
28
29	const int NO_POSITION = -1;
30	const int LEFT_BORDER = -2;
31	const int RIGHT_BORDER = -3;
32
33	struct softbase {
34	char base;
35	int origin; // position in source alignment
36	char last_gapchar; // last gap seen before base
37	int targetpos; // target position
38
39	softbase(char base_, int origin_, char last_gapchar_)
40	: base(base_)
41	, origin(origin_)
42	, last_gapchar(last_gapchar_)
43	, targetpos(NO_POSITION)
44	{
45	mg_assert(last_gapchar);
46	}
47
48	operator char () const { return base; }
49	};
50
51	typedef std::list<softbase> softbaselist;
52	typedef softbaselist::iterator softbaseiter;
53	typedef softbaselist::const_iterator const_softbaseiter;
54
55	class MG_remap : virtual Noncopyable {
56	int in_length;
57	int out_length;
58	int *remap_tab; // fixed mapping (targetPosition or NO_POSITION)
59
60	// soft-mapping:
61	int *soft_remap_tab; // soft-mapping (NO_POSITION, targetPosition, LEFT_BORDER or RIGHT_BORDER)
62
63	char *calc_softmapping(softbaselist& softbases, int start, int end, int &outlen);
64	int softmap_to(softbaselist& softbases, int start, int end, GBS_strstruct *outs);
65
66	bool have_softmapping() const { return soft_remap_tab; }
67	void create_softmapping();
68	void forget_softmapping() {
69	delete [] soft_remap_tab;
70	soft_remap_tab = NULL;
71	}
72
73	static int build_initial_mapping(const char iref, int ilen, const char *oref, int olen);
74	void merge_mapping(MG_remap &other, int& inconsistent, int& added);
75
76	public:
77
78	MG_remap()
79	: in_length(0)
80	, out_length(0)
81	, remap_tab(NULL)
82	, soft_remap_tab(NULL)
83	{}
84	~MG_remap() {
85	forget_softmapping();
86	delete [] remap_tab;
87	}
88
89	const char add_reference(const char in_reference, const char *out_reference); // returns only warnings
90	char remap(const char sequence); // returns 0 on error, else copy of sequence
91
92	char *readable_inconsistent_positions();
93
94	#if defined(DUMP_MAPPING)
95	static void dump(const int data, int len, const char comment, int dontShow) {
96	fflush(stdout);
97	fflush(stderr);
98	fputc('>', stdout);
99	int digits = log10(len)+2;
100	for (int pos = 0; pos<len; ++pos) {
101	if (data[pos] == dontShow) {
102	fprintf(stdout, "%*s", digits, "_");
103	}
104	else {
105	fprintf(stdout, "%*i", digits, data[pos]);
106	}
107	}
108	fprintf(stdout, " (%s)\n", comment);
109	fflush(stdout);
110	}
111	void dump_remap(const char *comment) { dump(remap_tab, in_length, comment, NO_POSITION); }
112	#endif // DUMP_MAPPING
113	};
114
115	int MG_remap::build_initial_mapping(const char iref, int ilen, const char *oref, int olen) {
116	int *remap = new int[ilen];
117
118	const char *spacers = "-. n";
119
120	int ipos = 0;
121	int opos = 0;
122
123	while (ipos<ilen && opos<olen) {
124	size_t ispaces = strspn(iref+ipos, spacers);
125	size_t ospaces = strspn(oref+opos, spacers);
126
127	while (ispaces && ipos<ilen) {
128	remap[ipos++] = NO_POSITION;
129	ispaces--;
130	}
131	opos += ospaces;
132	if (ipos<ilen && opos<olen) remap[ipos++] = opos++;
133	}
134	while (ipos<ilen) remap[ipos++] = NO_POSITION;
135
136	return remap;
137	}
138
139	#if defined(DEBUG)
140	// #define DUMP_INCONSISTENCIES
141	#endif // DEBUG
142
143
144	void MG_remap::merge_mapping(MG_remap &other, int& inconsistent, int& added) {
145	const int *primary = remap_tab;
146	int *secondary = other.remap_tab;
147	bool increased_len = false;
148
149	if (other.in_length>in_length) {
150	// re-use memory of bigger map
151	std::swap(other.in_length, in_length);
152	std::swap(other.remap_tab, remap_tab);
153
154	increased_len = true;
155	}
156	out_length = std::max(out_length, other.out_length); // remember biggest output length
157
158	int mixlen = std::min(in_length, other.in_length);
159
160	// eliminate inconsistant positions from secondary mapping (forward)
161	inconsistent = 0;
162	{
163	int max_target_pos = 0;
164	for (int pos = 0; pos<mixlen; ++pos) {
165	max_target_pos = std::max(max_target_pos, primary[pos]);
166	if (secondary[pos]<max_target_pos) {
167	if (secondary[pos] != NO_POSITION) {
168	#if defined(DUMP_INCONSISTENCIES)
169	fprintf(stderr, "merge-inconsistency: pos=%i primary[]=%i secondary[]=%i max_target_pos=%i\n",
170	pos, primary[pos], secondary[pos], max_target_pos);
171	#endif // DUMP_INCONSISTENCIES
172	inconsistent++;
173	}
174	secondary[pos] = NO_POSITION; // consistency error -> ignore position
175	}
176	}
177	}
178	// (backward)
179	{
180	int min_target_pos = out_length-1;
181	for (int pos = mixlen-1; pos >= 0; --pos) {
182	if (primary[pos] >= 0 && primary[pos]<min_target_pos) min_target_pos = primary[pos];
183	if (secondary[pos] > min_target_pos) {
184	#if defined(DUMP_INCONSISTENCIES)
185	fprintf(stderr, "merge-inconsistency: pos=%i primary[]=%i secondary[]=%i min_target_pos=%i\n",
186	pos, primary[pos], secondary[pos], min_target_pos);
187	#endif // DUMP_INCONSISTENCIES
188	inconsistent++;
189	secondary[pos] = NO_POSITION; // consistency error -> ignore position
190	}
191	}
192	}
193
194	// merge mappings
195	added = 0;
196	for (int pos = 0; pos < mixlen; ++pos) {
197	if (primary[pos] == NO_POSITION) {
198	remap_tab[pos] = secondary[pos];
199	added += (remap_tab[pos] != NO_POSITION);
200	}
201	else {
202	remap_tab[pos] = primary[pos];
203	}
204	// remap_tab[pos] = primary[pos] == NO_POSITION ? secondary[pos] : primary[pos];
205	mg_assert(remap_tab[pos]<out_length);
206	}
207
208	if (increased_len) { // count positions appended at end of sequence
209	for (int pos = other.in_length; pos<in_length; ++pos) {
210	added += (remap_tab[pos] != NO_POSITION);
211	}
212	}
213
214	// Note: copying the rest from larger mapping is not necessary
215	// (it's already there, because of memory-reuse)
216	}
217
218	const char MG_remap::add_reference(const char in_reference, const char *out_reference) {
219	// returns warning about inconsistencies/useless reference
220	const char *warning = NULL;
221
222	if (have_softmapping()) forget_softmapping();
223
224	if (!remap_tab) {
225	in_length = strlen(in_reference);
226	out_length = strlen(out_reference);
227	remap_tab = build_initial_mapping(in_reference, in_length, out_reference, out_length);
228	#if defined(DUMP_MAPPING)
229	dump_remap("initial");
230	#endif // DUMP_MAPPING
231	}
232	else {
233	MG_remap tmp;
234	tmp.add_reference(in_reference, out_reference);
235
236	int inconsistent, added;
237	merge_mapping(tmp, inconsistent, added);
238
239	if (inconsistent>0) warning = GBS_global_string("contains %i inconsistent positions", inconsistent);
240	if (added<1) {
241	const char *useless_warning = "doesn't add useful information";
242	if (warning) warning = GBS_global_string("%s and %s", warning, useless_warning);
243	else warning = useless_warning;
244	}
245
246	#if defined(DUMP_MAPPING)
247	dump_remap("merged");
248	#endif // DUMP_MAPPING
249	}
250
251	return warning;
252	}
253
254	void MG_remap::create_softmapping() {
255	soft_remap_tab = new int[in_length];
256
257	int last_fixed_position = NO_POSITION;
258	int pos;
259	for (pos = 0; pos<in_length && last_fixed_position == NO_POSITION; ++pos) {
260	if (remap_tab[pos] == NO_POSITION) {
261	soft_remap_tab[pos] = LEFT_BORDER;
262	}
263	else {
264	soft_remap_tab[pos] = NO_POSITION;
265	last_fixed_position = pos;
266	}
267	}
268	if (last_fixed_position != NO_POSITION) {
269	for ( ; pos<in_length; ++pos) {
270	if (remap_tab[pos] != NO_POSITION) {
271	int softstart = last_fixed_position+1;
272	int softsize = pos-softstart;
273
274	if (softsize>0) {
275	int target_softstart = remap_tab[last_fixed_position]+1;
276	int target_softsize = remap_tab[pos]-target_softstart;
277
278	double target_step;
279	if (softsize>1 && target_softsize>1) {
280	target_step = double(target_softsize-1)/(softsize-1);
281	}
282	else {
283	target_step = 0.0;
284	}
285
286	if (target_step >= 1.0 && target_softsize>softsize) {
287	// target range > source range -> split softmapping in the middle
288	int halfsoftsize = softsize/2;
289	int target_softpos = softstart;
290	int off;
291	for (off = 0; off<halfsoftsize; ++off) {
292	soft_remap_tab[softstart+off] = target_softpos++;
293	}
294	target_softpos += target_softsize-softsize;
295	for (; off<softsize; ++off) {
296	soft_remap_tab[softstart+off] = target_softpos++;
297	}
298	}
299	else {
300	double target_softpos = target_softstart;
301	for (int off = 0; off<softsize; ++off) {
302	soft_remap_tab[softstart+off] = int(target_softpos+0.5);
303	target_softpos += target_step;
304	}
305	}
306	}
307	last_fixed_position = pos;
308	soft_remap_tab[pos] = NO_POSITION;
309	}
310	}
311
312	for (--pos; pos>last_fixed_position; --pos) {
313	soft_remap_tab[pos] = RIGHT_BORDER;
314	}
315	}
316
317	#if defined(DUMP_MAPPING)
318	dump(soft_remap_tab, in_length, "softmap", -1);
319	#endif // DUMP_MAPPING
320	}
321
322	static void drop_dots(softbaselist& softbases, int excessive_positions) {
323	// drop consecutive dots
324	bool justseendot = false;
325
326	softbaseiter next = softbases.begin();
327	while (excessive_positions && next != softbases.end()) {
328	bool isdot = (next->base == '.');
329	if (isdot && justseendot) {
330	excessive_positions--;
331	next = softbases.erase(next);
332	}
333	else {
334	justseendot = isdot;
335	++next;
336	}
337	}
338
339	if (excessive_positions) {
340	// drop single dots
341	next = softbases.begin();
342	while (excessive_positions && next != softbases.end()) {
343	if (next->base == '.') {
344	next = softbases.erase(next);
345	}
346	else {
347	++next;
348	}
349	}
350	}
351	}
352
353	#if defined(DUMP_SOFTMAPPING)
354	static char *softbaselist_2_string(const softbaselist& softbases) {
355	GBS_strstruct out = GBS_stropen(softbases.size()100);
356
357	for (const_softbaseiter base = softbases.begin(); base != softbases.end(); ++base) {
358	const char *info = GBS_global_string(" %c'%c' %i->",
359	base->last_gapchar ? base->last_gapchar : ' ',
360	base->base,
361	base->origin);
362	GBS_strcat(out, info);
363	if (base->targetpos == NO_POSITION) {
364	GBS_strcat(out, "NP");
365	}
366	else {
367	GBS_intcat(out, base->targetpos);
368	}
369	}
370
371	return GBS_strclose(out);
372	}
373	#endif // DUMP_SOFTMAPPING
374
375	char *MG_remap::calc_softmapping(softbaselist& softbases, int start, int end, int& outlen) {
376	//! tries to map all bases(+dots) in 'softbases' into range [start, end]
377	//! @return heap-copy of range-content (may be oversized, if too many bases in list)
378
379	int wanted_size = end-start+1;
380	int listsize = softbases.size();
381
382	#if defined(DUMP_SOFTMAPPING)
383	char *sbl_initial = softbaselist_2_string(softbases);
384	char *sbl_exdots = NULL;
385	char *sbl_target = NULL;
386	char *sbl_exclash = NULL;
387	#endif // DUMP_SOFTMAPPING
388
389	if (listsize>wanted_size) {
390	int excessive_positions = listsize-wanted_size;
391	drop_dots(softbases, excessive_positions);
392	listsize = softbases.size();
393
394	#if defined(DUMP_SOFTMAPPING)
395	sbl_exdots = softbaselist_2_string(softbases);
396	#endif // DUMP_SOFTMAPPING
397	}
398
399	char *result = NULL;
400	if (listsize >= wanted_size) { // not or just enough space -> plain copy
401	result = (char*)malloc(listsize+1);
402	*std::copy(softbases.begin(), softbases.end(), result) = 0;
403	outlen = listsize;
404	}
405	else { // otherwise do soft-mapping
406	result = (char*)malloc(wanted_size+1);
407	mg_assert(listsize < wanted_size);
408
409	// calculate target positions and detect mapping conflicts
410	bool conflicts = false;
411	{
412	int lasttargetpos = NO_POSITION;
413	for (softbaseiter base = softbases.begin(); base != softbases.end(); ++base) {
414	// // int targetpos = calc_softpos(base->origin);
415	int targetpos = soft_remap_tab[base->origin];
416	if (targetpos == lasttargetpos) {
417	mg_assert(targetpos != NO_POSITION);
418	conflicts = true;
419	}
420	base->targetpos = lasttargetpos = targetpos;
421	}
422	}
423
424	#if defined(DUMP_SOFTMAPPING)
425	sbl_target = softbaselist_2_string(softbases);
426	#endif // // DUMP_SOFTMAPPING
427
428	if (conflicts) {
429	int nextpos = end+1;
430	for (softbaselist::reverse_iterator base = softbases.rbegin(); base != softbases.rend(); ++base) {
431	if (base->targetpos >= nextpos) {
432	base->targetpos = nextpos-1;
433	}
434	nextpos = base->targetpos;
435	mg_assert(base->targetpos >= start);
436	mg_assert(base->targetpos <= end);
437	}
438	mg_assert(nextpos >= start);
439
440	#if defined(DUMP_SOFTMAPPING)
441	sbl_exclash = softbaselist_2_string(softbases);
442	#endif // // DUMP_SOFTMAPPING
443	}
444
445	int idx = 0;
446	for (softbaseiter base = softbases.begin(); base != softbases.end(); ++base) {
447	int pos = base->targetpos - start;
448
449	if (idx<pos) {
450	char gapchar = base->last_gapchar;
451	while (idx<pos) result[idx++] = gapchar;
452	}
453	result[idx++] = base->base;
454	}
455	result[idx] = 0;
456	outlen = idx;
457	mg_assert(idx <= wanted_size);
458	}
459
460	#if defined(DUMP_SOFTMAPPING)
461	fflush(stdout);
462	fflush(stderr);
463	printf("initial:%s\n", sbl_initial);
464	if (sbl_exdots) printf("exdots :%s\n", sbl_exdots);
465	if (sbl_target) printf("target :%s\n", sbl_target);
466	if (sbl_exclash) printf("exclash:%s\n", sbl_exclash);
467	printf("calc_softmapping(%i, %i) -> \"%s\"\n", start, end, result);
468	fflush(stdout);
469	free(sbl_exclash);
470	free(sbl_target);
471	free(sbl_exdots);
472	free(sbl_initial);
473	#endif // DUMP_SOFTMAPPING
474
475	return result;
476	}
477
478	int MG_remap::softmap_to(softbaselist& softbases, int start, int end, GBS_strstruct *outs) {
479	int mappedlen;
480	char *softmapped = calc_softmapping(softbases, start, end, mappedlen);
481
482	GBS_strcat(outs, softmapped);
483	free(softmapped);
484	softbases.clear();
485
486	return mappedlen;
487	}
488
489	char MG_remap::remap(const char sequence) {
490	int slen = strlen(sequence);
491	int minlen = std::min(slen, in_length);
492	GBS_strstruct *outs = GBS_stropen(out_length+1);
493	int written = 0;
494	softbaselist softbases;
495	int pos;
496
497	if (!have_softmapping()) create_softmapping();
498
499	// remap left border
500	for (pos = 0; pos<in_length && soft_remap_tab[pos] == LEFT_BORDER; ++pos) {
501	char c = pos<slen ? sequence[pos] : '-';
502	if (c != '-' && c != '.') {
503	softbases.push_back(softbase(c, pos, '.'));
504	}
505	}
506	char last_gapchar = 0;
507	{
508	int bases = softbases.size();
509	if (bases) {
510	int fixed = remap_tab[pos];
511	mg_assert(fixed != NO_POSITION);
512
513	int bases_start = fixed-bases;
514	if (written<bases_start) {
515	GBS_chrncat(outs, '.', bases_start-written);
516	written = bases_start;
517	}
518
519	written += softmap_to(softbases, written, fixed-1, outs);
520	}
521	else {
522	last_gapchar = '.';
523	}
524	}
525
526
527	// remap center (fixed mapping and softmapping inbetween)
528	for (; pos<in_length; ++pos) {
529	char c = pos<slen ? sequence[pos] : '-';
530	int target_pos = remap_tab[pos];
531
532	if (target_pos == NO_POSITION) { // softmap
533	if (c == '-') {
534	if (!last_gapchar) last_gapchar = '-';
535	}
536	else {
537	if (!last_gapchar) last_gapchar = c == '.' ? '.' : '-';
538
539	softbases.push_back(softbase(c, pos, last_gapchar)); // remember bases for softmapping
540	if (c != '.') last_gapchar = 0;
541	}
542	}
543	else { // fixed mapping
544	if (!softbases.empty()) {
545	written += softmap_to(softbases, written, target_pos-1, outs);
546	}
547	if (written<target_pos) {
548	if (!last_gapchar) last_gapchar = c == '.' ? '.' : '-';
549
550	GBS_chrncat(outs, last_gapchar, target_pos-written);
551	written = target_pos;
552	}
553
554	if (c == '-') {
555	if (!last_gapchar) last_gapchar = '-';
556	GBS_chrcat(outs, last_gapchar); written++;
557	}
558	else {
559	GBS_chrcat(outs, c); written++;
560	if (c != '.') last_gapchar = 0;
561	}
562	}
563	}
564
565	// Spool leftover softbases.
566	// Happens when
567	// - sequence ends before last fixed position
568	// - sequence starts after last fixed position
569
570	if (!softbases.empty()) {
571	// softmap_to(softbases, written, written+softbases.size()-1, outs);
572	softmap_to(softbases, written, written, outs);
573	}
574	mg_assert(softbases.empty());
575
576	// copy overlength rest of sequence:
577	const char *gap_chars = "- .";
578	for (int i = minlen; i < slen; i++) {
579	int c = sequence[i];
580	if (!strchr(gap_chars, c)) GBS_chrcat(outs, c); // don't fill with gaps
581	}
582
583	// convert end of seq ('-' -> '.')
584	{
585	char *out = GBS_mempntr(outs);
586	int end = GBS_memoffset(outs);
587
588	for (int off = end-1; off >= 0; --off) {
589	if (out[off] == '-') out[off] = '.';
590	else if (out[off] != '.') break;
591	}
592	}
593
594	return GBS_strclose(outs);
595	}
596
597	// --------------------------------------------------------------------------------
598
599	static MG_remap MG_create_remap(GBDATA gb_left, GBDATA gb_right, const char reference_species_names, const char *alignment_name) {
600	MG_remap *rem = new MG_remap;
601	char *ref_list = strdup(reference_species_names);
602
603	for (char *tok = strtok(ref_list, " \n,;"); tok; tok = strtok(NULL, " \n,;")) {
604	bool is_SAI = strncmp(tok, "SAI:", 4) == 0;
605	GBDATA *gb_species_left = 0;
606	GBDATA *gb_species_right = 0;
607
608	if (is_SAI) {
609	gb_species_left = GBT_find_SAI(gb_left, tok+4);
610	gb_species_right = GBT_find_SAI(gb_right, tok+4);
611	}
612	else {
613	gb_species_left = GBT_find_species(gb_left, tok);
614	gb_species_right = GBT_find_species(gb_right, tok);
615	}
616
617	if (!gb_species_left \|\| !gb_species_right) {
618	aw_message(GBS_global_string("Warning: Couldn't find %s'%s' in %s DB.\nPlease read ADAPT ALIGNMENT section in help file!",
619	is_SAI ? "" : "species ",
620	tok,
621	gb_species_left ? "right" : "left"));
622	}
623	else {
624	// look for sequence/SAI "data"
625	GBDATA *gb_seq_left = GBT_read_sequence(gb_species_left, alignment_name);
626	GBDATA *gb_seq_right = GBT_read_sequence(gb_species_right, alignment_name);
627
628	if (gb_seq_left && gb_seq_right) {
629	GB_TYPES type_left = GB_read_type(gb_seq_left);
630	GB_TYPES type_right = GB_read_type(gb_seq_right);
631
632	if (type_left != type_right) {
633	aw_message(GBS_global_string("Warning: data type of '%s' differs in both databases", tok));
634	}
635	else {
636	const char *warning;
637	if (type_right == GB_STRING) {
638	warning = rem->add_reference(GB_read_char_pntr(gb_seq_left), GB_read_char_pntr(gb_seq_right));
639	}
640	else {
641	char *sleft = GB_read_as_string(gb_seq_left);
642	char *sright = GB_read_as_string(gb_seq_right);
643
644	warning = rem->add_reference(sleft, sright);
645
646	free(sleft);
647	free(sright);
648	}
649
650	if (warning) aw_message(GBS_global_string("Warning: '%s' %s", tok, warning));
651	}
652	}
653	}
654	}
655	free(ref_list);
656
657	return rem;
658	}
659
660	// --------------------------------------------------------------------------------
661
662	MG_remaps::MG_remaps(GBDATA gb_left, GBDATA gb_right, bool enable, const char *reference_species_names)
663	: n_remaps(0),
664	remaps(NULL)
665	{
666	if (enable) { // otherwise nothing will be remapped!
667	GBT_get_alignment_names(alignment_names, gb_left);
668	for (n_remaps = 0; alignment_names[n_remaps]; n_remaps++) {} // count alignments
669
670	remaps = (MG_remap*)GB_calloc(sizeof(remaps), n_remaps);
671	for (int i = 0; i<n_remaps; ++i) {
672	remaps[i] = MG_create_remap(gb_left, gb_right, reference_species_names, alignment_names[i]);
673	}
674	}
675	}
676
677	MG_remaps::~MG_remaps() {
678	for (int i=0; i<n_remaps; i++) delete remaps[i];
679	free(remaps);
680	}
681
682	static GB_ERROR MG_transfer_sequence(MG_remap remap, GBDATA source_species, GBDATA destination_species, const char alignment_name) {
683	// align sequence after copy
684	GB_ERROR error = NULL;
685
686	if (remap) { // shall remap?
687	GBDATA *gb_seq_left = GBT_read_sequence(source_species, alignment_name);
688	GBDATA *gb_seq_right = GBT_read_sequence(destination_species, alignment_name);
689
690	if (gb_seq_left && gb_seq_right) { // if one DB hasn't sequence -> sequence was not copied
691	char *ls = GB_read_string(gb_seq_left);
692	char *rs = GB_read_string(gb_seq_right);
693
694	if (strcmp(ls, rs) == 0) { // if sequences are not identical -> sequence was not copied
695	free(rs);
696	rs = remap->remap(ls);
697
698	if (!rs) error = GB_await_error();
699	else {
700	long old_check = GBS_checksum(ls, 0, ".- ");
701	long new_check = GBS_checksum(rs, 0, ".- ");
702
703	if (old_check == new_check) {
704	error = GB_write_string(gb_seq_right, rs);
705	}
706	else {
707	error = GBS_global_string("Failed to adapt alignment of '%s' (checksum changed)", GBT_read_name(source_species));
708	}
709	}
710	}
711	free(rs);
712	free(ls);
713	}
714	}
715	return error;
716	}
717
718	GB_ERROR MG_transfer_all_alignments(MG_remaps remaps, GBDATA source_species, GBDATA *destination_species) {
719	GB_ERROR error = NULL;
720
721	if (remaps->remaps) {
722	for (int i=0; i<remaps->n_remaps && !error; i++) {
723	error = MG_transfer_sequence(remaps->remaps[i], source_species, destination_species, remaps->alignment_names[i]);
724	}
725	}
726
727	return error;
728	}
729
730	// --------------------------------------------------------------------------------
731
732	#ifdef UNIT_TESTS
733
734	#include <test_unit.h>
735
736	// #define VERBOOSE_REMAP_TESTS
737
738	#if defined(VERBOOSE_REMAP_TESTS)
739	#define DUMP_REMAP(id, comment, from, to) fprintf(stderr, "%s %s '%s' -> '%s'\n", id, comment, from, to)
740	#define DUMP_REMAP_STR(str) fputs(str, stderr)
741	#else
742	#define DUMP_REMAP(id, comment, from, to)
743	#define DUMP_REMAP_STR(str)
744	#endif
745
746	#define TEST_REMAP(id, remapper, ASS_EQ, seqold, expected) \
747	char *calculated = remapper.remap(seqold); \
748	DUMP_REMAP(id, "want", seqold, expected); \
749	DUMP_REMAP(id, "calc", seqold, calculated); \
750	ASS_EQ(expected, calculated); \
751	DUMP_REMAP_STR("\n"); \
752	free(calculated);
753
754	#define TEST_REMAP1REF_INT(id, ASS_EQ, refold, refnew, seqold, expected) \
755	do { \
756	MG_remap remapper; \
757	DUMP_REMAP(id, "ref ", refold, refnew); \
758	remapper.add_reference(refold, refnew); \
759	TEST_REMAP(id, remapper, ASS_EQ, seqold, expected); \
760	} while(0)
761
762	#define TEST_REMAP2REFS_INT(id, ASS_EQ, refold1, refnew1, refold2, refnew2, seqold, expected) \
763	do { \
764	MG_remap remapper; \
765	DUMP_REMAP(id, "ref1", refold1, refnew1); \
766	DUMP_REMAP(id, "ref2", refold2, refnew2); \
767	remapper.add_reference(refold1, refnew1); \
768	remapper.add_reference(refold2, refnew2); \
769	TEST_REMAP(id, remapper, ASS_EQ, seqold, expected); \
770	} while(0)
771
772
773	#define TEST_REMAP1REF(id, ro, rn, seqold, expected) \
774	TEST_REMAP1REF_INT(id, TEST_EXPECT_EQUAL, ro, rn, seqold, expected)
775
776	#define TEST_REMAP1REF__BROKEN(id, ro, rn, seqold, expected) \
777	TEST_REMAP1REF_INT(id, TEST_EXPECT_EQUAL__BROKEN, ro, rn, seqold, expected)
778
779	#define TEST_REMAP2REFS(id, ro1, rn1, ro2, rn2, seqold, expected) \
780	TEST_REMAP2REFS_INT(id, TEST_EXPECT_EQUAL, ro1, rn1, ro2, rn2, seqold, expected)
781
782	#define TEST_REMAP2REFS__BROKEN(id, ro1, rn1, ro2, rn2, seqold, expected) \
783	TEST_REMAP2REFS_INT(id, TEST_EXPECT_EQUAL__BROKEN, ro1, rn1, ro2, rn2, seqold, expected)
784
785	#define TEST_REMAP1REF_FWDREV(id, ro, rn, so, sn) \
786	TEST_REMAP1REF(id "(fwd)", ro, rn, so, sn); \
787	TEST_REMAP1REF(id "(rev)", rn, ro, sn, so)
788
789	#define TEST_REMAP1REF_ALLDIR(id, ro, rn, so, sn) \
790	TEST_REMAP1REF_FWDREV(id "/ref=ref", ro, rn, so, sn); \
791	TEST_REMAP1REF_FWDREV(id "/ref=src", so, sn, ro, rn)
792
793	#define TEST_REMAP2REFS_FWDREV(id, ro1, rn1, ro2, rn2, so, sn) \
794	TEST_REMAP2REFS(id "(fwd)", ro1, rn1, ro2, rn2, so, sn); \
795	TEST_REMAP2REFS(id "(rev)", rn1, ro1, rn2, ro2, sn, so)
796
797	#define TEST_REMAP2REFS_ALLDIR(id, ro1, rn1, ro2, rn2, so, sn) \
798	TEST_REMAP2REFS_FWDREV(id "/ref=r1+r2 ", ro1, rn1, ro2, rn2, so, sn); \
799	TEST_REMAP2REFS_FWDREV(id "/ref=r1+src", ro1, rn1, so, sn, ro2, rn2); \
800	TEST_REMAP2REFS_FWDREV(id "/ref=r2+src", ro2, rn2, so, sn, ro1, rn1)
801
802
803	void TEST_remapping() {
804	// ----------------------------------------
805	// remap with 1 reference
806
807	TEST_REMAP1REF_ALLDIR("simple gap",
808	"ACGT", "AC--GT",
809	"TGCA", "TG--CA");
810
811	TEST_REMAP1REF_FWDREV("dotgap",
812	"ACGT", "AC..GT", // dots in reference do not get propagated
813	"TGCA", "TG--CA");
814
815	TEST_REMAP1REF("unused position (1)",
816	"A-C-G-T", "A--C--G--T",
817	"A---G-T", "A-----G--T");
818	TEST_REMAP1REF("unused position (2)",
819	"A-C-G-T", "A--C--G--T",
820	"Az-z-zT", "Az--z--z-T");
821
822	TEST_REMAP1REF("empty (1)",
823	"A-C-G", "A--C--G",
824	"", ".......");
825	TEST_REMAP1REF("empty (2)",
826	"...A-C-G...", "...A--C--G...",
827	"", "..........");
828
829	TEST_REMAP1REF("leadonly",
830	"...A-C-G...", "...A--C--G...",
831	"XX", ".XX.......");
832	TEST_REMAP1REF("trailonly",
833	"...A-C-G...", "...A--C--G...",
834	".........XX", "..........XX");
835	TEST_REMAP1REF__BROKEN("lead+trail",
836	"...A-C-G...", "...A--C--G...",
837	"XX.......XX", ".XX-------XX");
838
839	TEST_REMAP1REF("enforce leading dots (1)",
840	"--A-T", "------A--T",
841	"--T-A", "......T--A"); // leading gaps shall always be dots
842	TEST_REMAP1REF("enforce leading dots (2)",
843	"---A-T", "------A--T",
844	"-.-T-A", "......T--A"); // leading gaps shall always be dots
845	TEST_REMAP1REF("enforce leading dots (3)",
846	"---A-T", "------A--T",
847	"...T-A", "......T--A"); // leading gaps shall always be dots
848	TEST_REMAP1REF("enforce leading dots (4)",
849	"-..--A-T", "--.---A--T",
850	"-.-.-T-A", "......T--A"); // leading gaps shall always be dots
851	TEST_REMAP1REF("enforce leading dots (5)",
852	"---A-T-", "---A----T",
853	"-----A-", "........A"); // leading gaps shall always be dots
854	TEST_REMAP1REF("enforce leading dots (6)",
855	"---A-T-", "---A----T",
856	".....A-", "........A"); // leading gaps shall always be dots
857
858	TEST_REMAP1REF("no trailing gaps",
859	"A-T", "A--T---",
860	"T-A", "T--A");
861
862	TEST_REMAP1REF("should expand full-dot-gaps (1)",
863	"AC-GT", "AC--GT",
864	"TG.CA", "TG..CA"); // if gap was present and only contained dots -> new gap should only contain dots as well
865	TEST_REMAP1REF("should expand full-dot-gaps (2)",
866	"AC--GT", "AC----GT",
867	"TG..CA", "TG....CA"); // if gap was present and only contained dots -> new gap should only contain dots as well
868
869	TEST_REMAP1REF("keep 'missing bases' (1)",
870	"AC---GT", "AC---GT",
871	"TG-.-CA", "TG-.-CA"); // missing bases should not be deleted
872
873	TEST_REMAP1REF("keep 'missing bases' (2)",
874	"AC-D-GT", "AC-D-GT",
875	"TG-.-CA", "TG-.-CA"); // missing bases should not be deleted
876
877	// ----------------------------------------
878	// remap with 2 references
879	TEST_REMAP2REFS_ALLDIR("simple 3-way",
880	"A-CA-C-T", "AC-A---CT",
881	"A---GC-T", "A----G-CT",
882	"T-GAC--A", "TG-A-C--A");
883
884	TEST_REMAP2REFS("undefpos-nogap(2U)",
885	"---A", "------A",
886	"C---", "C------",
887	"GUUG", "GU---UG");
888	TEST_REMAP2REFS("undefpos-nogap(3U)",
889	"C----", "C------",
890	"----A", "------A",
891	"GUUUG", "GU--UUG");
892	TEST_REMAP2REFS("undefpos-nogap(4U)",
893	"-----A", "------A",
894	"C-----", "C------",
895	"GUUUUG", "GUU-UUG");
896	TEST_REMAP2REFS("undefpos-nogap(4U2)",
897	"-----A", "-------A",
898	"C-----", "C-------",
899	"GUUUUG", "GUU--UUG");
900
901	TEST_REMAP2REFS("undefpos1-gapleft",
902	"----A", "------A",
903	"C----", "C------",
904	"G-UUG", "G---UUG");
905	TEST_REMAP2REFS("undefpos1-gapright",
906	"----A", "------A",
907	"C----", "C------",
908	"GUU-G", "GU--U-G"); // behavior changed (prior: "GUU---G")
909
910	// test non-full sequences
911	TEST_REMAP2REFS("missing ali-pos (ref1-source)",
912	"C", "C--", // in-seq missing 1 ali pos
913	"-A", "--A",
914	"GG", "G-G");
915
916	TEST_REMAP2REFS("missing ali-pos (ref2-source)",
917	"-A", "--A",
918	"C", "C--", // in-seq missing 1 ali pos
919	"GG", "G-G");
920
921	TEST_REMAP2REFS("missing ali-pos (ref1-target)",
922	"C-", "C", // out-seq missing 2 ali pos
923	"-A", "--A",
924	"GG", "G-G");
925	TEST_REMAP2REFS("missing ali-pos (ref2-target)",
926	"-A", "--A",
927	"C-", "C", // out-seq missing 2 ali pos
928	"GG", "G-G");
929
930	TEST_REMAP2REFS("missing ali-pos (seq-source)",
931	"C--", "C--",
932	"-A-", "--A",
933	"GG", "G-G"); // in-seq missing 1 ali pos
934
935	// --------------------
936	// test (too) small gaps
937
938	TEST_REMAP1REF("gap gets too small (1)",
939	"A---T---A", "A--T--A",
940	"AGGGT---A", "AGGGT-A");
941
942	TEST_REMAP1REF("gap gets too small (2)",
943	"A---T---A", "A--T--A",
944	"AGGGT--GA", "AGGGTGA");
945
946	TEST_REMAP1REF("gap gets too small (3)",
947	"A---T---A", "A--T--A",
948	"AGGGT-G-A", "AGGGTGA");
949
950	TEST_REMAP1REF("gap gets too small (4)",
951	"A---T---A", "A--T--A",
952	"AGGGTG--A", "AGGGTGA");
953
954	TEST_REMAP1REF("gap tightens to fit (1)",
955	"A---T---A", "A--T--A",
956	"AGG-T---A", "AGGT--A");
957
958	TEST_REMAP1REF("gap tightens to fit (2)",
959	"A---T---A", "A--T--A",
960	"AGC-T---A", "AGCT--A");
961	TEST_REMAP1REF("gap tightens to fit (2)",
962	"A---T---A", "A--T--A",
963	"A-CGT---A", "ACGT--A");
964
965	TEST_REMAP1REF("gap with softmapping conflict (1)",
966	"A-------A", "A----A",
967	"A-CGT---A", "ACGT-A");
968	TEST_REMAP1REF("gap with softmapping conflict (2)",
969	"A-------A", "A----A",
970	"A--CGT--A", "ACGT-A");
971	TEST_REMAP1REF("gap with softmapping conflict (3)",
972	"A-------A", "A----A",
973	"A---CGT-A", "A-CGTA");
974	TEST_REMAP1REF("gap with softmapping conflict (4)",
975	"A-------A", "A----A",
976	"A----CGTA", "A-CGTA");
977	TEST_REMAP1REF("gap with softmapping conflict (5)",
978	"A-------A", "A----A",
979	"AC----GTA", "AC-GTA");
980
981	TEST_REMAP1REF("drop missing bases to avoid misalignment (1)",
982	"A---T", "A--T",
983	"AG.GT", "AGGT");
984	TEST_REMAP1REF("drop missing bases to avoid misalignment (2)",
985	"A----T", "A--T",
986	"AG..GT", "AGGT");
987	TEST_REMAP1REF("dont drop missing bases if fixed map",
988	"A-C-T", "A-CT",
989	"AG.GT", "AG.GT");
990	TEST_REMAP1REF("dont drop gaps if fixed map",
991	"A-C-T", "A-CT",
992	"AG-GT", "AG-GT");
993
994	// --------------------
995
996	TEST_REMAP1REF("append rest of seq (1)",
997	"AT" , "A--T",
998	"AGG---T...A...", "A--GGTA");
999	TEST_REMAP1REF("append rest of seq (2)",
1000	"AT." , "A--T--.",
1001	"AGG---T...A...", "A--GGTA");
1002	TEST_REMAP1REF("append rest of seq (3)",
1003	"AT--" , "A--T---",
1004	"AGG---T...A...", "A--GGTA");
1005	TEST_REMAP1REF("append rest of seq (4)",
1006	"ATC" , "A--T--C",
1007	"AGG---T...A...", "A--G--GTA");
1008	TEST_REMAP1REF("append rest of seq (4)",
1009	"AT-C" , "A--T--C",
1010	"AGG---T...A...", "A--GG--TA");
1011
1012	// --------------------
1013
1014	TEST_REMAP2REFS("impossible references",
1015	"AC-TA", "A---CT--A", // impossible references
1016	"A-GTA", "AG---T--A",
1017	"TGCAT", "T---GCA-T"); // solve by insertion (partial misalignment)
1018
1019	TEST_REMAP2REFS("impossible references",
1020	"AC-TA", "A--C-T--A", // impossible references
1021	"A-GTA", "AG---T--A",
1022	"TGCAT", "T--GCA--T"); // solve by insertion (partial misalignment)
1023	}
1024
1025	#endif

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

Download in other formats: