source: tags/arb-6.0/NTREE/AP_pos_var_pars.cxx

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