source: tags/ms_r16q3/PHYLO/PH_filt.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: 18.8 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    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
131    for (i=0; i<256; i++) {
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    {
203        case DONT_COUNT:
204            for (i=0; rest_chars[i]; i++) mask[(unsigned char)rest_chars[i]] = false;
205            break;
206
207        case SKIP_COLUMN_IF_MAX:
208            strcat(delete_when_max, "X");
209            strcat(get_maximum_from, "X");
210            mapRestToX = true;
211            break;
212
213        case SKIP_COLUMN_IF_OCCUR:
214            reference_table[(unsigned char)'X'] = num_all_chars;
215            mapRestToX = true;
216            break;
217
218        case COUNT_DONT_USE_MAX: // use like another valid base/acid while not maximal
219            // do nothing: don't get maximum of this charcater
220            // but use character ( true in mask )
221            break;
222
223        case TREAT_AS_REGULAR:
224            strcat(get_maximum_from, upper_rest_chars); // add uppercase rest chars to maximas
225            // lowercase rest chars are handled together with normal lowercase chars (see below)
226            break;
227
228        case TREAT_AS_UPPERCASE:
229            ph_assert(0); break;  // illegal value!
230    }
231
232    if (mapRestToX) {
233        // map all rest_chars to 'X'
234        for (i=0; rest_chars[i]; i++) {
235            reference_table[(unsigned char)rest_chars[i]] = reference_table[(unsigned char)'X'];
236        }
237    }
238
239    switch (filter_lower) { // 'acgtu' in column
240        case DONT_COUNT:
241            for (i=0; low_chars[i]; i++) mask[(unsigned char)low_chars[i]] = false;
242            break;
243
244        case SKIP_COLUMN_IF_MAX:
245            // count all low_chars to 'a'
246            for (i=0; low_chars[i]; i++) reference_table[(unsigned char)low_chars[i]] = reference_table[(unsigned char)'a'];
247            strcat(delete_when_max, "a");
248            strcat(get_maximum_from, "a");
249            break;
250
251        case SKIP_COLUMN_IF_OCCUR:
252            for (i=0; low_chars[i]; i++) reference_table[(unsigned char)low_chars[i]] = num_all_chars;
253            break;
254
255        case COUNT_DONT_USE_MAX:  // use like another valid base/acid while not maximal
256            // do nothing: don't get maximum of this charcater
257            // but use character ( true in mask )
258            break;
259
260        case TREAT_AS_UPPERCASE: // use like corresponding uppercase characters
261            for (i=0; low_chars[i]; i++)     reference_table[(unsigned char)low_chars[i]]      = reference_table[toupper(low_chars[i])];
262            for (i=0; low_rest_chars[i]; i++) reference_table[(unsigned char)low_rest_chars[i]] = reference_table[toupper(low_rest_chars[i])];
263            break;
264
265        case TREAT_AS_REGULAR:
266            ph_assert(0); break;  // illegal value!
267    }
268
269    GB_ERROR error = NULL;
270    if (PHDATA::ROOT->nentries) {
271        arb_progress progress("Counting", PHDATA::ROOT->nentries);
272        // counting routine
273        for (i=0; i<long(PHDATA::ROOT->nentries) && !error; i++) {
274            sequence_buffer = (unsigned char*)GB_read_char_pntr(PHDATA::ROOT->hash_elements[i]->gb_species_data_ptr);
275            long send = stopcol;
276            long slen = GB_read_string_count(PHDATA::ROOT->hash_elements[i]->gb_species_data_ptr);
277            if (slen< send) send = slen;
278            for (j=startcol; j<send; j++) {
279                if (mask[sequence_buffer[j]]) {
280                    chars_counted[j-startcol][reference_table[sequence_buffer[j]]]++;
281                }
282                else {
283                    chars_counted[j-startcol][num_all_chars+1]++;
284                }
285            }
286            progress.inc_and_check_user_abort(error);
287        }
288    }
289    if (!error) {
290        // calculate similarity
291        arb_progress progress("Calculate similarity", len);
292        for (i=0; i<len && !error; i++) {
293            if (chars_counted[i][num_all_chars]==0) // else: forget whole column
294            {
295                max=0; max_char=' ';
296                for (j=0; get_maximum_from[j]!='\0'; j++) {
297                    if (max<chars_counted[i][reference_table[(unsigned char)get_maximum_from[j]]]) {
298                        max_char = get_maximum_from[j];
299                        max      = chars_counted[i][reference_table[(unsigned char)max_char]];
300                    }
301                }
302                if ((max!=0) && !strchr(delete_when_max, max_char)) {
303                    // delete SKIP_COLUMN_IF_MAX classes for counting
304                    for (j=0; delete_when_max[j]!='\0'; j++) {
305                        chars_counted[i][num_all_chars+1] += chars_counted[i][reference_table[(unsigned char)delete_when_max[j]]];
306                        chars_counted[i][reference_table[(unsigned char)delete_when_max[j]]]=0;
307                    }
308                    mline[i+startcol] = (max/
309                                         ((float) PHDATA::ROOT->nentries -
310                                          (float) chars_counted[i][num_all_chars+1]))*100.0;
311                    // (maximum in column / number of counted positions) * 100
312                }
313            }
314            progress.inc_and_check_user_abort(error);
315        }
316    }
317
318    for (i=0; i<len; i++) {
319        free(chars_counted[i]);
320    }
321    free(chars_counted);
322
323    if (!error) {
324        char *filt  = ARB_calloc<char>(PHDATA::ROOT->get_seq_len()+1);
325        for (i=0; i<PHDATA::ROOT->get_seq_len(); i++) {
326            filt[i] = minhom<=mline[i] && maxhom>=mline[i]  ? '1' : '0';
327        }
328        filt[i] = '\0';
329        aw_root->awar(AWAR_PHYLO_FILTER_FILTER)->write_string(filt);
330        free(filt);
331
332        return mline;
333    }
334    else {
335        free(mline);
336        return 0;
337    }
338}
339
340static void update_on_config_change_cb(AW_root *aw_root) {
341    if (aw_root->awar(AWAR_PHYLO_FILTER_AUTOCALC)->read_int()) {
342        ph_calc_filter_cb();
343    }
344    display_status_cb();
345    expose_cb();
346}
347
348static void correct_startstop_cb(AW_root *aw_root, bool start_changed) {
349    AW_awar *awar_startcol = aw_root->awar(AWAR_PHYLO_FILTER_STARTCOL);
350    AW_awar *awar_stopcol  = aw_root->awar(AWAR_PHYLO_FILTER_STOPCOL);
351
352    int startcol = awar_startcol->read_int();
353    int stopcol  = awar_stopcol->read_int();
354
355    if (startcol>stopcol) {
356        if (start_changed) awar_stopcol ->write_int(startcol);
357        else               awar_startcol->write_int(stopcol);
358    }
359}
360static void correct_minmaxhom_cb(AW_root *aw_root, bool min_changed) {
361    AW_awar *awar_minhom = aw_root->awar(AWAR_PHYLO_FILTER_MINHOM);
362    AW_awar *awar_maxhom = aw_root->awar(AWAR_PHYLO_FILTER_MAXHOM);
363
364    int minhom = awar_minhom->read_int();
365    int maxhom = awar_maxhom->read_int();
366
367    if (minhom>maxhom) {
368        if (min_changed) awar_maxhom->write_int(minhom);
369        else             awar_minhom->write_int(maxhom);
370    }
371}
372
373void PH_create_filter_variables(AW_root *aw_root, AW_default default_file, GBDATA *gb_main) {
374    // filter awars
375    long alilength;
376    {
377        GB_transaction ta(gb_main);
378        char *aliname = GBT_get_default_alignment(gb_main);
379        alilength     = GBT_get_alignment_len(gb_main, aliname);
380        free(aliname);
381    }
382
383    RootCallback update_on_config_change = makeRootCallback(update_on_config_change_cb);
384
385    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));
386    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));
387    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));
388    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));
389
390    aw_root->awar_int(AWAR_PHYLO_FILTER_DOT,   DONT_COUNT, default_file)->add_callback(update_on_config_change); // '.' in column
391    aw_root->awar_int(AWAR_PHYLO_FILTER_MINUS, DONT_COUNT, default_file)->add_callback(update_on_config_change); // '-' in column
392    aw_root->awar_int(AWAR_PHYLO_FILTER_AMBIG, DONT_COUNT, default_file)->add_callback(update_on_config_change); // 'MNY....' in column
393    aw_root->awar_int(AWAR_PHYLO_FILTER_LOWER, DONT_COUNT, default_file)->add_callback(update_on_config_change); // 'acgtu' in column
394
395    aw_root->awar_int(AWAR_PHYLO_FILTER_AUTOCALC, 0, default_file)->add_callback(update_on_config_change); // auto-recalculate?
396}
397
398static AWT_config_mapping_def phyl_filter_config_mapping[] = {
399    { AWAR_PHYLO_FILTER_STARTCOL, "startcol" },
400    { AWAR_PHYLO_FILTER_STOPCOL,  "stopcol" },
401    { AWAR_PHYLO_FILTER_MINHOM,   "minhom" },
402    { AWAR_PHYLO_FILTER_MAXHOM,   "maxhom" },
403    { AWAR_PHYLO_FILTER_DOT,      "filtdot" },
404    { AWAR_PHYLO_FILTER_MINUS,    "filtminus" },
405    { AWAR_PHYLO_FILTER_AMBIG,    "filtambig" },
406    { AWAR_PHYLO_FILTER_LOWER,    "filtlower" },
407
408    { 0, 0 }
409};
410
411AW_window *PH_create_filter_window(AW_root *aw_root) {
412    AW_window_simple *aws = new AW_window_simple;
413    aws->init(aw_root, "PHYL_FILTER", "PHYL FILTER");
414    aws->load_xfig("phylo/filter.fig");
415    aws->button_length(10);
416
417    aws->at("close");
418    aws->callback(AW_POPDOWN);
419    aws->create_button("CLOSE", "CLOSE", "C");
420
421    const int SCALERWIDTH = 200;
422    aws->auto_space(5, 5);
423
424    aws->at("startcol");
425    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_STARTCOL, 6, SCALERWIDTH);
426
427    aws->at("stopcol");
428    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_STOPCOL, 6, SCALERWIDTH);
429
430    aws->at("minhom");
431    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_MINHOM, 3, SCALERWIDTH);
432
433    aws->at("maxhom");
434    aws->create_input_field_with_scaler(AWAR_PHYLO_FILTER_MAXHOM, 3, SCALERWIDTH);
435
436    aws->label_length(20);
437
438    aws->at("point_opts");
439    aws->label("'.'");
440    aws->create_option_menu(AWAR_PHYLO_FILTER_DOT, true);
441    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
442    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
443    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
444    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
445    aws->update_option_menu();
446
447    aws->at("minus_opts");
448    aws->label("'-'");
449    aws->create_option_menu(AWAR_PHYLO_FILTER_MINUS, true);
450    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
451    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
452    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
453    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
454    aws->update_option_menu();
455
456    aws->at("rest_opts");
457    aws->label("ambiguity codes");
458    aws->create_option_menu(AWAR_PHYLO_FILTER_AMBIG, true);
459    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
460    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
461    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
462    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
463    aws->insert_option(filter_text[TREAT_AS_REGULAR],     "0", TREAT_AS_REGULAR);
464    aws->update_option_menu();
465
466    aws->at("lower_opts");
467    aws->label("lowercase chars");
468    aws->create_option_menu(AWAR_PHYLO_FILTER_LOWER, true);
469    aws->insert_option(filter_text[DONT_COUNT],           "0", DONT_COUNT);
470    aws->insert_option(filter_text[SKIP_COLUMN_IF_MAX],   "0", SKIP_COLUMN_IF_MAX);
471    aws->insert_option(filter_text[SKIP_COLUMN_IF_OCCUR], "0", SKIP_COLUMN_IF_OCCUR);
472    aws->insert_option(filter_text[COUNT_DONT_USE_MAX],   "0", COUNT_DONT_USE_MAX);
473    aws->insert_option(filter_text[TREAT_AS_UPPERCASE],   "0", TREAT_AS_UPPERCASE);
474    aws->update_option_menu();
475
476    aws->at("config");
477    AWT_insert_config_manager(aws, AW_ROOT_DEFAULT, "phylfilter", phyl_filter_config_mapping);
478
479    aws->at("auto");
480    aws->label("Auto recalculate?");
481    aws->create_toggle(AWAR_PHYLO_FILTER_AUTOCALC);
482
483    return aws;
484}
485
Note: See TracBrowser for help on using the repository browser.