source: tags/ms_r16q3/SL/CONSENSUS/chartable.cxx

Last change on this file was 15176, checked in by westram, 8 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 41.2 KB
Line 
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
23using 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
33static long bases_table_size = 0L;
34#endif
35
36void 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
58void 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
80void 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}
111void 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}
142void 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
195int 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}
243int 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
293SepBaseFreq::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
301SepBaseFreq::~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
310void 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
355int 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
378char *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
394void 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
564bool              BaseFrequencies::initialized       = false;
565Ambiguity         BaseFrequencies::ambiguity_table[MAX_AMBIGUITY_CODES];
566uint8_t           BaseFrequencies::unitsPerSequence  = 1;
567unsigned char     BaseFrequencies::char_to_index_tab[MAXCHARTABLE];
568bool              BaseFrequencies::is_gap[MAXCHARTABLE];
569unsigned char     BaseFrequencies::upper_index_chars[MAX_USED_BASES_TABLES];
570int               BaseFrequencies::used_bases_tables = 0;
571GB_alignment_type BaseFrequencies::ali_type          = GB_AT_RNA;
572
573inline 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
580void 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
587void 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
698BaseFrequencies::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
712void 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
722void 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
748double 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
797BaseFrequencies::~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
807const 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 &range;
836        }
837    }
838    return NULL;
839}
840void 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}
847void 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
866void BaseFrequencies::sub(const BaseFrequencies& other)
867{
868    if (other.ignore) return;
869    sub(other, 0, other.size()-1);
870}
871
872void 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
890void 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
906const 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 &range;
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
994void BaseFrequencies::incr_short(unsigned char c, int offset) {
995    INC_DEC(inc_short);
996}
997
998void BaseFrequencies::decr_short(unsigned char c, int offset) {
999    INC_DEC(dec_short);
1000}
1001
1002void BaseFrequencies::incr_long(unsigned char c, int offset) {
1003    INC_DEC(inc_long);
1004}
1005
1006void BaseFrequencies::decr_long(unsigned char c, int offset) {
1007    INC_DEC(dec_long);
1008}
1009
1010#undef INC_DEC
1011
1012void 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
1049void 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
1084void 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
1120void 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
1132bool 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
1143bool 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
1185void 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
1204void 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
1309void 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
Note: See TracBrowser for help on using the repository browser.