source: branches/port5/PHYLO/PH_filt.cxx

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