source: branches/stable/NALIGNER/ali_solution.cxx

Last change on this file was 6385, checked in by westram, 10 years ago
  • removed trailing whitespace
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.5 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : ali_solution.cxx                                  //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "ali_solution.hxx"
12
13
14// ----------------
15//      ALI_MAP
16
17ALI_MAP::ALI_MAP(unsigned long first_seq, unsigned long last_seq,
18                 unsigned long first_ref, unsigned long last_ref)
19{
20    unsigned long l;
21
22    insert_counter = 0;
23    first_seq_base = first_seq;
24    last_seq_base = last_seq;
25    first_ref_base = first_ref;
26    last_ref_base = last_ref;
27
28    mapping   = (long **)          CALLOC((unsigned int) (last_seq_base - first_seq_base + 1), sizeof(long));
29    inserted  = (unsigned char **) CALLOC((unsigned int) ((last_seq_base - first_seq_base)/8) + 1, sizeof(unsigned char));
30    undefined = (unsigned char **) CALLOC((unsigned int) ((last_seq_base - first_seq_base)/8) + 1, sizeof(unsigned char));
31
32    if (mapping == 0 || inserted == 0 || undefined == 0) {
33        ali_fatal_error("Out of memory");
34    }
35
36    for (l = 0; l < (last_seq_base - first_seq_base)/8 + 1; l++)
37        (*undefined)[l] = 0xff;
38}
39
40ALI_MAP::ALI_MAP(ALI_MAP *map)
41{
42    unsigned long l;
43
44    first_seq_base = map->first_seq_base;
45    last_seq_base = map->last_seq_base;
46    first_ref_base = map->first_ref_base;
47    last_ref_base = map->last_ref_base;
48
49    mapping   = (long **)          CALLOC((unsigned int) (last_seq_base - first_seq_base + 1), sizeof(long));
50    inserted  = (unsigned char **) CALLOC((unsigned int) ((last_seq_base - first_seq_base) / 8) + 1, sizeof(unsigned char));
51    undefined = (unsigned char **) CALLOC((unsigned int) ((last_seq_base - first_seq_base) / 8) + 1, sizeof(unsigned char));
52
53    if (mapping == 0 || inserted == 0 || undefined == 0) {
54        ali_fatal_error("Out of memory");
55    }
56
57    for (l = 0; l < last_seq_base - first_seq_base + 1; l++) {
58        if (l < (last_seq_base - first_seq_base)/8 + 1) {
59            (*inserted)[l] = (*map->inserted)[l];
60            (*undefined)[l] = (*map->undefined)[l];
61        }
62        (*mapping)[l] = (*map->mapping)[l];
63    }
64}
65
66int ALI_MAP::is_konsistent()
67{
68    unsigned long i;
69
70    for (i = 1; i <= last_seq_base - first_seq_base; i++)
71        if ((*mapping)[i-1] >= (*mapping)[i])
72            return 0;
73
74    return 1;
75}
76
77ALI_SEQUENCE *ALI_MAP::sequence(ALI_NORM_SEQUENCE *ref_seq)
78{
79    int ins_counter;
80    int begin_flag = 0, undefined_flag;
81    unsigned long map_pos, seq_pos;
82    unsigned char *seq, *seq_buffer;
83
84    seq_buffer = (unsigned char *) CALLOC((unsigned int)
85                                          (last_ref_base - first_ref_base + insert_counter + 1),
86                                          sizeof(unsigned char));
87
88    if (seq_buffer == 0)
89        ali_fatal_error("Out of memory");
90
91    seq = seq_buffer;
92    seq_pos = 0;
93    undefined_flag = 0;
94    ins_counter = 0;
95    for (map_pos = 0; map_pos <= last_seq_base - first_seq_base; map_pos++) {
96        if (!undefined_flag)
97            begin_flag = ref_seq->is_begin(first_seq_base + map_pos);
98        if (!(((*undefined)[map_pos/8]>>(7-(map_pos%8))) & 0x01)) {
99            for (; seq_pos < (unsigned long)((*mapping)[map_pos] + ins_counter); seq_pos++) {
100                if (begin_flag)
101                    *seq++ = ALI_DOT_CODE;
102                else
103                    *seq++ = ALI_GAP_CODE;
104            }
105            *seq++ = ref_seq->base(first_seq_base + map_pos);
106            seq_pos++;
107            undefined_flag = 0;
108        }
109        else {
110            undefined_flag = 1;
111        }
112        if ((*inserted)[map_pos/8]>>(7-(map_pos%8)) & 0x01)
113            ins_counter++;
114    }
115
116    begin_flag = ref_seq->is_begin(first_seq_base + map_pos);
117    for (; seq_pos <= last_ref_base - first_ref_base + ins_counter; seq_pos++) {
118        if (begin_flag)
119            *seq++ = ALI_DOT_CODE;
120        else
121            *seq++ = ALI_GAP_CODE;
122    }
123
124    return new ALI_SEQUENCE(ref_seq->name(), seq_buffer,
125                            last_ref_base - first_ref_base + insert_counter + 1);
126}
127
128ALI_SEQUENCE *ALI_MAP::sequence_without_inserts(ALI_NORM_SEQUENCE *ref_seq)
129{
130    int begin_flag = 0, undefined_flag;
131    unsigned long map_pos, seq_pos;
132    unsigned char *seq, *seq_buffer;
133
134    seq_buffer = (unsigned char *) CALLOC((unsigned int)
135                                          (last_ref_base - first_ref_base + 1), sizeof(unsigned char));
136    if (seq_buffer == 0)
137        ali_fatal_error("Out of memory");
138
139    seq = seq_buffer;
140    seq_pos = 0;
141    undefined_flag = 0;
142    for (map_pos = 0; map_pos <= last_seq_base - first_seq_base; map_pos++) {
143        if (!undefined_flag)
144            begin_flag = ref_seq->is_begin(first_seq_base + map_pos);
145        if (!((*undefined)[map_pos/8]>>(7-(map_pos%8)) & 0x01) &&
146            !((*inserted)[map_pos/8]>>(7-(map_pos%8)) & 0x01)) {
147            for (; seq_pos < (unsigned long)((*mapping)[map_pos]); seq_pos++) {
148                if (begin_flag)
149                    *seq++ = ALI_DOT_CODE;
150                else
151                    *seq++ = ALI_GAP_CODE;
152            }
153            *seq++ = ref_seq->base(first_seq_base + map_pos);
154            seq_pos++;
155            undefined_flag = 0;
156        }
157        else
158            undefined_flag = 1;
159    }
160
161    begin_flag = ref_seq->is_begin(first_seq_base + map_pos);
162    for (; seq_pos <= last_ref_base - first_ref_base; seq_pos++) {
163        if (begin_flag)
164            *seq++ = ALI_DOT_CODE;
165        else
166            *seq++ = ALI_GAP_CODE;
167    }
168
169    return new ALI_SEQUENCE(ref_seq->name(), seq_buffer,
170                            last_ref_base - first_ref_base + 1);
171}
172
173ALI_MAP *ALI_MAP::inverse_without_inserts()
174{
175    unsigned long map_pos;
176    ALI_MAP *inv_map;
177
178    inv_map = new ALI_MAP(first_ref_base, last_ref_base,
179                          first_seq_base, last_seq_base);
180
181    for (map_pos = 0; map_pos <= last_seq_base - first_seq_base; map_pos++) {
182        if (!((*undefined)[map_pos/8]>>(7-(map_pos%8)) & 0x01) &&
183            !((*inserted)[map_pos/8]>>(7-(map_pos%8)) & 0x01)) {
184            inv_map->set(first_ref_base + (*mapping)[map_pos],
185                         map_pos);
186        }
187    }
188
189    return inv_map;
190}
191
192char *ALI_MAP::insert_marker()
193{
194    int ins_counter;
195    unsigned long map_pos, seq_pos;
196    char *seq, *seq_buffer;
197
198    seq_buffer = (char *) CALLOC((last_ref_base - first_ref_base + insert_counter + 2),
199                                  sizeof(char));
200
201    seq = seq_buffer;
202    seq_pos = 0;
203    ins_counter = 0;
204    for (map_pos = 0; map_pos <= last_seq_base - first_seq_base; map_pos++) {
205        if (!(((*undefined)[map_pos/8]>>(7-(map_pos%8))) & 0x01)) {
206            for (; seq_pos < (unsigned long)((*mapping)[map_pos] + ins_counter); seq_pos++) {
207                *seq++ = '.';
208            }
209            if ((*inserted)[map_pos/8]>>(7-(map_pos%8)) & 0x01)
210                *seq++ = 'X';
211            else
212                *seq++ = '.';
213            seq_pos++;
214        }
215        if ((*inserted)[map_pos/8]>>(7-(map_pos%8)) & 0x01)
216            ins_counter++;
217    }
218
219    for (; seq_pos <= last_ref_base - first_ref_base + ins_counter; seq_pos++) {
220        *seq++ = '.';
221    }
222
223    *seq = '\0';
224
225    return seq_buffer;
226}
227
228// -------------------------
229//      ALI_SUB_SOLUTION
230
231ALI_SUB_SOLUTION::ALI_SUB_SOLUTION(ALI_SUB_SOLUTION *solution)
232{
233    ALI_MAP *map;
234    ALI_TLIST<ALI_MAP *> *list;
235
236    profile = solution->profile;
237
238    if (!solution->map_list.is_empty()) {
239        list = &solution->map_list;
240        map = new ALI_MAP(list->first());
241        map_list.append_end(map);
242        while (list->is_next()) {
243            map = new ALI_MAP(list->next());
244            map_list.append_end(map);
245        }
246    }
247}
248
249ALI_SUB_SOLUTION::~ALI_SUB_SOLUTION()
250{
251    ALI_MAP *map;
252
253    if (!map_list.is_empty()) {
254        map = map_list.first();
255        delete map;
256        while (map_list.is_next()) {
257            map = map_list.next();
258            delete map;
259        }
260    }
261}
262
263int ALI_SUB_SOLUTION::free_area(
264                                unsigned long *start, unsigned long *end,
265                                unsigned long *start_ref, unsigned long *end_ref,
266                                unsigned long area_number)
267{
268    ALI_MAP *map;
269    unsigned long last_of_prev, last_of_prev_ref;
270    unsigned long area_number_akt;
271
272    if (map_list.is_empty()) {
273        if (area_number > 0)
274            return 0;
275        *start = 0;
276        *end = profile->sequence_length() - 1;
277        *start_ref = 0;
278        *end_ref = profile->length() - 1;
279        return 1;
280    }
281
282    area_number_akt = 0;
283    map = map_list.first();
284    if (map->first_base() > 0 &&
285        map->first_reference_base() > 0) {
286        if (area_number_akt == area_number) {
287            *start = 0;
288            *end = map->first_base() - 1;
289            *start_ref = 0;
290            *end_ref = map->first_reference_base() - 1;
291            return 1;
292        }
293        area_number_akt++;
294    }
295
296    last_of_prev = map->last_base();
297    last_of_prev_ref = map->last_reference_base();
298    while (map_list.is_next()) {
299        map = map_list.next();
300        if (map->first_base() > last_of_prev + 1) {
301            if (area_number_akt == area_number) {
302                *start = last_of_prev + 1;
303                *end = map->first_base() - 1;
304                *start_ref = last_of_prev_ref + 1;
305                *end_ref = map->first_reference_base() - 1;
306                return 1;
307            }
308            area_number_akt++;
309        }
310        last_of_prev = map->last_base();
311        last_of_prev_ref = map->last_reference_base();
312    }
313
314    if (map->last_base() < profile->sequence_length() - 1 &&
315        area_number_akt == area_number) {
316        *start = map->last_base() + 1;
317        *end = profile->sequence_length() - 1;
318        *start_ref = map->last_reference_base() + 1;
319        *end_ref = profile->length() - 1;
320        return 1;
321    }
322
323    return 0;
324}
325
326unsigned long ALI_SUB_SOLUTION::number_of_free_areas()
327{
328    ALI_MAP *map;
329    unsigned long last_of_prev;
330    unsigned long counter;
331
332    if (map_list.is_empty())
333        return 1;
334
335    counter = 0;
336    map = map_list.first();
337    if (map->first_base() > 0 && map->first_reference_base() > 0)
338        counter++;
339
340    last_of_prev = map->last_base();
341    while (map_list.is_next()) {
342        map = map_list.next();
343        if (map->first_base() > last_of_prev + 1)
344            counter++;
345        last_of_prev = map->last_base();
346    }
347
348    if (map->last_base() < profile->sequence_length() - 1)
349        counter++;
350
351    return counter;
352}
353
354
355int ALI_SUB_SOLUTION::is_konsistent(ALI_MAP *in_map)
356{
357    ALI_MAP *map;
358    unsigned long last_of_prev, last_of_prev_ref;
359
360    if (map_list.is_empty()) {
361        if (in_map->last_base() < profile->sequence_length() &&
362            in_map->last_reference_base() < profile->length())
363            return 1;
364        return 0;
365    }
366
367    map = map_list.first();
368    if (in_map->last_base() < map->first_base() &&
369        in_map->last_reference_base() < map->first_reference_base()) {
370        return 1;
371    }
372
373    last_of_prev = map->last_base();
374    last_of_prev_ref = map->last_reference_base();
375    while (map_list.is_next()) {
376        map = map_list.next();
377        if (last_of_prev < in_map->first_base() &&
378            in_map->last_base() < map->first_base() &&
379            last_of_prev_ref < in_map->first_reference_base() &&
380            in_map->last_reference_base() < map->first_reference_base()) {
381            return 1;
382        }
383        last_of_prev = map->last_base();
384        last_of_prev_ref = map->last_reference_base();
385    }
386
387    if (map->last_base() < in_map->first_base() &&
388        in_map->last_base() < profile->sequence_length() &&
389        map->last_reference_base() < in_map->first_base() &&
390        in_map->last_base() < profile->length()) {
391        return 1;
392    }
393
394    return 0;
395}
396
397int ALI_SUB_SOLUTION::insert(ALI_MAP *in_map)
398{
399    ALI_MAP *map;
400    unsigned long last_of_prev, last_of_prev_ref;
401
402    if (map_list.is_empty()) {
403        if (in_map->last_base() < profile->sequence_length() &&
404            in_map->last_reference_base() < profile->length()) {
405            map_list.append_end(in_map);
406            return 1;
407        }
408        return 0;
409    }
410
411    map = map_list.first();
412    if (in_map->last_base() < map->first_base() &&
413        in_map->last_reference_base() < map->first_reference_base()) {
414        map_list.insert_bevor(in_map);
415        return 1;
416    }
417
418    last_of_prev = map->last_base();
419    last_of_prev_ref = map->last_reference_base();
420    while (map_list.is_next()) {
421        map = map_list.next();
422        if (last_of_prev < in_map->first_base() &&
423            in_map->last_base() < map->first_base() &&
424            last_of_prev_ref < in_map->first_reference_base() &&
425            in_map->last_reference_base() < map->first_reference_base()) {
426            map_list.insert_bevor(in_map);
427            return 1;
428        }
429        last_of_prev = map->last_base();
430        last_of_prev_ref = map->last_reference_base();
431    }
432
433    if (map->last_base() < in_map->first_base() &&
434        in_map->last_base() < profile->sequence_length() &&
435        map->last_reference_base() < in_map->first_reference_base() &&
436        in_map->last_reference_base() < profile->length()) {
437        map_list.append_end(in_map);
438        return 1;
439    }
440
441    return 0;
442}
443
444int ALI_SUB_SOLUTION::delete_map(ALI_MAP *del_map)
445{
446    ALI_MAP *map;
447
448    if (map_list.is_empty())
449        return 0;
450
451    map = map_list.first();
452    if (map == del_map) {
453        map_list.delete_element();
454        return 1;
455    }
456
457    while (map_list.is_next()) {
458        map = map_list.next();
459        if (map == del_map) {
460            map_list.delete_element();
461            return 1;
462        }
463    }
464
465    return 0;
466}
467
468ALI_MAP *ALI_SUB_SOLUTION::make_one_map()
469{
470    ALI_MAP *map, *new_map;
471    unsigned long i;
472    unsigned long last_pos;
473    unsigned long first_base_of_first, first_reference_of_first;
474    unsigned long last_base_of_last, last_reference_of_last;
475
476    // check if maps are closed
477    if (map_list.is_empty())
478        return 0;
479
480    map = map_list.first();
481    first_base_of_first = map->first_base();
482    first_reference_of_first = map->first_reference_base();
483    last_base_of_last = map->last_base();
484    last_reference_of_last = map->last_reference_base();
485
486    while (map_list.is_next()) {
487        map = map_list.next();
488        if (last_base_of_last != map->first_base() - 1 ||
489            last_reference_of_last != map->first_reference_base() - 1)
490            ali_fatal_error("maps are not compact",
491                            "ALI_SUB_SOLUTION::make_one_map()");
492        last_base_of_last = map->last_base();
493        last_reference_of_last = map->last_reference_base();
494    }
495
496    new_map = new ALI_MAP(first_base_of_first, last_base_of_last,
497                          first_reference_of_first, last_reference_of_last);
498
499    map = map_list.first();
500    do {
501        last_pos = 0;
502        for (i = map->first_base(); i <= map->last_base(); i++) {
503            if (map->is_undefined(i))
504                ali_fatal_error("Unexpected value",
505                                "ALI_SUB_SOLUTION::make_one_map()");
506            if ((unsigned long)(map->position(i)) < last_pos)
507                ali_fatal_error("Inconsistent positions",
508                                "ALI_SUB_SOLUTION::make_one_map()");
509            last_pos = map->position(i);
510
511            if (map->is_inserted(i))
512                new_map->set(i, map->first_reference_base() +
513                             map->position(i) - first_reference_of_first, 1);
514            else
515                new_map->set(i, map->first_reference_base() +
516                             map->position(i) - first_reference_of_first, 0);
517        }
518
519        if (map_list.is_next())
520            map = map_list.next();
521        else
522            map = 0;
523    } while (map != 0);
524
525    return new_map;
526}
527
528void ALI_SUB_SOLUTION::print()
529{
530    ALI_MAP *map;
531
532    printf("ALI_SUB_SOLUTION:\n");
533    if (!map_list.is_empty()) {
534        map = map_list.first();
535        printf("(%ld to %ld) -> (%ld to %ld)\n",
536               map->first_base(), map->last_base(),
537               map->first_reference_base(), map->last_reference_base());
538        while (map_list.is_next()) {
539            map = map_list.next();
540            printf("(%ld to %ld) -> (%ld to %ld)\n",
541                   map->first_base(), map->last_base(),
542                   map->first_reference_base(), map->last_reference_base());
543        }
544    }
545}
546
Note: See TracBrowser for help on using the repository browser.