source: tags/arb_5.3/NTREE/AP_pos_var_pars.cxx

Last change on this file was 5825, checked in by westram, 15 years ago
  • replaced GB_entry/GB_create/GB_search+GB_read_XXX by GBT_read_XXX
  • improved error handling in many functions
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.2 KB
Line 
1#include <stdlib.h>
2#include <stdio.h>
3#include <ctype.h>
4#include <string.h>
5// #include <malloc.h>
6#include <memory.h>
7#include <math.h>
8#include <arbdb.h>
9#include <arbdbt.h>
10#include <aw_root.hxx>
11#include <aw_device.hxx>
12#include <aw_window.hxx>
13#include <awt.hxx>
14#include <awt_tree.hxx>
15#include <awt_sel_boxes.hxx>
16
17#include "ap_pos_var_pars.hxx"
18
19#define ap_assert(cond) arb_assert(cond)
20
21extern GBDATA *GLOBAL_gb_main;
22
23
24AP_pos_var::AP_pos_var(GBDATA *gb_maini,char *ali_namei, long ali_leni, int isdna, char *tree_namei) {
25    memset((char  *)this, 0, sizeof(AP_pos_var) ) ;
26    this->gb_main   = gb_maini;
27    this->ali_name  = strdup(ali_namei);
28    this->is_dna    = isdna;
29    this->ali_len   = ali_leni;
30    this->tree_name = strdup(tree_namei);
31}
32
33AP_pos_var::~AP_pos_var() {
34    free(ali_name);
35    free(tree_name);
36    free(transitions);
37    free(transversions);
38    int i;
39    for (i=0;i<256;i++) free(frequencies[i]);
40}
41
42long    AP_pos_var::getsize(GBT_TREE *tree){
43    if (!tree) return 0;
44    if (tree->is_leaf) return 1;
45    return getsize(tree->leftson) + getsize(tree->rightson) + 1;
46}
47
48const char *AP_pos_var::parsimony(GBT_TREE *tree, GB_UINT4 *bases, GB_UINT4 *tbases){
49    GB_ERROR error = 0;
50    timer ++;
51    long i;
52    long l,r;
53
54    if (tree->is_leaf) {
55        unsigned char *sequence;
56        long           seq_len = ali_len;
57        if (!tree->gb_node) return 0; // zombie
58
59        GBDATA *gb_data = GBT_read_sequence(tree->gb_node,ali_name);
60        if (!gb_data) return 0; // no sequence
61        if (GB_read_string_count(gb_data) < seq_len)
62            seq_len     = GB_read_string_count(gb_data);
63        sequence        = (unsigned char*)GB_read_char_pntr(gb_data);
64
65        for (i = 0; i< seq_len; i++ ) {
66            l = char_2_freq[sequence[i]];
67            if (l) {
68                ap_assert(frequencies[l]);
69                frequencies[l][i]++;
70            }
71        }
72
73        if (bases){
74            for (i = 0; i< seq_len; i++ ) bases[i] = char_2_transition[sequence[i]];
75        }
76        if (tbases){
77            for (i = 0; i< seq_len; i++ ) tbases[i] = char_2_transversion[sequence[i]];
78        }
79        return 0;
80    }
81    if (aw_status(timer/(double)treesize)) return "Operation aborted";
82
83    GB_UINT4 *ls = (GB_UINT4 *)calloc(sizeof(GB_UINT4),(int)ali_len);
84    GB_UINT4 *rs = (GB_UINT4 *)calloc(sizeof(GB_UINT4),(int)ali_len);
85    GB_UINT4 *lts = (GB_UINT4 *)calloc(sizeof(GB_UINT4),(int)ali_len);
86    GB_UINT4 *rts = (GB_UINT4 *)calloc(sizeof(GB_UINT4),(int)ali_len);
87
88    if (!error) error = this->parsimony(tree->leftson,ls,lts);
89    if (!error) error = this->parsimony(tree->rightson,rs,rts);
90    if (!error){
91        for (i=0; i< ali_len; i++ ) {
92            l = ls[i];
93            r = rs[i];
94            if ( l & r ) {
95                if (bases) bases[i] = l&r;
96            }else{
97                transitions[i] ++;
98                if (bases) bases[i] = l|r;
99            }
100            l = lts[i];
101            r = rts[i];
102            if ( l & r ) {
103                if (tbases) tbases[i] = l&r;
104            }else{
105                transversions[i] ++;
106                if (tbases) tbases[i] = l|r;
107            }
108        }
109    }
110
111    free(lts);
112    free(rts);
113
114    free(ls);
115    free(rs);
116    return error;
117}
118
119
120// Calculate the positional variability: control procedure
121GB_ERROR AP_pos_var::retrieve( GBT_TREE *tree){
122    GB_ERROR error = 0;
123    int i;
124
125    if (is_dna) {
126        long           base;
127        unsigned char *char_2_bitstring;
128        char_2_freq[(unsigned char)'a'] = 'A';
129        char_2_freq[(unsigned char)'A'] = 'A';
130        char_2_freq[(unsigned char)'c'] = 'C';
131        char_2_freq[(unsigned char)'C'] = 'C';
132        char_2_freq[(unsigned char)'g'] = 'G';
133        char_2_freq[(unsigned char)'G'] = 'G';
134        char_2_freq[(unsigned char)'t'] = 'U';
135        char_2_freq[(unsigned char)'T'] = 'U';
136        char_2_freq[(unsigned char)'u'] = 'U';
137        char_2_freq[(unsigned char)'U'] = 'U';
138        char_2_bitstring = (unsigned char *)AP_create_dna_to_ap_bases();
139        for (i=0;i<256;i++) {
140            int j;
141            if (i=='-') j = '.'; else j = i;
142            base = char_2_transition[i] = char_2_bitstring[j];
143            char_2_transversion[i] = 0;
144            if (base & (AP_A | AP_G) ) char_2_transversion[i] = 1;
145            if (base & (AP_C | AP_T) ) char_2_transversion[i] |= 2;
146        }
147        delete [] char_2_bitstring;
148    }
149    else {
150        long            base;
151        AWT_translator *translator       = AWT_get_user_translator(gb_main);
152        const long     *char_2_bitstring = translator->Pro2Bitset();
153 
154        for (i=0;i<256;i++){
155            char_2_transversion[i] = 0;
156            base = char_2_transition[i] = char_2_bitstring[i];
157            if (base) char_2_freq[i] = toupper(i);
158        }
159    }
160    this->treesize = this->getsize(tree);
161    this->timer = 0;
162
163    for (i=0;i<256;i++) {
164        int j;
165        if ( (j = char_2_freq[i]) ) {
166            if (!frequencies[j]) {
167                frequencies[j] = (GB_UINT4 *)calloc(sizeof(GB_UINT4),(int)ali_len);
168            }
169        }
170    }
171
172    transitions = (GB_UINT4 *)calloc(sizeof(GB_UINT4),(int)ali_len);
173    transversions = (GB_UINT4 *)calloc(sizeof(GB_UINT4),(int)ali_len);
174
175    error = this->parsimony(tree);
176
177    return error;
178}
179
180GB_ERROR AP_pos_var::delete_old_sai(const char *sai_name) {
181    GBDATA *gb_extended = GBT_find_SAI(gb_main, sai_name);
182    if (gb_extended) {
183        GBDATA *gb_ali = GB_search(gb_extended, ali_name, GB_FIND);
184        if (gb_ali) {
185            return GB_delete(gb_ali);
186        }
187    }
188    return NULL; // sai/ali did not exist
189}
190
191GB_ERROR AP_pos_var::save_sai(const char *sai_name) {
192    GB_ERROR  error       = 0;
193    GBDATA   *gb_extended = GBT_find_or_create_SAI(gb_main, sai_name);
194
195    if (!gb_extended) error = GB_await_error();
196    else {
197        GBDATA *gb_ali     = GB_search(gb_extended, ali_name, GB_DB);
198        if (!gb_ali) error = GB_await_error();
199        else {
200            const char *description = 
201                GBS_global_string("PVP: Positional Variability by Parsimony: tree '%s' ntaxa %li",
202                                  tree_name, treesize/2);
203           
204            error = GBT_write_string(gb_ali, "_TYPE", description);
205        }
206
207        if (!error) {
208            char *data = (char*)calloc(1,(int)ali_len+1);
209            int  *sum  = (int*)calloc(sizeof(int), (int)ali_len);
210           
211            for (int j=0; j<256 && !error; j++) {                   // get sum of frequencies
212                if (frequencies[j]) {
213                    for (int i=0; i<ali_len; i++) {
214                        sum[i] += frequencies[j][i];
215                    }
216
217                    if (j >= 'A' && j <= 'Z') {
218                        GBDATA *gb_freq     = GB_search(gb_ali, GBS_global_string("FREQUENCIES/N%c", j), GB_INTS);
219                        if (!gb_freq) error = GB_await_error();
220                        else    error       = GB_write_ints(gb_freq, frequencies[j], ali_len);
221                    }
222                }
223            }
224
225            if (!error) {
226                GBDATA *gb_transi     = GB_search(gb_ali, "FREQUENCIES/TRANSITIONS", GB_INTS);
227                if (!gb_transi) error = GB_await_error();
228                else    error         = GB_write_ints(gb_transi,transitions,ali_len);
229            }
230            if (!error) {
231                GBDATA *gb_transv     = GB_search(gb_ali, "FREQUENCIES/TRANSVERSIONS", GB_INTS);
232                if (!gb_transv) error = GB_await_error();
233                error                 = GB_write_ints(gb_transv,transversions,ali_len);
234            }
235
236            if (!error) {
237                int    max_categ = 0;
238                double logbase   = sqrt(2.0);
239                double lnlogbase = log(logbase);
240                double b         = .75;
241                double max_rate  = 1.0;
242
243                for (int i=0; i<ali_len; i++) {
244                    if (sum[i] * 10 <= treesize) {
245                        data[i] = '.';
246                        continue;
247                    }
248                    if (transitions[i] == 0) {
249                        data[i] = '-';
250                        continue;
251                    }
252                    double rate = transitions[i]/ (double)sum[i];
253                    if (rate >= b * .95) {
254                        rate = b * .95;
255                    }
256                    rate = -b * log(1-rate/b);
257                    if (rate > max_rate) rate = max_rate;
258                    rate /= max_rate;       // scaled  1.0 == fast rate
259                    // ~0.0 slow rate
260                    double dcat = -log(rate)/lnlogbase;
261                    int icat = (int)dcat;
262                    if (icat > 35) icat = 35;
263                    if (icat >= max_categ) max_categ = icat +1;
264                    data[i] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"[icat];
265                }
266
267                error = GBT_write_string(gb_ali, "data", data);
268
269                if (!error) {
270                    // Generate Categories
271                    GBS_strstruct *out = GBS_stropen(1000);
272                    for (int i = 0; i<max_categ; i++) {
273                        GBS_floatcat(out, pow(1.0/logbase, i) );
274                        GBS_chrcat(out,' ');
275                    }
276
277                    error = GBT_write_string(gb_ali, "_CATEGORIES", GBS_mempntr(out));
278                    GBS_strforget(out);
279                }
280            }
281
282            free(sum);
283            free(data);
284        }
285    }
286
287    return error;
288}
289
290
291
292// Calculate the positional variability: window interface
293void AP_calc_pos_var_pars(AW_window *aww) {
294    AW_root        *root  = aww->get_root();
295    GB_transaction  dummy(GLOBAL_gb_main);
296    char           *tree_name;
297    GB_ERROR        error = 0;
298
299    aw_openstatus("Calculating positional variability");
300    aw_status("Loading Tree");
301    GBT_TREE *tree;
302    {           // get tree
303        tree_name = root->awar(AWAR_PVP_TREE)->read_string();
304        tree = GBT_read_tree(GLOBAL_gb_main,tree_name,sizeof(GBT_TREE));
305        if (!tree) {
306            error = "Please select a valid tree";
307        }
308        else {
309            GBT_link_tree(tree,GLOBAL_gb_main, GB_TRUE, 0, 0);
310        }
311    }
312
313    if (!error) {
314        aw_status("Counting Mutations");
315
316        char *ali_name = GBT_get_default_alignment(GLOBAL_gb_main);
317        long  ali_len  = GBT_get_alignment_len(GLOBAL_gb_main,ali_name);
318
319        if (ali_len <=0) {
320            error = "Please select a valid alignment";
321        }
322        else {
323            GB_alignment_type  at       = GBT_get_alignment_type(GLOBAL_gb_main, ali_name);
324            int                isdna    = at==GB_AT_DNA || at==GB_AT_RNA;
325            char              *sai_name = root->awar(AWAR_PVP_SAI)->read_string();
326            AP_pos_var         pv(GLOBAL_gb_main, ali_name, ali_len, isdna, tree_name);
327
328            error             = pv.delete_old_sai(sai_name);
329            if (!error) error = pv.retrieve( tree);
330            if (!error) error = pv.save_sai(sai_name);
331
332            free(sai_name);
333        }
334        free(ali_name);
335    }
336
337    if (tree) GBT_delete_tree(tree);
338    free(tree_name);
339
340    aw_closestatus();
341
342    if (error) aw_message(error);
343    return;
344}
345
346AW_window *AP_open_pos_var_pars_window( AW_root *root ){
347
348    GB_transaction dummy(GLOBAL_gb_main);
349
350    AW_window_simple *aws = new AW_window_simple;
351    aws->init( root, "CSP_BY_PARSIMONY", "Conservation Profile: Parsimony Method");
352    aws->load_xfig("cpro/parsimony.fig");
353
354    root->awar_string(AWAR_PVP_SAI, "POS_VAR_BY_PARSIMONY",AW_ROOT_DEFAULT);
355    char *largest_tree = GBT_find_largest_tree(GLOBAL_gb_main);
356    root->awar_string(AWAR_PVP_TREE, "tree_full",AW_ROOT_DEFAULT);
357    root->awar(AWAR_PVP_TREE)->write_string(largest_tree);
358    freeset(largest_tree, 0);
359
360    aws->at("close");aws->callback((AW_CB0)AW_POPDOWN);
361    aws->create_button("CLOSE","CLOSE","C");
362
363    aws->at("help");aws->callback(AW_POPUP_HELP,(AW_CL)"pos_var_pars.hlp");
364    aws->create_button("HELP","HELP","H");
365
366    aws->at("name");
367    aws->create_input_field(AWAR_PVP_SAI);
368
369    aws->at("box");
370    awt_create_selection_list_on_extendeds(GLOBAL_gb_main,aws,AWAR_PVP_SAI);
371
372    aws->at("trees");
373    awt_create_selection_list_on_trees(GLOBAL_gb_main,aws,AWAR_PVP_TREE);
374
375    aws->at("go");
376    aws->highlight();
377    aws->callback(AP_calc_pos_var_pars);
378    aws->create_button("GO","GO");
379
380    return (AW_window *)aws;
381}
Note: See TracBrowser for help on using the repository browser.