source: tags/terminal-6.x/PHYLO/PH_filt.cxx

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