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 |
---|