source: tags/initial/NTREE/AP_pos_var_pars.cxx

Last change on this file was 2, checked in by oldcode, 23 years ago

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