Context Navigation

source: branches/profile/MERGE/MG_adapt_ali.cxx

Visit:

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

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

Download in other formats: