Context Navigation

source: tags/ms_r16q4/SL/GUI_ALIVIEW/ColumnStat.cxx

Visit:

Last change on this file was 14474, checked in by westram, 9 years ago
use typed callbacks for SAI-selectionbox-filter
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 14.8 KB

Line
1	// =============================================================== //
2	// //
3	// File : ColumnStat.cxx //
4	// Purpose : //
5	// //
6	// Institute of Microbiology (Technical University Munich) //
7	// http://www.arb-home.de/ //
8	// //
9	// =============================================================== //
10
11	#include "ColumnStat.hxx"
12	#include "awt_sel_boxes.hxx"
13	#include "ga_local.h"
14	#include <AP_filter.hxx>
15	#include <BI_helix.hxx>
16
17	#include <aw_root.hxx>
18	#include <aw_awar.hxx>
19	#include <aw_msg.hxx>
20	#include <aw_select.hxx>
21
22	#include <arbdbt.h>
23	#include <arb_strbuf.h>
24
25	#define AWAR_COLSTAT_NAME "/name=/name"
26	#define AWAR_COLSTAT_ALIGNMENT "/name=/alignment"
27	#define AWAR_COLSTAT_SMOOTH "/name=/smooth"
28	#define AWAR_COLSTAT_ENABLE_HELIX "/name=/enable_helix"
29
30	void ColumnStat::refresh_sai_selection_list() {
31	GB_transaction ta(gb_main);
32
33	freeset(alignment_name, awr->awar(awar_alignment)->read_string());
34	freeset(type_path, GBS_string_eval(alignment_name, "=1/_TYPE", 0));
35
36	if (saisel) saisel->refresh();
37	}
38
39	static void refresh_columnstat_selection(AW_root, ColumnStat column_stat) {
40	column_stat->refresh_sai_selection_list();
41	}
42
43	ColumnStat::ColumnStat(GBDATA gb_maini, AW_root awri, const char awar_template, AW_awar awar_used_alignment) {
44	/* awar_template == ".../name"
45	* -> generated ".../alignment"
46	* ".../smooth"
47	* ".../enable_helix"
48	*/
49
50	memset((char *)this, 0, sizeof(ColumnStat));
51
52	gb_main = gb_maini;
53	awr = awri;
54
55	awar_name = GBS_string_eval(awar_template, AWAR_COLSTAT_NAME, 0);
56	awar_alignment = GBS_string_eval(awar_template, AWAR_COLSTAT_ALIGNMENT, 0);
57	awar_smooth = GBS_string_eval(awar_template, AWAR_COLSTAT_SMOOTH, 0);
58	awar_enable_helix = GBS_string_eval(awar_template, AWAR_COLSTAT_ENABLE_HELIX, 0);
59
60	ga_assert(strcmp(awar_name, awar_alignment) != 0); // awar_template must end with (or contain) "/name"
61
62	awr->awar_string(awar_name, "");
63	awr->awar_string(awar_alignment)->map(awar_used_alignment);
64	awr->awar_int(awar_smooth);
65	awr->awar_int(awar_enable_helix, 1);
66
67	awr->awar(awar_alignment)->add_callback(makeRootCallback(refresh_columnstat_selection, this));
68	refresh_sai_selection_list(); // same as refresh_columnstat_selection(this)
69	}
70
71	void ColumnStat::forget() {
72	delete [] weights; weights = NULL;
73	delete [] rates; rates = NULL;
74	delete [] ttratio; ttratio = NULL;
75	delete [] is_helix; is_helix = NULL;
76	delete [] mut_sum; mut_sum = NULL;
77	delete [] freq_sum; freq_sum = NULL;
78	delete desc; desc = NULL;
79
80	for (int i=0; i<256; i++) {
81	delete [] frequency[i];
82	frequency[i] = NULL;
83	}
84	seq_len = 0; // mark invalid
85	}
86
87	ColumnStat::~ColumnStat() {
88	forget();
89	delete awar_name;
90	delete awar_alignment;
91	delete awar_smooth;
92	delete awar_enable_helix;
93	}
94
95	GB_ERROR ColumnStat::calculate(AP_filter *filter) {
96	forget(); // delete previously calculated stats
97
98	GB_transaction ta(gb_main);
99	GB_ERROR error = filter ? filter->is_invalid() : NULL;
100	size_t alignment_length = 0;
101
102	if (!error) {
103	long alignment_length_l = GBT_get_alignment_len(gb_main, alignment_name);
104	if (alignment_length_l <= 1) {
105	error = GBS_global_string("Unknown size for alignment '%s'", alignment_name);
106	}
107	else {
108	alignment_length = alignment_length_l;
109	if (filter && filter->get_length() != alignment_length) {
110	error = GBS_global_string("Incompatible filter (filter=%zu bp, alignment=%zu bp)",
111	filter->get_length(), alignment_length);
112	}
113	}
114	}
115	seq_len = 0;
116
117	if (!error) {
118	char *sai_name = awr->awar(awar_name)->read_string();
119	GBDATA *gb_sai = GBT_find_SAI(gb_main, sai_name);
120
121	if (!gb_sai) {
122	if (strcmp(sai_name, "") != 0) { // no error if SAI name is empty
123	error = GBS_global_string("Column statistic: no such SAI: '%s'", sai_name);
124	}
125	}
126
127	if (gb_sai) {
128	GBDATA *gb_ali = 0;
129	GBDATA *gb_freqs = 0;
130	if (!error) {
131	gb_ali = GB_entry(gb_sai, alignment_name);
132	if (!gb_ali) error = GBS_global_string("SAI '%s' does not exist in alignment '%s'", sai_name, alignment_name);
133	}
134	if (!error) {
135	gb_freqs = GB_entry(gb_ali, "FREQUENCIES");
136	if (!gb_freqs) error = GBS_global_string("Column statistic '%s' is invalid (has no FREQUENCIES)", sai_name);
137	}
138
139	if (!error) {
140	seq_len = filter ? filter->get_filtered_length() : alignment_length;
141
142	delete [] weights; weights = new GB_UINT4[seq_len];
143	delete [] rates; rates = new float[seq_len];
144	delete [] ttratio; ttratio = new float[seq_len];
145	delete [] is_helix; is_helix = new bool[seq_len];
146	delete [] mut_sum; mut_sum = new GB_UINT4[seq_len];
147	delete [] freq_sum; freq_sum = new GB_UINT4[seq_len];
148
149	delete desc; desc = 0;
150
151	size_t i;
152	size_t j;
153	for (i=0; i<256; i++) { delete frequency[i]; frequency[i]=0; }
154
155	long use_helix = awr->awar(awar_enable_helix)->read_int();
156
157	for (i=0; i<seq_len; i++) {
158	is_helix[i] = false;
159	weights[i] = 1;
160	}
161
162	if (!error && use_helix) { // calculate weights and helix filter
163	BI_helix helix;
164	error = helix.init(this->gb_main, alignment_name);
165	if (error) {
166	aw_message(error);
167	error = 0;
168	goto no_helix;
169	}
170	error = 0;
171	for (j=i=0; i<alignment_length; i++) {
172	if (!filter \|\| filter->use_position(i)) {
173	if (helix.pairtype(i) == HELIX_PAIR) {
174	is_helix[j] = true;
175	weights[j] = 1;
176	}
177	else {
178	is_helix[j] = false;
179	weights[j] = 2;
180	}
181	j++;
182	}
183	}
184	}
185	no_helix :
186
187	for (i=0; i<seq_len; i++) freq_sum[i] = 0;
188
189	GBDATA *gb_freq;
190	GB_UINT4 *freqi[256];
191	for (i=0; i<256; i++) freqi[i] = 0;
192	int wf; // ********* read the frequency statistic
193	for (gb_freq = GB_child(gb_freqs); gb_freq; gb_freq = GB_nextChild(gb_freq)) {
194	char *key = GB_read_key(gb_freq);
195	if (key[0] == 'N' && key[1] && !key[2]) {
196	wf = key[1];
197	freqi[wf] = GB_read_ints(gb_freq);
198	}
199	free(key);
200	}
201
202	GB_UINT4 *minmut = 0;
203	GB_UINT4 *transver = 0;
204	GBDATA *gb_minmut = GB_entry(gb_freqs, "TRANSITIONS");
205	GBDATA *gb_transver = GB_entry(gb_freqs, "TRANSVERSIONS");
206
207	if (!gb_minmut) error = GBS_global_string("Column statistic '%s' is invalid (has no FREQUENCIES/TRANSITIONS)", sai_name);
208	if (!gb_transver) error = GBS_global_string("Column statistic '%s' is invalid (has no FREQUENCIES/TRANSVERSIONS)", sai_name);
209
210	if (gb_minmut) {
211	minmut = GB_read_ints(gb_minmut);
212	if (!minmut) error = GBS_global_string("Column statistic '%s' is invalid (error reading FREQUENCIES/TRANSITIONS: %s)", sai_name, GB_await_error());
213	}
214	if (gb_transver) {
215	transver = GB_read_ints(gb_transver);
216	if (!transver) error = GBS_global_string("Column statistic '%s' is invalid (error reading FREQUENCIES/TRANSVERSIONS: %s)", sai_name, GB_await_error());
217	}
218
219	if (!error) {
220	unsigned long max_freq_sum = 0;
221	for (wf = 0; wf<256; wf++) { // ********* calculate sum of mutations
222	if (!freqi[wf]) continue; // ********* calculate nbases for each col.
223	for (j=i=0; i<alignment_length; i++) {
224	if (filter && !filter->use_position(i)) continue;
225	freq_sum[j] += freqi[wf][i];
226	mut_sum[j] = minmut[i];
227	j++;
228	}
229	}
230	for (i=0; i<seq_len; i++)
231	if (freq_sum[i] > max_freq_sum)
232	max_freq_sum = freq_sum[i];
233	unsigned long min_freq_sum = DIST_MIN_SEQ(max_freq_sum);
234
235	for (wf = 0; wf<256; wf++) {
236	if (!freqi[wf]) continue;
237	frequency[wf] = new float[seq_len];
238	for (j=i=0; i<alignment_length; i++) {
239	if (filter && !filter->use_position(i)) continue;
240	if (freq_sum[j] > min_freq_sum) {
241	frequency[wf][j] = freqi[wf][i]/(float)freq_sum[j];
242	}
243	else {
244	frequency[wf][j] = 0;
245	weights[j] = 0;
246	}
247	j++;
248	}
249	}
250
251	for (j=i=0; i<alignment_length; i++) { // ******* calculate rates
252	if (filter && !filter->use_position(i)) continue;
253	if (!weights[j]) {
254	rates[j] = 1.0;
255	ttratio[j] = 0.5;
256	j++;
257	continue;
258	}
259	rates[j] = (mut_sum[j] / (float)freq_sum[j]);
260	if (transver[i] > 0) {
261	ttratio[j] = (minmut[i] - transver[i]) / (float)transver[i];
262	}
263	else {
264	ttratio[j] = 2.0;
265	}
266	if (ttratio[j] < 0.05) ttratio[j] = 0.05;
267	if (ttratio[j] > 5.0) ttratio[j] = 5.0;
268	j++;
269	}
270
271	// ****** normalize rates
272
273	double sum_rates = 0;
274	for (i=0; i<seq_len; i++) sum_rates += rates[i];
275	sum_rates /= seq_len;
276	for (i=0; i<seq_len; i++) rates[i] /= sum_rates;
277	}
278
279	free(transver);
280	free(minmut);
281	for (i=0; i<256; i++) free(freqi[i]);
282	}
283	}
284	free(sai_name);
285	}
286
287	return error;
288	}
289
290	void ColumnStat::weight_by_inverseRates() const {
291	for (size_t i = 0; i<seq_len; ++i) {
292	if (rates[i]>0.0000001) {
293	weights[i] *= (int)(2.0 / rates[i]);
294	}
295	}
296	}
297
298
299
300	void ColumnStat::print() {
301	if (seq_len) {
302	double sum_rates[2] = { 0.0, 0.0 };
303	double sum_tt[2] = { 0.0, 0.0 };
304	long count[2] = { 0, 0 };
305
306	for (unsigned j=0; j<seq_len; j++) {
307	if (weights[j]) {
308	fputc(".#"[is_helix[j]], stdout);
309	printf("%u: minmut %5i freqs %5i rates %5f tt %5f ",
310	j, mut_sum[j], freq_sum[j], rates[j], ttratio[j]);
311	for (int wf = 0; wf<256; wf++) {
312	if (frequency[wf]) printf("%c:%5f ", wf, frequency[wf][j]);
313	}
314	sum_rates[is_helix[j]] += rates[j];
315	sum_tt[is_helix[j]] += ttratio[j];
316	count[is_helix[j]]++;
317	printf("w: %i\n", weights[j]);
318	}
319	}
320	printf("Helical Rates %5f Non Hel. Rates %5f\n",
321	sum_rates[1]/count[1], sum_rates[0]/count[0]);
322	printf("Helical TT %5f Non Hel. TT %5f\n",
323	sum_tt[1]/count[1], sum_tt[0]/count[0]);
324	}
325	}
326
327	static char filter_columnstat_SAIs(GBDATA gb_extended, ColumnStat *column_stat) {
328	// returns NULL for non-columnstat SAIs
329	char *result = NULL;
330	if (column_stat->has_valid_alignment()) {
331	GBDATA *gb_type = GB_search(gb_extended, column_stat->get_type_path(), GB_FIND);
332
333	if (gb_type) {
334	const char *type = GB_read_char_pntr(gb_type);
335
336	if (GBS_string_matches(type, "PV?:*", GB_MIND_CASE)) {
337	GBS_strstruct *strstruct = GBS_stropen(100);
338
339	GBS_strcat(strstruct, GBT_read_name(gb_extended));
340	GBS_strcat(strstruct, ": <");
341	GBS_strcat(strstruct, type);
342	GBS_strcat(strstruct, ">");
343
344	result = GBS_strclose(strstruct);
345	}
346	}
347	}
348	return result;
349	}
350
351	void ColumnStat::create_sai_selection_list(AW_window *aww) {
352	GB_transaction ta(gb_main);
353	saisel = awt_create_SAI_selection_list(gb_main, aww, awar_name, true, makeSaiSelectionlistFilterCallback(filter_columnstat_SAIs, this));
354	}
355
356	void COLSTAT_create_selection_list(AW_window aws, ColumnStat column_stat) {
357	column_stat->create_sai_selection_list(aws);
358	}
359
360	AW_window COLSTAT_create_selection_window(AW_root aw_root, ColumnStat *column_stat) {
361	GB_transaction ta(column_stat->get_gb_main());
362	AW_window_simple *aws = new AW_window_simple;
363	aws->init(aw_root, "SELECT_CSP", "Select Column Statistic");
364	aws->load_xfig("awt/col_statistic.fig");
365
366	aws->at("close");
367	aws->callback(AW_POPDOWN);
368	aws->create_button("CLOSE", "CLOSE", "C");
369
370	aws->at("help"); aws->callback(makeHelpCallback("awt_csp.hlp"));
371	aws->create_button("HELP", "HELP", "H");
372
373	aws->at("box");
374	COLSTAT_create_selection_list(aws, column_stat);
375
376	aws->at("smooth");
377	aws->create_toggle_field(column_stat->get_awar_smooth());
378	aws->insert_toggle("Calculate each column (nearly) independently", "D", 0);
379	aws->insert_toggle("Smooth parameter estimates a little", "M", 1);
380	aws->insert_toggle("Smooth parameter estimates across columns", "S", 2);
381	aws->update_toggle_field();
382
383	aws->at("helix");
384	aws->label("Use Helix Information (SAI 'HELIX')");
385	aws->create_toggle(column_stat->get_awar_enable_helix());
386	return (AW_window *)aws;
387	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: