1 | // =============================================================== // |
---|
2 | // // |
---|
3 | // File : ED4_mini_classes.cxx // |
---|
4 | // Purpose : // |
---|
5 | // // |
---|
6 | // Institute of Microbiology (Technical University Munich) // |
---|
7 | // http://www.arb-home.de/ // |
---|
8 | // // |
---|
9 | // =============================================================== // |
---|
10 | |
---|
11 | |
---|
12 | #include "ed4_class.hxx" |
---|
13 | #include "ed4_edit_string.hxx" |
---|
14 | #include "ed4_awars.hxx" |
---|
15 | #include "ed4_tools.hxx" |
---|
16 | |
---|
17 | #include <aw_awar.hxx> |
---|
18 | #include <aw_root.hxx> |
---|
19 | #include <iupac.h> |
---|
20 | |
---|
21 | #include <cctype> |
---|
22 | |
---|
23 | // ------------------------ |
---|
24 | // ED4_bases_table |
---|
25 | |
---|
26 | #ifdef DEBUG |
---|
27 | // # define COUNT_BASES_TABLE_SIZE |
---|
28 | #endif |
---|
29 | |
---|
30 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
31 | static long bases_table_size = 0L; |
---|
32 | #endif |
---|
33 | |
---|
34 | void ED4_bases_table::init(int length) |
---|
35 | { |
---|
36 | if (length) { |
---|
37 | if (no_of_bases.shortTable) { |
---|
38 | delete no_of_bases.shortTable; |
---|
39 | no_of_bases.shortTable = 0; |
---|
40 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
41 | bases_table_size -= (no_of_entries+1)*table_entry_size; |
---|
42 | #endif |
---|
43 | } |
---|
44 | |
---|
45 | no_of_entries = length; |
---|
46 | int new_size = (no_of_entries+1)*table_entry_size; |
---|
47 | no_of_bases.shortTable = (unsigned char*)new char[new_size]; |
---|
48 | memset(no_of_bases.shortTable, 0, new_size); |
---|
49 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
50 | bases_table_size += new_size; |
---|
51 | printf("bases_table_size = %li\n", bases_table_size); |
---|
52 | #endif |
---|
53 | } |
---|
54 | } |
---|
55 | |
---|
56 | void ED4_bases_table::expand_table_entry_size() { // converts short table into long table |
---|
57 | e4_assert(table_entry_size==SHORT_TABLE_ELEM_SIZE); |
---|
58 | |
---|
59 | int *new_table = new int[no_of_entries+1]; |
---|
60 | int i; |
---|
61 | |
---|
62 | for (i=0; i<=no_of_entries; i++) { |
---|
63 | new_table[i] = no_of_bases.shortTable[i]; |
---|
64 | } |
---|
65 | delete [] no_of_bases.shortTable; |
---|
66 | |
---|
67 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
68 | bases_table_size += (no_of_entries+1)*(LONG_TABLE_ELEM_SIZE-SHORT_TABLE_ELEM_SIZE); |
---|
69 | printf("bases_table_size = %li\n", bases_table_size); |
---|
70 | #endif |
---|
71 | |
---|
72 | no_of_bases.longTable = new_table; |
---|
73 | table_entry_size = LONG_TABLE_ELEM_SIZE; |
---|
74 | } |
---|
75 | |
---|
76 | void ED4_bases_table::add(const ED4_bases_table& other, int start, int end) |
---|
77 | { |
---|
78 | e4_assert(no_of_entries==other.no_of_entries); |
---|
79 | e4_assert(start>=0); |
---|
80 | e4_assert(end<no_of_entries); |
---|
81 | e4_assert(start<=end); |
---|
82 | |
---|
83 | int i; |
---|
84 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
85 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
86 | for (i=start; i<=end; i++) { |
---|
87 | set_elem_short(i, get_elem_short(i)+other.get_elem_short(i)); |
---|
88 | } |
---|
89 | } |
---|
90 | else { |
---|
91 | for (i=start; i<=end; i++) { |
---|
92 | set_elem_short(i, get_elem_short(i)+other.get_elem_long(i)); |
---|
93 | } |
---|
94 | } |
---|
95 | } |
---|
96 | else { |
---|
97 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
98 | for (i=start; i<=end; i++) { |
---|
99 | set_elem_long(i, get_elem_long(i)+other.get_elem_short(i)); |
---|
100 | } |
---|
101 | } |
---|
102 | else { |
---|
103 | for (i=start; i<=end; i++) { |
---|
104 | set_elem_long(i, get_elem_long(i)+other.get_elem_long(i)); |
---|
105 | } |
---|
106 | } |
---|
107 | } |
---|
108 | } |
---|
109 | void ED4_bases_table::sub(const ED4_bases_table& other, int start, int end) |
---|
110 | { |
---|
111 | e4_assert(no_of_entries==other.no_of_entries); |
---|
112 | e4_assert(start>=0); |
---|
113 | e4_assert(end<no_of_entries); |
---|
114 | e4_assert(start<=end); |
---|
115 | |
---|
116 | int i; |
---|
117 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
118 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
119 | for (i=start; i<=end; i++) { |
---|
120 | set_elem_short(i, get_elem_short(i)-other.get_elem_short(i)); |
---|
121 | } |
---|
122 | } |
---|
123 | else { |
---|
124 | for (i=start; i<=end; i++) { |
---|
125 | set_elem_short(i, get_elem_short(i)-other.get_elem_long(i)); |
---|
126 | } |
---|
127 | } |
---|
128 | } |
---|
129 | else { |
---|
130 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
131 | for (i=start; i<=end; i++) { |
---|
132 | set_elem_long(i, get_elem_long(i)-other.get_elem_short(i)); |
---|
133 | } |
---|
134 | } |
---|
135 | else { |
---|
136 | for (i=start; i<=end; i++) { |
---|
137 | set_elem_long(i, get_elem_long(i)-other.get_elem_long(i)); |
---|
138 | } |
---|
139 | } |
---|
140 | } |
---|
141 | } |
---|
142 | void ED4_bases_table::sub_and_add(const ED4_bases_table& Sub, const ED4_bases_table& Add, PosRange range) |
---|
143 | { |
---|
144 | e4_assert(no_of_entries==Sub.no_of_entries); |
---|
145 | e4_assert(no_of_entries==Add.no_of_entries); |
---|
146 | |
---|
147 | int start = range.start(); |
---|
148 | int end = range.end(); |
---|
149 | |
---|
150 | e4_assert(start>=0); |
---|
151 | e4_assert(end<no_of_entries); |
---|
152 | e4_assert(start<=end); |
---|
153 | |
---|
154 | int i; |
---|
155 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
156 | if (Sub.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
157 | if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
158 | for (i=start; i<=end; i++) { |
---|
159 | set_elem_short(i, get_elem_short(i)-Sub.get_elem_short(i)+Add.get_elem_short(i)); |
---|
160 | } |
---|
161 | } |
---|
162 | else { |
---|
163 | for (i=start; i<=end; i++) { |
---|
164 | set_elem_short(i, get_elem_short(i)-Sub.get_elem_short(i)+Add.get_elem_long(i)); |
---|
165 | } |
---|
166 | } |
---|
167 | } |
---|
168 | else { |
---|
169 | if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
170 | for (i=start; i<=end; i++) { |
---|
171 | set_elem_short(i, get_elem_short(i)-Sub.get_elem_long(i)+Add.get_elem_short(i)); |
---|
172 | } |
---|
173 | } |
---|
174 | else { |
---|
175 | for (i=start; i<=end; i++) { |
---|
176 | set_elem_short(i, get_elem_short(i)-Sub.get_elem_long(i)+Add.get_elem_long(i)); |
---|
177 | } |
---|
178 | } |
---|
179 | } |
---|
180 | } |
---|
181 | else { |
---|
182 | if (Sub.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
183 | if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
184 | for (i=start; i<=end; i++) { |
---|
185 | set_elem_long(i, get_elem_long(i)-Sub.get_elem_short(i)+Add.get_elem_short(i)); |
---|
186 | } |
---|
187 | } |
---|
188 | else { |
---|
189 | for (i=start; i<=end; i++) { |
---|
190 | set_elem_long(i, get_elem_long(i)-Sub.get_elem_short(i)+Add.get_elem_long(i)); |
---|
191 | } |
---|
192 | } |
---|
193 | } |
---|
194 | else { |
---|
195 | if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
196 | for (i=start; i<=end; i++) { |
---|
197 | set_elem_long(i, get_elem_long(i)-Sub.get_elem_long(i)+Add.get_elem_short(i)); |
---|
198 | } |
---|
199 | } |
---|
200 | else { |
---|
201 | for (i=start; i<=end; i++) { |
---|
202 | set_elem_long(i, get_elem_long(i)-Sub.get_elem_long(i)+Add.get_elem_long(i)); |
---|
203 | } |
---|
204 | } |
---|
205 | } |
---|
206 | } |
---|
207 | } |
---|
208 | |
---|
209 | int ED4_bases_table::firstDifference(const ED4_bases_table& other, int start, int end, int *firstDifferentPos) const |
---|
210 | { |
---|
211 | int i; |
---|
212 | int result = 0; |
---|
213 | |
---|
214 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
215 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
216 | for (i=start; i<=end; i++) { |
---|
217 | if (get_elem_short(i)!=other.get_elem_short(i)) { |
---|
218 | *firstDifferentPos = i; |
---|
219 | result = 1; |
---|
220 | break; |
---|
221 | } |
---|
222 | } |
---|
223 | } |
---|
224 | else { |
---|
225 | for (i=start; i<=end; i++) { |
---|
226 | if (get_elem_short(i)!=other.get_elem_long(i)) { |
---|
227 | *firstDifferentPos = i; |
---|
228 | result = 1; |
---|
229 | break; |
---|
230 | } |
---|
231 | } |
---|
232 | } |
---|
233 | } |
---|
234 | else { |
---|
235 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
236 | for (i=start; i<=end; i++) { |
---|
237 | if (get_elem_long(i)!=other.get_elem_short(i)) { |
---|
238 | *firstDifferentPos = i; |
---|
239 | result = 1; |
---|
240 | break; |
---|
241 | } |
---|
242 | } |
---|
243 | } |
---|
244 | else { |
---|
245 | for (i=start; i<=end; i++) { |
---|
246 | if (get_elem_long(i)!=other.get_elem_long(i)) { |
---|
247 | *firstDifferentPos = i; |
---|
248 | result = 1; |
---|
249 | break; |
---|
250 | } |
---|
251 | } |
---|
252 | } |
---|
253 | } |
---|
254 | |
---|
255 | return result; |
---|
256 | } |
---|
257 | int ED4_bases_table::lastDifference(const ED4_bases_table& other, int start, int end, int *lastDifferentPos) const |
---|
258 | { |
---|
259 | int i; |
---|
260 | int result = 0; |
---|
261 | |
---|
262 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
263 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
264 | for (i=end; i>=start; i--) { |
---|
265 | if (get_elem_short(i)!=other.get_elem_short(i)) { |
---|
266 | *lastDifferentPos = i; |
---|
267 | result = 1; |
---|
268 | break; |
---|
269 | } |
---|
270 | } |
---|
271 | } |
---|
272 | else { |
---|
273 | for (i=end; i>=start; i--) { |
---|
274 | if (get_elem_short(i)!=other.get_elem_long(i)) { |
---|
275 | *lastDifferentPos = i; |
---|
276 | result = 1; |
---|
277 | break; |
---|
278 | } |
---|
279 | } |
---|
280 | } |
---|
281 | } |
---|
282 | else { |
---|
283 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
284 | for (i=end; i>=start; i--) { |
---|
285 | if (get_elem_long(i)!=other.get_elem_short(i)) { |
---|
286 | *lastDifferentPos = i; |
---|
287 | result = 1; |
---|
288 | break; |
---|
289 | } |
---|
290 | } |
---|
291 | } |
---|
292 | else { |
---|
293 | for (i=end; i>=start; i--) { |
---|
294 | if (get_elem_long(i)!=other.get_elem_long(i)) { |
---|
295 | *lastDifferentPos = i; |
---|
296 | result = 1; |
---|
297 | break; |
---|
298 | } |
---|
299 | } |
---|
300 | } |
---|
301 | } |
---|
302 | |
---|
303 | return result; |
---|
304 | } |
---|
305 | |
---|
306 | |
---|
307 | ED4_bases_table::ED4_bases_table(int maxseqlength) |
---|
308 | { |
---|
309 | memset((char*)this, 0, sizeof(*this)); |
---|
310 | table_entry_size = SHORT_TABLE_ELEM_SIZE; |
---|
311 | if (maxseqlength) init(maxseqlength); |
---|
312 | } |
---|
313 | |
---|
314 | ED4_bases_table::~ED4_bases_table() |
---|
315 | { |
---|
316 | delete [] no_of_bases.shortTable; |
---|
317 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
318 | bases_table_size -= (no_of_entries+1)*table_entry_size; |
---|
319 | printf("bases_table_size = %li\n", bases_table_size); |
---|
320 | #endif |
---|
321 | } |
---|
322 | |
---|
323 | void ED4_bases_table::change_table_length(int new_length, int default_entry) |
---|
324 | { |
---|
325 | if (new_length!=no_of_entries) { |
---|
326 | int min_length = new_length<no_of_entries ? new_length : no_of_entries; |
---|
327 | int growth = new_length-no_of_entries; |
---|
328 | |
---|
329 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
330 | unsigned char *new_table = new unsigned char[new_length+1]; |
---|
331 | |
---|
332 | e4_assert(default_entry>=0 && default_entry<256); |
---|
333 | |
---|
334 | memcpy(new_table, no_of_bases.shortTable, min_length*sizeof(*new_table)); |
---|
335 | new_table[new_length] = no_of_bases.shortTable[no_of_entries]; |
---|
336 | if (growth>0) { |
---|
337 | for (int e=no_of_entries; e<new_length; ++e) new_table[e] = default_entry; |
---|
338 | } |
---|
339 | |
---|
340 | delete [] no_of_bases.shortTable; |
---|
341 | no_of_bases.shortTable = new_table; |
---|
342 | } |
---|
343 | else { |
---|
344 | int *new_table = new int[new_length+1]; |
---|
345 | |
---|
346 | memcpy(new_table, no_of_bases.longTable, min_length*sizeof(*new_table)); |
---|
347 | new_table[new_length] = no_of_bases.longTable[no_of_entries]; |
---|
348 | if (growth>0) { |
---|
349 | for (int e=no_of_entries; e<new_length; ++e) new_table[e] = default_entry; |
---|
350 | } |
---|
351 | |
---|
352 | delete [] no_of_bases.longTable; |
---|
353 | no_of_bases.longTable = new_table; |
---|
354 | } |
---|
355 | |
---|
356 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
357 | bases_table_size += growth*table_entry_size; |
---|
358 | printf("bases_table_size = %li\n", bases_table_size); |
---|
359 | #endif |
---|
360 | no_of_entries = new_length; |
---|
361 | } |
---|
362 | } |
---|
363 | |
---|
364 | #if defined(TEST_CHAR_TABLE_INTEGRITY) || defined(ASSERTION_USED) |
---|
365 | |
---|
366 | int ED4_bases_table::empty() const |
---|
367 | { |
---|
368 | int i; |
---|
369 | |
---|
370 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
371 | for (i=0; i<no_of_entries; i++) { |
---|
372 | if (get_elem_short(i)) { |
---|
373 | return 0; |
---|
374 | } |
---|
375 | } |
---|
376 | } |
---|
377 | else { |
---|
378 | for (i=0; i<no_of_entries; i++) { |
---|
379 | if (get_elem_long(i)) { |
---|
380 | return 0; |
---|
381 | } |
---|
382 | } |
---|
383 | } |
---|
384 | |
---|
385 | return 1; |
---|
386 | } |
---|
387 | #endif // ASSERTION_USED |
---|
388 | |
---|
389 | // ------------------------ |
---|
390 | // Build consensus |
---|
391 | |
---|
392 | // we make static copies of the awars to avoid performance breakdown (BK_up_to_date is changed by callback ED4_consensus_definition_changed) |
---|
393 | |
---|
394 | struct ConsensusBuildParams { |
---|
395 | int countgaps; |
---|
396 | int gapbound; |
---|
397 | int group; |
---|
398 | int considbound; |
---|
399 | int upper; |
---|
400 | int lower; |
---|
401 | |
---|
402 | ConsensusBuildParams(AW_root *awr) { |
---|
403 | countgaps = awr->awar(ED4_AWAR_CONSENSUS_COUNTGAPS)->read_int(); |
---|
404 | gapbound = awr->awar(ED4_AWAR_CONSENSUS_GAPBOUND)->read_int(); |
---|
405 | group = awr->awar(ED4_AWAR_CONSENSUS_GROUP)->read_int(); |
---|
406 | considbound = awr->awar(ED4_AWAR_CONSENSUS_CONSIDBOUND)->read_int(); |
---|
407 | upper = awr->awar(ED4_AWAR_CONSENSUS_UPPER)->read_int(); |
---|
408 | lower = awr->awar(ED4_AWAR_CONSENSUS_LOWER)->read_int(); |
---|
409 | } |
---|
410 | #if defined(UNIT_TESTS) // UT_DIFF |
---|
411 | ConsensusBuildParams() { |
---|
412 | // uses awar defaults |
---|
413 | countgaps = 1; |
---|
414 | gapbound = 60; |
---|
415 | group = 1; |
---|
416 | considbound = 30; |
---|
417 | upper = 95; |
---|
418 | lower = 70; |
---|
419 | } |
---|
420 | #endif |
---|
421 | }; |
---|
422 | |
---|
423 | static ConsensusBuildParams *BK = NULL; // @@@ make member of ED4_char_table ? |
---|
424 | |
---|
425 | void ED4_consensus_definition_changed(AW_root*) { |
---|
426 | delete BK; BK = 0; // invalidate |
---|
427 | ED4_ROOT->request_refresh_for_consensus_terminals(); |
---|
428 | } |
---|
429 | |
---|
430 | static ARB_ERROR toggle_consensus_display(ED4_base *base, AW_CL show) { |
---|
431 | if (base->is_consensus_manager()) { |
---|
432 | ED4_manager *consensus_man = base->to_manager(); |
---|
433 | ED4_spacer_terminal *spacer = consensus_man->parent->get_defined_level(ED4_L_SPACER)->to_spacer_terminal(); |
---|
434 | |
---|
435 | if (show) { |
---|
436 | consensus_man->unhide_children(); |
---|
437 | spacer->extension.size[HEIGHT] = SPACERHEIGHT; |
---|
438 | } |
---|
439 | else { |
---|
440 | consensus_man->hide_children(); |
---|
441 | |
---|
442 | ED4_group_manager *group_man = consensus_man->get_parent(ED4_L_GROUP)->to_group_manager(); |
---|
443 | spacer->extension.size[HEIGHT] = (group_man->dynamic_prop&ED4_P_IS_FOLDED) ? SPACERNOCONSENSUSHEIGHT : SPACERHEIGHT; |
---|
444 | } |
---|
445 | } |
---|
446 | |
---|
447 | return NULL; |
---|
448 | } |
---|
449 | |
---|
450 | void ED4_consensus_display_changed(AW_root *root) { |
---|
451 | int show = root->awar(ED4_AWAR_CONSENSUS_SHOW)->read_int(); |
---|
452 | ED4_ROOT->root_group_man->route_down_hierarchy(toggle_consensus_display, show).expect_no_error(); |
---|
453 | } |
---|
454 | |
---|
455 | char *ED4_char_table::build_consensus_string(PosRange r) const { |
---|
456 | ExplicitRange range(r, size()); |
---|
457 | long entries = range.size(); |
---|
458 | char *new_buf = (char*)malloc(entries+1); |
---|
459 | |
---|
460 | build_consensus_string_to(new_buf, range); |
---|
461 | new_buf[entries] = 0; |
---|
462 | |
---|
463 | return new_buf; |
---|
464 | } |
---|
465 | |
---|
466 | void ED4_char_table::build_consensus_string_to(char *consensus_string, ExplicitRange range) const { |
---|
467 | // 'consensus_string' has to be a buffer of size 'range.size()+1' |
---|
468 | // Note : Always check that consensus behavior is identical to that used in CON_evaluatestatistic() |
---|
469 | |
---|
470 | if (!BK) BK = new ConsensusBuildParams(ED4_ROOT->aw_root); |
---|
471 | |
---|
472 | e4_assert(consensus_string); |
---|
473 | e4_assert(range.end()<size()); |
---|
474 | |
---|
475 | #define PERCENT(part, all) ((100*(part))/(all)) |
---|
476 | #define MAX_BASES_TABLES 41 // 25 |
---|
477 | |
---|
478 | e4_assert(used_bases_tables<=MAX_BASES_TABLES); // this is correct for DNA/RNA -> build_consensus_string() must be corrected for AMI/PRO |
---|
479 | |
---|
480 | const int left_idx = range.start(); |
---|
481 | const int right_idx = range.end(); |
---|
482 | |
---|
483 | if (sequences) { |
---|
484 | for (int i=left_idx; i<=right_idx; i++) { |
---|
485 | int o = i-left_idx; |
---|
486 | int bases = 0; // count of all bases together |
---|
487 | int base[MAX_BASES_TABLES]; // count of specific base |
---|
488 | int j; |
---|
489 | int max_base = -1; // maximum count of specific base |
---|
490 | int max_j = -1; // index of this base |
---|
491 | |
---|
492 | for (j=0; j<used_bases_tables; j++) { |
---|
493 | base[j] = linear_table(j)[i]; |
---|
494 | if (!ADPP_IS_ALIGN_CHARACTER(index_to_upperChar(j))) { |
---|
495 | bases += base[j]; |
---|
496 | if (base[j]>max_base) { // search for most used base |
---|
497 | max_base = base[j]; |
---|
498 | max_j = j; |
---|
499 | } |
---|
500 | } |
---|
501 | } |
---|
502 | |
---|
503 | int gaps = sequences-bases; // count of gaps |
---|
504 | |
---|
505 | // What to do with gaps? |
---|
506 | |
---|
507 | if (gaps == sequences) { |
---|
508 | consensus_string[o] = '='; |
---|
509 | } |
---|
510 | else if (BK->countgaps && PERCENT(gaps, sequences)>=BK->gapbound) { |
---|
511 | consensus_string[o] = '-'; |
---|
512 | } |
---|
513 | else { |
---|
514 | // Simplify using IUPAC : |
---|
515 | |
---|
516 | char kchar; // character for consensus |
---|
517 | int kcount; // count this character |
---|
518 | |
---|
519 | if (BK->group) { // group -> iupac |
---|
520 | if (ali_type == GB_AT_RNA || ali_type == GB_AT_DNA) { |
---|
521 | int no_of_bases = 0; // count of different bases used to create iupac |
---|
522 | char used_bases[MAX_BASES_TABLES+1]; // string containing those bases |
---|
523 | |
---|
524 | kcount = 0; |
---|
525 | for (j=0; j<used_bases_tables; j++) { |
---|
526 | int bchar = index_to_upperChar(j); |
---|
527 | |
---|
528 | if (!ADPP_IS_ALIGN_CHARACTER(bchar)) { |
---|
529 | if (PERCENT(base[j], sequences)>=BK->considbound) { |
---|
530 | used_bases[no_of_bases++] = index_to_upperChar(j); |
---|
531 | kcount += base[j]; |
---|
532 | } |
---|
533 | } |
---|
534 | } |
---|
535 | used_bases[no_of_bases] = 0; |
---|
536 | kchar = ED4_encode_iupac(used_bases, ali_type); |
---|
537 | } |
---|
538 | else { |
---|
539 | e4_assert(ali_type == GB_AT_AA); |
---|
540 | |
---|
541 | const int amino_groups = iupac::AA_GROUP_COUNT; |
---|
542 | int group_count[amino_groups]; |
---|
543 | |
---|
544 | for (j=0; j<amino_groups; j++) { |
---|
545 | group_count[j] = 0; |
---|
546 | } |
---|
547 | for (j=0; j<used_bases_tables; j++) { |
---|
548 | unsigned char bchar = index_to_upperChar(j); |
---|
549 | |
---|
550 | if (!ADPP_IS_ALIGN_CHARACTER(bchar)) { |
---|
551 | if (PERCENT(base[j], sequences)>=BK->considbound) { |
---|
552 | group_count[iupac::get_amino_group_for(bchar)] += base[j]; |
---|
553 | } |
---|
554 | } |
---|
555 | } |
---|
556 | |
---|
557 | kcount = 0; |
---|
558 | int bestGroup = 0; |
---|
559 | for (j=0; j<amino_groups; j++) { |
---|
560 | if (group_count[j]>kcount) { |
---|
561 | bestGroup = j; |
---|
562 | kcount = group_count[j]; |
---|
563 | } |
---|
564 | } |
---|
565 | |
---|
566 | kchar = iupac::get_amino_consensus_char(iupac::Amino_Group(bestGroup)); |
---|
567 | } |
---|
568 | } |
---|
569 | else { |
---|
570 | e4_assert(max_base); // expect at least one base to occur |
---|
571 | e4_assert(max_j >= 0); |
---|
572 | |
---|
573 | kchar = index_to_upperChar(max_j); |
---|
574 | kcount = max_base; |
---|
575 | e4_assert(kchar); |
---|
576 | } |
---|
577 | |
---|
578 | // show as upper or lower case ? |
---|
579 | |
---|
580 | int percent = PERCENT(kcount, sequences); |
---|
581 | |
---|
582 | if (percent>=BK->upper) { |
---|
583 | consensus_string[o] = kchar; |
---|
584 | } |
---|
585 | else if (percent>=BK->lower) { |
---|
586 | consensus_string[o] = tolower(kchar); |
---|
587 | } |
---|
588 | else { |
---|
589 | consensus_string[o] = '.'; |
---|
590 | } |
---|
591 | } |
---|
592 | e4_assert(consensus_string[o]); |
---|
593 | } |
---|
594 | } |
---|
595 | else { |
---|
596 | memset(consensus_string, '?', right_idx-left_idx+1); |
---|
597 | } |
---|
598 | #undef PERCENT |
---|
599 | } |
---|
600 | |
---|
601 | // ----------------------- |
---|
602 | // ED4_char_table |
---|
603 | |
---|
604 | bool ED4_char_table::initialized = false; |
---|
605 | unsigned char ED4_char_table::char_to_index_tab[MAXCHARTABLE]; |
---|
606 | unsigned char *ED4_char_table::upper_index_chars = 0; |
---|
607 | unsigned char *ED4_char_table::lower_index_chars = 0; |
---|
608 | int ED4_char_table::used_bases_tables = 0; |
---|
609 | GB_alignment_type ED4_char_table::ali_type = GB_AT_RNA; |
---|
610 | |
---|
611 | inline void ED4_char_table::set_char_to_index(unsigned char c, int index) |
---|
612 | { |
---|
613 | e4_assert(index>=0 && index<used_bases_tables); |
---|
614 | char_to_index_tab[upper_index_chars[index] = toupper(c)] = index; |
---|
615 | char_to_index_tab[lower_index_chars[index] = tolower(c)] = index; |
---|
616 | } |
---|
617 | unsigned char ED4_char_table::index_to_upperChar(int index) const |
---|
618 | { |
---|
619 | return upper_index_chars[index]; |
---|
620 | } |
---|
621 | unsigned char ED4_char_table::index_to_lowerChar(int index) const |
---|
622 | { |
---|
623 | return lower_index_chars[index]; |
---|
624 | } |
---|
625 | |
---|
626 | void ED4_char_table::expand_tables() { |
---|
627 | int i; |
---|
628 | for (i=0; i<used_bases_tables; i++) { |
---|
629 | linear_table(i).expand_table_entry_size(); |
---|
630 | } |
---|
631 | } |
---|
632 | |
---|
633 | void ED4_char_table::initial_setup(const char *gap_chars, GB_alignment_type ali_type_) { |
---|
634 | const char *groups = 0; |
---|
635 | used_bases_tables = strlen(gap_chars); |
---|
636 | ali_type = ali_type_; |
---|
637 | |
---|
638 | switch (ali_type) { |
---|
639 | case GB_AT_RNA: |
---|
640 | groups = "A,C,G,TU,MRWSYKVHDBN,"; |
---|
641 | break; |
---|
642 | case GB_AT_DNA: |
---|
643 | groups = "A,C,G,UT,MRWSYKVHDBN,"; |
---|
644 | break; |
---|
645 | case GB_AT_AA: |
---|
646 | groups = "P,A,G,S,T,Q,N,E,D,B,Z,H,K,R,L,I,V,M,F,Y,W,C,X"; // @@@ DRY (create 'groups' from AA_GROUP_...) |
---|
647 | break; |
---|
648 | default: |
---|
649 | e4_assert(0); |
---|
650 | break; |
---|
651 | } |
---|
652 | |
---|
653 | e4_assert(groups); |
---|
654 | |
---|
655 | int i; |
---|
656 | |
---|
657 | for (i=0; groups[i]; i++) { |
---|
658 | if (groups[i]==',') { |
---|
659 | used_bases_tables++; |
---|
660 | } |
---|
661 | } |
---|
662 | |
---|
663 | lower_index_chars = new unsigned char[used_bases_tables]; |
---|
664 | upper_index_chars = new unsigned char[used_bases_tables]; |
---|
665 | |
---|
666 | int idx = 0; |
---|
667 | unsigned char init_val = used_bases_tables-1; // map all unknown stuff to last table |
---|
668 | memset(char_to_index_tab, init_val, MAXCHARTABLE*sizeof(*char_to_index_tab)); |
---|
669 | |
---|
670 | for (i=0; groups[i]; i++) { |
---|
671 | if (groups[i]==',') { |
---|
672 | idx++; |
---|
673 | } |
---|
674 | else { |
---|
675 | e4_assert(isupper(groups[i])); |
---|
676 | set_char_to_index(groups[i], idx); |
---|
677 | } |
---|
678 | } |
---|
679 | |
---|
680 | const char *align_string_ptr = gap_chars; |
---|
681 | while (1) { |
---|
682 | char c = *align_string_ptr++; |
---|
683 | if (!c) break; |
---|
684 | set_char_to_index(c, idx++); |
---|
685 | } |
---|
686 | |
---|
687 | e4_assert(idx==used_bases_tables); |
---|
688 | initialized = true; |
---|
689 | } |
---|
690 | |
---|
691 | ED4_char_table::ED4_char_table(int maxseqlength) |
---|
692 | { |
---|
693 | memset((char*)this, 0, sizeof(*this)); |
---|
694 | |
---|
695 | if (!initialized) { |
---|
696 | char *align_string = ED4_ROOT->aw_root->awar_string(ED4_AWAR_GAP_CHARS)->read_string(); |
---|
697 | initial_setup(align_string, ED4_ROOT->alignment_type); |
---|
698 | free(align_string); |
---|
699 | } |
---|
700 | |
---|
701 | bases_table = new ED4_bases_table_ptr[used_bases_tables]; |
---|
702 | |
---|
703 | int i; |
---|
704 | for (i=0; i<used_bases_tables; i++) { |
---|
705 | bases_table[i] = new ED4_bases_table(maxseqlength); |
---|
706 | } |
---|
707 | |
---|
708 | sequences = 0; |
---|
709 | } |
---|
710 | |
---|
711 | void ED4_char_table::init(int maxseqlength) |
---|
712 | { |
---|
713 | int i; |
---|
714 | for (i=0; i<used_bases_tables; i++) { |
---|
715 | linear_table(i).init(maxseqlength); |
---|
716 | } |
---|
717 | |
---|
718 | sequences = 0; |
---|
719 | } |
---|
720 | |
---|
721 | void ED4_char_table::bases_and_gaps_at(int column, int *bases, int *gaps) const |
---|
722 | { |
---|
723 | int i, |
---|
724 | b = 0, |
---|
725 | g = 0; |
---|
726 | |
---|
727 | for (i=0; i<used_bases_tables; i++) { |
---|
728 | char c = upper_index_chars[i]; |
---|
729 | |
---|
730 | if (ADPP_IS_ALIGN_CHARACTER(c)) { |
---|
731 | g += table(c)[column]; |
---|
732 | } |
---|
733 | else { |
---|
734 | b += table(c)[column]; |
---|
735 | } |
---|
736 | } |
---|
737 | |
---|
738 | if (bases) *bases = b; |
---|
739 | if (gaps) *gaps = g; |
---|
740 | } |
---|
741 | |
---|
742 | ED4_char_table::~ED4_char_table() |
---|
743 | { |
---|
744 | int i; |
---|
745 | for (i=0; i<used_bases_tables; i++) { |
---|
746 | delete bases_table[i]; |
---|
747 | } |
---|
748 | |
---|
749 | delete [] bases_table; |
---|
750 | } |
---|
751 | |
---|
752 | const PosRange *ED4_char_table::changed_range(const ED4_char_table& other) const |
---|
753 | { |
---|
754 | int i; |
---|
755 | int Size = size(); |
---|
756 | int start = Size-1; |
---|
757 | int end = 0; |
---|
758 | |
---|
759 | e4_assert(Size==other.size()); |
---|
760 | for (i=0; i<used_bases_tables; i++) { |
---|
761 | if (linear_table(i).firstDifference(other.linear_table(i), 0, start, &start)) { |
---|
762 | if (end<start) { |
---|
763 | end = start; |
---|
764 | } |
---|
765 | linear_table(i).lastDifference(other.linear_table(i), end+1, Size-1, &end); |
---|
766 | |
---|
767 | for (; i<used_bases_tables; i++) { |
---|
768 | linear_table(i).firstDifference(other.linear_table(i), 0, start-1, &start); |
---|
769 | if (end<start) { |
---|
770 | end = start; |
---|
771 | } |
---|
772 | linear_table(i).lastDifference(other.linear_table(i), end+1, Size-1, &end); |
---|
773 | } |
---|
774 | |
---|
775 | e4_assert(start<=end); |
---|
776 | |
---|
777 | static PosRange range; |
---|
778 | range = PosRange(start, end); |
---|
779 | |
---|
780 | return ⦥ |
---|
781 | } |
---|
782 | } |
---|
783 | return NULL; |
---|
784 | } |
---|
785 | void ED4_char_table::add(const ED4_char_table& other) |
---|
786 | { |
---|
787 | if (other.ignore) return; |
---|
788 | if (other.size()==0) return; |
---|
789 | prepare_to_add_elements(other.sequences); |
---|
790 | add(other, 0, other.size()-1); |
---|
791 | } |
---|
792 | void ED4_char_table::add(const ED4_char_table& other, int start, int end) |
---|
793 | { |
---|
794 | if (other.ignore) return; |
---|
795 | |
---|
796 | test(); |
---|
797 | other.test(); |
---|
798 | |
---|
799 | e4_assert(start<=end); |
---|
800 | |
---|
801 | int i; |
---|
802 | for (i=0; i<used_bases_tables; i++) { |
---|
803 | linear_table(i).add(other.linear_table(i), start, end); |
---|
804 | } |
---|
805 | |
---|
806 | sequences += other.sequences; |
---|
807 | |
---|
808 | test(); |
---|
809 | } |
---|
810 | |
---|
811 | void ED4_char_table::sub(const ED4_char_table& other) |
---|
812 | { |
---|
813 | if (other.ignore) return; |
---|
814 | sub(other, 0, other.size()-1); |
---|
815 | } |
---|
816 | |
---|
817 | void ED4_char_table::sub(const ED4_char_table& other, int start, int end) |
---|
818 | { |
---|
819 | if (other.ignore) return; |
---|
820 | |
---|
821 | test(); |
---|
822 | other.test(); |
---|
823 | e4_assert(start<=end); |
---|
824 | |
---|
825 | int i; |
---|
826 | for (i=0; i<used_bases_tables; i++) { |
---|
827 | linear_table(i).sub(other.linear_table(i), start, end); |
---|
828 | } |
---|
829 | |
---|
830 | sequences -= other.sequences; |
---|
831 | |
---|
832 | test(); |
---|
833 | } |
---|
834 | |
---|
835 | void ED4_char_table::sub_and_add(const ED4_char_table& Sub, const ED4_char_table& Add) { |
---|
836 | test(); |
---|
837 | |
---|
838 | if (Sub.ignore) { |
---|
839 | e4_assert(Add.ignore); |
---|
840 | return; |
---|
841 | } |
---|
842 | |
---|
843 | Sub.test(); |
---|
844 | Add.test(); |
---|
845 | |
---|
846 | const PosRange *range = Sub.changed_range(Add); |
---|
847 | if (range) { |
---|
848 | prepare_to_add_elements(Add.added_sequences()-Sub.added_sequences()); |
---|
849 | sub_and_add(Sub, Add, *range); |
---|
850 | } |
---|
851 | |
---|
852 | test(); |
---|
853 | } |
---|
854 | |
---|
855 | void ED4_char_table::sub_and_add(const ED4_char_table& Sub, const ED4_char_table& Add, PosRange range) { |
---|
856 | test(); |
---|
857 | Sub.test(); |
---|
858 | Add.test(); |
---|
859 | |
---|
860 | e4_assert(!Sub.ignore && !Add.ignore); |
---|
861 | e4_assert(range.is_part()); |
---|
862 | |
---|
863 | int i; |
---|
864 | for (i=0; i<used_bases_tables; i++) { |
---|
865 | linear_table(i).sub_and_add(Sub.linear_table(i), Add.linear_table(i), range); |
---|
866 | } |
---|
867 | |
---|
868 | test(); |
---|
869 | } |
---|
870 | |
---|
871 | const PosRange *ED4_char_table::changed_range(const char *string1, const char *string2, int min_len) |
---|
872 | { |
---|
873 | const unsigned long *range1 = (const unsigned long *)string1; |
---|
874 | const unsigned long *range2 = (const unsigned long *)string2; |
---|
875 | const int step = sizeof(*range1); |
---|
876 | |
---|
877 | int l = min_len/step; |
---|
878 | int r = min_len%step; |
---|
879 | int i; |
---|
880 | int j; |
---|
881 | int k; |
---|
882 | |
---|
883 | int start = -1, end = -1; |
---|
884 | |
---|
885 | for (i=0; i<l; i++) { // search wordwise (it's faster) |
---|
886 | if (range1[i] != range2[i]) { |
---|
887 | k = i*step; |
---|
888 | for (j=0; j<step; j++) { |
---|
889 | if (string1[k+j] != string2[k+j]) { |
---|
890 | start = end = k+j; |
---|
891 | break; |
---|
892 | } |
---|
893 | } |
---|
894 | |
---|
895 | break; |
---|
896 | } |
---|
897 | } |
---|
898 | |
---|
899 | k = l*step; |
---|
900 | if (i==l) { // no difference found in word -> look at rest |
---|
901 | for (j=0; j<r; j++) { |
---|
902 | if (string1[k+j] != string2[k+j]) { |
---|
903 | start = end = k+j; |
---|
904 | break; |
---|
905 | } |
---|
906 | } |
---|
907 | |
---|
908 | if (j==r) { // the strings are equal |
---|
909 | return 0; |
---|
910 | } |
---|
911 | } |
---|
912 | |
---|
913 | e4_assert(start != -1 && end != -1); |
---|
914 | |
---|
915 | for (j=r-1; j>=0; j--) { // search rest backwards |
---|
916 | if (string1[k+j] != string2[k+j]) { |
---|
917 | end = k+j; |
---|
918 | break; |
---|
919 | } |
---|
920 | } |
---|
921 | |
---|
922 | if (j==-1) { // not found in rest -> search words backwards |
---|
923 | int m = start/step; |
---|
924 | for (i=l-1; i>=m; i--) { |
---|
925 | if (range1[i] != range2[i]) { |
---|
926 | k = i*step; |
---|
927 | for (j=step-1; j>=0; j--) { |
---|
928 | if (string1[k+j] != string2[k+j]) { |
---|
929 | end = k+j; |
---|
930 | break; |
---|
931 | } |
---|
932 | } |
---|
933 | break; |
---|
934 | } |
---|
935 | } |
---|
936 | } |
---|
937 | |
---|
938 | e4_assert(start<=end); |
---|
939 | |
---|
940 | static PosRange range; |
---|
941 | range = PosRange(start, end); |
---|
942 | |
---|
943 | return ⦥ |
---|
944 | } |
---|
945 | |
---|
946 | void ED4_char_table::add(const char *scan_string, int len) |
---|
947 | { |
---|
948 | test(); |
---|
949 | |
---|
950 | prepare_to_add_elements(1); |
---|
951 | |
---|
952 | int i; |
---|
953 | int sz = size(); |
---|
954 | |
---|
955 | if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) { |
---|
956 | for (i=0; i<len; i++) { |
---|
957 | unsigned char c = scan_string[i]; |
---|
958 | if (!c) break; |
---|
959 | table(c).inc_short(i); |
---|
960 | } |
---|
961 | ED4_bases_table& t = table('.'); |
---|
962 | for (; i<sz; i++) { |
---|
963 | t.inc_short(i); |
---|
964 | } |
---|
965 | } |
---|
966 | else { |
---|
967 | for (i=0; i<len; i++) { |
---|
968 | unsigned char c = scan_string[i]; |
---|
969 | if (!c) break; |
---|
970 | table(c).inc_long(i); |
---|
971 | } |
---|
972 | ED4_bases_table& t = table('.'); |
---|
973 | for (; i<sz; i++) { |
---|
974 | t.inc_long(i); |
---|
975 | } |
---|
976 | } |
---|
977 | |
---|
978 | sequences++; |
---|
979 | |
---|
980 | test(); |
---|
981 | } |
---|
982 | |
---|
983 | void ED4_char_table::sub(const char *scan_string, int len) |
---|
984 | { |
---|
985 | test(); |
---|
986 | |
---|
987 | int i; |
---|
988 | int sz = size(); |
---|
989 | |
---|
990 | if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) { |
---|
991 | for (i=0; i<len; i++) { |
---|
992 | unsigned char c = scan_string[i]; |
---|
993 | if (!c) break; |
---|
994 | table(c).dec_short(i); |
---|
995 | } |
---|
996 | ED4_bases_table& t = table('.'); |
---|
997 | for (; i<sz; i++) { |
---|
998 | t.dec_short(i); |
---|
999 | } |
---|
1000 | } |
---|
1001 | else { |
---|
1002 | for (i=0; i<len; i++) { |
---|
1003 | unsigned char c = scan_string[i]; |
---|
1004 | if (!c) break; |
---|
1005 | table(c).dec_long(i); |
---|
1006 | } |
---|
1007 | ED4_bases_table& t = table('.'); |
---|
1008 | for (; i<sz; i++) { |
---|
1009 | t.dec_long(i); |
---|
1010 | } |
---|
1011 | } |
---|
1012 | |
---|
1013 | sequences--; |
---|
1014 | |
---|
1015 | test(); |
---|
1016 | } |
---|
1017 | |
---|
1018 | void ED4_char_table::sub_and_add(const char *old_string, const char *new_string, PosRange range) { |
---|
1019 | test(); |
---|
1020 | int start = range.start(); |
---|
1021 | int end = range.end(); |
---|
1022 | e4_assert(start<=end); |
---|
1023 | |
---|
1024 | int i; |
---|
1025 | ED4_bases_table& t = table('.'); |
---|
1026 | |
---|
1027 | if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) { |
---|
1028 | for (i=start; i<=end; i++) { |
---|
1029 | unsigned char o = old_string[i]; |
---|
1030 | unsigned char n = new_string[i]; |
---|
1031 | |
---|
1032 | e4_assert(o); // @@@ if these never fail, some code below is superfluous |
---|
1033 | e4_assert(n); |
---|
1034 | |
---|
1035 | if (!o) { |
---|
1036 | for (; n && i<=end; i++) { |
---|
1037 | n = new_string[i]; |
---|
1038 | table(n).inc_short(i); |
---|
1039 | t.dec_short(i); |
---|
1040 | } |
---|
1041 | } |
---|
1042 | else if (!n) { |
---|
1043 | for (; o && i<=end; i++) { |
---|
1044 | o = old_string[i]; |
---|
1045 | table(o).dec_short(i); |
---|
1046 | t.inc_short(i); |
---|
1047 | } |
---|
1048 | |
---|
1049 | } |
---|
1050 | else if (o!=n) { |
---|
1051 | table(o).dec_short(i); |
---|
1052 | table(n).inc_short(i); |
---|
1053 | |
---|
1054 | } |
---|
1055 | } |
---|
1056 | } |
---|
1057 | else { |
---|
1058 | for (i=start; i<=end; i++) { |
---|
1059 | unsigned char o = old_string[i]; |
---|
1060 | unsigned char n = new_string[i]; |
---|
1061 | |
---|
1062 | e4_assert(o); // @@@ if these never fail, some code below is superfluous |
---|
1063 | e4_assert(n); |
---|
1064 | |
---|
1065 | if (!o) { |
---|
1066 | for (; n && i<=end; i++) { |
---|
1067 | n = new_string[i]; |
---|
1068 | table(n).inc_long(i); |
---|
1069 | t.dec_long(i); |
---|
1070 | } |
---|
1071 | } |
---|
1072 | else if (!n) { |
---|
1073 | for (; o && i<=end; i++) { |
---|
1074 | o = old_string[i]; |
---|
1075 | table(o).dec_long(i); |
---|
1076 | t.inc_long(i); |
---|
1077 | } |
---|
1078 | |
---|
1079 | } |
---|
1080 | else if (o!=n) { |
---|
1081 | table(o).dec_long(i); |
---|
1082 | table(n).inc_long(i); |
---|
1083 | } |
---|
1084 | } |
---|
1085 | } |
---|
1086 | |
---|
1087 | test(); |
---|
1088 | } |
---|
1089 | |
---|
1090 | void ED4_char_table::change_table_length(int new_length) { |
---|
1091 | for (int c=0; c<used_bases_tables; c++) { |
---|
1092 | linear_table(c).change_table_length(new_length, 0); |
---|
1093 | } |
---|
1094 | test(); |
---|
1095 | } |
---|
1096 | |
---|
1097 | // -------------- |
---|
1098 | // tests: |
---|
1099 | |
---|
1100 | #if defined(TEST_CHAR_TABLE_INTEGRITY) || defined(ASSERTION_USED) |
---|
1101 | |
---|
1102 | bool ED4_char_table::empty() const |
---|
1103 | { |
---|
1104 | int i; |
---|
1105 | for (i=0; i<used_bases_tables; i++) { |
---|
1106 | if (!linear_table(i).empty()) { |
---|
1107 | return false; |
---|
1108 | } |
---|
1109 | } |
---|
1110 | return true; |
---|
1111 | } |
---|
1112 | |
---|
1113 | bool ED4_char_table::ok() const |
---|
1114 | { |
---|
1115 | if (empty()) return true; |
---|
1116 | if (is_ignored()) return true; |
---|
1117 | |
---|
1118 | if (sequences<0) { |
---|
1119 | fprintf(stderr, "Negative # of sequences (%i) in ED4_char_table\n", sequences); |
---|
1120 | return false; |
---|
1121 | } |
---|
1122 | |
---|
1123 | int i; |
---|
1124 | for (i=0; i<MAXSEQUENCECHARACTERLENGTH; i++) { |
---|
1125 | int bases; |
---|
1126 | int gaps; |
---|
1127 | |
---|
1128 | bases_and_gaps_at(i, &bases, &gaps); |
---|
1129 | |
---|
1130 | if (!bases && !gaps) { // this occurs after insert column |
---|
1131 | int j; |
---|
1132 | for (j=i+1; j<MAXSEQUENCECHARACTERLENGTH; j++) { // test if we find any bases till end of table !!! |
---|
1133 | |
---|
1134 | bases_and_gaps_at(j, &bases, 0); |
---|
1135 | if (bases) { // bases found -> empty position was illegal |
---|
1136 | fprintf(stderr, "Zero in ED4_char_table at column %i\n", i); |
---|
1137 | return false; |
---|
1138 | } |
---|
1139 | } |
---|
1140 | break; |
---|
1141 | } |
---|
1142 | else { |
---|
1143 | if ((bases+gaps)!=sequences) { |
---|
1144 | fprintf(stderr, "bases+gaps should be equal to # of sequences at column %i\n", i); |
---|
1145 | return false; |
---|
1146 | } |
---|
1147 | } |
---|
1148 | } |
---|
1149 | |
---|
1150 | return true; |
---|
1151 | } |
---|
1152 | |
---|
1153 | #endif // TEST_CHAR_TABLE_INTEGRITY || ASSERTION_USED |
---|
1154 | |
---|
1155 | #if defined(TEST_CHAR_TABLE_INTEGRITY) |
---|
1156 | |
---|
1157 | void ED4_char_table::test() const { |
---|
1158 | |
---|
1159 | if (!ok()) { |
---|
1160 | GBK_terminate("ED4_char_table::test() failed"); |
---|
1161 | } |
---|
1162 | } |
---|
1163 | |
---|
1164 | #endif // TEST_CHAR_TABLE_INTEGRITY |
---|
1165 | |
---|
1166 | |
---|
1167 | // -------------------------------------------------------------------------------- |
---|
1168 | |
---|
1169 | #ifdef UNIT_TESTS |
---|
1170 | #ifndef TEST_UNIT_H |
---|
1171 | #include <test_unit.h> |
---|
1172 | #endif |
---|
1173 | |
---|
1174 | void TEST_char_table() { |
---|
1175 | const char alphabeth[] = "ACGTN-"; |
---|
1176 | const int alphabeth_size = strlen(alphabeth); |
---|
1177 | |
---|
1178 | const int seqlen = 30; |
---|
1179 | char seq[seqlen+1]; |
---|
1180 | |
---|
1181 | const char *gapChars = ".-"; // @@@ doesnt make any difference which gaps are used - why? |
---|
1182 | // const char *gapChars = "-"; |
---|
1183 | // const char *gapChars = "."; |
---|
1184 | // const char *gapChars = "A"; // makes me fail |
---|
1185 | ED4_char_table::initial_setup(gapChars, GB_AT_RNA); |
---|
1186 | ED4_init_is_align_character(gapChars); |
---|
1187 | |
---|
1188 | TEST_REJECT(BK); |
---|
1189 | BK = new ConsensusBuildParams; |
---|
1190 | |
---|
1191 | // BK->lower = 70; BK->upper = 95; BK->gapbound = 60; // defaults from awars |
---|
1192 | BK->lower = 40; BK->upper = 70; BK->gapbound = 40; |
---|
1193 | |
---|
1194 | srand(100); |
---|
1195 | for (int loop = 0; loop<5; ++loop) { |
---|
1196 | unsigned seed = rand(); |
---|
1197 | |
---|
1198 | srand(seed); |
---|
1199 | |
---|
1200 | ED4_char_table tab(seqlen); |
---|
1201 | |
---|
1202 | // build some seq |
---|
1203 | for (int c = 0; c<seqlen; ++c) seq[c] = alphabeth[rand()%alphabeth_size]; |
---|
1204 | seq[seqlen] = 0; |
---|
1205 | |
---|
1206 | TEST_EXPECT_EQUAL(strlen(seq), size_t(seqlen)); |
---|
1207 | |
---|
1208 | const int add_count = 300; |
---|
1209 | char *added_seqs[add_count]; |
---|
1210 | |
---|
1211 | // add seq multiple times |
---|
1212 | for (int a = 0; a<add_count; ++a) { |
---|
1213 | tab.add(seq, seqlen); |
---|
1214 | added_seqs[a] = strdup(seq); |
---|
1215 | |
---|
1216 | // modify 1 character in seq: |
---|
1217 | int sidx = rand()%seqlen; |
---|
1218 | int aidx = rand()%alphabeth_size; |
---|
1219 | seq[sidx] = alphabeth[aidx]; |
---|
1220 | } |
---|
1221 | |
---|
1222 | // build consensi (just check regression) |
---|
1223 | { |
---|
1224 | char *consensus = tab.build_consensus_string(); |
---|
1225 | switch (seed) { |
---|
1226 | case 677741240: TEST_EXPECT_EQUAL(consensus, ".-s-NW..aWu.NnWYa.R.mgcNK.c..."); break; |
---|
1227 | case 721151648: TEST_EXPECT_EQUAL(consensus, "a.nn..K..-gU.RW-SNcau.WNNacn.u"); break; |
---|
1228 | case 345295160: TEST_EXPECT_EQUAL(consensus, "..-g...MSn...guc.n.u.R.n.-Ng.k"); break; |
---|
1229 | case 346389111: TEST_EXPECT_EQUAL(consensus, ".unAn...gN.kc-cS.Raun...Sa-gY."); break; |
---|
1230 | case 367171911: TEST_EXPECT_EQUAL(consensus, "na.Na.nu.c-.-NU.aYgn-nng-.Wa.M"); break; |
---|
1231 | default: TEST_EXPECT_EQUAL(consensus, "undef"); |
---|
1232 | } |
---|
1233 | free(consensus); |
---|
1234 | } |
---|
1235 | |
---|
1236 | for (int a = 0; a<add_count; a++) { |
---|
1237 | tab.sub(added_seqs[a], seqlen); |
---|
1238 | freenull(added_seqs[a]); |
---|
1239 | } |
---|
1240 | { |
---|
1241 | char *consensus = tab.build_consensus_string(); |
---|
1242 | TEST_EXPECT_EQUAL(consensus, "??????????????????????????????"); // check tab is empty |
---|
1243 | free(consensus); |
---|
1244 | } |
---|
1245 | } |
---|
1246 | |
---|
1247 | delete BK; |
---|
1248 | BK = 0; |
---|
1249 | } |
---|
1250 | |
---|
1251 | #endif // UNIT_TESTS |
---|
1252 | |
---|
1253 | // -------------------------------------------------------------------------------- |
---|