source: tags/ms_r16q2/PHYLO/PH_filt.cxx

Last change on this file was 14622, checked in by westram, 6 years ago
  • reintegrates 'ui' into 'trunk'
    • implements #656
    • adds instant update for
      • next neighbour search
      • branch analysis (mark functions)
      • arb-phylo (filter setup)
  • adds: log:branches/ui@14588:14621
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.9 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 <cctype>
19#include <awt_config_manager.hxx>
20
21static long PH_timer() {
22    static long time = 0;
23    return ++time;
24}
25
26PH_filter::PH_filter()
27{
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 0;
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 = 0;
62
63    if (!PHDATA::ROOT) return 0;                    // 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    free(mline);
105    mline = (float *) calloc((int) PHDATA::ROOT->get_seq_len(), sizeof(float));
106
107
108    const long startcol = aw_root->awar(AWAR_PHYLO_FILTER_STARTCOL)->read_int();
109    const long stopcol  = aw_root->awar(AWAR_PHYLO_FILTER_STOPCOL)->read_int();
110    const long minhom   = aw_root->awar(AWAR_PHYLO_FILTER_MINHOM)->read_int();
111    const long maxhom   = aw_root->awar(AWAR_PHYLO_FILTER_MAXHOM)->read_int();
112
113    const FilterMode filter_dot   = FilterMode(aw_root->awar(AWAR_PHYLO_FILTER_DOT)  ->read_int()); // '.' in column
114    const FilterMode filter_minus = FilterMode(aw_root->awar(AWAR_PHYLO_FILTER_MINUS)->read_int()); // '-' in column
115    const FilterMode filter_ambig = FilterMode(aw_root->awar(AWAR_PHYLO_FILTER_AMBIG)->read_int()); // 'MNY....' in column
116    const FilterMode filter_lower = FilterMode(aw_root->awar(AWAR_PHYLO_FILTER_LOWER)->read_int()); // 'acgtu' in column
117
118    delete_when_max[0] = '\0';
119
120    long len = stopcol - startcol;
121
122    // chars_counted[column][index] counts the occurrences of single characters per column
123    // index = num_all_chars   -> count chars which act as column stopper ( = forget whole column if char occurs)
124    // index = num_all_chars+1 -> count masked characters ( = don't count)
125
126    chars_counted=(long **) calloc((int) (len), sizeof(long *));
127    for (i=0; i<len; i++) {
128        chars_counted[i]=(long *) calloc((int) num_all_chars+2, sizeof(long));
129        for (j=0; j<num_all_chars+2; j++) chars_counted[i][j]=0;
130    }
131
132    for (i=0; i<PHDATA::ROOT->get_seq_len(); i++) mline[i]=-1.0; // all columns invalid
133    for (i=0; i<256; i++) {
134        mask[i]            = false;
135        reference_table[i] = num_all_chars;         // invalid and synonym characters
136    }
137
138    // set valid characters
139    for (i=0; i<num_all_chars; i++) {
140        mask[(unsigned char)all_chars[i]]            = true;
141        reference_table[(unsigned char)all_chars[i]] = i;
142    }
143
144    // rna or dna sequence: set synonyms
145    if (isNUC) {
146        reference_table[(unsigned char)'U'] = reference_table[(unsigned char)'T']; // T=U
147        reference_table[(unsigned char)'u'] = reference_table[(unsigned char)'t'];
148        reference_table[(unsigned char)'N'] = reference_table[(unsigned char)'X'];
149        reference_table[(unsigned char)'n'] = reference_table[(unsigned char)'x'];
150    }
151
152    // set mappings according to options
153    // be careful the elements of rest and low are mapped to 'X' and 'a'
154    switch (filter_dot) {     // '.' in column
155        case DONT_COUNT:
156            mask[(unsigned char)'.'] = false;
157            break;
158
159        case SKIP_COLUMN_IF_MAX:
160            strcat(delete_when_max, ".");
161            strcat(get_maximum_from, ".");
162            break;
163
164        case SKIP_COLUMN_IF_OCCUR:
165            reference_table[(unsigned char)'.']=num_all_chars;    // map to invalid position
166            break;
167
168        case COUNT_DONT_USE_MAX: // use like another valid base/acid while not maximal
169            // do nothing: don't get maximum of this charcater
170            // but use character ( true in mask )
171            break;
172
173        case TREAT_AS_UPPERCASE:
174        case TREAT_AS_REGULAR:
175            ph_assert(0); break;  // illegal value!
176    }
177
178    switch (filter_minus) {     // '-' in column
179        case DONT_COUNT:
180            mask[(unsigned char)'-'] = false;
181            break;
182
183        case SKIP_COLUMN_IF_MAX:
184            strcat(delete_when_max, "-");
185            strcat(get_maximum_from, "-");
186            break;
187
188        case SKIP_COLUMN_IF_OCCUR:
189            reference_table[(unsigned char)'-']=num_all_chars;
190            break;
191
192        case COUNT_DONT_USE_MAX: // use like another valid base/acid while not maximal
193            // do nothing: don't get maximum of this charcater
194            // but use character ( true in mask )
195            break;
196
197        case TREAT_AS_UPPERCASE:
198        case TREAT_AS_REGULAR:
199            ph_assert(0); break;  // illegal value!
200    }
201    // 'MNY....' in column
202    bool mapRestToX = false;
203    switch (filter_ambig) // all rest characters counted to 'X' (see below)
204    {
205        case DONT_COUNT:
206            for (i=0; rest_chars[i]; i++) mask[(unsigned char)rest_chars[i]] = false;
207            break;
208
209        case SKIP_COLUMN_IF_MAX:
210            strcat(delete_when_max, "X");
211            strcat(get_maximum_from, "X");
212            mapRestToX = true;
213            break;
214
215        case SKIP_COLUMN_IF_OCCUR:
216            reference_table[(unsigned char)'X'] = num_all_chars;
217            mapRestToX = true;
218            break;
219
220        case COUNT_DONT_USE_MAX: // use like another valid base/acid while not maximal
221            // do nothing: don't get maximum of this charcater
222            // but use character ( true in mask )
223            break;
224
225        case TREAT_AS_REGULAR:
226            strcat(get_maximum_from, upper_rest_chars); // add uppercase rest chars to maximas
227            // lowercase rest chars are handled together with normal lowercase chars (see below)
228            break;
229
230        case TREAT_AS_UPPERCASE:
231            ph_assert(0); break;  // illegal value!
232    }
233
234    if (mapRestToX) {
235        // map all rest_chars to 'X'
236        for (i=0; rest_chars[i]; i++) {
237            reference_table[(unsigned char)rest_chars[i]] = reference_table[(unsigned char)'X'];
238        }
239    }
240
241    switch (filter_lower) { // 'acgtu' in column
242        case DONT_COUNT:
243            for (i=0; low_chars[i]; i++) mask[(unsigned char)low_chars[i]] = false;
244            break;
245
246        case SKIP_COLUMN_IF_MAX:
247            // count all low_chars to 'a'
248            for (i=0; low_chars[i]; i++) reference_table[(unsigned char)low_chars[i]] = reference_table[(unsigned char)'a'];
249            strcat(delete_when_max, "a");
250            strcat(get_maximum_from, "a");
251            break;
252
253        case SKIP_COLUMN_IF_OCCUR:
254            for (i=0; low_chars[i]; i++) reference_table[(unsigned char)low_chars[i]] = num_all_chars;
255            break;
256
257        case COUNT_DONT_USE_MAX:  // use like another valid base/acid while not maximal
258            // do nothing: don't get maximum of this charcater
259            // but use character ( true in mask )
260            break;
261
262        case TREAT_AS_UPPERCASE: // use like corresponding uppercase characters
263            for (i=0; low_chars[i]; i++)     reference_table[(unsigned char)low_chars[i]]      = reference_table[toupper(low_chars[i])];
264            for (i=0; low_rest_chars[i]; i++) reference_table[(unsigned char)low_rest_chars[i]] = reference_table[toupper(low_rest_chars[i])];
265            break;
266
267        case TREAT_AS_REGULAR:
268            ph_assert(0); break;  // illegal value!
269    }
270
271    GB_ERROR error = NULL;
272    if (PHDATA::ROOT->nentries) {
273        arb_progress progress("Counting", PHDATA::ROOT->nentries);
274        // counting routine
275        for (i=0; i<long(PHDATA::ROOT->nentries) && !error; i++) {
276            sequence_buffer = (unsigned char*)GB_read_char_pntr(PHDATA::ROOT->hash_elements[i]->gb_species_data_ptr);
277            long send = stopcol;
278            long slen = GB_read_string_count(PHDATA::ROOT->hash_elements[i]->gb_species_data_ptr);
279            if (slen< send) send = slen;
280            for (j=startcol; j<send; j++) {
281                if (mask[sequence_buffer[j]]) {
282                    chars_counted[j-startcol][reference_table[sequence_buffer[j]]]++;
283                }
284                else {
285                    chars_counted[j-startcol][num_all_chars+1]++;
286                }
287            }
288            progress.inc_and_check_user_abort(error);
289        }
290    }
291    if (!error) {
292        // calculate similarity
293        arb_progress progress("Calculate similarity", len);
294        for (i=0; i<len && !error; i++) {
295            if (chars_counted[i][num_all_chars]==0) // else: forget whole column
296            {
297                max=0; max_char=' ';
298                for (j=0; get_maximum_from[j]!='\0'; j++) {
299                    if (max<chars_counted[i][reference_table[(unsigned char)get_maximum_from[j]]]) {
300                        max_char = get_maximum_from[j];
301                        max      = chars_counted[i][reference_table[(unsigned char)max_char]];
302                    }
303                }
304                if ((max!=0) && !strchr(delete_when_max, max_char)) {
305                    // delete SKIP_COLUMN_IF_MAX classes for counting
306                    for (j=0; delete_when_max[j]!='\0'; j++) {
307                        chars_counted[i][num_all_chars+1] += chars_counted[i][reference_table[(unsigned char)delete_when_max[j]]];
308                        chars_counted[i][reference_table[(unsigned char)delete_when_max[j]]]=0;
309                    }
310                    mline[i+startcol] = (max/
311                                         ((float) PHDATA::ROOT->nentries -
312                                          (float) chars_counted[i][num_all_chars+1]))*100.0;
313                    // (maximum in column / number of counted positions) * 100
314                }
315            }
316            progress.inc_and_check_user_abort(error);
317        }
318    }
319
320    for (i=0; i<len; i++) {
321        free(chars_counted[i]);
322    }
323    free(chars_counted);
324
325    if (!error) {
326        char *filt  = (char *)calloc((int) PHDATA::ROOT->get_seq_len()+1, sizeof(char));
327        for (i=0; i<PHDATA::ROOT->get_seq_len(); i++) {
328            filt[i] = minhom<=mline[i] && maxhom>=mline[i]  ? '1' : '0';
329        }
330        filt[i] = '\0';
331        aw_root->awar(AWAR_PHYLO_FILTER_FILTER)->write_string(filt);
332        free(filt);
333
334        return mline;
335    }
336    else {
337        free(mline);
338        return 0;
339    }
340}
341
342static void update_on_config_change_cb(AW_root *aw_root) {
343    if (aw_root->awar(AWAR_PHYLO_FILTER_AUTOCALC)->read_int()) {
344        ph_calc_filter_cb();
345    }
346    display_status_cb();
347    expose_cb();
348}
349
350static void correct_startstop_cb(AW_root *aw_root, bool start_changed) {
351    AW_awar *awar_startcol = aw_root->awar(AWAR_PHYLO_FILTER_STARTCOL);
352    AW_awar *awar_stopcol  = aw_root->awar(AWAR_PHYLO_FILTER_STOPCOL);
353
354    int startcol = awar_startcol->read_int();
355    int stopcol  = awar_stopcol->read_int();
356
357    if (startcol>stopcol) {
358        if (start_changed) awar_stopcol ->write_int(startcol);
359        else               awar_startcol->write_int(stopcol);
360    }
361}
362static void correct_minmaxhom_cb(AW_root *aw_root, bool min_changed) {
363    AW_awar *awar_minhom = aw_root->awar(AWAR_PHYLO_FILTER_MINHOM);
364    AW_awar *awar_maxhom = aw_root->awar(AWAR_PHYLO_FILTER_MAXHOM);
365
366    int minhom = awar_minhom->read_int();
367    int maxhom = awar_maxhom->read_int();
368
369    if (minhom>maxhom) {
370        if (min_changed) awar_maxhom->write_int(minhom);
371        else             awar_minhom->write_int(maxhom);
372    }
373}
374
375void PH_create_filter_variables(AW_root *aw_root, AW_default default_file, GBDATA *gb_main) {
376    // filter awars
377    long alilength;
378    {
379        GB_transaction ta(gb_main);
380        char *aliname = GBT_get_default_alignment(gb_main);
381        alilength     = GBT_get_alignment_len(gb_main, aliname);
382        free(aliname);
383    }
384
385    RootCallback update_on_config_change = makeRootCallback(update_on_config_change_cb);
386
387    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));
388    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));
389    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));
390    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));
391
392    aw_root->awar_int(AWAR_PHYLO_FILTER_DOT,   DONT_COUNT, default_file)->add_callback(update_on_config_change); // '.' in column
393    aw_root->awar_int(AWAR_PHYLO_FILTER_MINUS, DONT_COUNT, default_file)->add_callback(update_on_config_change); // '-' in column
394    aw_root->awar_int(AWAR_PHYLO_FILTER_AMBIG, DONT_COUNT, default_file)->add_callback(update_on_config_change); // 'MNY....' in column
395    aw_root->awar_int(AWAR_PHYLO_FILTER_LOWER, DONT_COUNT, default_file)->add_callback(update_on_config_change); // 'acgtu' in column
396
397    aw_root->awar_int(AWAR_PHYLO_FILTER_AUTOCALC, 0, default_file)->add_callback(update_on_config_change); // auto-recalculate?
398}
399
400static AWT_config_mapping_def phyl_filter_config_mapping[] = {
401    { AWAR_PHYLO_FILTER_STARTCOL, "startcol" },
402    { AWAR_PHYLO_FILTER_STOPCOL,  "stopcol" },
403    { AWAR_PHYLO_FILTER_MINHOM,   "minhom" },
404    { AWAR_PHYLO_FILTER_MAXHOM,   "maxhom" },
405    { AWAR_PHYLO_FILTER_DOT,      "filtdot" },
406    { AWAR_PHYLO_FILTER_MINUS,    "filtminus" },
407    { AWAR_PHYLO_FILTER_AMBIG,    "filtambig" },
408    { AWAR_PHYLO_FILTER_LOWER,    "filtlower" },
409
410    { 0, 0 }
411};
412
413AW_window *PH_create_filter_window(AW_root *aw_root) {
414    AW_window_simple *aws = new AW_window_simple;
415    aws->init(aw_root, "PHYL_FILTER", "PHYL FILTER");
416    aws->load_xfig("phylo/filter.fig");
417    aws->button_length(10);
418
419    aws->at("close");
420    aws->callback(AW_POPDOWN);
421    aws->create_button("CLOSE", "CLOSE", "C");
422
423    const int SCALERWIDTH = 200;
424    aws->auto_space(5, 5);
425
426    aws->at("startcol");
427    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_STARTCOL, 6, SCALERWIDTH);
428
429    aws->at("stopcol");
430    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_STOPCOL, 6, SCALERWIDTH);
431
432    aws->at("minhom");
433    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_MINHOM, 3, SCALERWIDTH);
434
435    aws->at("maxhom");
436    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_MAXHOM, 3, SCALERWIDTH);
437
438    aws->label_length(20);
439
440    aws->at("point_opts");
441    aws->label("'.'");
442    aws->create_option_menu(AWAR_PHYLO_FILTER_DOT, true);
443    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
444    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
445    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
446    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
447    aws->update_option_menu();
448
449    aws->at("minus_opts");
450    aws->label("'-'");
451    aws->create_option_menu(AWAR_PHYLO_FILTER_MINUS, true);
452    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
453    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
454    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
455    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
456    aws->update_option_menu();
457
458    aws->at("rest_opts");
459    aws->label("ambiguity codes");
460    aws->create_option_menu(AWAR_PHYLO_FILTER_AMBIG, true);
461    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
462    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
463    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
464    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
465    aws->insert_option(filter_text[TREAT_AS_REGULAR],     "0", TREAT_AS_REGULAR);
466    aws->update_option_menu();
467
468    aws->at("lower_opts");
469    aws->label("lowercase chars");
470    aws->create_option_menu(AWAR_PHYLO_FILTER_LOWER, true);
471    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
472    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
473    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
474    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
475    aws->insert_option(filter_text[TREAT_AS_UPPERCASE],   "0", TREAT_AS_UPPERCASE);
476    aws->update_option_menu();
477
478    aws->at("config");
479    AWT_insert_config_manager(aws, AW_ROOT_DEFAULT, "phylfilter", phyl_filter_config_mapping);
480
481    aws->at("auto");
482    aws->label("Auto recalculate?");
483    aws->create_toggle(AWAR_PHYLO_FILTER_AUTOCALC);
484
485    return aws;
486}
487
Note: See TracBrowser for help on using the repository browser.