1 | // ================================================================= // |
---|
2 | // // |
---|
3 | // File : chartable.cxx // |
---|
4 | // Purpose : count characters of multiple sequences // |
---|
5 | // // |
---|
6 | // Institute of Microbiology (Technical University Munich) // |
---|
7 | // http://www.arb-home.de/ // |
---|
8 | // // |
---|
9 | // ================================================================= // |
---|
10 | |
---|
11 | #include "chartable.h" |
---|
12 | |
---|
13 | #include <arb_msg.h> |
---|
14 | #include <arb_mem.h> |
---|
15 | #include <iupac.h> |
---|
16 | #include <arb_defs.h> |
---|
17 | |
---|
18 | #include <cctype> |
---|
19 | |
---|
20 | #define CONSENSUS_AWAR_SOURCE CAS_INTERNAL // not awar-constructable here |
---|
21 | #include <consensus.h> |
---|
22 | |
---|
23 | using namespace chartable; |
---|
24 | |
---|
25 | // ------------------------ |
---|
26 | // SepBaseFreq |
---|
27 | |
---|
28 | #ifdef DEBUG |
---|
29 | // # define COUNT_BASES_TABLE_SIZE |
---|
30 | #endif |
---|
31 | |
---|
32 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
33 | static long bases_table_size = 0L; |
---|
34 | #endif |
---|
35 | |
---|
36 | void SepBaseFreq::init(int length) |
---|
37 | { |
---|
38 | if (length) { |
---|
39 | if (no_of_bases.shortTable) { |
---|
40 | delete no_of_bases.shortTable; |
---|
41 | no_of_bases.shortTable = 0; |
---|
42 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
43 | bases_table_size -= (no_of_entries+1)*table_entry_size; |
---|
44 | #endif |
---|
45 | } |
---|
46 | |
---|
47 | no_of_entries = length; |
---|
48 | int new_size = (no_of_entries+1)*table_entry_size; |
---|
49 | no_of_bases.shortTable = (unsigned char*)new char[new_size]; |
---|
50 | memset(no_of_bases.shortTable, 0, new_size); |
---|
51 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
52 | bases_table_size += new_size; |
---|
53 | printf("bases_table_size = %li\n", bases_table_size); |
---|
54 | #endif |
---|
55 | } |
---|
56 | } |
---|
57 | |
---|
58 | void SepBaseFreq::expand_table_entry_size() { // converts short table into long table |
---|
59 | ct_assert(table_entry_size==SHORT_TABLE_ELEM_SIZE); |
---|
60 | |
---|
61 | int *new_table = new int[no_of_entries+1]; |
---|
62 | int i; |
---|
63 | |
---|
64 | for (i=0; i<=no_of_entries; i++) { |
---|
65 | new_table[i] = no_of_bases.shortTable[i]; |
---|
66 | } |
---|
67 | delete [] no_of_bases.shortTable; |
---|
68 | |
---|
69 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
70 | bases_table_size += (no_of_entries+1)*(LONG_TABLE_ELEM_SIZE-SHORT_TABLE_ELEM_SIZE); |
---|
71 | printf("bases_table_size = %li\n", bases_table_size); |
---|
72 | #endif |
---|
73 | |
---|
74 | no_of_bases.longTable = new_table; |
---|
75 | table_entry_size = LONG_TABLE_ELEM_SIZE; |
---|
76 | } |
---|
77 | |
---|
78 | #define INVALID_TABLE_OPERATION() GBK_terminatef("SepBaseFreq: invalid operation at %i", __LINE__) |
---|
79 | |
---|
80 | void SepBaseFreq::add(const SepBaseFreq& other, int start, int end) |
---|
81 | { |
---|
82 | ct_assert(no_of_entries==other.no_of_entries); |
---|
83 | ct_assert(start>=0); |
---|
84 | ct_assert(end<no_of_entries); |
---|
85 | ct_assert(start<=end); |
---|
86 | |
---|
87 | int i; |
---|
88 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
89 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
90 | for (i=start; i<=end; i++) { |
---|
91 | set_elem_short(i, get_elem_short(i)+other.get_elem_short(i)); |
---|
92 | } |
---|
93 | } |
---|
94 | else { |
---|
95 | INVALID_TABLE_OPERATION(); // cannot add long to short |
---|
96 | } |
---|
97 | } |
---|
98 | else { |
---|
99 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
100 | for (i=start; i<=end; i++) { // LOOP_VECTORIZED |
---|
101 | set_elem_long(i, get_elem_long(i)+other.get_elem_short(i)); |
---|
102 | } |
---|
103 | } |
---|
104 | else { |
---|
105 | for (i=start; i<=end; i++) { // LOOP_VECTORIZED |
---|
106 | set_elem_long(i, get_elem_long(i)+other.get_elem_long(i)); |
---|
107 | } |
---|
108 | } |
---|
109 | } |
---|
110 | } |
---|
111 | void SepBaseFreq::sub(const SepBaseFreq& other, int start, int end) |
---|
112 | { |
---|
113 | ct_assert(no_of_entries==other.no_of_entries); |
---|
114 | ct_assert(start>=0); |
---|
115 | ct_assert(end<no_of_entries); |
---|
116 | ct_assert(start<=end); |
---|
117 | |
---|
118 | int i; |
---|
119 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
120 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
121 | for (i=start; i<=end; i++) { |
---|
122 | set_elem_short(i, get_elem_short(i)-other.get_elem_short(i)); |
---|
123 | } |
---|
124 | } |
---|
125 | else { |
---|
126 | INVALID_TABLE_OPERATION(); // cannot sub long from short |
---|
127 | } |
---|
128 | } |
---|
129 | else { |
---|
130 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
131 | for (i=start; i<=end; i++) { // LOOP_VECTORIZED |
---|
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++) { // LOOP_VECTORIZED |
---|
137 | set_elem_long(i, get_elem_long(i)-other.get_elem_long(i)); |
---|
138 | } |
---|
139 | } |
---|
140 | } |
---|
141 | } |
---|
142 | void SepBaseFreq::sub_and_add(const SepBaseFreq& Sub, const SepBaseFreq& Add, PosRange range) |
---|
143 | { |
---|
144 | ct_assert(no_of_entries==Sub.no_of_entries); |
---|
145 | ct_assert(no_of_entries==Add.no_of_entries); |
---|
146 | |
---|
147 | int start = range.start(); |
---|
148 | int end = range.end(); |
---|
149 | |
---|
150 | ct_assert(start>=0); |
---|
151 | ct_assert(end<no_of_entries); |
---|
152 | ct_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 | Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) |
---|
158 | { |
---|
159 | for (i=start; i<=end; i++) { |
---|
160 | set_elem_short(i, get_elem_short(i)-Sub.get_elem_short(i)+Add.get_elem_short(i)); |
---|
161 | } |
---|
162 | } |
---|
163 | else { |
---|
164 | INVALID_TABLE_OPERATION(); // cannot add or sub long to/from short |
---|
165 | } |
---|
166 | } |
---|
167 | else { |
---|
168 | if (Sub.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
169 | if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
170 | for (i=start; i<=end; i++) { // LOOP_VECTORIZED |
---|
171 | set_elem_long(i, get_elem_long(i)-Sub.get_elem_short(i)+Add.get_elem_short(i)); |
---|
172 | } |
---|
173 | } |
---|
174 | else { |
---|
175 | for (i=start; i<=end; i++) { // LOOP_VECTORIZED |
---|
176 | set_elem_long(i, get_elem_long(i)-Sub.get_elem_short(i)+Add.get_elem_long(i)); |
---|
177 | } |
---|
178 | } |
---|
179 | } |
---|
180 | else { |
---|
181 | if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
182 | for (i=start; i<=end; i++) { // LOOP_VECTORIZED |
---|
183 | set_elem_long(i, get_elem_long(i)-Sub.get_elem_long(i)+Add.get_elem_short(i)); |
---|
184 | } |
---|
185 | } |
---|
186 | else { |
---|
187 | for (i=start; i<=end; i++) { // LOOP_VECTORIZED |
---|
188 | set_elem_long(i, get_elem_long(i)-Sub.get_elem_long(i)+Add.get_elem_long(i)); |
---|
189 | } |
---|
190 | } |
---|
191 | } |
---|
192 | } |
---|
193 | } |
---|
194 | |
---|
195 | int SepBaseFreq::firstDifference(const SepBaseFreq& other, int start, int end, int *firstDifferentPos) const |
---|
196 | { |
---|
197 | int i; |
---|
198 | int result = 0; |
---|
199 | |
---|
200 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
201 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
202 | for (i=start; i<=end; i++) { |
---|
203 | if (get_elem_short(i)!=other.get_elem_short(i)) { |
---|
204 | *firstDifferentPos = i; |
---|
205 | result = 1; |
---|
206 | break; |
---|
207 | } |
---|
208 | } |
---|
209 | } |
---|
210 | else { |
---|
211 | for (i=start; i<=end; i++) { |
---|
212 | if (get_elem_short(i)!=other.get_elem_long(i)) { |
---|
213 | *firstDifferentPos = i; |
---|
214 | result = 1; |
---|
215 | break; |
---|
216 | } |
---|
217 | } |
---|
218 | } |
---|
219 | } |
---|
220 | else { |
---|
221 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
222 | for (i=start; i<=end; i++) { |
---|
223 | if (get_elem_long(i)!=other.get_elem_short(i)) { |
---|
224 | *firstDifferentPos = i; |
---|
225 | result = 1; |
---|
226 | break; |
---|
227 | } |
---|
228 | } |
---|
229 | } |
---|
230 | else { |
---|
231 | for (i=start; i<=end; i++) { |
---|
232 | if (get_elem_long(i)!=other.get_elem_long(i)) { |
---|
233 | *firstDifferentPos = i; |
---|
234 | result = 1; |
---|
235 | break; |
---|
236 | } |
---|
237 | } |
---|
238 | } |
---|
239 | } |
---|
240 | |
---|
241 | return result; |
---|
242 | } |
---|
243 | int SepBaseFreq::lastDifference(const SepBaseFreq& other, int start, int end, int *lastDifferentPos) const |
---|
244 | { |
---|
245 | int i; |
---|
246 | int result = 0; |
---|
247 | |
---|
248 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
249 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
250 | for (i=end; i>=start; i--) { |
---|
251 | if (get_elem_short(i)!=other.get_elem_short(i)) { |
---|
252 | *lastDifferentPos = i; |
---|
253 | result = 1; |
---|
254 | break; |
---|
255 | } |
---|
256 | } |
---|
257 | } |
---|
258 | else { |
---|
259 | for (i=end; i>=start; i--) { |
---|
260 | if (get_elem_short(i)!=other.get_elem_long(i)) { |
---|
261 | *lastDifferentPos = i; |
---|
262 | result = 1; |
---|
263 | break; |
---|
264 | } |
---|
265 | } |
---|
266 | } |
---|
267 | } |
---|
268 | else { |
---|
269 | if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
270 | for (i=end; i>=start; i--) { |
---|
271 | if (get_elem_long(i)!=other.get_elem_short(i)) { |
---|
272 | *lastDifferentPos = i; |
---|
273 | result = 1; |
---|
274 | break; |
---|
275 | } |
---|
276 | } |
---|
277 | } |
---|
278 | else { |
---|
279 | for (i=end; i>=start; i--) { |
---|
280 | if (get_elem_long(i)!=other.get_elem_long(i)) { |
---|
281 | *lastDifferentPos = i; |
---|
282 | result = 1; |
---|
283 | break; |
---|
284 | } |
---|
285 | } |
---|
286 | } |
---|
287 | } |
---|
288 | |
---|
289 | return result; |
---|
290 | } |
---|
291 | |
---|
292 | |
---|
293 | SepBaseFreq::SepBaseFreq(int maxseqlength) |
---|
294 | : table_entry_size(SHORT_TABLE_ELEM_SIZE), |
---|
295 | no_of_entries(0) |
---|
296 | { |
---|
297 | no_of_bases.shortTable = NULL; |
---|
298 | if (maxseqlength) init(maxseqlength); |
---|
299 | } |
---|
300 | |
---|
301 | SepBaseFreq::~SepBaseFreq() |
---|
302 | { |
---|
303 | delete [] no_of_bases.shortTable; |
---|
304 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
305 | bases_table_size -= (no_of_entries+1)*table_entry_size; |
---|
306 | printf("bases_table_size = %li\n", bases_table_size); |
---|
307 | #endif |
---|
308 | } |
---|
309 | |
---|
310 | void SepBaseFreq::change_table_length(int new_length, int default_entry) |
---|
311 | { |
---|
312 | if (new_length!=no_of_entries) { |
---|
313 | int min_length = new_length<no_of_entries ? new_length : no_of_entries; |
---|
314 | int growth = new_length-no_of_entries; |
---|
315 | |
---|
316 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
317 | unsigned char *new_table = new unsigned char[new_length+1]; |
---|
318 | |
---|
319 | ct_assert(default_entry>=0 && default_entry<256); |
---|
320 | |
---|
321 | memcpy(new_table, no_of_bases.shortTable, min_length*sizeof(*new_table)); |
---|
322 | new_table[new_length] = no_of_bases.shortTable[no_of_entries]; |
---|
323 | if (growth>0) { |
---|
324 | for (int e=no_of_entries; e<new_length; ++e) new_table[e] = default_entry; |
---|
325 | } |
---|
326 | |
---|
327 | delete [] no_of_bases.shortTable; |
---|
328 | no_of_bases.shortTable = new_table; |
---|
329 | } |
---|
330 | else { |
---|
331 | int *new_table = new int[new_length+1]; |
---|
332 | |
---|
333 | memcpy(new_table, no_of_bases.longTable, min_length*sizeof(*new_table)); |
---|
334 | new_table[new_length] = no_of_bases.longTable[no_of_entries]; |
---|
335 | if (growth>0) { |
---|
336 | for (int e=no_of_entries; e<new_length; ++e) { // LOOP_VECTORIZED |
---|
337 | new_table[e] = default_entry; |
---|
338 | } |
---|
339 | } |
---|
340 | |
---|
341 | delete [] no_of_bases.longTable; |
---|
342 | no_of_bases.longTable = new_table; |
---|
343 | } |
---|
344 | |
---|
345 | #ifdef COUNT_BASES_TABLE_SIZE |
---|
346 | bases_table_size += growth*table_entry_size; |
---|
347 | printf("bases_table_size = %li\n", bases_table_size); |
---|
348 | #endif |
---|
349 | no_of_entries = new_length; |
---|
350 | } |
---|
351 | } |
---|
352 | |
---|
353 | #if defined(TEST_CHAR_TABLE_INTEGRITY) || defined(ASSERTION_USED) |
---|
354 | |
---|
355 | int SepBaseFreq::empty() const |
---|
356 | { |
---|
357 | int i; |
---|
358 | |
---|
359 | if (table_entry_size==SHORT_TABLE_ELEM_SIZE) { |
---|
360 | for (i=0; i<no_of_entries; i++) { |
---|
361 | if (get_elem_short(i)) { |
---|
362 | return 0; |
---|
363 | } |
---|
364 | } |
---|
365 | } |
---|
366 | else { |
---|
367 | for (i=0; i<no_of_entries; i++) { |
---|
368 | if (get_elem_long(i)) { |
---|
369 | return 0; |
---|
370 | } |
---|
371 | } |
---|
372 | } |
---|
373 | |
---|
374 | return 1; |
---|
375 | } |
---|
376 | #endif // ASSERTION_USED |
---|
377 | |
---|
378 | char *BaseFrequencies::build_consensus_string(PosRange r, const ConsensusBuildParams& cbp) const { |
---|
379 | ExplicitRange range(r, size()); |
---|
380 | |
---|
381 | long entries = range.size(); |
---|
382 | char *new_buf = ARB_alloc<char>(entries+1); |
---|
383 | |
---|
384 | build_consensus_string_to(new_buf, range, cbp); |
---|
385 | new_buf[entries] = 0; |
---|
386 | |
---|
387 | return new_buf; |
---|
388 | } |
---|
389 | |
---|
390 | #if defined(DEBUG) |
---|
391 | // #define DEBUG_CONSENSUS |
---|
392 | #endif |
---|
393 | |
---|
394 | void BaseFrequencies::build_consensus_string_to(char *consensus_string, ExplicitRange range, const ConsensusBuildParams& BK) const { |
---|
395 | // 'consensus_string' has to be a buffer of size 'range.size()+1' |
---|
396 | // Note : Always check that consensus behavior is identical to that used in CON_evaluatestatistic() |
---|
397 | |
---|
398 | ct_assert(consensus_string); |
---|
399 | ct_assert(range.end()<size()); |
---|
400 | |
---|
401 | #define PERCENT(part, all) ((100*(part))/(all)) |
---|
402 | #define MAX_BASES_TABLES 41 // 25 |
---|
403 | |
---|
404 | ct_assert(used_bases_tables<=MAX_BASES_TABLES); // this is correct for DNA/RNA -> build_consensus_string() must be corrected for AMI/PRO |
---|
405 | |
---|
406 | const int left_idx = range.start(); |
---|
407 | const int right_idx = range.end(); |
---|
408 | |
---|
409 | #if defined(DEBUG_CONSENSUS) |
---|
410 | #define DUMPINT(var) do { if (dumpcol) fprintf(stderr, "%s=%i ", #var, var); } while(0) |
---|
411 | #else |
---|
412 | #define DUMPINT(var) |
---|
413 | #endif |
---|
414 | |
---|
415 | if (sequenceUnits) { |
---|
416 | for (int i=left_idx; i<=right_idx; i++) { |
---|
417 | #if defined(DEBUG_CONSENSUS) |
---|
418 | int column = i+1; |
---|
419 | bool dumpcol = (column == 21); |
---|
420 | DUMPINT(column); |
---|
421 | #endif |
---|
422 | |
---|
423 | const int o = i-left_idx; // output offset |
---|
424 | |
---|
425 | int bases = 0; // count of all bases together |
---|
426 | int base[MAX_BASES_TABLES]; // count of specific base |
---|
427 | int max_base = -1; // maximum count of specific base |
---|
428 | int max_base_idx = -1; // index of this base |
---|
429 | |
---|
430 | for (int j=0; j<used_bases_tables; j++) { |
---|
431 | base[j] = linear_table(j)[i]; |
---|
432 | if (!isGap(index_to_upperChar(j))) { |
---|
433 | bases += base[j]; |
---|
434 | if (base[j]>max_base) { // search for most used base |
---|
435 | max_base = base[j]; |
---|
436 | max_base_idx = j; |
---|
437 | } |
---|
438 | } |
---|
439 | } |
---|
440 | |
---|
441 | int gaps = sequenceUnits-bases; // count of gaps |
---|
442 | |
---|
443 | DUMPINT(sequenceUnits); |
---|
444 | DUMPINT(gaps); |
---|
445 | DUMPINT(bases); |
---|
446 | |
---|
447 | // What to do with gaps? |
---|
448 | |
---|
449 | if (gaps == sequenceUnits) { |
---|
450 | consensus_string[o] = '='; |
---|
451 | } |
---|
452 | else if (BK.countgaps && PERCENT(gaps, sequenceUnits) >= BK.gapbound) { |
---|
453 | DUMPINT(PERCENT(gaps,sequenceUnits)); |
---|
454 | DUMPINT(BK.gapbound); |
---|
455 | consensus_string[o] = '-'; |
---|
456 | } |
---|
457 | else { |
---|
458 | char kchar = 0; // character for consensus |
---|
459 | int kcount = 0; // count this character |
---|
460 | |
---|
461 | if (BK.group) { // simplify using IUPAC |
---|
462 | if (ali_type == GB_AT_RNA || ali_type == GB_AT_DNA) { |
---|
463 | int no_of_bases = 0; // count of different bases used to create iupac |
---|
464 | char used_bases[MAX_BASES_TABLES+1]; // string containing those bases |
---|
465 | |
---|
466 | for (int j=0; j<used_bases_tables; j++) { |
---|
467 | int bchar = index_to_upperChar(j); |
---|
468 | |
---|
469 | if (!isGap(bchar)) { |
---|
470 | if (base[j]>0 && PERCENT(base[j],bases) >= BK.considbound) { |
---|
471 | #if defined(DEBUG_CONSENSUS) |
---|
472 | if (!kcount) DUMPINT(BK.considbound); |
---|
473 | #endif |
---|
474 | used_bases[no_of_bases++] = index_to_upperChar(j); |
---|
475 | kcount += base[j]; |
---|
476 | |
---|
477 | DUMPINT(base[j]); |
---|
478 | DUMPINT(PERCENT(base[j],bases)); |
---|
479 | DUMPINT(kcount); |
---|
480 | } |
---|
481 | } |
---|
482 | } |
---|
483 | used_bases[no_of_bases] = 0; |
---|
484 | kchar = iupac::encode(used_bases, ali_type); |
---|
485 | } |
---|
486 | else { |
---|
487 | ct_assert(ali_type == GB_AT_AA); |
---|
488 | |
---|
489 | const int amino_groups = iupac::AA_GROUP_COUNT; |
---|
490 | int group_count[amino_groups]; |
---|
491 | |
---|
492 | for (int j=0; j<amino_groups; j++) { |
---|
493 | group_count[j] = 0; |
---|
494 | } |
---|
495 | for (int j=0; j<used_bases_tables; j++) { |
---|
496 | unsigned char bchar = index_to_upperChar(j); |
---|
497 | |
---|
498 | if (!isGap(bchar)) { |
---|
499 | group_count[iupac::get_amino_group_for(bchar)] += base[j]; |
---|
500 | } |
---|
501 | } |
---|
502 | |
---|
503 | int bestGroup = 0; |
---|
504 | for (int j=0; j<amino_groups; j++) { |
---|
505 | if (group_count[j]>kcount) { |
---|
506 | bestGroup = j; |
---|
507 | kcount = group_count[j]; |
---|
508 | } |
---|
509 | } |
---|
510 | |
---|
511 | ct_assert(kcount>0); |
---|
512 | if (PERCENT(kcount, bases) >= BK.considbound) { |
---|
513 | kchar = iupac::get_amino_consensus_char(iupac::Amino_Group(bestGroup)); |
---|
514 | } |
---|
515 | else { |
---|
516 | kchar = 'X'; |
---|
517 | kcount = bases; |
---|
518 | } |
---|
519 | } |
---|
520 | } |
---|
521 | if (!kcount) { // IUPAC grouping is either off OR didnt consider any bases |
---|
522 | ct_assert(max_base); // expect at least one base to occur |
---|
523 | ct_assert(max_base_idx >= 0); |
---|
524 | |
---|
525 | kchar = index_to_upperChar(max_base_idx); |
---|
526 | kcount = max_base; |
---|
527 | } |
---|
528 | |
---|
529 | ct_assert(kchar); |
---|
530 | ct_assert(kcount); |
---|
531 | ct_assert(kcount<=bases); |
---|
532 | |
---|
533 | // show as upper or lower case ? |
---|
534 | int percent = PERCENT(kcount, BK.countgaps ? sequenceUnits : bases); |
---|
535 | DUMPINT(percent); |
---|
536 | DUMPINT(BK.upper); |
---|
537 | DUMPINT(BK.lower); |
---|
538 | if (percent>=BK.upper) { |
---|
539 | consensus_string[o] = kchar; |
---|
540 | } |
---|
541 | else if (percent>=BK.lower) { |
---|
542 | consensus_string[o] = tolower(kchar); |
---|
543 | } |
---|
544 | else { |
---|
545 | consensus_string[o] = '.'; |
---|
546 | } |
---|
547 | } |
---|
548 | ct_assert(consensus_string[o]); |
---|
549 | |
---|
550 | #if defined(DEBUG_CONSENSUS) |
---|
551 | if (dumpcol) fprintf(stderr, "result='%c'\n", consensus_string[o]); |
---|
552 | #endif |
---|
553 | } |
---|
554 | } |
---|
555 | else { |
---|
556 | memset(consensus_string, '?', right_idx-left_idx+1); |
---|
557 | } |
---|
558 | #undef PERCENT |
---|
559 | } |
---|
560 | |
---|
561 | // ----------------------- |
---|
562 | // BaseFrequencies |
---|
563 | |
---|
564 | bool BaseFrequencies::initialized = false; |
---|
565 | Ambiguity BaseFrequencies::ambiguity_table[MAX_AMBIGUITY_CODES]; |
---|
566 | uint8_t BaseFrequencies::unitsPerSequence = 1; |
---|
567 | unsigned char BaseFrequencies::char_to_index_tab[MAXCHARTABLE]; |
---|
568 | bool BaseFrequencies::is_gap[MAXCHARTABLE]; |
---|
569 | unsigned char BaseFrequencies::upper_index_chars[MAX_USED_BASES_TABLES]; |
---|
570 | int BaseFrequencies::used_bases_tables = 0; |
---|
571 | GB_alignment_type BaseFrequencies::ali_type = GB_AT_RNA; |
---|
572 | |
---|
573 | inline void BaseFrequencies::set_char_to_index(unsigned char c, int index) |
---|
574 | { |
---|
575 | ct_assert(index>=0 && index<used_bases_tables); |
---|
576 | char_to_index_tab[upper_index_chars[index] = toupper(c)] = index; |
---|
577 | char_to_index_tab[tolower(c)] = index; |
---|
578 | } |
---|
579 | |
---|
580 | void BaseFrequencies::expand_tables() { |
---|
581 | int i; |
---|
582 | for (i=0; i<used_bases_tables; i++) { |
---|
583 | linear_table(i).expand_table_entry_size(); |
---|
584 | } |
---|
585 | } |
---|
586 | |
---|
587 | void BaseFrequencies::setup(const char *gap_chars, GB_alignment_type ali_type_) { |
---|
588 | const char *groups = 0; |
---|
589 | const char *ambiguous = 0; |
---|
590 | |
---|
591 | ali_type = ali_type_; |
---|
592 | |
---|
593 | switch (ali_type) { |
---|
594 | case GB_AT_RNA: |
---|
595 | groups = "A,C,G,TU,"; // Note: terminal ',' defines a group for unknown characters |
---|
596 | ambiguous = "MRWSYKVHDBN"; |
---|
597 | unitsPerSequence = 12; |
---|
598 | break; |
---|
599 | |
---|
600 | case GB_AT_DNA: |
---|
601 | groups = "A,C,G,UT,"; |
---|
602 | ambiguous = "MRWSYKVHDBN"; |
---|
603 | unitsPerSequence = 12; |
---|
604 | break; |
---|
605 | |
---|
606 | case GB_AT_AA: |
---|
607 | 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_...) |
---|
608 | unitsPerSequence = 1; |
---|
609 | break; |
---|
610 | |
---|
611 | default: |
---|
612 | ct_assert(0); |
---|
613 | break; |
---|
614 | } |
---|
615 | |
---|
616 | ct_assert(groups); |
---|
617 | #if defined(ASSERTION_USED) |
---|
618 | if (ali_type == GB_AT_AA) { |
---|
619 | ct_assert(!ambiguous); |
---|
620 | // ct_assert(strlen(ambiguous) <= MAX_AMBIGUITY_CODES); // @@@ later |
---|
621 | } |
---|
622 | else { |
---|
623 | ct_assert(strlen(ambiguous) == MAX_AMBIGUITY_CODES); |
---|
624 | } |
---|
625 | #endif |
---|
626 | |
---|
627 | for (int i = 0; i<MAXCHARTABLE; ++i) is_gap[i] = false; |
---|
628 | |
---|
629 | int gapcharcount = 0; |
---|
630 | char *unique_gap_chars = strdup(gap_chars); |
---|
631 | for (int i = 0; gap_chars[i]; ++i) { |
---|
632 | if (!is_gap[safeCharIndex(gap_chars[i])]) { // ignore duplicate occurrences in 'gap_chars' |
---|
633 | is_gap[safeCharIndex(gap_chars[i])] = true; |
---|
634 | unique_gap_chars[gapcharcount++] = gap_chars[i]; |
---|
635 | } |
---|
636 | } |
---|
637 | |
---|
638 | used_bases_tables = gapcharcount; |
---|
639 | |
---|
640 | for (int i=0; groups[i]; i++) { |
---|
641 | if (groups[i]==',') { |
---|
642 | used_bases_tables++; |
---|
643 | } |
---|
644 | } |
---|
645 | |
---|
646 | ct_assert(used_bases_tables<=MAX_USED_BASES_TABLES); |
---|
647 | |
---|
648 | int idx = 0; |
---|
649 | unsigned char init_val = used_bases_tables-1; // map all unknown stuff to last table |
---|
650 | memset(char_to_index_tab, init_val, MAXCHARTABLE*sizeof(*char_to_index_tab)); |
---|
651 | |
---|
652 | for (int i=0; groups[i]; i++) { |
---|
653 | if (groups[i]==',') { |
---|
654 | idx++; |
---|
655 | } |
---|
656 | else { |
---|
657 | ct_assert(isupper(groups[i])); |
---|
658 | set_char_to_index(groups[i], idx); |
---|
659 | } |
---|
660 | } |
---|
661 | |
---|
662 | const char *align_string_ptr = unique_gap_chars; |
---|
663 | while (1) { |
---|
664 | char c = *align_string_ptr++; |
---|
665 | if (!c) break; |
---|
666 | set_char_to_index(c, idx++); |
---|
667 | } |
---|
668 | |
---|
669 | if (ambiguous) { |
---|
670 | ct_assert(ali_type == GB_AT_DNA || ali_type == GB_AT_RNA); // @@@ amino ambiguities not impl |
---|
671 | |
---|
672 | const uint8_t indices2increment[MAX_TARGET_INDICES+1] = { 0, 0, 6, 4, 3 }; |
---|
673 | for (int i = 0; ambiguous[i]; ++i) { |
---|
674 | const char *contained = iupac::decode(ambiguous[i], ali_type, false); |
---|
675 | |
---|
676 | Ambiguity& amb = ambiguity_table[i]; |
---|
677 | |
---|
678 | amb.indices = strlen(contained); |
---|
679 | ct_assert(amb.indices>=2 && amb.indices<=MAX_TARGET_INDICES); |
---|
680 | amb.increment = indices2increment[amb.indices]; |
---|
681 | |
---|
682 | for (int j = 0; j<amb.indices; ++j) { |
---|
683 | int cidx = char_to_index_tab[safeCharIndex(contained[j])]; |
---|
684 | ct_assert(cidx<MAX_INDEX_TABLES); // has to be non-ambiguous |
---|
685 | amb.index[j] = cidx; |
---|
686 | } |
---|
687 | |
---|
688 | char_to_index_tab[toupper(ambiguous[i])] = |
---|
689 | char_to_index_tab[tolower(ambiguous[i])] = i+MAX_INDEX_TABLES; |
---|
690 | } |
---|
691 | } |
---|
692 | |
---|
693 | free(unique_gap_chars); |
---|
694 | ct_assert(idx==used_bases_tables); |
---|
695 | initialized = true; |
---|
696 | } |
---|
697 | |
---|
698 | BaseFrequencies::BaseFrequencies(int maxseqlength) |
---|
699 | : ignore(0) |
---|
700 | { |
---|
701 | ct_assert(initialized); |
---|
702 | |
---|
703 | bases_table = new SepBaseFreqPtr[used_bases_tables]; |
---|
704 | |
---|
705 | for (int i=0; i<used_bases_tables; i++) { |
---|
706 | bases_table[i] = new SepBaseFreq(maxseqlength); |
---|
707 | } |
---|
708 | |
---|
709 | sequenceUnits = 0; |
---|
710 | } |
---|
711 | |
---|
712 | void BaseFrequencies::init(int maxseqlength) |
---|
713 | { |
---|
714 | int i; |
---|
715 | for (i=0; i<used_bases_tables; i++) { |
---|
716 | linear_table(i).init(maxseqlength); |
---|
717 | } |
---|
718 | |
---|
719 | sequenceUnits = 0; |
---|
720 | } |
---|
721 | |
---|
722 | void BaseFrequencies::bases_and_gaps_at(int column, int *bases, int *gaps) const |
---|
723 | { |
---|
724 | int b = 0; |
---|
725 | int g = 0; |
---|
726 | |
---|
727 | for (int i=0; i<used_bases_tables; i++) { |
---|
728 | char c = upper_index_chars[i]; |
---|
729 | |
---|
730 | if (isGap(c)) { |
---|
731 | g += table(c)[column]; |
---|
732 | } |
---|
733 | else { |
---|
734 | b += table(c)[column]; |
---|
735 | } |
---|
736 | } |
---|
737 | |
---|
738 | if (bases) { |
---|
739 | *bases = b/unitsPerSequence; |
---|
740 | ct_assert((b%unitsPerSequence) == 0); // could happen if an ambiguous code contains a gap |
---|
741 | } |
---|
742 | if (gaps) { |
---|
743 | *gaps = g/unitsPerSequence; |
---|
744 | ct_assert((g%unitsPerSequence) == 0); // could happen if an ambiguous code contains a gap |
---|
745 | } |
---|
746 | } |
---|
747 | |
---|
748 | double BaseFrequencies::max_frequency_at(int column, bool ignore_gaps) const { |
---|
749 | int gaps = 0; |
---|
750 | int minus = 0; |
---|
751 | int maxb = 0; |
---|
752 | |
---|
753 | for (int i=0; i<used_bases_tables; i++) { |
---|
754 | char c = upper_index_chars[i]; |
---|
755 | |
---|
756 | if (isGap(c)) { |
---|
757 | if (c == '-') { |
---|
758 | minus += table(c)[column]; |
---|
759 | } |
---|
760 | else { |
---|
761 | gaps += table(c)[column]; |
---|
762 | } |
---|
763 | } |
---|
764 | else { |
---|
765 | maxb = std::max(maxb, table(c)[column]); |
---|
766 | } |
---|
767 | } |
---|
768 | |
---|
769 | if (ignore_gaps) { |
---|
770 | gaps += minus; |
---|
771 | } |
---|
772 | else { |
---|
773 | maxb = std::max(maxb, minus); |
---|
774 | } |
---|
775 | |
---|
776 | double mfreq; |
---|
777 | if (sequenceUnits>gaps) { |
---|
778 | mfreq = maxb/double(sequenceUnits-gaps); |
---|
779 | } |
---|
780 | else { // only gaps (but no '-') |
---|
781 | ct_assert(sequenceUnits == gaps); |
---|
782 | mfreq = 0.0; |
---|
783 | } |
---|
784 | |
---|
785 | #if defined(ASSERTION_USED) |
---|
786 | if (mfreq != 0.0) { |
---|
787 | if (ali_type != GB_AT_AA) { |
---|
788 | ct_assert(mfreq>=0.2); // wrong for amino acids |
---|
789 | } |
---|
790 | ct_assert(mfreq<=1.0); |
---|
791 | } |
---|
792 | #endif |
---|
793 | |
---|
794 | return mfreq; |
---|
795 | } |
---|
796 | |
---|
797 | BaseFrequencies::~BaseFrequencies() |
---|
798 | { |
---|
799 | int i; |
---|
800 | for (i=0; i<used_bases_tables; i++) { |
---|
801 | delete bases_table[i]; |
---|
802 | } |
---|
803 | |
---|
804 | delete [] bases_table; |
---|
805 | } |
---|
806 | |
---|
807 | const PosRange *BaseFrequencies::changed_range(const BaseFrequencies& other) const |
---|
808 | { |
---|
809 | int i; |
---|
810 | int Size = size(); |
---|
811 | int start = Size-1; |
---|
812 | int end = 0; |
---|
813 | |
---|
814 | ct_assert(Size==other.size()); |
---|
815 | for (i=0; i<used_bases_tables; i++) { |
---|
816 | if (linear_table(i).firstDifference(other.linear_table(i), 0, start, &start)) { |
---|
817 | if (end<start) { |
---|
818 | end = start; |
---|
819 | } |
---|
820 | linear_table(i).lastDifference(other.linear_table(i), end+1, Size-1, &end); |
---|
821 | |
---|
822 | for (; i<used_bases_tables; i++) { |
---|
823 | linear_table(i).firstDifference(other.linear_table(i), 0, start-1, &start); |
---|
824 | if (end<start) { |
---|
825 | end = start; |
---|
826 | } |
---|
827 | linear_table(i).lastDifference(other.linear_table(i), end+1, Size-1, &end); |
---|
828 | } |
---|
829 | |
---|
830 | ct_assert(start<=end); |
---|
831 | |
---|
832 | static PosRange range; |
---|
833 | range = PosRange(start, end); |
---|
834 | |
---|
835 | return ⦥ |
---|
836 | } |
---|
837 | } |
---|
838 | return NULL; |
---|
839 | } |
---|
840 | void BaseFrequencies::add(const BaseFrequencies& other) |
---|
841 | { |
---|
842 | if (other.ignore) return; |
---|
843 | if (other.size()==0) return; |
---|
844 | prepare_to_add_elements(other.added_sequences()); |
---|
845 | add(other, 0, other.size()-1); |
---|
846 | } |
---|
847 | void BaseFrequencies::add(const BaseFrequencies& other, int start, int end) |
---|
848 | { |
---|
849 | if (other.ignore) return; |
---|
850 | |
---|
851 | test(); |
---|
852 | other.test(); |
---|
853 | |
---|
854 | ct_assert(start<=end); |
---|
855 | |
---|
856 | int i; |
---|
857 | for (i=0; i<used_bases_tables; i++) { |
---|
858 | linear_table(i).add(other.linear_table(i), start, end); |
---|
859 | } |
---|
860 | |
---|
861 | sequenceUnits += other.sequenceUnits; |
---|
862 | |
---|
863 | test(); |
---|
864 | } |
---|
865 | |
---|
866 | void BaseFrequencies::sub(const BaseFrequencies& other) |
---|
867 | { |
---|
868 | if (other.ignore) return; |
---|
869 | sub(other, 0, other.size()-1); |
---|
870 | } |
---|
871 | |
---|
872 | void BaseFrequencies::sub(const BaseFrequencies& other, int start, int end) |
---|
873 | { |
---|
874 | if (other.ignore) return; |
---|
875 | |
---|
876 | test(); |
---|
877 | other.test(); |
---|
878 | ct_assert(start<=end); |
---|
879 | |
---|
880 | int i; |
---|
881 | for (i=0; i<used_bases_tables; i++) { |
---|
882 | linear_table(i).sub(other.linear_table(i), start, end); |
---|
883 | } |
---|
884 | |
---|
885 | sequenceUnits -= other.sequenceUnits; |
---|
886 | |
---|
887 | test(); |
---|
888 | } |
---|
889 | |
---|
890 | void BaseFrequencies::sub_and_add(const BaseFrequencies& Sub, const BaseFrequencies& Add, PosRange range) { |
---|
891 | test(); |
---|
892 | Sub.test(); |
---|
893 | Add.test(); |
---|
894 | |
---|
895 | ct_assert(!Sub.ignore && !Add.ignore); |
---|
896 | ct_assert(range.is_part()); |
---|
897 | |
---|
898 | int i; |
---|
899 | for (i=0; i<used_bases_tables; i++) { |
---|
900 | linear_table(i).sub_and_add(Sub.linear_table(i), Add.linear_table(i), range); |
---|
901 | } |
---|
902 | |
---|
903 | test(); |
---|
904 | } |
---|
905 | |
---|
906 | const PosRange *BaseFrequencies::changed_range(const char *string1, const char *string2, int min_len) |
---|
907 | { |
---|
908 | const unsigned long *range1 = (const unsigned long *)string1; |
---|
909 | const unsigned long *range2 = (const unsigned long *)string2; |
---|
910 | const int step = sizeof(*range1); |
---|
911 | |
---|
912 | int l = min_len/step; |
---|
913 | int r = min_len%step; |
---|
914 | int i; |
---|
915 | int j; |
---|
916 | int k; |
---|
917 | |
---|
918 | int start = -1, end = -1; |
---|
919 | |
---|
920 | for (i=0; i<l; i++) { // search wordwise (it's faster) |
---|
921 | if (range1[i] != range2[i]) { |
---|
922 | k = i*step; |
---|
923 | for (j=0; j<step; j++) { |
---|
924 | if (string1[k+j] != string2[k+j]) { |
---|
925 | start = end = k+j; |
---|
926 | break; |
---|
927 | } |
---|
928 | } |
---|
929 | |
---|
930 | break; |
---|
931 | } |
---|
932 | } |
---|
933 | |
---|
934 | k = l*step; |
---|
935 | if (i==l) { // no difference found in word -> look at rest |
---|
936 | for (j=0; j<r; j++) { |
---|
937 | if (string1[k+j] != string2[k+j]) { |
---|
938 | start = end = k+j; |
---|
939 | break; |
---|
940 | } |
---|
941 | } |
---|
942 | |
---|
943 | if (j==r) { // the strings are equal |
---|
944 | return 0; |
---|
945 | } |
---|
946 | } |
---|
947 | |
---|
948 | ct_assert(start != -1 && end != -1); |
---|
949 | |
---|
950 | for (j=r-1; j>=0; j--) { // search rest backwards |
---|
951 | if (string1[k+j] != string2[k+j]) { |
---|
952 | end = k+j; |
---|
953 | break; |
---|
954 | } |
---|
955 | } |
---|
956 | |
---|
957 | if (j==-1) { // not found in rest -> search words backwards |
---|
958 | int m = start/step; |
---|
959 | for (i=l-1; i>=m; i--) { |
---|
960 | if (range1[i] != range2[i]) { |
---|
961 | k = i*step; |
---|
962 | for (j=step-1; j>=0; j--) { |
---|
963 | if (string1[k+j] != string2[k+j]) { |
---|
964 | end = k+j; |
---|
965 | break; |
---|
966 | } |
---|
967 | } |
---|
968 | break; |
---|
969 | } |
---|
970 | } |
---|
971 | } |
---|
972 | |
---|
973 | ct_assert(start<=end); |
---|
974 | |
---|
975 | static PosRange range; |
---|
976 | range = PosRange(start, end); |
---|
977 | |
---|
978 | return ⦥ |
---|
979 | } |
---|
980 | |
---|
981 | #define INC_DEC(incdec) do { \ |
---|
982 | int idx = char_to_index_tab[c]; \ |
---|
983 | if (idx<MAX_INDEX_TABLES) { \ |
---|
984 | table(c).incdec(offset, unitsPerSequence); \ |
---|
985 | } \ |
---|
986 | else { \ |
---|
987 | Ambiguity& amb = ambiguity_table[idx-MAX_INDEX_TABLES]; \ |
---|
988 | for (int i = 0; i<amb.indices; ++i) { \ |
---|
989 | linear_table(amb.index[i]).incdec(offset, amb.increment); \ |
---|
990 | } \ |
---|
991 | } \ |
---|
992 | } while(0) |
---|
993 | |
---|
994 | void BaseFrequencies::incr_short(unsigned char c, int offset) { |
---|
995 | INC_DEC(inc_short); |
---|
996 | } |
---|
997 | |
---|
998 | void BaseFrequencies::decr_short(unsigned char c, int offset) { |
---|
999 | INC_DEC(dec_short); |
---|
1000 | } |
---|
1001 | |
---|
1002 | void BaseFrequencies::incr_long(unsigned char c, int offset) { |
---|
1003 | INC_DEC(inc_long); |
---|
1004 | } |
---|
1005 | |
---|
1006 | void BaseFrequencies::decr_long(unsigned char c, int offset) { |
---|
1007 | INC_DEC(dec_long); |
---|
1008 | } |
---|
1009 | |
---|
1010 | #undef INC_DEC |
---|
1011 | |
---|
1012 | void BaseFrequencies::add(const char *scan_string, int len) |
---|
1013 | { |
---|
1014 | test(); |
---|
1015 | |
---|
1016 | prepare_to_add_elements(1); |
---|
1017 | |
---|
1018 | int i; |
---|
1019 | int sz = size(); |
---|
1020 | |
---|
1021 | if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) { |
---|
1022 | for (i=0; i<len; i++) { |
---|
1023 | unsigned char c = scan_string[i]; |
---|
1024 | if (!c) break; |
---|
1025 | incr_short(c, i); |
---|
1026 | } |
---|
1027 | SepBaseFreq& t = table('.'); |
---|
1028 | for (; i<sz; i++) { |
---|
1029 | t.inc_short(i, unitsPerSequence); |
---|
1030 | } |
---|
1031 | } |
---|
1032 | else { |
---|
1033 | for (i=0; i<len; i++) { |
---|
1034 | unsigned char c = scan_string[i]; |
---|
1035 | if (!c) break; |
---|
1036 | incr_long(c, i); |
---|
1037 | } |
---|
1038 | SepBaseFreq& t = table('.'); |
---|
1039 | for (; i<sz; i++) { // LOOP_VECTORIZED |
---|
1040 | t.inc_long(i, unitsPerSequence); |
---|
1041 | } |
---|
1042 | } |
---|
1043 | |
---|
1044 | sequenceUnits += unitsPerSequence; |
---|
1045 | |
---|
1046 | test(); |
---|
1047 | } |
---|
1048 | |
---|
1049 | void BaseFrequencies::sub(const char *scan_string, int len) |
---|
1050 | { |
---|
1051 | test(); |
---|
1052 | |
---|
1053 | int i; |
---|
1054 | int sz = size(); |
---|
1055 | |
---|
1056 | if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) { |
---|
1057 | for (i=0; i<len; i++) { |
---|
1058 | unsigned char c = scan_string[i]; |
---|
1059 | if (!c) break; |
---|
1060 | decr_short(c, i); |
---|
1061 | } |
---|
1062 | SepBaseFreq& t = table('.'); |
---|
1063 | for (; i<sz; i++) { |
---|
1064 | t.dec_short(i, unitsPerSequence); |
---|
1065 | } |
---|
1066 | } |
---|
1067 | else { |
---|
1068 | for (i=0; i<len; i++) { |
---|
1069 | unsigned char c = scan_string[i]; |
---|
1070 | if (!c) break; |
---|
1071 | decr_long(c, i); |
---|
1072 | } |
---|
1073 | SepBaseFreq& t = table('.'); |
---|
1074 | for (; i<sz; i++) { // LOOP_VECTORIZED |
---|
1075 | t.dec_long(i, unitsPerSequence); |
---|
1076 | } |
---|
1077 | } |
---|
1078 | |
---|
1079 | sequenceUnits -= unitsPerSequence; |
---|
1080 | |
---|
1081 | test(); |
---|
1082 | } |
---|
1083 | |
---|
1084 | void BaseFrequencies::sub_and_add(const char *old_string, const char *new_string, PosRange range) { |
---|
1085 | test(); |
---|
1086 | int start = range.start(); |
---|
1087 | int end = range.end(); |
---|
1088 | ct_assert(start<=end); |
---|
1089 | |
---|
1090 | if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) { |
---|
1091 | for (int i=start; i<=end; i++) { |
---|
1092 | unsigned char o = old_string[i]; |
---|
1093 | unsigned char n = new_string[i]; |
---|
1094 | |
---|
1095 | ct_assert(o && n); // passed string have to be defined in range |
---|
1096 | |
---|
1097 | if (o!=n) { |
---|
1098 | decr_short(o, i); |
---|
1099 | incr_short(n, i); |
---|
1100 | } |
---|
1101 | } |
---|
1102 | } |
---|
1103 | else { |
---|
1104 | for (int i=start; i<=end; i++) { |
---|
1105 | unsigned char o = old_string[i]; |
---|
1106 | unsigned char n = new_string[i]; |
---|
1107 | |
---|
1108 | ct_assert(o && n); // passed string have to be defined in range |
---|
1109 | |
---|
1110 | if (o!=n) { |
---|
1111 | decr_long(o, i); |
---|
1112 | incr_long(n, i); |
---|
1113 | } |
---|
1114 | } |
---|
1115 | } |
---|
1116 | |
---|
1117 | test(); |
---|
1118 | } |
---|
1119 | |
---|
1120 | void BaseFrequencies::change_table_length(int new_length) { |
---|
1121 | for (int c=0; c<used_bases_tables; c++) { |
---|
1122 | linear_table(c).change_table_length(new_length, 0); |
---|
1123 | } |
---|
1124 | test(); |
---|
1125 | } |
---|
1126 | |
---|
1127 | // -------------- |
---|
1128 | // tests: |
---|
1129 | |
---|
1130 | #if defined(TEST_CHAR_TABLE_INTEGRITY) || defined(ASSERTION_USED) |
---|
1131 | |
---|
1132 | bool BaseFrequencies::empty() const |
---|
1133 | { |
---|
1134 | int i; |
---|
1135 | for (i=0; i<used_bases_tables; i++) { |
---|
1136 | if (!linear_table(i).empty()) { |
---|
1137 | return false; |
---|
1138 | } |
---|
1139 | } |
---|
1140 | return true; |
---|
1141 | } |
---|
1142 | |
---|
1143 | bool BaseFrequencies::ok() const |
---|
1144 | { |
---|
1145 | if (empty()) return true; |
---|
1146 | if (is_ignored()) return true; |
---|
1147 | |
---|
1148 | if (sequenceUnits<0) { |
---|
1149 | fprintf(stderr, "Negative # of sequences (%i) in BaseFrequencies\n", added_sequences()); |
---|
1150 | return false; |
---|
1151 | } |
---|
1152 | |
---|
1153 | const int table_size = size(); |
---|
1154 | for (int i=0; i<table_size; i++) { |
---|
1155 | int bases, gaps; |
---|
1156 | bases_and_gaps_at(i, &bases, &gaps); |
---|
1157 | |
---|
1158 | if (!bases && !gaps) { // this occurs after insert column |
---|
1159 | for (int j=i+1; j<table_size; j++) { // test if we find any bases till end of table !!! |
---|
1160 | |
---|
1161 | bases_and_gaps_at(j, &bases, 0); |
---|
1162 | if (bases) { // bases found -> empty position was illegal |
---|
1163 | fprintf(stderr, "Zero in BaseFrequencies at column %i\n", i); |
---|
1164 | return false; |
---|
1165 | } |
---|
1166 | } |
---|
1167 | break; |
---|
1168 | } |
---|
1169 | else { |
---|
1170 | if ((bases+gaps)!=added_sequences()) { |
---|
1171 | fprintf(stderr, "bases+gaps should be equal to # of sequences at column %i (bases=%i, gaps=%i, added_sequences=%i)\n", |
---|
1172 | i, bases, gaps, added_sequences()); |
---|
1173 | return false; |
---|
1174 | } |
---|
1175 | } |
---|
1176 | } |
---|
1177 | |
---|
1178 | return true; |
---|
1179 | } |
---|
1180 | |
---|
1181 | #endif // TEST_CHAR_TABLE_INTEGRITY || ASSERTION_USED |
---|
1182 | |
---|
1183 | #if defined(TEST_CHAR_TABLE_INTEGRITY) |
---|
1184 | |
---|
1185 | void BaseFrequencies::test() const { |
---|
1186 | |
---|
1187 | if (!ok()) { |
---|
1188 | GBK_terminate("BaseFrequencies::test() failed"); |
---|
1189 | } |
---|
1190 | } |
---|
1191 | |
---|
1192 | #endif // TEST_CHAR_TABLE_INTEGRITY |
---|
1193 | |
---|
1194 | |
---|
1195 | // -------------------------------------------------------------------------------- |
---|
1196 | |
---|
1197 | #ifdef UNIT_TESTS |
---|
1198 | #ifndef TEST_UNIT_H |
---|
1199 | #include <test_unit.h> |
---|
1200 | #endif |
---|
1201 | |
---|
1202 | #define SETUP(gapChars,alitype) BaseFrequencies::setup(gapChars,alitype) |
---|
1203 | |
---|
1204 | void TEST_char_table() { |
---|
1205 | const char alphabeth[] = "ACGTN-"; |
---|
1206 | const int alphabeth_size = strlen(alphabeth); |
---|
1207 | |
---|
1208 | const int seqlen = 30; |
---|
1209 | char seq[seqlen+1]; |
---|
1210 | |
---|
1211 | const char *gapChars = ".-"; // @@@ doesnt make any difference which gaps are used - why? |
---|
1212 | // const char *gapChars = "-"; |
---|
1213 | // const char *gapChars = "."; |
---|
1214 | // const char *gapChars = "A"; // makes me fail |
---|
1215 | SETUP(gapChars, GB_AT_RNA); |
---|
1216 | |
---|
1217 | ConsensusBuildParams BK; |
---|
1218 | // BK.lower = 70; BK.upper = 95; BK.gapbound = 60; // defaults from awars |
---|
1219 | BK.lower = 40; BK.upper = 70; BK.gapbound = 40; |
---|
1220 | |
---|
1221 | srand(100); |
---|
1222 | for (int loop = 0; loop<5; ++loop) { |
---|
1223 | unsigned seed = rand(); |
---|
1224 | |
---|
1225 | srand(seed); |
---|
1226 | |
---|
1227 | BaseFrequencies tab(seqlen); |
---|
1228 | |
---|
1229 | // build some seq |
---|
1230 | for (int c = 0; c<seqlen; ++c) seq[c] = alphabeth[rand()%alphabeth_size]; |
---|
1231 | seq[seqlen] = 0; |
---|
1232 | |
---|
1233 | TEST_EXPECT_EQUAL(strlen(seq), size_t(seqlen)); |
---|
1234 | |
---|
1235 | const int add_count = 300; |
---|
1236 | char *added_seqs[add_count]; |
---|
1237 | |
---|
1238 | const char *s1 = "ACGTACGTAcgtacgtaCGTACGTACGTAC"; |
---|
1239 | const char *s2 = "MRWSYKVHDBNmrwsykvhdbnSYKVHDBN"; // use ambiguities |
---|
1240 | |
---|
1241 | // add seq multiple times |
---|
1242 | for (int a = 0; a<add_count; ++a) { |
---|
1243 | tab.add(seq, seqlen); |
---|
1244 | added_seqs[a] = strdup(seq); |
---|
1245 | |
---|
1246 | // modify 1 character in seq: |
---|
1247 | int sidx = rand()%seqlen; |
---|
1248 | int aidx = rand()%alphabeth_size; |
---|
1249 | seq[sidx] = alphabeth[aidx]; |
---|
1250 | |
---|
1251 | if (a == 16) { // with 15 sequences in tab |
---|
1252 | // test sub_and_add (short table version) |
---|
1253 | char *consensus = tab.build_consensus_string(BK); |
---|
1254 | |
---|
1255 | tab.add(s1, seqlen); |
---|
1256 | PosRange r(0, seqlen-1); |
---|
1257 | tab.sub_and_add(s1, s2, r); |
---|
1258 | tab.sub_and_add(s2, consensus, r); |
---|
1259 | tab.sub(consensus, seqlen); |
---|
1260 | |
---|
1261 | char *consensus2 = tab.build_consensus_string(BK); |
---|
1262 | TEST_EXPECT_EQUAL(consensus, consensus2); |
---|
1263 | |
---|
1264 | free(consensus2); |
---|
1265 | free(consensus); |
---|
1266 | } |
---|
1267 | } |
---|
1268 | |
---|
1269 | // build consensi (just check regression) |
---|
1270 | // Note: semantic tests for consensus are in ../../NTREE/AP_consensus.cxx@TEST_nucleotide_consensus_and_maxFrequency |
---|
1271 | { |
---|
1272 | char *consensus = tab.build_consensus_string(BK); |
---|
1273 | switch (seed) { |
---|
1274 | case 677741240: TEST_EXPECT_EQUAL(consensus, "k-s-cWs.aWu.a.WYa.R.mKcaK.c.ry"); break; |
---|
1275 | case 721151648: TEST_EXPECT_EQUAL(consensus, "awm...Ka.-gUyRW-Sacau.Wa.acy.W"); break; |
---|
1276 | case 345295160: TEST_EXPECT_EQUAL(consensus, "yy-gmwkMS....gucy..Y.Rrr.-gg.K"); break; |
---|
1277 | case 346389111: TEST_EXPECT_EQUAL(consensus, "suwA.yYwSurKc-cSmRaY.mg.Sa-kY."); break; |
---|
1278 | case 367171911: TEST_EXPECT_EQUAL(consensus, ".augaY.u.c-c-cUDaYg.-.cg-cWWsM"); break; |
---|
1279 | default: TEST_EXPECT_EQUAL(consensus, "undef"); |
---|
1280 | } |
---|
1281 | |
---|
1282 | // test sub_and_add() |
---|
1283 | |
---|
1284 | tab.add(s1, seqlen); |
---|
1285 | PosRange r(0, seqlen-1); |
---|
1286 | tab.sub_and_add(s1, s2, r); |
---|
1287 | tab.sub_and_add(s2, consensus, r); |
---|
1288 | tab.sub(consensus, seqlen); |
---|
1289 | |
---|
1290 | char *consensus2 = tab.build_consensus_string(BK); |
---|
1291 | TEST_EXPECT_EQUAL(consensus, consensus2); |
---|
1292 | |
---|
1293 | free(consensus2); |
---|
1294 | free(consensus); |
---|
1295 | } |
---|
1296 | |
---|
1297 | for (int a = 0; a<add_count; a++) { |
---|
1298 | tab.sub(added_seqs[a], seqlen); |
---|
1299 | freenull(added_seqs[a]); |
---|
1300 | } |
---|
1301 | { |
---|
1302 | char *consensus = tab.build_consensus_string(BK); |
---|
1303 | TEST_EXPECT_EQUAL(consensus, "??????????????????????????????"); // check tab is empty |
---|
1304 | free(consensus); |
---|
1305 | } |
---|
1306 | } |
---|
1307 | } |
---|
1308 | |
---|
1309 | void TEST_SepBaseFreq() { |
---|
1310 | const int LEN = 20; |
---|
1311 | const int OFF1 = 6; |
---|
1312 | const int OFF2 = 7; // adjacent to OFF1 |
---|
1313 | |
---|
1314 | SepBaseFreq short1(LEN), short2(LEN); |
---|
1315 | for (int i = 0; i<100; ++i) short1.inc_short(OFF1, 1); |
---|
1316 | for (int i = 0; i<70; ++i) short1.inc_short(OFF2, 1); |
---|
1317 | for (int i = 0; i<150; ++i) short2.inc_short(OFF1, 1); |
---|
1318 | for (int i = 0; i<80; ++i) short2.inc_short(OFF2, 1); |
---|
1319 | |
---|
1320 | SepBaseFreq long1(LEN), long2(LEN); |
---|
1321 | long1.expand_table_entry_size(); |
---|
1322 | long2.expand_table_entry_size(); |
---|
1323 | for (int i = 0; i<2000; ++i) long1.inc_long(OFF1, 1); |
---|
1324 | for (int i = 0; i<2500; ++i) long2.inc_long(OFF1, 1); |
---|
1325 | |
---|
1326 | SepBaseFreq shortsum(LEN); |
---|
1327 | shortsum.add(short1, 0, LEN-1); TEST_EXPECT_EQUAL(shortsum[OFF1], 100); TEST_EXPECT_EQUAL(shortsum[OFF2], 70); |
---|
1328 | shortsum.add(short2, 0, LEN-1); TEST_EXPECT_EQUAL(shortsum[OFF1], 250); TEST_EXPECT_EQUAL(shortsum[OFF2], 150); |
---|
1329 | shortsum.sub(short1, 0, LEN-1); TEST_EXPECT_EQUAL(shortsum[OFF1], 150); TEST_EXPECT_EQUAL(shortsum[OFF2], 80); |
---|
1330 | |
---|
1331 | shortsum.sub_and_add(short1, short2, PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(shortsum[OFF1], 200); TEST_EXPECT_EQUAL(shortsum[OFF2], 90); |
---|
1332 | |
---|
1333 | // shortsum.add(long1, 0, LEN-1); // invalid operation -> terminate |
---|
1334 | |
---|
1335 | SepBaseFreq longsum(LEN); |
---|
1336 | longsum.expand_table_entry_size(); |
---|
1337 | longsum.add(long1, 0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 2000); |
---|
1338 | longsum.add(long2, 0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 4500); |
---|
1339 | longsum.sub(long1, 0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 2500); |
---|
1340 | longsum.add(short1, 0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 2600); |
---|
1341 | longsum.sub(short2, 0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 2450); |
---|
1342 | |
---|
1343 | longsum.sub_and_add(long1, long2, PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(longsum[OFF1], 2950); |
---|
1344 | longsum.sub_and_add(short1, short2, PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(longsum[OFF1], 3000); |
---|
1345 | longsum.sub_and_add(long1, short2, PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(longsum[OFF1], 1150); |
---|
1346 | longsum.sub_and_add(short1, long2, PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(longsum[OFF1], 3550); |
---|
1347 | |
---|
1348 | int pos = -1; |
---|
1349 | TEST_EXPECT_EQUAL(short1.firstDifference(short2, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1); |
---|
1350 | TEST_EXPECT_EQUAL(short1.lastDifference (short2, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF2); |
---|
1351 | TEST_EXPECT_EQUAL(short1.firstDifference(long1, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1); |
---|
1352 | TEST_EXPECT_EQUAL(short1.lastDifference (long1, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF2); |
---|
1353 | TEST_EXPECT_EQUAL(long2.firstDifference (short2, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1); |
---|
1354 | TEST_EXPECT_EQUAL(long2.lastDifference (short2, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF2); |
---|
1355 | TEST_EXPECT_EQUAL(long2.firstDifference (long1, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1); |
---|
1356 | TEST_EXPECT_EQUAL(long2.lastDifference (long1, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1); |
---|
1357 | } |
---|
1358 | |
---|
1359 | #endif // UNIT_TESTS |
---|
1360 | |
---|
1361 | // -------------------------------------------------------------------------------- |
---|
1362 | |
---|
1363 | |
---|
1364 | |
---|