source: tags/svn.1.5.4/NTREE/AP_consensus.cxx

Last change on this file was 8313, checked in by westram, 14 years ago
  • removed dead code

(partly reverted by [8316])

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.9 KB
Line 
1/* -----------------------------------------------------------------
2 * Project:                       ARB
3 *
4 * Module:                        consensus [abbrev.:CON]
5 *
6 * Exported Classes:              x
7 *
8 * Global Functions:              x
9 *
10 * Global Variables:              AWARS:
11 *                  AW_STRING, "tmp/con/alignment"     : name of alignment
12 *                  AW_STRING, "tmp/con/which_species" : only marked species ?
13 *                  AW_STRING, "con/group"         : allow Sgrouping ?
14 *                  AW_STRING, "con/countgaps"
15 *                  AW_DOUBLE, "con/fupper"
16 *                  AW_INT, "con/lower"
17 *                  AW_INT, "con/gapbound"  : result is '-' if more than
18 *                                            'gapbound' per cent gaps
19 *                  AW_DOUBLE, "con/fconsidbound" : if frequency of character is
20 *                                     more than 'fconsidbound' per cent, this
21 *                                     character can contribute to groups
22 *                  AW_STRING, "tmp/con/name"   : save with this name
23 *                  AW_STRING, "con/result" : result has one or more lines ?
24 *
25 * Global Structures:             x
26 *
27 * Private Classes:               .
28 *
29 * Private Functions:             .
30 *
31 * Private Variables:             .
32 *
33 * Dependencies:       Needs consens.fig and CON_groups.fig
34 *
35 * Description: This module calculates the consensus of sequences of
36 *              bases or amino acids. The result can be one or more lines
37 *              of characters and it is written to the extended data of
38 *              the alignment.
39 *
40 * Integration Notes: To use this module the main function must have a
41 *                    callback to the function
42 *                    AW_window *AP_open_consensus_window( AW_root *aw_root)
43 *                    and the function void AP_create_consensus_var
44 *                    (AW_root *aw_root, AW_default aw_def) has to be called.
45 *
46 * -----------------------------------------------------------------
47 */
48
49#include <awt_sel_boxes.hxx>
50#include <aw_window.hxx>
51#include <arb_progress.h>
52#include <aw_root.hxx>
53#include <aw_msg.hxx>
54#include <aw_awar.hxx>
55#include <arbdbt.h>
56#include <arb_strbuf.h>
57
58// AISC_MKPT_PROMOTE:#ifndef AW_BASE_HXX
59// AISC_MKPT_PROMOTE:#include <aw_base.hxx>
60// AISC_MKPT_PROMOTE:#endif
61
62#define AWAR_MAX_FREQ   "tmp/CON_MAX_FREQ/"
63
64#define AWAR_MAX_FREQ_NO_GAPS AWAR_MAX_FREQ "no_gaps"
65#define AWAR_MAX_FREQ_SAI_NAME AWAR_MAX_FREQ "sai_name"
66
67
68extern GBDATA *GLOBAL_gb_main;
69enum {
70    BAS_GAP,
71    BAS_A,  BAS_C,  BAS_G,  BAS_T, BAS_N,
72
73    MAX_BASES,
74    MAX_AMINOS  = 27,
75    MAX_GROUPS  = 40
76};
77
78/* -----------------------------------------------------------------
79 * Function:                     CON_evaluatestatistic
80 *
81 * Arguments:                    int **statistic,char **groupflags,
82 *                               char *groupnames
83 * Returns:                      char *result
84 *
85 * Description:       This function creates one or more result strings out
86 *                    of table statistic. Memory for result is allocated
87 *                    and later freed in function CON_calculate_cb
88 *
89 * NOTE:              Usage of groupflags and groupnames see function
90 *                    CON_makegrouptable.
91 *
92 * Global Variables referenced:  .
93 *
94 * Global Variables modified:    x
95 *
96 * AWARs referenced:             .
97 *
98 * AWARs modified:               x
99 *
100 * Dependencies:                 Always check that behavior is identical to that
101 *                               of ED4_char_table::build_consensus_string()
102 * -----------------------------------------------------------------
103 */
104static void CON_evaluatestatistic(char   *&result, int **statistic, char **groupflags,
105                           char    *groupnames, int alignlength, double fupper, int lower,
106                           double   fconsidbound, int gapbound, int countgap, int numgroups)
107{
108    int row=0;
109    int j = 0;
110    int groupfr[MAX_GROUPS]; // frequency of group
111    int highestfr, highestgr;
112    long numentries;
113
114    arb_progress progress("calculating result", alignlength);
115
116    result=(char *)GB_calloc(alignlength+1, 1);
117
118    for (int column=0; column<alignlength; column++) {
119        numentries=0;
120        for (row=0; statistic[row]; row++) numentries+=statistic[row][column];
121        if (numentries==0) {
122            result[column]='.';
123        }
124        else {
125            if (numentries-statistic[0][column]==0) {
126                result[column]='='; // 100 per cent `-` -> `=`
127            }
128            else {
129                if (!countgap) {
130                    numentries -= statistic[0][column];
131                    statistic[0][column]=0;
132                }
133
134                if ((statistic[0][column]*100/numentries)>gapbound) {
135                    result[column]='-';
136                }
137                else {
138                    for (j=0; j<numgroups; j++) {
139                        groupfr[j]=0;
140                    }
141
142                    row=0;
143                    while (statistic[row]) {
144                        if ((statistic[row][column]*100.0)>=fconsidbound*numentries) {
145                            for (j=numgroups-1; j>=0; j--) {
146                                if (groupflags[j][row]) {
147                                    groupfr[j] += statistic[row][column];
148                                }
149                            }
150                        }
151                        row++;
152                    }
153
154                    highestfr=0;
155                    highestgr=0;
156                    for (j=0; j<numgroups; j++) {
157                        if (groupfr[j] > highestfr) {
158                            highestfr=groupfr[j];
159                            highestgr=j;
160                        }
161                    }
162
163                    if ((highestfr*100.0/numentries)>=fupper) {
164                        result[column]=groupnames[highestgr];
165                    }
166                    else if ((highestfr*100/numentries)>=lower) {
167                        char c=groupnames[highestgr];
168                        if (c>='A' && c<='Z') {
169                            c=c-'A'+'a';
170                        }
171                        result[column]=c;
172                    }
173                    else {
174                        result[column]='.';
175                    }
176                }
177            }
178        }
179        ++progress;
180    }
181}
182
183
184
185/* -----------------------------------------------------------------
186 * Function:                     CON_makegrouptable
187 *
188 * Arguments:                    .
189 *
190 * Returns:                      char **gf [groupflags],char *groupnames
191 *
192 * Description:       Creates table groupflags that is used in function
193 *                    CON_evaluatestatistic. E.g. gf[10][3]=1 means, that
194 *                    all characters c with convtable[c]==3 are members
195 *                    of group 10.
196 *                    Table groupnames is also created.
197 *                    E.g. c=groupnames[5] gives abbrev of 5th group.
198 *
199 * NOTE:                         .
200 *
201 * Global Variables referenced:  .
202 *
203 * Global Variables modified:    x
204 *
205 * AWARs referenced:             .
206 *
207 * AWARs modified:               x
208 *
209 * Dependencies:                 .
210 * -----------------------------------------------------------------
211 */
212static int CON_makegrouptable(char **gf, char *groupnames,
213                       int isamino, int groupallowed)
214{
215    for (int j=0; j<MAX_GROUPS; j++) {
216        gf[j]=(char *)GB_calloc(MAX_GROUPS, 1); }
217
218    if (!isamino)
219    {
220        strcpy(groupnames, "-ACGUMRWSYKVHDBN\0");
221        for (int i=1; i<MAX_BASES; i++) {
222            gf[i][i]=1; }
223        if (groupallowed)
224        {
225            gf[5][BAS_A]=1; gf[5][BAS_C]=1;
226            gf[6][BAS_A]=1; gf[6][BAS_G]=1;
227            gf[7][BAS_A]=1; gf[7][BAS_T]=1;
228            gf[8][BAS_C]=1; gf[8][BAS_G]=1;
229            gf[9][BAS_C]=1; gf[9][BAS_T]=1;
230            gf[10][BAS_G]=1; gf[10][BAS_T]=1;
231            gf[11][BAS_A]=1; gf[11][BAS_C]=1; gf[11][BAS_G]=1;
232            gf[12][BAS_A]=1; gf[12][BAS_C]=1; gf[12][BAS_T]=1;
233            gf[13][BAS_A]=1; gf[13][BAS_G]=1; gf[13][BAS_T]=1;
234            gf[14][BAS_C]=1; gf[14][BAS_G]=1; gf[14][BAS_T]=1;
235            gf[15][BAS_A]=1; gf[15][BAS_C]=1;
236            gf[15][BAS_G]=1; gf[15][BAS_T]=1;
237            return (16);
238        }
239        return (5);
240    }
241    else
242    {
243        strcpy(groupnames, "-ABCDEFGHI*KLMN.PQRST.VWXYZADHIF\0");
244        for (int i=1; i<MAX_AMINOS; i++) {
245            gf[i][i]=1; }
246        if (groupallowed)
247#define SC(x, P) gf[x][P-'A'+1] = 1
248        {
249            SC(27, 'P'); SC(27, 'A'); SC(27, 'G'); SC(27, 'S'); SC(27, 'T');
250            // PAGST
251            SC(28, 'Q'); SC(28, 'N'); SC(28, 'E'); SC(28, 'D'); SC(28, 'B');
252            SC(28, 'Z');   // QNEDBZ
253            SC(29, 'H'); SC(29, 'K'); SC(29, 'R');
254            // HKR
255            SC(30, 'L'); SC(30, 'I'); SC(30, 'V'); SC(30, 'M');
256            // LIVM
257            SC(31, 'F'); SC(31, 'Y'); SC(31, 'W');
258            // FYW
259            return (32);
260        }
261#undef SC
262        return (27);
263    }
264}
265
266
267/* -----------------------------------------------------------------
268 * Function:                     CON_makestatistic
269 *
270 * Arguments:                    int *convtable
271 *
272 * Returns:                      int **statistic
273 *
274 * Description:       This function fills the table statistic, that is used
275 *                    later by function CON_evaluatestatistic. The filling
276 *                    is done depending on table convtable, that is created
277 *                    in function CON_maketables. Memory for statistic is
278 *                    allocated also in function CON_maketables.
279 *
280 *
281 * NOTE:                         .
282 *
283 * Global Variables referenced:  .
284 *
285 * Global Variables modified:    x
286 *
287 * AWARs referenced:             .
288 *
289 * AWARs modified:               x
290 *
291 * Dependencies:                 .
292 * -----------------------------------------------------------------
293 */
294
295
296static long CON_makestatistic(int **statistic, int *convtable, char *align, int onlymarked)
297{
298    long maxalignlen=GBT_get_alignment_len(GLOBAL_gb_main, align);
299    GBDATA *gb_species, *alidata;
300    int i, nrofspecies=0;
301
302    if (onlymarked) {
303        nrofspecies = GBT_count_marked_species(GLOBAL_gb_main);
304    }
305    else {
306        nrofspecies = GBT_get_species_count(GLOBAL_gb_main);
307    }
308
309    arb_progress progress(nrofspecies);
310    progress.auto_subtitles("Examining sequence");
311
312    if (onlymarked) {
313        gb_species = GBT_first_marked_species(GLOBAL_gb_main);
314    }
315    else {
316        gb_species = GBT_first_species(GLOBAL_gb_main);
317    }
318
319    while (gb_species) {
320        if ((alidata=GBT_read_sequence(gb_species, align))) {
321            unsigned char        c;
322            const unsigned char *data = (const unsigned char *)GB_read_char_pntr(alidata);
323
324            i = 0;
325            while ((c=data[i])) {
326                if ((c=='-') || ((c>='a')&&(c<='z')) || ((c>='A')&&(c<='Z'))
327                    || (c=='*')) {
328                    if (i > maxalignlen) break;
329                    statistic[convtable[c]][i] += 1;
330                }
331                i++;
332            }
333        }
334        if (onlymarked) {
335            gb_species = GBT_next_marked_species(gb_species);
336        }
337        else {
338            gb_species = GBT_next_species(gb_species);
339        }
340        ++progress;
341    }
342    return (nrofspecies);
343}
344
345
346/* -----------------------------------------------------------------
347 * Function:          CON_maketables
348 *
349 * Arguments:         long maxalignlen,int isamino
350 *
351 * Returns:           return parameters: int *convtable,int **statistic
352 *
353 * Description:       Creates tables convtable and statistic, that are
354 *                    used by function CON_makestatistic. The memory
355 *                    allocated for both tables is freed in the
356 *                    function CON_calculate_cb.
357 *                    E.g. convtable['c']=k means that the character c
358 *                    is counted in table statistic in row k.
359 *
360 * NOTE:                         .
361 *
362 * Global Variables referenced:  .
363 *
364 * Global Variables modified:    x
365 *
366 * AWARs referenced:             .
367 *
368 * AWARs modified:               x
369 *
370 * Dependencies:                 .
371 * -----------------------------------------------------------------
372 */
373static void CON_maketables(int *convtable, int **statistic, long maxalignlen, int isamino)
374{
375    int i;
376    for (i=0; i<256; i++) { convtable[i]=0; }
377    if (!isamino) {
378        for (i='a'; i <= 'z'; i++) convtable[i] = BAS_N;
379        for (i='A'; i <= 'Z'; i++) convtable[i] = BAS_N;
380
381        convtable['a']=BAS_A;   convtable['A']=BAS_A;
382        convtable['c']=BAS_C;   convtable['C']=BAS_C;
383        convtable['g']=BAS_G;   convtable['G']=BAS_G;
384        convtable['t']=BAS_T;   convtable['T']=BAS_T;
385        convtable['u']=BAS_T;   convtable['U']=BAS_T;
386
387
388        for (i=0; i<MAX_BASES; i++) {
389            statistic[i]=(int*)GB_calloc((unsigned int)maxalignlen, sizeof(int));
390        }
391        statistic[MAX_BASES]=NULL;
392    }
393    else {
394        for (i=0; i<MAX_AMINOS-1; i++)
395        {
396            convtable['a'+i]=i+1;
397            convtable['A'+i]=i+1;
398        }
399        for (i=0; i<MAX_AMINOS; i++) {
400            statistic[i]=(int*)GB_calloc((unsigned int)maxalignlen, sizeof(int));
401        }
402        statistic[MAX_AMINOS]=NULL;
403        convtable['*']=10; // 'J'
404    }
405}
406
407// export results into database
408static GB_ERROR CON_export(char *savename, char *align, int **statistic, char *result, int *convtable, char *groupnames, int onlymarked, long nrofspecies, long maxalignlen, int countgaps, int gapbound, int groupallowed, double fconsidbound, double fupper, int lower, int resultiscomplex)
409{
410    GB_ERROR    err;
411    const char *off    = "off";
412    const char *on     = "on";
413    char       *buffer = (char *)GB_calloc(2000, sizeof(char));
414
415    GBDATA *gb_extended = GBT_find_or_create_SAI(GLOBAL_gb_main, savename);
416    GBDATA *gb_data     = GBT_add_data(gb_extended, align, "data", GB_STRING);
417    err                 = GB_write_string(gb_data, result);
418    GBDATA *gb_options  = GBT_add_data(gb_extended, align, "_TYPE", GB_STRING);
419
420    const char *allvsmarked    = "all";
421    if (onlymarked) allvsmarked = "marked";
422
423    const char *countgapsstring   = off;
424    if (countgaps) countgapsstring = on;
425
426    const char *simplifystring      = off;
427    if (groupallowed) simplifystring = on;
428
429    sprintf(buffer, "CON: [species: %s]  [number: %ld]  [count gaps: %s] "
430            "[threshold for gaps: %d]  [simplify: %s] "
431            "[threshold for character: %f]  [upper: %f]  [lower: %d]",
432            allvsmarked, nrofspecies, countgapsstring,
433            gapbound, simplifystring,
434            fconsidbound, fupper, lower);
435
436    err=GB_write_string(gb_options, buffer);
437
438    GBDATA *gb_names = GB_search(GB_get_father(gb_options), "_SPECIES", GB_FIND);
439    if (gb_names) GB_delete(gb_names); // delete old entry
440
441    if (nrofspecies<20) {
442        GBDATA        *gb_species;
443        GBS_strstruct *strstruct = GBS_stropen(1000);
444
445        if (onlymarked) gb_species = GBT_first_marked_species(GLOBAL_gb_main);
446        else gb_species            = GBT_first_species(GLOBAL_gb_main);
447
448        while (gb_species) {
449            if (GBT_read_sequence(gb_species, align)) {
450                GBDATA     *gb_speciesname = GB_search(gb_species, "name", GB_FIND);
451                const char *name           = GB_read_char_pntr(gb_speciesname);
452
453                GBS_strcat(strstruct, name);
454                GBS_chrcat(strstruct, ' ');
455            }
456            if (onlymarked) gb_species = GBT_next_marked_species(gb_species);
457            else gb_species            = GBT_next_species(gb_species);
458        }
459
460        char *allnames = GBS_strclose(strstruct);
461        err            = GBT_write_string(GB_get_father(gb_options), "_SPECIES", allnames);
462        free(allnames);
463    }
464
465    {
466        char buffer2[256];
467        sprintf(buffer2, "%s/FREQUENCIES", align);
468        GBDATA *gb_graph = GB_search(gb_extended, buffer2, GB_FIND);
469        if (gb_graph) GB_delete(gb_graph);  // delete old entry
470    }
471    // export additional information
472    if (resultiscomplex)
473    {
474        GBDATA *gb_graph = GBT_add_data(gb_extended, align, "FREQUENCIES", GB_DB);
475        char *charname=(char *)GB_calloc(5, sizeof(char));
476
477        float **additional=0;
478        // problem : aminos, especially '*' -> new order
479
480        int *allreadycounted=(int*)GB_calloc((unsigned int)256, sizeof(char));
481        int *neworder=(int*)GB_calloc((unsigned int)256, sizeof(int));
482        int k;
483        int numdiffchars=1;  // first additional row (nr. 0) is max-row
484        for (int c=0; c<256; c++)
485        {
486            if ((k=convtable[c]))
487            {
488                if (!(allreadycounted[k]))
489                {
490                    allreadycounted[k]=1;
491                    neworder[numdiffchars++]=k;
492                }
493            }
494        }
495
496        additional=(float**)GB_calloc((unsigned int)numdiffchars, sizeof(float*));
497        int group;
498        for (group=0; group<numdiffchars; group++)
499        {
500            additional[group]=(float*)GB_calloc((unsigned int)maxalignlen,
501                                                sizeof(float));
502        }
503
504        int *absolutrow=(int*)GB_calloc((unsigned int)maxalignlen, sizeof(int));
505        long col;
506        for (col=0; col<maxalignlen; col++)
507        {
508            int group2=1;
509            int colsum=0;
510            while (neworder[group2])
511            {
512                colsum+=statistic[neworder[group2++]][col];
513            }
514            if (countgaps) colsum+=statistic[0][col];
515            absolutrow[col]=colsum;
516        }
517
518        for (col=0; col<maxalignlen; col++)
519        {
520            int group2=1;
521            float highest=0, relative;
522            int diffchar;
523            if (absolutrow[col]) {
524                while ((diffchar=neworder[group2++])) {
525                    relative=(float)statistic[diffchar][col]
526                        /(float)absolutrow[col];
527                    if (relative>highest) highest=relative;
528                    additional[diffchar][col]=relative;
529                }
530                additional[0][col]=highest;
531            }
532            else {
533                additional[0][col]=0.0;
534            }
535        }
536
537        GBDATA *gb_relative=GB_search(gb_graph, "MAX", GB_FLOATS);
538        err=GB_write_floats(gb_relative, additional[0], maxalignlen);
539
540        for (group=1; group<numdiffchars; group++)
541        {
542            char ch = groupnames[neworder[group]];
543            if (ch <'A' || ch>'Z') continue;
544            sprintf(charname, "N%c", ch);
545            gb_relative=GB_search(gb_graph, charname, GB_FLOATS);
546            err=GB_write_floats(gb_relative, additional[group], maxalignlen);
547        }
548
549        free(charname);
550        free(neworder);
551        free(allreadycounted);
552        for (group=0; group<numdiffchars; group++) free(additional[group]);
553        free(additional);
554    }
555    free(buffer);
556    return err;
557}
558
559
560static void CON_cleartables(int **statistic, int isamino) {
561    int i;
562    int no_of_tables = isamino ? MAX_AMINOS : MAX_BASES;
563
564    for (i = 0; i<no_of_tables; ++i) {
565        free(statistic[i]);
566    }
567}
568
569/* -----------------------------------------------------------------
570 * Function:                     void CON_calculate_cb(AW_window *aw )
571 *
572 * Arguments:                    .
573 *
574 * Returns:                      .
575 *
576 * Description:                  main callback
577 *                               This function calculates the consensus.
578 *                               Function CON_makestatistic creates the
579 *                               statistic and function CON_evaluatestatistic
580 *                               evaluates the statistic. Then the result
581 *                               string(s) are written to the extended data of
582 *                               the alignment.
583 * NOTE:                         .
584 *
585 * Global Variables referenced:  .
586 *
587 * Global Variables modified:    x
588 *
589 * AWARs referenced:
590 *                  AW_STRING, "tmp/con/alignment"     : name of alignment
591 *                  AW_STRING, "tmp/con/which_species" : only marked species ?
592 *                  AW_STRING, "con/group"         : allow grouping ?
593 *                  AW_STRING, "con/countgaps"
594 *                  AW_DOUBLE, "con/fupper"
595 *                  AW_INT, "con/lower"
596 *                  AW_INT, "con/gapbound"  : result is '-' if more than
597 *                                            'gapbound' per cent gaps
598 *                  AW_DOUBLE, "con/fconsidbound" : if frequency of character is
599 *                                     more than 'considbound' per cent, this
600 *                                     character can contribute to groups
601 *                  AW_STRING, "tmp/intcon/name"   : save with this name
602 *                  AW_STRING, "con/result" : result has one or more lines ?
603 *
604 * AWARs modified:               x
605 *
606 * Dependencies:                 CON_maketables
607 *                               CON_makestatistic
608 *                               CON_makegrouptable
609 *                               CON_evaluatestatistic
610 * -----------------------------------------------------------------
611 */
612static void CON_calculate_cb(AW_window *aw)
613{
614    AW_root  *awr   = aw->get_root();
615    char     *align = awr->awar("tmp/con/alignment")->read_string();
616    GB_ERROR  error = 0;
617
618    GB_push_transaction(GLOBAL_gb_main);
619
620    long maxalignlen = GBT_get_alignment_len(GLOBAL_gb_main, align);
621    if (maxalignlen <= 0) error = GB_export_errorf("alignment '%s' doesn't exist", align);
622
623    if (!error) {
624        int isamino         = GBT_is_alignment_protein(GLOBAL_gb_main, align);
625        int onlymarked      = 1;
626        int resultiscomplex = 1;
627
628        {
629            char *marked        = awr->awar("tmp/con/which_species")->read_string();
630            char *complexresult = awr->awar("con/result")->read_string();
631
632            if (strcmp("marked", marked)        != 0) onlymarked = 0;
633            if (strcmp("complex", complexresult) != 0) resultiscomplex = 0;
634
635            free(complexresult);
636            free(marked);
637        }
638
639
640        // creating the table for characters and allocating memory for 'statistic'
641        int *statistic[MAX_AMINOS+1];
642        int  convtable[256];
643        CON_maketables(convtable, statistic, maxalignlen, isamino);
644
645        // filling the statistic table
646        arb_progress progress("Calculating consensus");
647
648        long   nrofspecies = CON_makestatistic(statistic, convtable, align, onlymarked);
649        double fupper      = awr->awar("con/fupper")->read_float();
650        int    lower       = (int)awr->awar("con/lower")->read_int();
651
652        if (fupper>100.0)   fupper = 100;
653        if (fupper<0)       fupper = 0;
654        if (lower<0)        lower  = 0;
655
656        if (lower>fupper) {
657            error = "fault: lower greater than upper";
658        }
659        else {
660
661            double fconsidbound                = awr->awar("con/fconsidbound")->read_float();
662            if (fconsidbound>100) fconsidbound = 100;
663            if (fconsidbound<0)   fconsidbound = 0;
664
665            int gapbound               = (int)awr->awar("con/gapbound")->read_int();
666            if (gapbound<0)   gapbound = 0;
667            if (gapbound>100) gapbound = 100;
668
669            int groupallowed, countgap;
670            {
671                char *group     = awr->awar("con/group")->read_string();
672                char *countgaps = awr->awar("con/countgaps")->read_string();
673
674                groupallowed = strcmp("on", group) == 0;
675                countgap     = strcmp("on", countgaps) == 0;
676
677                free(countgaps);
678                free(group);
679            }
680
681            // creating the table for groups
682            char *groupflags[40];
683            char  groupnames[40];
684            int   numgroups = CON_makegrouptable(groupflags, groupnames, isamino, groupallowed);
685
686            // calculate and export the result strings
687            char *result = 0;
688            CON_evaluatestatistic(result, statistic, groupflags, groupnames, (int)maxalignlen, fupper, lower, fconsidbound, gapbound, countgap, numgroups);
689
690            char *savename = awr->awar("tmp/con/name")->read_string();
691
692            error = CON_export(savename, align, statistic, result, convtable, groupnames,
693                               onlymarked, nrofspecies, maxalignlen, countgap, gapbound, groupallowed,
694                               fconsidbound, fupper, lower, resultiscomplex);
695
696            // freeing allocated memory
697            free(savename);
698            free(result);
699            for (int i=0; i<MAX_GROUPS; i++) free(groupflags[i]);
700        }
701
702        CON_cleartables(statistic, isamino);
703    }
704
705    free(align);
706    GB_end_transaction_show_error(GLOBAL_gb_main, error, aw_message);
707}
708
709void AP_create_consensus_var(AW_root *aw_root, AW_default aw_def)
710{
711    GB_transaction dummy(GLOBAL_gb_main);
712    {
713        char *defali = GBT_get_default_alignment(GLOBAL_gb_main);
714        aw_root->awar_string("tmp/con/alignment", defali, aw_def);
715        free(defali);
716    }
717    aw_root->awar_string("tmp/con/which_species", "marked", aw_def);
718    aw_root->awar_string("con/group", "off", aw_def);
719    aw_root->awar_string("con/countgaps", "on", aw_def);
720    aw_root->awar_float("con/fupper", 95, aw_def);
721    aw_root->awar_int("con/lower", 70, aw_def);
722    aw_root->awar_int("con/gapbound", 60, aw_def);
723    aw_root->awar_float("con/fconsidbound", 30, aw_def);
724    aw_root->awar_string("tmp/con/name", "CONSENSUS", aw_def);
725    aw_root->awar_string("con/result", "single line", aw_def);
726
727    aw_root->awar_string(AWAR_MAX_FREQ_SAI_NAME, "MAX_FREQUENCY", aw_def);
728    aw_root->awar_int(AWAR_MAX_FREQ_NO_GAPS, 1, aw_def);
729
730}
731
732// Open window to show IUPAC tables
733static AW_window * CON_showgroupswin_cb(AW_root *aw_root)
734{
735    AW_window_simple *aws = new AW_window_simple;
736    aws->init(aw_root, "SHOW_IUPAC", "Show IUPAC");
737    aws->load_xfig("consensus/groups.fig");
738    aws->button_length(7);
739
740    aws->at("ok"); aws->callback((AW_CB0)AW_POPDOWN);
741    aws->create_button("CLOSE", "CLOSE", "O");
742
743    return (AW_window *)aws;
744}
745
746AW_window *
747AP_open_con_expert_window(AW_root *aw_root)
748{
749    AW_window_simple *aws = new AW_window_simple;
750    aws->init(aw_root, "BUILD_CONSENSUS", "CONSENSUS OF SEQUENCES");
751    aws->load_xfig("consensus/expert.fig");
752    aws->button_length(6);
753
754    aws->at("cancel"); aws->callback((AW_CB0)AW_POPDOWN);
755    aws->create_button("CLOSE", "CLOSE", "C");
756
757    aws->at("help"); aws->callback(AW_POPUP_HELP, (AW_CL)"consensus.hlp");
758    aws->create_button("HELP", "HELP", "H");
759
760    aws->button_length(10);
761    aws->at("showgroups"); aws->callback(AW_POPUP, (AW_CL)CON_showgroupswin_cb, 0);
762    aws->create_button("SHOW_IUPAC", "show\nIUPAC...", "s");
763    aws->button_length(10);
764
765    aws->at("which_species");
766    aws->create_toggle_field("tmp/con/which_species", NULL, "");
767    aws->insert_toggle("all", "1", "all");
768    aws->insert_default_toggle("marked",   "1", "marked");
769    aws->update_toggle_field();
770
771    aws->at("which_alignment");
772    awt_create_selection_list_on_alignments(GLOBAL_gb_main, (AW_window *)aws, "tmp/con/alignment", "*=");
773
774    aws->button_length(15);
775
776    // activation of consensus calculation by button ...
777    aws->at("calculate"); aws->callback((AW_CB0)CON_calculate_cb);
778    aws->create_button("GO", "GO", "G");
779
780    aws->at("group");
781    aws->create_toggle_field("con/group", NULL, "");
782    aws->insert_toggle("on", "1", "on");
783    aws->insert_default_toggle("off", "1", "off");
784    aws->update_toggle_field();
785
786    aws->at("countgaps");
787    aws->create_toggle_field("con/countgaps", NULL, "");
788    aws->insert_toggle("on", "1", "on");
789    aws->insert_default_toggle("off", "1", "off");
790    aws->update_toggle_field();
791
792    aws->at("upper");
793    aws->create_input_field("con/fupper", 4);
794
795    aws->at("lower");
796    aws->create_input_field("con/lower", 4);
797
798    aws->at("considbound");
799    aws->create_input_field("con/fconsidbound", 4);
800
801    aws->at("gapbound");
802    aws->create_input_field("con/gapbound", 4);
803
804    aws->at("name");
805    aws->create_input_field("tmp/con/name", 10);
806
807    aws->at("result");
808    aws->create_toggle_field("con/result", NULL, "");
809    aws->insert_toggle("single line", "1", "single line");
810    aws->insert_default_toggle("complex", "1", "complex");
811    aws->update_toggle_field();
812
813    aws->at("save_box");
814    awt_create_selection_list_on_extendeds(GLOBAL_gb_main, aws, "tmp/con/name");
815
816    return aws;
817}
818
819/* -----------------------------------------------------------------
820 * Function:           CON_calc_max_freq_cb( AW_window *aw)
821 *
822 * Description:       Gets the maximum frequency for each columns.
823 *
824 * NOTE:
825 *
826 * -----------------------------------------------------------------
827 */
828
829static void CON_calc_max_freq_cb(AW_window *aw) {
830
831    AW_root *awr=aw->get_root();
832    long maxalignlen, i;
833    int *statistic[MAX_AMINOS+1];
834    int convtable[256];
835    int onlymarked=1, isamino=1;
836    long no_gaps;
837
838    GB_ERROR  error       = 0;
839    GB_push_transaction(GLOBAL_gb_main);
840
841    char *align = GBT_get_default_alignment(GLOBAL_gb_main);
842    maxalignlen = GBT_get_alignment_len(GLOBAL_gb_main, align);
843    no_gaps     = awr->awar(AWAR_MAX_FREQ_NO_GAPS)->read_int();
844
845    if (maxalignlen<=0) {
846        GB_pop_transaction(GLOBAL_gb_main);
847        aw_message("alignment doesn't exist!");
848        delete align;
849        return;
850    }
851    isamino = GBT_is_alignment_protein(GLOBAL_gb_main, align);
852
853    arb_progress progress("Calculating max. frequency");
854    long nrofspecies;
855
856    CON_maketables(convtable, statistic, maxalignlen, isamino);
857    nrofspecies = CON_makestatistic(statistic, convtable, align, onlymarked);
858
859    int   end            = MAX_BASES;
860    if (isamino) end     = MAX_AMINOS;
861    int   pos;
862    char *result         = new char[maxalignlen+1];
863    char *result2        = new char[maxalignlen+1];
864    result[maxalignlen]  = 0;
865    result2[maxalignlen] = 0;
866
867    for (pos = 0;  pos < maxalignlen; pos++) {
868        int sum                               = 0;
869        int max                               = 0;
870        for (i=0; i<end; i++) {
871            if (no_gaps && (i == convtable[(unsigned char)'-'])) continue;
872            sum                              += statistic[i][pos];
873            if (statistic[i][pos] > max) max  = statistic[i][pos];
874        }
875        if (sum == 0) {
876            result[pos] = '=';
877            result2[pos] = '=';
878        }
879        else {
880            result[pos] = "01234567890"[((10*max)/sum)%11];
881            result2[pos] = "01234567890"[((100*max)/sum)%10];
882        }
883    }
884
885    char   *savename    = awr->awar(AWAR_MAX_FREQ_SAI_NAME)->read_string();
886    GBDATA *gb_extended = GBT_find_or_create_SAI(GLOBAL_gb_main, savename);
887    free(savename);
888    GBDATA *gb_data     = GBT_add_data(gb_extended, align, "data", GB_STRING);
889    GBDATA *gb_data2    = GBT_add_data(gb_extended, align, "dat2", GB_STRING);
890    error               = GB_write_string(gb_data, result);
891    if (!error) error   = GB_write_string(gb_data2, result2);
892    GBDATA *gb_options  = GBT_add_data(gb_extended, align, "_TYPE", GB_STRING);
893    delete [] result;
894    delete [] result2;
895
896    if (!error) {
897        char    buffer[2000];
898        sprintf (buffer, "MFQ: [species: %li] [exclude gaps: %li]",
899                 nrofspecies, no_gaps);
900        error=GB_write_string(gb_options, buffer);
901    }
902
903    GB_pop_transaction(GLOBAL_gb_main);
904
905    CON_cleartables(statistic, isamino);
906    free(align);
907    if (error) aw_message(error);
908}
909
910AW_window *
911AP_open_max_freq_window(AW_root *aw_root)
912{
913    AW_window_simple *aws = new AW_window_simple;
914    aws->init(aw_root, "MAX_FREQUENCY", "MAX FREQUENCY");
915    aws->load_xfig("consensus/max_freq.fig");
916
917    GB_push_transaction(GLOBAL_gb_main);
918
919    aws->button_length(6);
920
921    aws->at("cancel"); aws->callback((AW_CB0)AW_POPDOWN);
922    aws->create_button("CLOSE", "CLOSE", "C");
923
924    aws->at("help"); aws->callback(AW_POPUP_HELP, (AW_CL)"max_freq.hlp");
925    aws->create_button("HELP", "HELP", "H");
926
927    // activation of consensus calculation by button ...
928    aws->at("go"); aws->callback((AW_CB0)CON_calc_max_freq_cb);
929    aws->create_button("GO", "GO", "C");
930
931    aws->at("save");
932    aws->create_input_field(AWAR_MAX_FREQ_SAI_NAME, 1);
933
934    aws->at("sai");
935    awt_create_selection_list_on_extendeds(GLOBAL_gb_main, aws, AWAR_MAX_FREQ_SAI_NAME);
936
937    aws->at("gaps");
938    aws->create_toggle(AWAR_MAX_FREQ_NO_GAPS);
939
940    GB_pop_transaction(GLOBAL_gb_main);
941
942    return (AW_window *)aws;
943}
944
Note: See TracBrowser for help on using the repository browser.