source: tags/ms_r18q1/PHYLO/PH_filt.cxx

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