1 | #include <stdio.h> |
---|
2 | #include <stdlib.h> |
---|
3 | #include <memory.h> |
---|
4 | #include <arbdb.h> |
---|
5 | #include <arbdbt.h> |
---|
6 | #include <aw_root.hxx> |
---|
7 | #include <aw_device.hxx> |
---|
8 | #include <aw_window.hxx> |
---|
9 | #include "awt.hxx" |
---|
10 | #include "awt_tree.hxx" |
---|
11 | #include "BI_helix.hxx" |
---|
12 | #include "awt_csp.hxx" |
---|
13 | #include "awt_sel_boxes.hxx" |
---|
14 | |
---|
15 | |
---|
16 | |
---|
17 | void awt_csp_rescan_sais(AW_root *awr, AW_CL csp_cd){ |
---|
18 | AWT_csp *csp = (AWT_csp *)csp_cd; |
---|
19 | GB_transaction dummy(csp->gb_main); |
---|
20 | free(csp->alignment_name); |
---|
21 | free(csp->type_path); |
---|
22 | csp->alignment_name = awr->awar(csp->awar_alignment)->read_string(); |
---|
23 | csp->type_path = GBS_string_eval(csp->alignment_name,"*=*1/_TYPE",0); |
---|
24 | if (csp->sai_sel_box_id) { |
---|
25 | awt_create_selection_list_on_extendeds_update(0,csp->sai_sel_box_id); |
---|
26 | } |
---|
27 | } |
---|
28 | |
---|
29 | AWT_csp::AWT_csp(GBDATA *gb_maini, AW_root *awri, const char *awar_template) { |
---|
30 | /* awar_template == ".../name" |
---|
31 | * -> generated ".../alignment" |
---|
32 | * ".../smooth" |
---|
33 | * ".../enable_helix" |
---|
34 | */ |
---|
35 | |
---|
36 | memset((char *)this,0,sizeof(AWT_csp)); |
---|
37 | |
---|
38 | this->gb_main = gb_maini; |
---|
39 | this->awr = awri; |
---|
40 | |
---|
41 | this->awar_name = GBS_string_eval(awar_template,AWT_CSP_AWAR_CSP_NAME,0); |
---|
42 | this->awar_alignment = GBS_string_eval(awar_template,AWT_CSP_AWAR_CSP_ALIGNMENT,0); |
---|
43 | this->awar_smooth = GBS_string_eval(awar_template,AWT_CSP_AWAR_CSP_SMOOTH,0); |
---|
44 | this->awar_enable_helix = GBS_string_eval(awar_template,AWT_CSP_AWAR_CSP_ENABLE_HELIX,0); |
---|
45 | |
---|
46 | awr->awar_string(awar_name, "NONE"); |
---|
47 | awr->awar_string(awar_alignment); |
---|
48 | awr->awar_int(awar_smooth); |
---|
49 | awr->awar_int(awar_enable_helix, 1); |
---|
50 | |
---|
51 | awr->awar(this->awar_alignment)->add_callback( awt_csp_rescan_sais, (AW_CL)this); |
---|
52 | awt_csp_rescan_sais(awr, (AW_CL)this); |
---|
53 | } |
---|
54 | |
---|
55 | void AWT_csp::exit(){ |
---|
56 | delete [] weights; weights = NULL; |
---|
57 | delete [] rates; rates = NULL; |
---|
58 | delete [] ttratio; ttratio = NULL; |
---|
59 | delete [] is_helix;is_helix= NULL; |
---|
60 | delete [] mut_sum; mut_sum = NULL; |
---|
61 | delete [] freq_sum;freq_sum= NULL; |
---|
62 | delete desc; desc = NULL; |
---|
63 | int i; |
---|
64 | for (i=0;i<256;i++){ |
---|
65 | delete [] frequency[i]; |
---|
66 | frequency[i] = NULL; |
---|
67 | } |
---|
68 | } |
---|
69 | AWT_csp::~AWT_csp(void){ |
---|
70 | this->exit(); |
---|
71 | delete awar_name; |
---|
72 | delete awar_alignment; |
---|
73 | delete awar_smooth; |
---|
74 | delete awar_enable_helix; |
---|
75 | } |
---|
76 | |
---|
77 | |
---|
78 | GB_ERROR AWT_csp::go(AP_filter *filter){ |
---|
79 | this->exit(); |
---|
80 | GB_transaction dummy(this->gb_main); |
---|
81 | long alignment_length = GBT_get_alignment_len(this->gb_main, alignment_name); |
---|
82 | GB_ERROR error = 0; |
---|
83 | if (alignment_length <= 1) |
---|
84 | return GB_export_errorf("Unknown Alignment Size: Name '%s'\n" |
---|
85 | " Select a Valid Alignment",alignment_name); |
---|
86 | if (filter && filter->filter_len != alignment_length) |
---|
87 | return GB_export_error( "Incompatible filter_len" ); |
---|
88 | seq_len = 0; |
---|
89 | |
---|
90 | char *sai_name = awr->awar(awar_name)->read_string(); |
---|
91 | GBDATA *gb_sai = GBT_find_SAI(this->gb_main, sai_name); |
---|
92 | if (!gb_sai) |
---|
93 | error= GB_export_error("Please select a valid Column Statistic"); |
---|
94 | GBDATA *gb_ali = 0; |
---|
95 | GBDATA *gb_freqs = 0; |
---|
96 | if (!error) { |
---|
97 | gb_ali = GB_entry(gb_sai,alignment_name); |
---|
98 | if (!gb_ali) error = GB_export_error("Please select a valid Column Statistic"); |
---|
99 | } |
---|
100 | if (!error) { |
---|
101 | gb_freqs = GB_entry(gb_ali,"FREQUENCIES"); |
---|
102 | if (!gb_ali) error = GB_export_error("Please select a valid Column Statistic"); |
---|
103 | } |
---|
104 | if (error) { |
---|
105 | free(sai_name); |
---|
106 | return error; |
---|
107 | } |
---|
108 | if (filter) seq_len = (size_t)filter->real_len; |
---|
109 | |
---|
110 | else seq_len = (size_t)alignment_length; |
---|
111 | |
---|
112 | unsigned long i; |
---|
113 | unsigned long j; |
---|
114 | delete [] weights; weights = new GB_UINT4[seq_len]; |
---|
115 | delete [] rates; rates = new float[seq_len]; |
---|
116 | delete [] ttratio; ttratio = new float[seq_len]; |
---|
117 | delete [] is_helix;is_helix = new unsigned char[seq_len]; |
---|
118 | delete [] mut_sum; mut_sum = new GB_UINT4[seq_len]; |
---|
119 | delete [] freq_sum;freq_sum = new GB_UINT4[seq_len]; |
---|
120 | delete desc; desc = 0; |
---|
121 | for (i=0;i<256;i++) { delete frequency[i];frequency[i]=0;} |
---|
122 | |
---|
123 | long use_helix = awr->awar(awar_enable_helix)->read_int(); |
---|
124 | |
---|
125 | for (j=i=0;i<seq_len;i++) { |
---|
126 | is_helix[i] = 0; |
---|
127 | weights[i] = 1; |
---|
128 | } |
---|
129 | |
---|
130 | if (!error && use_helix) { // calculate weights and helix filter |
---|
131 | BI_helix helix; |
---|
132 | error = helix.init(this->gb_main,alignment_name); |
---|
133 | if (error){ |
---|
134 | aw_message(error); |
---|
135 | error = 0; |
---|
136 | goto no_helix; |
---|
137 | } |
---|
138 | error = 0; |
---|
139 | for (j=i=0;i<(unsigned long)alignment_length;i++) { |
---|
140 | if (!filter || filter->filter_mask[i]) { |
---|
141 | if (helix.pairtype(i) == HELIX_PAIR) { |
---|
142 | is_helix[j] = 1; |
---|
143 | weights[j] = 1; |
---|
144 | } |
---|
145 | else{ |
---|
146 | is_helix[j] = 0; |
---|
147 | weights[j] = 2; |
---|
148 | } |
---|
149 | j++; |
---|
150 | } |
---|
151 | } |
---|
152 | } |
---|
153 | no_helix: |
---|
154 | |
---|
155 | for (i=0;i<seq_len;i++) freq_sum[i] = 0; |
---|
156 | |
---|
157 | GBDATA *gb_freq; |
---|
158 | GB_UINT4 *freqi[256]; |
---|
159 | for (i=0;i<256; i++) freqi[i] = 0; |
---|
160 | int wf; // ********* read the frequency statistic |
---|
161 | for (gb_freq = GB_child(gb_freqs); gb_freq; gb_freq = GB_nextChild(gb_freq)) { |
---|
162 | char *key = GB_read_key(gb_freq); |
---|
163 | if (key[0] == 'N' && key[1] && !key[2]) { |
---|
164 | wf = key[1]; |
---|
165 | freqi[wf] = GB_read_ints(gb_freq); |
---|
166 | } |
---|
167 | free(key); |
---|
168 | } |
---|
169 | |
---|
170 | GB_UINT4 *minmut = 0; |
---|
171 | GBDATA *gb_minmut = GB_entry(gb_freqs,"TRANSITIONS"); |
---|
172 | if (gb_minmut) minmut = GB_read_ints(gb_minmut); |
---|
173 | |
---|
174 | GB_UINT4 *transver = 0; |
---|
175 | GBDATA *gb_transver = GB_entry(gb_freqs,"TRANSVERSIONS"); |
---|
176 | if (gb_transver) transver = GB_read_ints(gb_transver); |
---|
177 | unsigned long max_freq_sum = 0; |
---|
178 | for (wf = 0; wf<256;wf++) { // ********* calculate sum of mutations |
---|
179 | if (!freqi[wf]) continue; // ********* calculate nbases for each col. |
---|
180 | for (j=i=0;i<(unsigned long)alignment_length;i++) { |
---|
181 | if (filter && !filter->filter_mask[i]) continue; |
---|
182 | freq_sum[j] += freqi[wf][i]; |
---|
183 | mut_sum[j] = minmut[i]; |
---|
184 | j++; |
---|
185 | } |
---|
186 | } |
---|
187 | for (i=0;i<seq_len;i++) |
---|
188 | if (freq_sum[i] > max_freq_sum) |
---|
189 | max_freq_sum = freq_sum[i]; |
---|
190 | unsigned long min_freq_sum = DIST_MIN_SEQ(max_freq_sum); |
---|
191 | |
---|
192 | for (wf = 0; wf<256;wf++) { |
---|
193 | if (!freqi[wf]) continue; |
---|
194 | frequency[wf] = new float[seq_len]; |
---|
195 | for (j=i=0;i<(unsigned long)alignment_length;i++) { |
---|
196 | if (filter && !filter->filter_mask[i]) continue; |
---|
197 | if (freq_sum[j] > min_freq_sum) { |
---|
198 | frequency[wf][j] = freqi[wf][i]/(float)freq_sum[j]; |
---|
199 | }else{ |
---|
200 | frequency[wf][j] = 0; |
---|
201 | weights[j] = 0; |
---|
202 | } |
---|
203 | j++; |
---|
204 | } |
---|
205 | } |
---|
206 | |
---|
207 | for (j=i=0;i<(unsigned long)alignment_length;i++) { // ******* calculate rates |
---|
208 | if (filter && !filter->filter_mask[i]) continue; |
---|
209 | if (!weights[j]) { |
---|
210 | rates[j] = 1.0; |
---|
211 | ttratio[j] = 0.5; |
---|
212 | j++; |
---|
213 | continue; |
---|
214 | } |
---|
215 | rates[j] = (mut_sum[j] / (float)freq_sum[j]); |
---|
216 | if (transver[i] > 0) { |
---|
217 | ttratio[j] = (minmut[i] - transver[i]) / (float)transver[i]; |
---|
218 | }else{ |
---|
219 | ttratio[j] = 2.0; |
---|
220 | } |
---|
221 | if (ttratio[j] < 0.05) ttratio[j] = 0.05; |
---|
222 | if (ttratio[j] > 5.0) ttratio[j] = 5.0; |
---|
223 | j++; |
---|
224 | } |
---|
225 | |
---|
226 | // ****** normalize rates |
---|
227 | |
---|
228 | double sum_rates = 0; |
---|
229 | for (i=0;i<seq_len;i++) sum_rates += rates[i]; |
---|
230 | sum_rates /= seq_len; |
---|
231 | for (i=0;i<seq_len;i++) rates[i] /= sum_rates; |
---|
232 | |
---|
233 | free(transver); |
---|
234 | free(minmut); |
---|
235 | for (i=0;i<256;i++) free(freqi[i]); |
---|
236 | free(sai_name); |
---|
237 | |
---|
238 | return error; |
---|
239 | } |
---|
240 | |
---|
241 | void AWT_csp::print(void) { |
---|
242 | unsigned int j; |
---|
243 | int wf; |
---|
244 | double sum_rates[2], sum_tt[2]; |
---|
245 | long count[2]; |
---|
246 | sum_rates[0] = sum_rates[1] = sum_tt[0] = sum_tt[1] = 0; |
---|
247 | count[0] = count[1] = 0; |
---|
248 | if (!seq_len) return; |
---|
249 | for (j=0;j<seq_len;j++){ |
---|
250 | if (!weights[j]) continue; |
---|
251 | if (is_helix[j]) { |
---|
252 | printf("#"); |
---|
253 | }else{ |
---|
254 | printf("."); |
---|
255 | } |
---|
256 | printf( "%i: minmut %5i freqs %5i rates %5f tt %5f ", |
---|
257 | j, mut_sum[j], freq_sum[j], rates[j], ttratio[j]); |
---|
258 | for (wf = 0;wf<256;wf++) { |
---|
259 | if (!frequency[wf]) continue; |
---|
260 | printf("%c:%5f ",wf,frequency[wf][j]); |
---|
261 | } |
---|
262 | sum_rates[is_helix[j]] += rates[j]; |
---|
263 | sum_tt[is_helix[j]] += ttratio[j]; |
---|
264 | count[is_helix[j]]++; |
---|
265 | printf("w: %i\n",weights[j]); |
---|
266 | } |
---|
267 | printf("Helical Rates %5f Non Hel. Rates %5f\n", |
---|
268 | sum_rates[1]/count[1], sum_rates[0]/count[0]); |
---|
269 | printf("Helical TT %5f Non Hel. TT %5f\n", |
---|
270 | sum_tt[1]/count[1], sum_tt[0]/count[0]); |
---|
271 | } |
---|
272 | |
---|
273 | char *awt_csp_sai_filter(GBDATA *gb_extended, AW_CL csp_cd) { |
---|
274 | AWT_csp *csp = (AWT_csp*)csp_cd; |
---|
275 | GBDATA *gb_type = GB_search(gb_extended, csp->type_path, GB_FIND); |
---|
276 | char *result = 0; |
---|
277 | |
---|
278 | if (gb_type) { |
---|
279 | const char *type = GB_read_char_pntr(gb_type); |
---|
280 | |
---|
281 | if (GBS_string_matches(type, "PV?:*", GB_MIND_CASE)) { |
---|
282 | GBS_strstruct *strstruct = GBS_stropen(100); |
---|
283 | |
---|
284 | GBS_strcat(strstruct, GBT_read_name(gb_extended)); |
---|
285 | GBS_strcat(strstruct, ": <"); |
---|
286 | GBS_strcat(strstruct, type); |
---|
287 | GBS_strcat(strstruct, ">"); |
---|
288 | |
---|
289 | result = GBS_strclose(strstruct); |
---|
290 | } |
---|
291 | } |
---|
292 | return result; |
---|
293 | } |
---|
294 | |
---|
295 | |
---|
296 | void create_selection_list_on_csp(AW_window *aws, AWT_csp *csp){ |
---|
297 | GB_transaction dummy(csp->gb_main); |
---|
298 | csp->sai_sel_box_id = awt_create_selection_list_on_extendeds(csp->gb_main, |
---|
299 | aws, csp->awar_name, awt_csp_sai_filter, (AW_CL)csp); |
---|
300 | } |
---|
301 | |
---|
302 | AW_window *create_csp_window(AW_root *aw_root, AWT_csp *csp){ |
---|
303 | GB_transaction dummy(csp->gb_main); |
---|
304 | AW_window_simple *aws = new AW_window_simple; |
---|
305 | aws->init( aw_root, "SELECT_CSP", "Select Column Statistic"); |
---|
306 | aws->load_xfig("awt/col_statistic.fig"); |
---|
307 | |
---|
308 | aws->at("close");aws->callback((AW_CB0)AW_POPDOWN); |
---|
309 | aws->create_button("CLOSE", "CLOSE","C"); |
---|
310 | |
---|
311 | aws->at("help");aws->callback(AW_POPUP_HELP,(AW_CL)"awt_csp.hlp"); |
---|
312 | aws->create_button("HELP", "HELP","H"); |
---|
313 | |
---|
314 | aws->at("box"); |
---|
315 | create_selection_list_on_csp(aws,csp); |
---|
316 | |
---|
317 | aws->at("smooth"); |
---|
318 | aws->create_toggle_field(csp->awar_smooth); |
---|
319 | aws->insert_toggle("Calculate each column (nearly) independently","D",0); |
---|
320 | aws->insert_toggle("Smooth parameter estimates a little","M",1); |
---|
321 | aws->insert_toggle("Smooth parameter estimates across columns","S",2); |
---|
322 | aws->update_toggle_field(); |
---|
323 | |
---|
324 | aws->at("helix"); |
---|
325 | aws->label("Use Helix Information (SAI 'HELIX')"); |
---|
326 | aws->create_toggle(csp->awar_enable_helix); |
---|
327 | return (AW_window *)aws; |
---|
328 | } |
---|