Context Navigation

source: tags/arb_5.0/DIST/DI_protdist.cxx

Visit:

Last change on this file was 6100, checked in by westram, 15 years ago
fix warning "format not a string literal and no format arguments" GB_export_error → GB_export_error/GB_export_errorf
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.5 KB

Line
1	#include <stdio.h>
2	#include <stdlib.h>
3	#include <memory.h>
4	// #include <malloc.h>
5	#include <string.h>
6
7	#include <math.h>
8
9	#include <arbdb.h>
10	#include <arbdbt.h>
11
12	#include <aw_root.hxx>
13	#include <aw_device.hxx>
14	#include <aw_window.hxx>
15	#include <aw_preset.hxx>
16	#include <awt.hxx>
17
18	#include <awt_tree.hxx>
19	#include "dist.hxx"
20	#include <awt_csp.hxx>
21
22	#include "di_matr.hxx"
23	#include "di_protdist.hxx"
24
25	#define epsilon 0.000001/* a small number */
26
27	double di_protdist::pameigs[20] = {
28	-0.022091252, -0.019297602, 0.000004760, -0.017477817,
29	-0.016575549, -0.015504543, -0.002112213, -0.002685727,
30	-0.002976402, -0.013440755, -0.012926992, -0.004293227,
31	-0.005356688, -0.011064786, -0.010480731, -0.008760449,
32	-0.007142318, -0.007381851, -0.007806557, -0.008127024
33	};
34
35	double di_protdist::pamprobs[20][20] = {
36	{
37	-0.01522976, -0.00746819, -0.13934468, 0.11755315, -0.00212101,
38	0.01558456, -0.07408235, -0.00322387, 0.01375826, 0.00448826,
39	0.00154174, 0.02013313, -0.00159183, -0.00069275, -0.00399898,
40	0.08414055, -0.01188178, -0.00029870, 0.00220371, 0.00042546
41	},
42	{
43	-0.07765582, -0.00712634, -0.03683209, -0.08065755, -0.00462872,
44	-0.03791039, 0.10642147, -0.00912185, 0.01436308, -0.00133243,
45	0.00166346, 0.00624657, -0.00003363, -0.00128729, -0.00690319,
46	0.17442028, -0.05373336, -0.00078751, -0.00038151, 0.01382718
47	},
48	{
49	-0.08810973, -0.04081786, -0.04066004, -0.04736004, -0.03275406,
50	-0.03761164, -0.05047487, -0.09086213, -0.03269598, -0.03558015,
51	-0.08407966, -0.07970977, -0.01504743, -0.04011920, -0.05182232,
52	-0.07026991, -0.05846931, -0.01016998, -0.03047472, -0.06280511
53	},
54	{
55	0.02513756, -0.00578333, 0.09865453, 0.01322314, -0.00310665,
56	0.05880899, -0.09252443, -0.02986539, -0.03127460, 0.01007539,
57	-0.00360119, -0.01995024, 0.00094940, -0.00145868, -0.01388816,
58	0.11358341, -0.12127513, -0.00054696, -0.00055627, 0.00417284
59	},
60	{
61	0.16517316, -0.00254742, -0.03318745, -0.01984173, 0.00031890,
62	-0.02817810, 0.02661678, -0.01761215, 0.01665112, 0.10513343,
63	-0.00545026, 0.01827470, -0.00207616, -0.00763758, -0.01322808,
64	-0.02202576, -0.07434204, 0.00020593, 0.00119979, -0.10827873
65	},
66	{
67	0.16088826, 0.00056313, -0.02579303, -0.00319655, 0.00037228,
68	-0.03193150, 0.01655305, -0.03028640, 0.01367746, -0.11248153,
69	0.00778371, 0.02675579, 0.00243718, 0.00895470, -0.01729803,
70	-0.02686964, -0.08262584, 0.00011794, -0.00225134, 0.09415650
71	},
72	{
73	-0.01739295, 0.00572017, -0.00712592, -0.01100922, -0.00870113,
74	-0.00663461, -0.01153857, -0.02248432, -0.00382264, -0.00358612,
75	-0.00139345, -0.00971460, -0.00133312, 0.01927783, -0.01053838,
76	-0.00911362, -0.01010908, 0.09417598, 0.01763850, -0.00955454
77	},
78	{
79	0.01728888, 0.01344211, 0.01200836, 0.01857259, -0.17088517,
80	0.01457592, 0.01997839, 0.02844884, 0.00839403, 0.00196862,
81	0.01391984, 0.03270465, 0.00347173, -0.01940984, 0.01233979,
82	0.00542887, 0.01008836, 0.00126491, -0.02863042, 0.00449764
83	},
84	{
85	-0.02881366, -0.02184155, -0.01566086, -0.02593764, -0.04050907,
86	-0.01539603, -0.02576729, -0.05089606, -0.00597430, 0.02181643,
87	0.09835597, -0.04040940, 0.00873512, 0.12139434, -0.02427882,
88	-0.02945238, -0.01566867, -0.01606503, 0.09475319, 0.02238670
89	},
90	{
91	0.04080274, -0.02869626, -0.05191093, -0.08435843, 0.00021141,
92	0.13043842, 0.00871530, 0.00496058, -0.02797641, -0.00636933,
93	0.02243277, 0.03640362, -0.05735517, 0.00196918, -0.02218934,
94	-0.00608972, 0.02872922, 0.00047619, 0.00151285, 0.00883489
95	},
96	{
97	-0.02623824, 0.00331152, 0.03640692, 0.04260231, -0.00038223,
98	-0.07480340, -0.01022492, -0.00426473, 0.01448116, 0.01456847,
99	0.05786680, 0.03368691, -0.10126924, -0.00147454, 0.01275395,
100	0.00017574, -0.01585206, -0.00015767, -0.00231848, 0.02310137
101	},
102	{
103	-0.00846258, -0.01508106, -0.01967505, -0.02772004, 0.01248253,
104	-0.01331243, -0.02569382, -0.04461524, -0.02207075, 0.04663443,
105	0.19347923, -0.02745691, 0.02288515, -0.04883849, -0.01084597,
106	-0.01947187, -0.00081675, 0.00516540, -0.07815919, 0.08035585
107	},
108	{
109	-0.06553111, 0.09756831, 0.00524326, -0.00885098, 0.00756653,
110	0.02783099, -0.00427042, -0.16680359, 0.03951331, -0.00490540,
111	0.01719610, 0.15018204, 0.00882722, -0.00423197, -0.01919217,
112	-0.02963619, -0.01831342, -0.00524338, 0.00011379, -0.02566864
113	},
114	{
115	-0.07494341, -0.11348850, 0.00241343, -0.00803016, 0.00492438,
116	0.00711909, -0.00829147, 0.05793337, 0.02734209, 0.02059759,
117	-0.02770280, 0.14128338, 0.01532479, 0.00364307, 0.05968116,
118	-0.06497960, -0.08113941, 0.00319445, -0.00104222, 0.03553497
119	},
120	{
121	0.05948223, -0.08959930, 0.03269977, -0.03272374, -0.00365667,
122	-0.03423294, -0.06418925, -0.05902138, 0.05746317, -0.02580596,
123	0.01259572, 0.05848832, 0.00672666, 0.00233355, -0.05145149,
124	0.07348503, 0.11427955, 0.00142592, -0.01030651, -0.04862799
125	},
126	{
127	-0.01606880, 0.05200845, -0.01212967, -0.06824429, -0.00234304,
128	0.01094203, -0.07375538, 0.08808629, 0.12394822, 0.02231351,
129	-0.03608265, -0.06978045, -0.00618360, 0.00274747, -0.01921876,
130	-0.01541969, -0.02223856, -0.00107603, -0.01251777, 0.05412534
131	},
132	{
133	0.01688843, 0.05784728, -0.02256966, -0.07072251, -0.00422551,
134	-0.06261233, -0.08502830, 0.08925346, -0.08529597, 0.01519343,
135	-0.05008258, 0.10931873, 0.00521033, 0.02593305, -0.00717855,
136	0.02291527, 0.02527388, -0.00266188, -0.00871160, 0.02708135
137	},
138	{
139	-0.04233344, 0.00076379, 0.01571257, 0.04003092, 0.00901468,
140	0.00670577, 0.03459487, 0.12420216, -0.00067366, -0.01515094,
141	0.05306642, 0.04338407, 0.00511287, 0.01036639, -0.17867462,
142	-0.02289440, -0.03213205, 0.00017924, -0.01187362, -0.03933874
143	},
144	{
145	0.01284817, -0.01685622, 0.00724363, 0.01687952, -0.00882070,
146	-0.00555957, 0.01676246, -0.05560456, -0.00966893, 0.06197684,
147	-0.09058758, 0.00880607, 0.00108629, -0.08308956, -0.08056832,
148	-0.00413297, 0.02973107, 0.00092948, 0.07010111, 0.13007418
149	},
150	{
151	0.00700223, -0.01347574, 0.00691332, 0.03122905, 0.00310308,
152	0.00946862, 0.03455040, -0.06712536, -0.00304506, 0.04267941,
153	-0.10422292, -0.01127831, -0.00549798, 0.11680505, -0.03352701,
154	-0.00084536, 0.01631369, 0.00095063, -0.09570217, 0.06480321
155	}
156	};
157
158	void di_protdist::cats(di_cattype wcat)
159	{
160	/* define categories of amino acids */
161	aas b;
162
163	/* fundamental subgroups */
164	cat[(long) cys - (long) ala] = 1;
165	cat[(long) met - (long) ala] = 2;
166	cat[(long) val - (long) ala] = 3;
167	cat[(long) leu - (long) ala] = 3;
168	cat[(long) ileu - (long) ala] = 3;
169	cat[(long) gly - (long) ala] = 4;
170	cat[0] = 4;
171	cat[(long) ser - (long) ala] = 4;
172	cat[(long) thr - (long) ala] = 4;
173	cat[(long) pro - (long) ala] = 5;
174	cat[(long) phe - (long) ala] = 6;
175	cat[(long) tyr - (long) ala] = 6;
176	cat[(long) trp - (long) ala] = 6;
177	cat[(long) glu - (long) ala] = 7;
178	cat[(long) gln - (long) ala] = 7;
179	cat[(long) asp - (long) ala] = 7;
180	cat[(long) asn - (long) ala] = 7;
181	cat[(long) lys - (long) ala] = 8;
182	cat[(long) arg - (long) ala] = 8;
183	cat[(long) his - (long) ala] = 8;
184	if (wcat == george) {
185	/*
186	* George, Hunt and Barker: sulfhydryl, small hydrophobic,
187	* small hydrophilic, aromatic, acid/acid-amide/hydrophilic,
188	* basic
189	*/
190	for (b = ala; (long) b <= (long) val; b = (aas) ((long) b + 1)) {
191	if (cat[(long) b - (long) ala] == 3)
192	cat[(long) b - (long) ala] = 2;
193	if (cat[(long) b - (long) ala] == 5)
194	cat[(long) b - (long) ala] = 4;
195	}
196	}
197	if (wcat == chemical) {
198	/*
199	* Conn and Stumpf: monoamino, aliphatic, heterocyclic,
200	* aromatic, dicarboxylic, basic
201	*/
202	for (b = ala; (long) b <= (long) val; b = (aas) ((long) b + 1)) {
203	if (cat[(long) b - (long) ala] == 2)
204	cat[(long) b - (long) ala] = 1;
205	if (cat[(long) b - (long) ala] == 4)
206	cat[(long) b - (long) ala] = 3;
207	}
208	}
209	/* Ben Hall's personal opinion */
210	if (wcat != hall)
211	return;
212	for (b = ala; (long) b <= (long) val; b = (aas) ((long) b + 1)) {
213	if (cat[(long) b - (long) ala] == 3)
214	cat[(long) b - (long) ala] = 2;
215	}
216	} /* cats */
217
218
219	void di_protdist::maketrans()
220	{
221	/*
222	* Make up transition probability matrix from code and category
223	* tables
224	*/
225	long i, j, k, m, n, s;
226	double x, sum = 0;
227	long sub[3], newsub[3];
228	double f[4], g[4];
229	aas b1, b2;
230	double TEMP, TEMP1, TEMP2, TEMP3;
231
232	for (i = 0; i <= 19; i++) {
233	pi[i] = 0.0;
234	for (j = 0; j <= 19; j++)
235	prob[i][j] = 0.0;
236	}
237	f[0] = freqt;
238	f[1] = freqc;
239	f[2] = freqa;
240	f[3] = freqg;
241	g[0] = freqc + freqt;
242	g[1] = freqc + freqt;
243	g[2] = freqa + freqg;
244	g[3] = freqa + freqg;
245	TEMP = f[0];
246	TEMP1 = f[1];
247	TEMP2 = f[2];
248	TEMP3 = f[3];
249	fracchange = xi * (2 * f[0] * f[1] / g[0] + 2 * f[2] * f[3] / g[2]) +
250	xv * (1 - TEMP * TEMP - TEMP1 * TEMP1 - TEMP2 * TEMP2 - TEMP3 * TEMP3);
251	for (i = 0; i <= 3; i++) {
252	for (j = 0; j <= 3; j++) {
253	for (k = 0; k <= 3; k++) {
254	if (trans[i][j][k] != stop)
255	sum += f[i] * f[j] * f[k];
256	}
257	}
258	}
259	for (i = 0; i <= 3; i++) {
260	sub[0] = i + 1;
261	for (j = 0; j <= 3; j++) {
262	sub[1] = j + 1;
263	for (k = 0; k <= 3; k++) {
264	sub[2] = k + 1;
265	b1 = trans[i][j][k];
266	for (m = 0; m <= 2; m++) {
267	s = sub[m];
268	for (n = 1; n <= 4; n++) {
269	memcpy(newsub, sub, sizeof(long) * 3L);
270	newsub[m] = n;
271	x = f[i] * f[j] * f[k] / (3.0 * sum);
272	if (((s == 1 \|\| s == 2) && (n == 3 \|\| n == 4)) \|\|
273	((n == 1 \|\| n == 2) && (s == 3 \|\| s == 4)))
274	x = xv f[n - 1];
275	else
276	x = xi f[n - 1] / g[n - 1] + xv * f[n - 1];
277	b2 = trans[newsub[0] - 1][newsub[1] - 1][newsub[2] - 1];
278	if (b1 != stop) {
279	pi[b1] += x;
280	if (b2 != stop) {
281	if (cat[b1] != cat[b2]) {
282	prob[b1][b2] += x * ease;
283	prob[b1][b1] += x * (1.0 - ease);
284	} else
285	prob[b1][b2] += x;
286	} else
287	prob[b1][b1] += x;
288	}
289	}
290	}
291	}
292	}
293	}
294	for (i = 0; i <= 19; i++)
295	prob[i][i] -= pi[i];
296	for (i = 0; i <= 19; i++) {
297	for (j = 0; j <= 19; j++)
298	prob[i][j] /= sqrt(pi[i] * pi[j]);
299	}
300	/* computes pi^(1/2)Bpi^(-1/2) */
301	} /* maketrans */
302
303	void di_protdist::code()
304	{
305	/* make up table of the code 0 = u, 1 = c, 2 = a, 3 = g */
306
307	trans[0][0][0] = phe;
308	trans[0][0][1] = phe;
309	trans[0][0][2] = leu;
310	trans[0][0][3] = leu;
311	trans[0][1][0] = ser;
312	trans[0][1][1] = ser;
313	trans[0][1][2] = ser;
314	trans[0][1][3] = ser;
315	trans[0][2][0] = tyr;
316	trans[0][2][1] = tyr;
317	trans[0][2][2] = stop;
318	trans[0][2][3] = stop;
319	trans[0][3][0] = cys;
320	trans[0][3][1] = cys;
321	trans[0][3][2] = stop;
322	trans[0][3][3] = trp;
323	trans[1][0][0] = leu;
324	trans[1][0][1] = leu;
325	trans[1][0][2] = leu;
326	trans[1][0][3] = leu;
327	trans[1][1][0] = pro;
328	trans[1][1][1] = pro;
329	trans[1][1][2] = pro;
330	trans[1][1][3] = pro;
331	trans[1][2][0] = his;
332	trans[1][2][1] = his;
333	trans[1][2][2] = gln;
334	trans[1][2][3] = gln;
335	trans[1][3][0] = arg;
336	trans[1][3][1] = arg;
337	trans[1][3][2] = arg;
338	trans[1][3][3] = arg;
339	trans[2][0][0] = ileu;
340	trans[2][0][1] = ileu;
341	trans[2][0][2] = ileu;
342	trans[2][0][3] = met;
343	trans[2][1][0] = thr;
344	trans[2][1][1] = thr;
345	trans[2][1][2] = thr;
346	trans[2][1][3] = thr;
347	trans[2][2][0] = asn;
348	trans[2][2][1] = asn;
349	trans[2][2][2] = lys;
350	trans[2][2][3] = lys;
351	trans[2][3][0] = ser;
352	trans[2][3][1] = ser;
353	trans[2][3][2] = arg;
354	trans[2][3][3] = arg;
355	trans[3][0][0] = val;
356	trans[3][0][1] = val;
357	trans[3][0][2] = val;
358	trans[3][0][3] = val;
359	trans[3][1][0] = ala;
360	trans[3][1][1] = ala;
361	trans[3][1][2] = ala;
362	trans[3][1][3] = ala;
363	trans[3][2][0] = asp;
364	trans[3][2][1] = asp;
365	trans[3][2][2] = glu;
366	trans[3][2][3] = glu;
367	trans[3][3][0] = gly;
368	trans[3][3][1] = gly;
369	trans[3][3][2] = gly;
370	trans[3][3][3] = gly;
371	if (whichcode == mito)
372	trans[0][3][2] = trp;
373	if (whichcode == vertmito) {
374	trans[0][3][2] = trp;
375	trans[2][3][2] = stop;
376	trans[2][3][3] = stop;
377	trans[2][0][2] = met;
378	}
379	if (whichcode == flymito) {
380	trans[0][3][2] = trp;
381	trans[2][0][2] = met;
382	trans[2][3][2] = ser;
383	}
384	if (whichcode == yeastmito) {
385	trans[0][3][2] = trp;
386	trans[1][0][2] = thr;
387	trans[2][0][2] = met;
388	}
389	} /* code */
390
391	void di_protdist::transition()
392	{
393	/* calculations related to transition-transversion ratio */
394	double aa, bb, freqr, freqy, freqgr, freqty;
395
396	freqr = freqa + freqg;
397	freqy = freqc + freqt;
398	freqgr = freqg / freqr;
399	freqty = freqt / freqy;
400	aa = ttratio * freqr * freqy - freqa * freqg - freqc * freqt;
401	bb = freqa * freqgr + freqc * freqty;
402	xi = aa / (aa + bb);
403	xv = 1.0 - xi;
404	if (xi <= 0.0 && xi >= -epsilon)
405	xi = 0.0;
406	if (xi < 0.0) {
407	printf("THIS TRANSITION-TRANSVERSION RATIO IS IMPOSSIBLE WITH");
408	printf(" THESE BASE FREQUENCIES\n");
409	exit(-1);
410	}
411	} /* transition */
412
413	void di_protdist::givens(di_aa_matrix a,long i,long j,long n,double ctheta,double stheta,GB_BOOL left)
414	{
415	/* Givens transform at i,j for 1..n with angle theta */
416	long k;
417	double d;
418
419	for (k = 0; k < n; k++) {
420	if (left) {
421	d = ctheta * a[i - 1][k] + stheta * a[j - 1][k];
422	a[j - 1][k] = ctheta * a[j - 1][k] - stheta * a[i - 1][k];
423	a[i - 1][k] = d;
424	} else {
425	d = ctheta * a[k][i - 1] + stheta * a[k][j - 1];
426	a[k][j - 1] = ctheta * a[k][j - 1] - stheta * a[k][i - 1];
427	a[k][i - 1] = d;
428	}
429	}
430	} /* givens */
431
432	void di_protdist::coeffs(double x,double y,double c,double s,double accuracy)
433	{
434	/* compute cosine and sine of theta */
435	double root;
436
437	root = sqrt(x * x + y * y);
438	if (root < accuracy) {
439	*c = 1.0;
440	*s = 0.0;
441	} else {
442	*c = x / root;
443	*s = y / root;
444	}
445	} /* coeffs */
446
447	void di_protdist::tridiag(di_aa_matrix a,long n,double accuracy)
448	{
449	/* Givens tridiagonalization */
450	long i, j;
451	double s, c;
452
453	for (i = 2; i < n; i++) {
454	for (j = i + 1; j <= n; j++) {
455	coeffs(a[i - 2][i - 1], a[i - 2][j - 1], &c, &s, accuracy);
456	givens(a, i, j, n, c, s, GB_TRUE);
457	givens(a, i, j, n, c, s, GB_FALSE);
458	givens(eigvecs, i, j, n, c, s, GB_TRUE);
459	}
460	}
461	} /* tridiag */
462
463	void di_protdist::shiftqr(di_aa_matrix a, long n, double accuracy)
464	{
465	/* QR eigenvalue-finder */
466	long i, j;
467	double approx, s, c, d, TEMP, TEMP1;
468
469	for (i = n; i >= 2; i--) {
470	do {
471	TEMP = a[i - 2][i - 2] - a[i - 1][i - 1];
472	TEMP1 = a[i - 1][i - 2];
473	d = sqrt(TEMP * TEMP + TEMP1 * TEMP1);
474	approx = a[i - 2][i - 2] + a[i - 1][i - 1];
475	if (a[i - 1][i - 1] < a[i - 2][i - 2])
476	approx = (approx - d) / 2.0;
477	else
478	approx = (approx + d) / 2.0;
479	for (j = 0; j < i; j++)
480	a[j][j] -= approx;
481	for (j = 1; j < i; j++) {
482	coeffs(a[j - 1][j - 1], a[j][j - 1], &c, &s, accuracy);
483	givens(a, j, j + 1, i, c, s, GB_TRUE);
484	givens(a, j, j + 1, i, c, s, GB_FALSE);
485	givens(eigvecs, j, j + 1, n, c, s, GB_TRUE);
486	}
487	for (j = 0; j < i; j++)
488	a[j][j] += approx;
489	} while (fabs(a[i - 1][i - 2]) > accuracy);
490	}
491	} /* shiftqr */
492
493
494	void di_protdist::qreigen(di_aa_matrix proba,long n)
495	{
496	/* QR eigenvector/eigenvalue method for symmetric matrix */
497	double accuracy;
498	long i, j;
499
500	accuracy = 1.0e-6;
501	for (i = 0; i < n; i++) {
502	for (j = 0; j < n; j++)
503	eigvecs[i][j] = 0.0;
504	eigvecs[i][i] = 1.0;
505	}
506	tridiag(proba, n, accuracy);
507	shiftqr(proba, n, accuracy);
508	for (i = 0; i < n; i++)
509	eig[i] = proba[i][i];
510	for (i = 0; i <= 19; i++) {
511	for (j = 0; j <= 19; j++)
512	proba[i][j] = sqrt(pi[j]) * eigvecs[i][j];
513	}
514	/* proba[i][j] is the value of U' times pi^(1/2) */
515	} /* qreigen */
516
517
518	void di_protdist::pameigen()
519	{
520	/* eigenanalysis for PAM matrix, precomputed */
521	memcpy(prob, pamprobs, sizeof(pamprobs));
522	memcpy(eig, pameigs, sizeof(pameigs));
523	fracchange = 0.01;
524	} /* pameigen */
525
526	void di_protdist::build_exptteig(double tt){
527	int m;
528	for (m = 0; m <= 19; m++) {
529	exptteig[m] = exp(tt * eig[m]);
530	}
531	}
532
533	void di_protdist::predict(double /tt/, long nb1,long nb2)
534	{
535	/* make contribution to prediction of this aa pair */
536	long m;
537	double q;
538	double TEMP;
539	for (m = 0; m <= 19; m++) {
540	q = prob[m][nb1] * prob[m][nb2] * exptteig[m];
541	p += q;
542	TEMP = eig[m];
543	dp += TEMP * q;
544	d2p += TEMP * TEMP * q;
545	}
546	} /* predict */
547
548	void di_protdist::build_predikt_table(int pos){
549	int b1, b2;
550	double tt = pos_2_tt(pos);
551	build_exptteig(tt);
552	akt_slopes = slopes[pos] = (di_paa_matrix *) calloc(sizeof(di_paa_matrix), 1);
553	akt_curves = curves[pos] = (di_paa_matrix *) calloc(sizeof(di_paa_matrix), 1);
554	akt_infs = infs[pos] = (di_bool_matrix *) calloc(sizeof(di_bool_matrix), 1);
555
556	for (b1 = ala; b1 < di_max_paa; b1++) {
557	for (b2 = ala; b2 <= b1; b2++) {
558	if (b1 != stop && b1 != del && b1 != quest && b1 != unk &&
559	b2 != stop && b2 != del && b2 != quest && b2 != unk) {
560	p = 0.0;
561	dp = 0.0;
562	d2p = 0.0;
563	if (b1 != asx && b1 != glx && b2 != asx && b2 != glx){
564	predict(tt, b1, b2);
565	} else {
566	if (b1 == asx) {
567	if (b2 == asx) {
568	predict(tt, 2L, 2L);
569	predict(tt, 2L, 3L);
570	predict(tt, 3L, 2L);
571	predict(tt, 3L, 3L);
572	} else {
573	if (b2 == glx) {
574	predict(tt, 2L, 5L);
575	predict(tt, 2L, 6L);
576	predict(tt, 3L, 5L);
577	predict(tt, 3L, 6L);
578	} else {
579	predict(tt, 2L, b2);
580	predict(tt, 3L, b2);
581	}
582	}
583	} else {
584	if (b1 == glx) {
585	if (b2 == asx) {
586	predict(tt, 5L, 2L);
587	predict(tt, 5L, 3L);
588	predict(tt, 6L, 2L);
589	predict(tt, 6L, 3L);
590	} else {
591	if (b2 == glx) {
592	predict(tt, 5L, 5L);
593	predict(tt, 5L, 6L);
594	predict(tt, 6L, 5L);
595	predict(tt, 6L, 6L);
596	} else {
597	predict(tt, 5L, b2);
598	predict(tt, 6L, b2);
599	}
600	}
601	} else {
602	if (b2 == asx) {
603	predict(tt, b1, 2L);
604	predict(tt, b1, 3L);
605	predict(tt, b1, 2L);
606	predict(tt, b1, 3L);
607	} else if (b2 == glx) {
608	predict(tt, b1, 5L);
609	predict(tt, b1, 6L);
610	predict(tt, b1, 5L);
611	predict(tt, b1, 6L);
612	}
613	}
614	}
615	}
616	if (p > 0.0){
617	akt_slopes[0][b1][b2] = dp / p;
618	akt_curves[0][b1][b2] = d2p / p - dp * dp / (p * p);
619	akt_infs[0][b1][b2] = 0;
620	akt_slopes[0][b2][b1] = akt_slopes[0][b1][b2];
621	akt_curves[0][b2][b1] = akt_curves[0][b1][b2];
622	akt_infs[0][b2][b1] = 0;
623	}else{
624	akt_infs[0][b1][b2] = 1;
625	akt_infs[0][b2][b1] = 1;
626	}
627	}//if
628	}// b2
629
630	} //for b1
631	}
632
633	int di_protdist::tt_2_pos(double tt) {
634	int pos = (int)(tt * fracchange * di_resolution);
635	if (pos >= di_resolution * di_max_dist )
636	pos = di_resolution * di_max_dist - 1;
637	if (pos < 0)
638	pos = 0;
639	return pos;
640	}
641
642	double di_protdist::pos_2_tt(int pos) {
643	double tt = pos / (fracchange * di_resolution);
644	return tt+epsilon;
645	}
646
647	void di_protdist::build_akt_predikt(double tt)
648	{
649	/* take an aktual slope from the hash table, else calculate a new one */
650	int pos = tt_2_pos(tt);
651	if (!slopes[pos]){
652	build_predikt_table(pos);
653	}
654	akt_slopes = slopes[pos];
655	akt_curves = curves[pos];
656	akt_infs = infs[pos];
657	return;
658
659	}
660
661	const char *di_protdist::makedists()
662	{
663	/* compute the distances */
664	long i, j, k, m, n, iterations;
665	double delta, slope, curv;
666	int b1=0, b2=0;
667	double tt=0;
668	int pos;
669
670	for (i = 0; i < spp; i++) {
671	matrix->set(i,i,0.0);
672	{
673	double gauge = (double)i/(double)spp;
674	if (aw_status(gauge*gauge)) return "Aborted";
675	}
676	{
677	/* move all unknown characters to del */
678	ap_pro *seq1 = entries[i]->sequence_protein->sequence;
679	for (k = 0; k <chars ; k++) {
680	b1 = seq1[k];
681	if (b1 <=val) continue;
682	if (b1 == asx \|\| b1 == glx) continue;
683	seq1[k] = del;
684	}
685	}
686
687	for (j = 0; j < i ; j++) {
688	if (whichcat > kimura ) {
689	if (whichcat == pam)
690	tt = 10.0;
691	else
692	tt = 1.0;
693	delta = tt / 2.0;
694	iterations = 0;
695	do {
696	slope = 0.0;
697	curv = 0.0;
698	pos = tt_2_pos(tt);
699	tt = pos_2_tt(pos);
700	build_akt_predikt(tt);
701	ap_pro *seq1 = entries[i]->sequence_protein->sequence;
702	ap_pro *seq2 = entries[j]->sequence_protein->sequence;
703	for (k = chars; k >0; k--) {
704	b1 = *(seq1++);
705	b2 = *(seq2++);
706	if (predict_infinity(b1,b2)){
707	break;
708	}
709	slope += predict_slope(b1,b2);
710	curv += predict_curve(b1,b2);
711	}
712	iterations++;
713	if (!predict_infinity(b1,b2)) {
714	if (curv < 0.0) {
715	tt -= slope / curv;
716	if (tt > 10000.0) {
717	aw_message(GB_export_errorf("INFINITE DISTANCE BETWEEN SPECIES %ld AND %ld; -1.0 WAS WRITTEN\n", i, j));
718	tt = -1.0 / fracchange;
719	break;
720	}
721	int npos = tt_2_pos(tt);
722	int d = npos - pos; if (d<0) d=-d;
723	if (d<=1){ // cannot optimize
724	break;
725	}
726
727	} else {
728	if ((slope > 0.0 && delta < 0.0) \|\| (slope < 0.0 && delta > 0.0))
729	delta /= -2;
730	if (tt + delta < 0 && tt<= epsilon) {
731	break;
732	}
733	tt += delta;
734	}
735	} else {
736	delta /= -2;
737	tt += delta;
738	if (tt < 0) tt = 0;
739	}
740	} while (iterations < 20);
741	} else { // cat < kimura
742	m = 0;
743	n = 0;
744	ap_pro *seq1 = entries[i]->sequence_protein->sequence;
745	ap_pro *seq2 = entries[j]->sequence_protein->sequence;
746	for (k = chars; k >0; k--) {
747	b1 = *(seq1++);
748	b2 = *(seq2++);
749	if (b1 <= val && b2 <= val) {
750	if (b1 == b2) m++;
751	n++;
752	}
753	}
754	if (n < 5) { // no info
755	tt = -1.0;
756	}else{
757	switch (whichcat) {
758	case kimura:
759	{
760	double rel = 1 - (double) m / n;
761	double drel = 1.0 - rel - 0.2 * rel * rel;
762	if (drel < 0.0) {
763	aw_message(GB_export_errorf("DISTANCE BETWEEN SEQUENCES %3ld AND %3ld IS TOO LARGE FOR KIMURA FORMULA", i, j));
764	tt = -1.0;
765	}else{
766	tt = -log(drel);
767	}
768	}
769	break;
770	case none:
771	tt = (n-m)/(double)n;
772	break;
773	case similarity:
774	tt = m/(double)n;
775	break;
776	default:
777	di_assert(0);
778	break;
779	}
780	}
781	}
782	matrix->set(i,j,fracchange * tt);
783	}
784	}
785	return 0;
786	} /* makedists */
787
788
789	void di_protdist::clean_slopes(){
790	int i;
791	if (slopes) {
792	for (i=0;i<di_resolution*di_max_dist;i++) {
793	delete slopes[i]; slopes[i] = 0;
794	delete curves[i]; curves[i] = 0;
795	delete infs[i]; infs[i] = 0;
796	}
797	}
798	akt_slopes = 0;
799	akt_curves = 0;
800	akt_infs = 0;
801	}
802
803	di_protdist::~di_protdist(){
804	clean_slopes();
805	}
806
807	di_protdist::di_protdist(di_codetype codei, di_cattype cati, long nentries, DI_ENTRY *entriesi, long seq_len, AP_smatrix matrixi){
808	memset((char *)this,0,sizeof(di_protdist));
809	entries = entriesi;
810	matrix = matrixi;
811	freqa = .25;
812	freqc = .25;
813	freqg = .25;
814	freqt = .25;
815	ttratio = 2.0;
816	ease = 0.457;
817	spp = nentries;
818	chars = seq_len;
819	transition();
820	whichcode = codei;
821	whichcat = cati;
822	switch(cati){
823	case none:
824	case similarity:
825	case kimura:
826	fracchange = 1.0;
827	break;
828	case pam:
829	code();
830	pameigen();
831	break;
832	default:
833	code();
834	maketrans();
835	qreigen(prob, 20L);
836	break;
837	}
838	}

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

Download in other formats: