source: branches/ali/PHYLO/PH_filt.cxx

Last change on this file was 19206, checked in by westram, 3 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.0 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : PH_filt.cxx                                       //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "phylo.hxx"
12#include "phwin.hxx"
13#include "PH_display.hxx"
14#include <arbdbt.h>
15#include <aw_awar.hxx>
16#include <arb_progress.h>
17#include <aw_root.hxx>
18#include <aw_msg.hxx>
19#include <cctype>
20#include <awt_config_manager.hxx>
21
22static long PH_timer() {
23    static long time = 0;
24    return ++time;
25}
26
27PH_filter::PH_filter() {
28    memset ((char *)this, 0, sizeof(PH_filter));
29}
30
31char *PH_filter::init(long size) {
32    delete [] filter;
33    filter = new char[size];
34    real_len = filter_len = size;
35    for (int i = 0; i < size; i++) {
36        filter[i] = 1;
37    }
38    update = PH_timer();
39    return NULp;
40}
41
42PH_filter::~PH_filter() {
43    delete [] filter;
44}
45
46inline void strlwr(char *s) {
47    for (int i = 0; s[i]; ++i) {
48        s[i] = tolower(s[i]);
49    }
50}
51
52float *PH_filter::calculate_column_homology() {
53    long           i, j, max, num_all_chars;
54    bool           mask[256];
55    char           delete_when_max[100], get_maximum_from[100], all_chars[100], max_char;
56    long           reference_table[256], **chars_counted;
57    char           real_chars[100], low_chars[100];
58    char           rest_chars[40], upper_rest_chars[20], low_rest_chars[20];
59    unsigned char *sequence_buffer;
60    AW_root       *aw_root;
61    float         *mline = NULp;
62
63    if (!PHDATA::ROOT) return NULp; // nothing loaded yet
64
65    arb_progress filt_progress("Calculating filter");
66
67    GB_transaction ta(PHDATA::ROOT->get_gb_main());
68
69    bool isNUC = true; // rna oder dna sequence : nur zum testen und Entwicklung
70    if (GBT_is_alignment_protein(PHDATA::ROOT->get_gb_main(), PHDATA::ROOT->use)) {
71        isNUC = false;
72    }
73
74    aw_root=PH_used_windows::windowList->phylo_main_window->get_root();
75
76    if (isNUC) {
77        strcpy(real_chars, "ACGTU");
78        strcpy(upper_rest_chars, "MRWSYKVHDBXN");
79    }
80    else {
81        strcpy(real_chars, "ABCDEFGHIKLMNPQRSTVWYZ");
82        strcpy(upper_rest_chars, "X");
83    }
84
85
86    strcpy(low_chars, real_chars); strlwr(low_chars); // build lower chars
87    { // create low_rest_chars and rest_chars
88        strcpy(low_rest_chars, upper_rest_chars);
89        strlwr(low_rest_chars);
90        strcpy(rest_chars, upper_rest_chars);
91        strcat(rest_chars, low_rest_chars);
92    }
93
94    strcpy(all_chars, real_chars);
95    strcat(all_chars, low_chars);
96    strcat(all_chars, rest_chars);
97
98    strcpy(get_maximum_from, real_chars); // get maximum from markerline from these characters
99    strcat(all_chars, ".-");
100
101    num_all_chars = strlen(all_chars);
102
103    // initialize variables
104    freeset(mline, ARB_calloc<typeof(*mline)>(PHDATA::ROOT->get_seq_len()));
105
106    const long startcol = aw_root->awar(AWAR_PHYLO_FILTER_STARTCOL)->read_int();
107    const long stopcol  = aw_root->awar(AWAR_PHYLO_FILTER_STOPCOL)->read_int();
108    const long minhom   = aw_root->awar(AWAR_PHYLO_FILTER_MINHOM)->read_int();
109    const long maxhom   = aw_root->awar(AWAR_PHYLO_FILTER_MAXHOM)->read_int();
110
111    const FilterMode filter_dot   = FilterMode(aw_root->awar(AWAR_PHYLO_FILTER_DOT)  ->read_int()); // '.' in column
112    const FilterMode filter_minus = FilterMode(aw_root->awar(AWAR_PHYLO_FILTER_MINUS)->read_int()); // '-' in column
113    const FilterMode filter_ambig = FilterMode(aw_root->awar(AWAR_PHYLO_FILTER_AMBIG)->read_int()); // 'MNY....' in column
114    const FilterMode filter_lower = FilterMode(aw_root->awar(AWAR_PHYLO_FILTER_LOWER)->read_int()); // 'acgtu' in column
115
116    delete_when_max[0] = '\0';
117
118    long len = stopcol - startcol;
119
120    // chars_counted[column][index] counts the occurrences of single characters per column
121    // index = num_all_chars   -> count chars which act as column stopper ( = forget whole column if char occurs)
122    // index = num_all_chars+1 -> count masked characters ( = don't count)
123
124    ARB_calloc(chars_counted, len);
125    for (i=0; i<len; i++) {
126        ARB_calloc(chars_counted[i], num_all_chars+2);
127        for (j=0; j<num_all_chars+2; j++) chars_counted[i][j] = 0;
128    }
129
130    for (i=0; i<PHDATA::ROOT->get_seq_len(); i++) mline[i]=-1.0; // all columns invalid; LOOP_VECTORIZED
131    for (i=0; i<256; i++) { // LOOP_VECTORIZED
132        mask[i]            = false;
133        reference_table[i] = num_all_chars;         // invalid and synonym characters
134    }
135
136    // set valid characters
137    for (i=0; i<num_all_chars; i++) {
138        mask[(unsigned char)all_chars[i]]            = true;
139        reference_table[(unsigned char)all_chars[i]] = i;
140    }
141
142    // rna or dna sequence: set synonyms
143    if (isNUC) {
144        reference_table[(unsigned char)'U'] = reference_table[(unsigned char)'T']; // T=U
145        reference_table[(unsigned char)'u'] = reference_table[(unsigned char)'t'];
146        reference_table[(unsigned char)'N'] = reference_table[(unsigned char)'X'];
147        reference_table[(unsigned char)'n'] = reference_table[(unsigned char)'x'];
148    }
149
150    // set mappings according to options
151    // be careful the elements of rest and low are mapped to 'X' and 'a'
152    switch (filter_dot) {     // '.' in column
153        case DONT_COUNT:
154            mask[(unsigned char)'.'] = false;
155            break;
156
157        case SKIP_COLUMN_IF_MAX:
158            strcat(delete_when_max, ".");
159            strcat(get_maximum_from, ".");
160            break;
161
162        case SKIP_COLUMN_IF_OCCUR:
163            reference_table[(unsigned char)'.']=num_all_chars;    // map to invalid position
164            break;
165
166        case COUNT_DONT_USE_MAX: // use like another valid base/acid while not maximal
167            // do nothing: don't get maximum of this charcater
168            // but use character ( true in mask )
169            break;
170
171        case TREAT_AS_UPPERCASE:
172        case TREAT_AS_REGULAR:
173            ph_assert(0); break;  // illegal value!
174    }
175
176    switch (filter_minus) {     // '-' in column
177        case DONT_COUNT:
178            mask[(unsigned char)'-'] = false;
179            break;
180
181        case SKIP_COLUMN_IF_MAX:
182            strcat(delete_when_max, "-");
183            strcat(get_maximum_from, "-");
184            break;
185
186        case SKIP_COLUMN_IF_OCCUR:
187            reference_table[(unsigned char)'-']=num_all_chars;
188            break;
189
190        case COUNT_DONT_USE_MAX: // use like another valid base/acid while not maximal
191            // do nothing: don't get maximum of this charcater
192            // but use character ( true in mask )
193            break;
194
195        case TREAT_AS_UPPERCASE:
196        case TREAT_AS_REGULAR:
197            ph_assert(0); break;  // illegal value!
198    }
199    // 'MNY....' in column
200    bool mapRestToX = false;
201    switch (filter_ambig) { // all rest characters counted to 'X' (see below)
202        case DONT_COUNT:
203            for (i=0; rest_chars[i]; i++) mask[(unsigned char)rest_chars[i]] = false;
204            break;
205
206        case SKIP_COLUMN_IF_MAX:
207            strcat(delete_when_max, "X");
208            strcat(get_maximum_from, "X");
209            mapRestToX = true;
210            break;
211
212        case SKIP_COLUMN_IF_OCCUR:
213            reference_table[(unsigned char)'X'] = num_all_chars;
214            mapRestToX = true;
215            break;
216
217        case COUNT_DONT_USE_MAX: // use like another valid base/acid while not maximal
218            // do nothing: don't get maximum of this charcater
219            // but use character ( true in mask )
220            break;
221
222        case TREAT_AS_REGULAR:
223            strcat(get_maximum_from, upper_rest_chars); // add uppercase rest chars to maximas
224            // lowercase rest chars are handled together with normal lowercase chars (see below)
225            break;
226
227        case TREAT_AS_UPPERCASE:
228            ph_assert(0); break;  // illegal value!
229    }
230
231    if (mapRestToX) {
232        // map all rest_chars to 'X'
233        for (i=0; rest_chars[i]; i++) {
234            reference_table[(unsigned char)rest_chars[i]] = reference_table[(unsigned char)'X'];
235        }
236    }
237
238    switch (filter_lower) { // 'acgtu' in column
239        case DONT_COUNT:
240            for (i=0; low_chars[i]; i++) mask[(unsigned char)low_chars[i]] = false;
241            break;
242
243        case SKIP_COLUMN_IF_MAX:
244            // count all low_chars to 'a'
245            for (i=0; low_chars[i]; i++) reference_table[(unsigned char)low_chars[i]] = reference_table[(unsigned char)'a'];
246            strcat(delete_when_max, "a");
247            strcat(get_maximum_from, "a");
248            break;
249
250        case SKIP_COLUMN_IF_OCCUR:
251            for (i=0; low_chars[i]; i++) reference_table[(unsigned char)low_chars[i]] = num_all_chars;
252            break;
253
254        case COUNT_DONT_USE_MAX:  // use like another valid base/acid while not maximal
255            // do nothing: don't get maximum of this charcater
256            // but use character ( true in mask )
257            break;
258
259        case TREAT_AS_UPPERCASE: // use like corresponding uppercase characters
260            for (i=0; low_chars[i]; i++)     reference_table[(unsigned char)low_chars[i]]      = reference_table[toupper(low_chars[i])];
261            for (i=0; low_rest_chars[i]; i++) reference_table[(unsigned char)low_rest_chars[i]] = reference_table[toupper(low_rest_chars[i])];
262            break;
263
264        case TREAT_AS_REGULAR:
265            ph_assert(0); break;  // illegal value!
266    }
267
268    GB_ERROR error = NULp;
269    if (PHDATA::ROOT->nentries) {
270        arb_progress progress("Counting", long(PHDATA::ROOT->nentries));
271        // counting routine
272        for (i=0; i<long(PHDATA::ROOT->nentries) && !error; i++) {
273            sequence_buffer = (unsigned char*)GB_read_char_pntr(PHDATA::ROOT->hash_elements[i]->gb_species_data_ptr);
274            long send = stopcol;
275            long slen = GB_read_string_count(PHDATA::ROOT->hash_elements[i]->gb_species_data_ptr);
276            if (slen< send) send = slen;
277            for (j=startcol; j<send; j++) {
278                if (mask[sequence_buffer[j]]) {
279                    chars_counted[j-startcol][reference_table[sequence_buffer[j]]]++;
280                }
281                else {
282                    chars_counted[j-startcol][num_all_chars+1]++;
283                }
284            }
285            progress.inc_and_check_user_abort(error);
286        }
287    }
288    if (!error) {
289        // calculate similarity
290        arb_progress progress("Calculate similarity", len);
291        for (i=0; i<len && !error; i++) {
292            if (chars_counted[i][num_all_chars]==0) { // else: forget whole column
293                max=0; max_char=' ';
294                for (j=0; get_maximum_from[j]!='\0'; j++) {
295                    if (max<chars_counted[i][reference_table[(unsigned char)get_maximum_from[j]]]) {
296                        max_char = get_maximum_from[j];
297                        max      = chars_counted[i][reference_table[(unsigned char)max_char]];
298                    }
299                }
300                if ((max!=0) && !strchr(delete_when_max, max_char)) {
301                    // delete SKIP_COLUMN_IF_MAX classes for counting
302                    for (j=0; delete_when_max[j]!='\0'; j++) {
303                        chars_counted[i][num_all_chars+1] += chars_counted[i][reference_table[(unsigned char)delete_when_max[j]]];
304                        chars_counted[i][reference_table[(unsigned char)delete_when_max[j]]]=0;
305                    }
306                    mline[i+startcol] = (max/
307                                         ((float) PHDATA::ROOT->nentries -
308                                          (float) chars_counted[i][num_all_chars+1]))*100.0;
309                    // (maximum in column / number of counted positions) * 100
310                }
311            }
312            progress.inc_and_check_user_abort(error);
313        }
314    }
315
316    for (i=0; i<len; i++) {
317        free(chars_counted[i]);
318    }
319    free(chars_counted);
320
321    if (!error) {
322        char *filt  = ARB_calloc<char>(PHDATA::ROOT->get_seq_len()+1);
323        for (i=0; i<PHDATA::ROOT->get_seq_len(); i++) {
324            filt[i] = minhom<=mline[i] && maxhom>=mline[i]  ? '1' : '0';
325        }
326        filt[i] = '\0';
327        aw_root->awar(AWAR_PHYLO_FILTER_FILTER)->write_string(filt);
328        free(filt);
329
330        return mline;
331    }
332    else {
333        free(mline);
334        return NULp;
335    }
336}
337
338static void update_on_config_change_cb(AW_root *aw_root) {
339    if (aw_root->awar(AWAR_PHYLO_FILTER_AUTOCALC)->read_int()) {
340        ph_calc_filter_cb();
341    }
342    display_status_cb();
343    expose_cb();
344}
345
346static void correct_startstop_cb(AW_root *aw_root, bool start_changed) {
347    AW_awar *awar_startcol = aw_root->awar(AWAR_PHYLO_FILTER_STARTCOL);
348    AW_awar *awar_stopcol  = aw_root->awar(AWAR_PHYLO_FILTER_STOPCOL);
349
350    int startcol = awar_startcol->read_int();
351    int stopcol  = awar_stopcol->read_int();
352
353    if (startcol>stopcol) {
354        if (start_changed) awar_stopcol ->write_int(startcol);
355        else               awar_startcol->write_int(stopcol);
356    }
357}
358static void correct_minmaxhom_cb(AW_root *aw_root, bool min_changed) {
359    AW_awar *awar_minhom = aw_root->awar(AWAR_PHYLO_FILTER_MINHOM);
360    AW_awar *awar_maxhom = aw_root->awar(AWAR_PHYLO_FILTER_MAXHOM);
361
362    int minhom = awar_minhom->read_int();
363    int maxhom = awar_maxhom->read_int();
364
365    if (minhom>maxhom) {
366        if (min_changed) awar_maxhom->write_int(minhom);
367        else             awar_minhom->write_int(maxhom);
368    }
369}
370
371void PH_create_filter_variables(AW_root *aw_root, AW_default default_file, GBDATA *gb_main) {
372    // filter awars
373    long alilength;
374    {
375        GB_transaction  ta(gb_main);
376        char *aliname = GBT_get_default_alignment(gb_main);
377        if (!aliname) {
378            aw_message_if(GB_await_error());
379            alilength = 2; // otherwise set_minmax below aborts
380        }
381        else {
382            alilength = GBT_get_alignment_len(gb_main, aliname);
383            ph_assert(alilength>0);
384            free(aliname);
385        }
386    }
387
388    RootCallback update_on_config_change = makeRootCallback(update_on_config_change_cb);
389
390    aw_root->awar_int(AWAR_PHYLO_FILTER_STARTCOL, 0,           default_file)->set_minmax(0, alilength-1)->add_callback(update_on_config_change)->add_callback(makeRootCallback(correct_startstop_cb, true));
391    aw_root->awar_int(AWAR_PHYLO_FILTER_STOPCOL,  alilength-1, default_file)->set_minmax(0, alilength-1)->add_callback(update_on_config_change)->add_callback(makeRootCallback(correct_startstop_cb, false));
392    aw_root->awar_int(AWAR_PHYLO_FILTER_MINHOM,   0,           default_file)->set_minmax(0, 100)        ->add_callback(update_on_config_change)->add_callback(makeRootCallback(correct_minmaxhom_cb, true));
393    aw_root->awar_int(AWAR_PHYLO_FILTER_MAXHOM,   100,         default_file)->set_minmax(0, 100)        ->add_callback(update_on_config_change)->add_callback(makeRootCallback(correct_minmaxhom_cb, false));
394
395    aw_root->awar_int(AWAR_PHYLO_FILTER_DOT,   DONT_COUNT, default_file)->add_callback(update_on_config_change); // '.' in column
396    aw_root->awar_int(AWAR_PHYLO_FILTER_MINUS, DONT_COUNT, default_file)->add_callback(update_on_config_change); // '-' in column
397    aw_root->awar_int(AWAR_PHYLO_FILTER_AMBIG, DONT_COUNT, default_file)->add_callback(update_on_config_change); // 'MNY....' in column
398    aw_root->awar_int(AWAR_PHYLO_FILTER_LOWER, DONT_COUNT, default_file)->add_callback(update_on_config_change); // 'acgtu' in column
399
400    aw_root->awar_int(AWAR_PHYLO_FILTER_AUTOCALC, 0, default_file)->add_callback(update_on_config_change); // auto-recalculate?
401}
402
403static AWT_config_mapping_def phyl_filter_config_mapping[] = {
404    { AWAR_PHYLO_FILTER_STARTCOL, "startcol" },
405    { AWAR_PHYLO_FILTER_STOPCOL,  "stopcol" },
406    { AWAR_PHYLO_FILTER_MINHOM,   "minhom" },
407    { AWAR_PHYLO_FILTER_MAXHOM,   "maxhom" },
408    { AWAR_PHYLO_FILTER_DOT,      "filtdot" },
409    { AWAR_PHYLO_FILTER_MINUS,    "filtminus" },
410    { AWAR_PHYLO_FILTER_AMBIG,    "filtambig" },
411    { AWAR_PHYLO_FILTER_LOWER,    "filtlower" },
412
413    { NULp, NULp }
414};
415
416AW_window *PH_create_filter_window(AW_root *aw_root) {
417    AW_window_simple *aws = new AW_window_simple;
418    aws->init(aw_root, "PHYL_FILTER", "PHYL FILTER");
419    aws->load_xfig("phylo/filter.fig");
420    aws->button_length(10);
421
422    aws->at("close");
423    aws->callback(AW_POPDOWN);
424    aws->create_button("CLOSE", "CLOSE", "C");
425
426    const int SCALERWIDTH = 200;
427    aws->auto_space(5, 5);
428
429    aws->at("startcol");
430    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_STARTCOL, 6, SCALERWIDTH);
431
432    aws->at("stopcol");
433    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_STOPCOL, 6, SCALERWIDTH);
434
435    aws->at("minhom");
436    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_MINHOM, 3, SCALERWIDTH);
437
438    aws->at("maxhom");
439    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_MAXHOM, 3, SCALERWIDTH);
440
441    aws->label_length(20);
442
443    aws->at("point_opts");
444    aws->label("'.'");
445    aws->create_option_menu(AWAR_PHYLO_FILTER_DOT);
446    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
447    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
448    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
449    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
450    aws->update_option_menu();
451
452    aws->at("minus_opts");
453    aws->label("'-'");
454    aws->create_option_menu(AWAR_PHYLO_FILTER_MINUS);
455    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
456    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
457    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
458    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
459    aws->update_option_menu();
460
461    aws->at("rest_opts");
462    aws->label("ambiguity codes");
463    aws->create_option_menu(AWAR_PHYLO_FILTER_AMBIG);
464    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
465    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
466    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
467    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
468    aws->insert_option(filter_text[TREAT_AS_REGULAR],     "0", TREAT_AS_REGULAR);
469    aws->update_option_menu();
470
471    aws->at("lower_opts");
472    aws->label("lowercase chars");
473    aws->create_option_menu(AWAR_PHYLO_FILTER_LOWER);
474    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
475    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
476    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
477    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
478    aws->insert_option(filter_text[TREAT_AS_UPPERCASE],   "0", TREAT_AS_UPPERCASE);
479    aws->update_option_menu();
480
481    aws->at("config");
482    AWT_insert_config_manager(aws, AW_ROOT_DEFAULT, "phylfilter", phyl_filter_config_mapping);
483
484    aws->at("auto");
485    aws->label("Auto recalculate?");
486    aws->create_toggle(AWAR_PHYLO_FILTER_AUTOCALC);
487
488    return aws;
489}
490
Note: See TracBrowser for help on using the repository browser.