Context Navigation

constants.c

Visit:

Last change on this file was 10371, checked in by aboeckma, 11 years ago
updated mafft version. Added extensions (no svn ignore, yet)
File size: 33.6 KB

Line
1	#include "mltaln.h"
2	#include "miyata.h"
3	#include "miyata5.h"
4	#include "DNA.h"
5
6	#include "JTT.c"
7	#include "blosum.c"
8
9	#define DEBUG 0
10	#define TEST 0
11
12	#define NORMALIZE1 1
13
14
15	static int shishagonyuu( double in )
16	{
17	int out;
18	if ( in > 0.0 ) out = ( (int)( in + 0.5 ) );
19	else if( in == 0.0 ) out = ( 0 );
20	else if( in < 0.0 ) out = ( (int)( in - 0.5 ) );
21	else out = 0;
22	return( out );
23	}
24
25	static void ambiguousscore( int *amino_n, int n_dis[26][26] )
26	{
27	int i;
28	for( i=0; i<26; i++ )
29	{
30	n_dis[i][amino_n['r']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['g']][i] ) );
31	n_dis[i][amino_n['y']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['c']][i] + n_dis[amino_n['t']][i] ) );
32	n_dis[i][amino_n['k']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['g']][i] + n_dis[amino_n['t']][i] ) );
33	n_dis[i][amino_n['m']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['c']][i] ) );
34	n_dis[i][amino_n['s']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['g']][i] + n_dis[amino_n['c']][i] ) );
35	n_dis[i][amino_n['w']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['t']][i] ) );
36	n_dis[i][amino_n['b']] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['c']][i] + n_dis[amino_n['g']][i] + n_dis[amino_n['t']][i] ) );
37	n_dis[i][amino_n['d']] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['g']][i] + n_dis[amino_n['t']][i] ) );
38	n_dis[i][amino_n['h']] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['c']][i] + n_dis[amino_n['t']][i] ) );
39	n_dis[i][amino_n['v']] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['c']][i] + n_dis[amino_n['g']][i] ) );
40
41	n_dis[amino_n['r']][i] = n_dis[i][amino_n['r']];
42	n_dis[amino_n['y']][i] = n_dis[i][amino_n['y']];
43	n_dis[amino_n['k']][i] = n_dis[i][amino_n['k']];
44	n_dis[amino_n['m']][i] = n_dis[i][amino_n['m']];
45	n_dis[amino_n['s']][i] = n_dis[i][amino_n['s']];
46	n_dis[amino_n['w']][i] = n_dis[i][amino_n['w']];
47	n_dis[amino_n['b']][i] = n_dis[i][amino_n['b']];
48	n_dis[amino_n['d']][i] = n_dis[i][amino_n['d']];
49	n_dis[amino_n['h']][i] = n_dis[i][amino_n['h']];
50	n_dis[amino_n['v']][i] = n_dis[i][amino_n['v']];
51	}
52
53	i = amino_n['r']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['g']][amino_n['g']] ) );
54	i = amino_n['y']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['c']][amino_n['c']] + n_dis[amino_n['t']][amino_n['t']] ) );
55	i = amino_n['k']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['g']][amino_n['g']] + n_dis[amino_n['t']][amino_n['t']] ) );
56	i = amino_n['m']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['c']][amino_n['c']] ) );
57	i = amino_n['s']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['g']][amino_n['g']] + n_dis[amino_n['c']][amino_n['c']] ) );
58	i = amino_n['w']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['t']][amino_n['t']] ) );
59	i = amino_n['b']; n_dis[i][i] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['c']][amino_n['c']] + n_dis[amino_n['g']][amino_n['g']] + n_dis[amino_n['t']][amino_n['t']] ) );
60	i = amino_n['d']; n_dis[i][i] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['g']][amino_n['g']] + n_dis[amino_n['t']][amino_n['t']] ) );
61	i = amino_n['h']; n_dis[i][i] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['c']][amino_n['c']] + n_dis[amino_n['t']][amino_n['t']] ) );
62	i = amino_n['v']; n_dis[i][i] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['c']][amino_n['c']] + n_dis[amino_n['g']][amino_n['g']] ) );
63	}
64
65
66	static void calcfreq_nuc( int nseq, char *seq, double datafreq )
67	{
68	int i, j, l;
69	int aan;
70	double total;
71	for( i=0; i<4; i++ )
72	datafreq[i] = 0.0;
73	total = 0.0;
74	for( i=0; i<nseq; i++ )
75	{
76	l = strlen( seq[i] );
77	for( j=0; j<l; j++ )
78	{
79	aan = amino_n[(int)seq[i][j]];
80	if( aan == 4 ) aan = 3;
81	if( aan >= 0 && aan < 4 )
82	{
83	datafreq[aan] += 1.0;
84	total += 1.0;
85	}
86	}
87	}
88	total = 0.0; for( i=0; i<4; i++ ) total += datafreq[i];
89	for( i=0; i<4; i++ ) datafreq[i] /= (double)total;
90	for( i=0; i<4; i++ ) if( datafreq[i] < 0.0001 ) datafreq[i] = 0.0001;
91
92
93	total = 0.0; for( i=0; i<4; i++ ) total += datafreq[i];
94	// fprintf( stderr, "total = %f\n", total );
95	for( i=0; i<4; i++ ) datafreq[i] /= (double)total;
96
97	#if 0
98	fprintf( stderr, "\ndatafreq = " );
99	for( i=0; i<4; i++ )
100	fprintf( stderr, "%10.5f ", datafreq[i] );
101	fprintf( stderr, "\n" );
102	exit( 1 );
103	#endif
104	}
105
106	static void calcfreq( int nseq, char *seq, double datafreq )
107	{
108	int i, j, l;
109	int aan;
110	double total;
111	for( i=0; i<20; i++ )
112	datafreq[i] = 0.0;
113	total = 0.0;
114	for( i=0; i<nseq; i++ )
115	{
116	l = strlen( seq[i] );
117	for( j=0; j<l; j++ )
118	{
119	aan = amino_n[(int)seq[i][j]];
120	if( aan >= 0 && aan < 20 )
121	{
122	datafreq[aan] += 1.0;
123	total += 1.0;
124	}
125	}
126	}
127	total = 0.0; for( i=0; i<20; i++ ) total += datafreq[i];
128	for( i=0; i<20; i++ ) datafreq[i] /= (double)total;
129	for( i=0; i<20; i++ ) if( datafreq[i] < 0.0001 ) datafreq[i] = 0.0001;
130
131	fprintf( stderr, "datafreq = \n" );
132	for( i=0; i<20; i++ )
133	fprintf( stderr, "%f\n", datafreq[i] );
134
135	total = 0.0; for( i=0; i<20; i++ ) total += datafreq[i];
136	fprintf( stderr, "total = %f\n", total );
137	for( i=0; i<20; i++ ) datafreq[i] /= (double)total;
138	}
139
140	static void generatenuc1pam( double *pam1, int kimuraR, double freq )
141	{
142	int i, j;
143	double R[4][4], mut[4], total, tmp;
144
145	R[0][0] = 0.0; R[0][1] = kimuraR; R[0][2] = 1.0; R[0][3] = 1.0;
146	R[1][0] = kimuraR; R[1][1] = 0.0; R[1][2] = 1.0; R[1][3] = 1.0;
147	R[2][0] = 1.0; R[2][1] = 1.0; R[2][2] = 0.0; R[2][3] = kimuraR;
148	R[3][0] = 1.0; R[3][1] = 1.0; R[3][2] = kimuraR; R[3][3] = 0.0;
149
150
151	total = 0.0;
152	for( i=0; i<4; i++ )
153	{
154	tmp = 0.0;
155	for( j=0; j<4; j++ ) tmp += R[i][j] * freq[j];
156	mut[i] = tmp;
157	total += tmp * freq[i];
158	}
159	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
160	{
161	if( i != j ) pam1[i][j] = 0.01 / total * R[i][j] * freq[j];
162	else pam1[i][j] = 1.0 - 0.01 / total * mut[i];
163	}
164	}
165
166
167	void constants( int nseq, char **seq )
168	{
169	int i, j, x;
170	// double tmp;
171
172	if( dorp == 'd' ) /* DNA */
173	{
174	int k, m;
175	double average;
176	double **pamx = AllocateDoubleMtx( 11,11 );
177	double **pam1 = AllocateDoubleMtx( 4, 4 );
178	double *freq = AllocateDoubleVec( 4 );
179
180
181	scoremtx = -1;
182	if( RNAppenalty == NOTSPECIFIED ) RNAppenalty = DEFAULTRNAGOP_N;
183	if( RNAppenalty_ex == NOTSPECIFIED ) RNAppenalty_ex = DEFAULTRNAGEP_N;
184	if( ppenalty == NOTSPECIFIED ) ppenalty = DEFAULTGOP_N;
185	if( ppenalty_OP == NOTSPECIFIED ) ppenalty_OP = DEFAULTGOP_N;
186	if( ppenalty_ex == NOTSPECIFIED ) ppenalty_ex = DEFAULTGEP_N;
187	if( ppenalty_EX == NOTSPECIFIED ) ppenalty_EX = DEFAULTGEP_N;
188	if( poffset == NOTSPECIFIED ) poffset = DEFAULTOFS_N;
189	if( RNApthr == NOTSPECIFIED ) RNApthr = DEFAULTRNATHR_N;
190	if( pamN == NOTSPECIFIED ) pamN = DEFAULTPAMN;
191	if( kimuraR == NOTSPECIFIED ) kimuraR = 2;
192
193	RNApenalty = (int)( 3 * 600.0 / 1000.0 * RNAppenalty + 0.5 );
194	RNApenalty_ex = (int)( 3 * 600.0 / 1000.0 * RNAppenalty_ex + 0.5 );
195	// fprintf( stderr, "DEFAULTRNAGOP_N = %d\n", DEFAULTRNAGOP_N );
196	// fprintf( stderr, "RNAppenalty = %d\n", RNAppenalty );
197	// fprintf( stderr, "RNApenalty = %d\n", RNApenalty );
198
199
200	RNAthr = (int)( 3 * 600.0 / 1000.0 * RNApthr + 0.5 );
201	penalty = (int)( 3 * 600.0 / 1000.0 * ppenalty + 0.5);
202	penalty_OP = (int)( 3 * 600.0 / 1000.0 * ppenalty_OP + 0.5);
203	penalty_ex = (int)( 3 * 600.0 / 1000.0 * ppenalty_ex + 0.5);
204	penalty_EX = (int)( 3 * 600.0 / 1000.0 * ppenalty_EX + 0.5);
205	offset = (int)( 3 * 600.0 / 1000.0 * poffset + 0.5);
206	offsetFFT = (int)( 3 * 600.0 / 1000.0 * (-0) + 0.5);
207	offsetLN = (int)( 3 * 600.0 / 1000.0 * 100 + 0.5);
208	penaltyLN = (int)( 3 * 600.0 / 1000.0 * -2000 + 0.5);
209	penalty_exLN = (int)( 3 * 600.0 / 1000.0 * -100 + 0.5);
210	sprintf( modelname, "%s%d (%d), %6.3f (%6.3f), %6.3f (%6.3f)", rnakozo?"RNA":"DNA", pamN, kimuraR,
211	-(double)ppenalty0.001, -(double)ppenalty0.003, -(double)poffset0.001, -(double)poffset0.003 );
212
213	for( i=0; i<26; i++ ) amino[i] = locaminon[i];
214	for( i=0; i<0x80; i++ ) amino_n[i] = -1;
215	for( i=0; i<26; i++ ) amino_n[(int)amino[i]] = i;
216	if( fmodel == 1 )
217	{
218	calcfreq_nuc( nseq, seq, freq );
219	fprintf( stderr, "a, freq[0] = %f\n", freq[0] );
220	fprintf( stderr, "g, freq[1] = %f\n", freq[1] );
221	fprintf( stderr, "c, freq[2] = %f\n", freq[2] );
222	fprintf( stderr, "t, freq[3] = %f\n", freq[3] );
223	}
224	else
225	{
226	freq[0] = 0.25;
227	freq[1] = 0.25;
228	freq[2] = 0.25;
229	freq[3] = 0.25;
230	}
231
232
233	if( kimuraR == 9999 )
234	{
235	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
236	pamx[i][j] = (double)locn_disn[i][j];
237	#if NORMALIZE1
238	average = 0.0;
239	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
240	average += pamx[i][j];
241	average /= 16.0;
242
243	if( disp )
244	fprintf( stderr, "average = %f\n", average );
245
246	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
247	pamx[i][j] -= average;
248
249	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
250	pamx[i][j] *= 600.0 / average;
251
252	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
253	pamx[i][j] -= offset;
254	#endif
255	}
256	else
257	{
258	#if 0
259	double f = 0.99;
260	double s = (double)kimuraR / ( 2 + kimuraR ) * 0.01;
261	double v = (double)1 / ( 2 + kimuraR ) * 0.01;
262	pam1[0][0] = f; pam1[0][1] = s; pam1[0][2] = v; pam1[0][3] = v;
263	pam1[1][0] = s; pam1[1][1] = f; pam1[1][2] = v; pam1[1][3] = v;
264	pam1[2][0] = v; pam1[2][1] = v; pam1[2][2] = f; pam1[2][3] = s;
265	pam1[3][0] = v; pam1[3][1] = v; pam1[3][2] = s; pam1[3][3] = f;
266	#else
267	generatenuc1pam( pam1, kimuraR, freq );
268	#endif
269
270	fprintf( stderr, "generating %dPAM scoring matrix for nucleotides ... ", pamN );
271
272	if( disp )
273	{
274	fprintf( stderr, " TPM \n" );
275	for( i=0; i<4; i++ )
276	{
277	for( j=0; j<4; j++ )
278	fprintf( stderr, "%+#6.10f", pam1[i][j] );
279	fprintf( stderr, "\n" );
280	}
281	fprintf( stderr, "\n" );
282	}
283
284
285	MtxuntDouble( pamx, 4 );
286	for( x=0; x < pamN; x++ ) MtxmltDouble( pamx, pam1, 4 );
287
288	if( disp )
289	{
290	fprintf( stderr, " TPM \n" );
291	for( i=0; i<4; i++ )
292	{
293	for( j=0; j<4; j++ )
294	fprintf( stderr, "%+#6.10f", pamx[i][j] );
295	fprintf( stderr, "\n" );
296	}
297	fprintf( stderr, "\n" );
298	}
299
300	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
301	pamx[i][j] /= freq[j];
302	// pamx[i][j] /= 0.25;
303
304	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
305	{
306	if( pamx[i][j] == 0.0 )
307	{
308	fprintf( stderr, "WARNING: pamx[i][j] = 0.0 ?\n" );
309	pamx[i][j] = 0.00001; /* by J. Thompson */
310	}
311	pamx[i][j] = log10( pamx[i][j] ) * 1000.0;
312	}
313
314	if( disp )
315	{
316	fprintf( stderr, " after log\n" );
317	for( i=0; i<4; i++ )
318	{
319	for( j=0; j<4; j++ )
320	fprintf( stderr, "%+10.6f ", pamx[i][j] );
321	fprintf( stderr, "\n" );
322	}
323	fprintf( stderr, "\n" );
324	}
325
326
327	// ?????
328
329	average = 0.0;
330	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
331	average += pamx[i][j] * freq[i] * freq[j];
332	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
333	pamx[i][j] -= average;
334
335	average = 0.0;
336	for( i=0; i<4; i++ )
337	average += pamx[i][i] * 1.0 / 4.0;
338
339	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
340	pamx[i][j] *= 600.0 / average;
341
342
343	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
344	pamx[i][j] -= offset;
345
346	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
347	pamx[i][j] = shishagonyuu( pamx[i][j] );
348
349	if( disp )
350	{
351	fprintf( stderr, " after shishagonyuu\n" );
352	for( i=0; i<4; i++ )
353	{
354	for( j=0; j<4; j++ )
355	fprintf( stderr, "%+#6.10f", pamx[i][j] );
356	fprintf( stderr, "\n" );
357	}
358	fprintf( stderr, "\n" );
359	}
360	fprintf( stderr, "done\n" );
361	}
362
363	for( i=0; i<5; i++ )
364	{
365	pamx[4][i] = pamx[3][i];
366	pamx[i][4] = pamx[i][3];
367	}
368
369	for( i=5; i<10; i++ ) for( j=5; j<10; j++ )
370	{
371	pamx[i][j] = pamx[i-5][j-5];
372	}
373
374	if( disp )
375	{
376	fprintf( stderr, " before dis\n" );
377	for( i=0; i<4; i++ )
378	{
379	for( j=0; j<4; j++ )
380	fprintf( stderr, "%+#6.10f", pamx[i][j] );
381	fprintf( stderr, "\n" );
382	}
383	fprintf( stderr, "\n" );
384	}
385
386	if( disp )
387	{
388	fprintf( stderr, " score matrix \n" );
389	for( i=0; i<4; i++ )
390	{
391	for( j=0; j<4; j++ )
392	fprintf( stderr, "%+#6.10f", pamx[i][j] );
393	fprintf( stderr, "\n" );
394	}
395	fprintf( stderr, "\n" );
396	exit( 1 );
397	}
398
399	for( i=0; i<26; i++ ) amino[i] = locaminon[i];
400	for( i=0; i<26; i++ ) amino_grp[(int)amino[i]] = locgrpn[i];
401	for( i=0; i<26; i++ ) for( j=0; j<26; j++ ) n_dis[i][j] = 0;
402	for( i=0; i<10; i++ ) for( j=0; j<10; j++ ) n_dis[i][j] = shishagonyuu( pamx[i][j] );
403
404	ambiguousscore( amino_n, n_dis );
405
406	if( disp )
407	{
408	fprintf( stderr, " score matrix \n" );
409	for( i=0; i<26; i++ )
410	{
411	for( j=0; j<26; j++ )
412	fprintf( stderr, "%+6d", n_dis[i][j] );
413	fprintf( stderr, "\n" );
414	}
415	fprintf( stderr, "\n" );
416	fprintf( stderr, "penalty = %d, penalty_ex = %d\n", penalty, penalty_ex );
417	}
418
419	// RIBOSUM
420	#if 1
421	average = 0.0;
422	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
423	average += ribosum4[i][j] * freq[i] * freq[j];
424	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
425	ribosum4[i][j] -= average;
426
427	average = 0.0;
428	for( i=0; i<4; i++ ) for( j=0; j<4; j++ ) for( k=0; k<4; k++ ) for( m=0; m<4; m++ )
429	{
430	// if( i%4==0&&j%4==3 \|\| i%4==3&&j%4==0 \|\| i%4==1&&j%4==2 \|\| i%4==2&&j%4==1 \|\| i%4==1&&j%4==3 \|\| i%4==3&&j%4==1 )
431	// if( k%4==0&&m%4==3 \|\| k%4==3&&m%4==0 \|\| k%4==1&&m%4==2 \|\| k%4==2&&m%4==1 \|\| k%4==1&&m%4==3 \|\| k%4==3&&m%4==1 )
432	average += ribosum16[i4+j][k4+m] * freq[i] * freq[j] * freq[k] * freq[m];
433	}
434	for( i=0; i<16; i++ ) for( j=0; j<16; j++ )
435	ribosum16[i][j] -= average;
436
437	average = 0.0;
438	for( i=0; i<4; i++ )
439	average += ribosum4[i][i] * freq[i];
440	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
441	ribosum4[i][j] *= 600.0 / average;
442
443	average = 0.0;
444	average += ribosum16[04+3][04+3] * freq[0] * freq[3]; // AU
445	average += ribosum16[34+0][34+0] * freq[3] * freq[0]; // UA
446	average += ribosum16[14+2][14+2] * freq[1] * freq[2]; // CG
447	average += ribosum16[24+1][24+1] * freq[2] * freq[1]; // GC
448	average += ribosum16[14+3][14+3] * freq[1] * freq[3]; // GU
449	average += ribosum16[34+1][34+1] * freq[3] * freq[1]; // UG
450	for( i=0; i<16; i++ ) for( j=0; j<16; j++ )
451	ribosum16[i][j] *= 600.0 / average;
452
453
454	#if 1
455	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
456	ribosum4[i][j] -= offset; /* extending gap cost ?????*/
457	for( i=0; i<16; i++ ) for( j=0; j<16; j++ )
458	ribosum16[i][j] -= offset; /* extending gap cost ?????*/
459	#endif
460
461	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
462	ribosum4[i][j] = shishagonyuu( ribosum4[i][j] );
463	for( i=0; i<16; i++ ) for( j=0; j<16; j++ )
464	ribosum16[i][j] = shishagonyuu( ribosum16[i][j] );
465
466	if( disp )
467	{
468	fprintf( stderr, "ribosum after shishagonyuu\n" );
469	for( i=0; i<4; i++ )
470	{
471	for( j=0; j<4; j++ )
472	fprintf( stderr, "%+#6.10f", ribosum4[i][j] );
473	fprintf( stderr, "\n" );
474	}
475	fprintf( stderr, "\n" );
476	fprintf( stderr, "ribosum16 after shishagonyuu\n" );
477	for( i=0; i<16; i++ )
478	{
479	for( j=0; j<16; j++ )
480	fprintf( stderr, "%+#7.0f", ribosum16[i][j] );
481	fprintf( stderr, "\n" );
482	}
483	fprintf( stderr, "\n" );
484	}
485	fprintf( stderr, "done\n" );
486
487	#if 1
488	for( i=0; i<37; i++ ) for( j=0; j<37; j++ ) ribosumdis[i][j] = 0.0; //iru
489	for( m=0; m<9; m++ ) for( i=0; i<4; i++ ) // loop
490	for( k=0; k<9; k++ ) for( j=0; j<4; j++ ) ribosumdis[m4+i][k4+j] = ribosum4[i][j]; // loop-loop
491	// for( k=0; k<9; k++ ) for( j=0; j<4; j++ ) ribosumdis[m4+i][k4+j] = n_dis[i][j]; // loop-loop
492
493	for( i=0; i<16; i++ ) for( j=0; j<16; j++ ) ribosumdis[i+4][j+4] = ribosum16[i][j]; // stem5-stem5
494	for( i=0; i<16; i++ ) for( j=0; j<16; j++ ) ribosumdis[i+20][j+20] = ribosum16[i][j]; // stem5-stem5
495	#else // do not use ribosum
496	for( i=0; i<37; i++ ) for( j=0; j<37; j++ ) ribosumdis[i][j] = 0.0; //iru
497	for( m=0; m<9; m++ ) for( i=0; i<4; i++ ) // loop
498	for( k=0; k<9; k++ ) for( j=0; j<4; j++ ) ribosumdis[m4+i][k4+j] = n_dis[i][j]; // loop-loop
499	#endif
500
501	if( disp )
502	{
503	fprintf( stderr, "ribosumdis\n" );
504	for( i=0; i<37; i++ )
505	{
506	for( j=0; j<37; j++ )
507	fprintf( stderr, "%+5d", ribosumdis[i][j] );
508	fprintf( stderr, "\n" );
509	}
510	fprintf( stderr, "\n" );
511	}
512	fprintf( stderr, "done\n" );
513	#endif
514
515	FreeDoubleMtx( pam1 );
516	FreeDoubleMtx( pamx );
517	free( freq );
518
519	}
520	else if( dorp == 'p' && scoremtx == 1 ) /* Blosum */
521	{
522	double *freq;
523	double *freq1;
524	double *datafreq;
525	double average;
526	// double tmp;
527	double **n_distmp;
528
529	n_distmp = AllocateDoubleMtx( 20, 20 );
530	datafreq = AllocateDoubleVec( 20 );
531	freq = AllocateDoubleVec( 20 );
532
533	if( ppenalty == NOTSPECIFIED ) ppenalty = DEFAULTGOP_B;
534	if( ppenalty_OP == NOTSPECIFIED ) ppenalty_OP = DEFAULTGOP_B;
535	if( ppenalty_ex == NOTSPECIFIED ) ppenalty_ex = DEFAULTGEP_B;
536	if( ppenalty_EX == NOTSPECIFIED ) ppenalty_EX = DEFAULTGEP_B;
537	if( poffset == NOTSPECIFIED ) poffset = DEFAULTOFS_B;
538	if( pamN == NOTSPECIFIED ) pamN = 0;
539	if( kimuraR == NOTSPECIFIED ) kimuraR = 1;
540	penalty = (int)( 600.0 / 1000.0 * ppenalty + 0.5 );
541	penalty_OP = (int)( 600.0 / 1000.0 * ppenalty_OP + 0.5 );
542	penalty_ex = (int)( 600.0 / 1000.0 * ppenalty_ex + 0.5 );
543	penalty_EX = (int)( 600.0 / 1000.0 * ppenalty_EX + 0.5 );
544	offset = (int)( 600.0 / 1000.0 * poffset + 0.5 );
545	offsetFFT = (int)( 600.0 / 1000.0 * (-0) + 0.5);
546	offsetLN = (int)( 600.0 / 1000.0 * 100 + 0.5);
547	penaltyLN = (int)( 600.0 / 1000.0 * -2000 + 0.5);
548	penalty_exLN = (int)( 600.0 / 1000.0 * -100 + 0.5);
549
550	BLOSUMmtx( nblosum, n_distmp, freq, amino, amino_grp );
551	if( nblosum == -1 )
552	sprintf( modelname, "User-defined, %6.3f, %+6.3f, %+6.3f", -(double)ppenalty/1000, -(double)poffset/1000, -(double)ppenalty_ex/1000 );
553	else
554	sprintf( modelname, "BLOSUM%d, %6.3f, %+6.3f, %+6.3f", nblosum, -(double)ppenalty/1000, -(double)poffset/1000, -(double)ppenalty_ex/1000 );
555	#if 0
556	for( i=0; i<26; i++ ) amino[i] = locaminod[i];
557	for( i=0; i<26; i++ ) amino_grp[(int)amino[i]] = locgrpd[i];
558	for( i=0; i<0x80; i++ ) amino_n[i] = 0;
559	for( i=0; i<26; i++ ) amino_n[(int)amino[i]] = i;
560	#endif
561	for( i=0; i<0x80; i++ )amino_n[i] = -1;
562	for( i=0; i<26; i++) amino_n[(int)amino[i]] = i;
563	if( fmodel == 1 )
564	{
565	calcfreq( nseq, seq, datafreq );
566	freq1 = datafreq;
567	}
568	else
569	freq1 = freq;
570	#if TEST
571	fprintf( stderr, "raw scoreing matrix : \n" );
572	for( i=0; i<20; i++ )
573	{
574	for( j=0; j<20; j++ )
575	{
576	fprintf( stdout, "%6.2f", n_distmp[i][j] );
577	}
578	fprintf( stdout, "\n" );
579	}
580	#endif
581	if( fmodel == -1 )
582	average = 0.0;
583	else
584	{
585	for( i=0; i<20; i++ )
586	#if TEST
587	fprintf( stdout, "freq[%c] = %f, datafreq[%c] = %f, freq1[] = %f\n", amino[i], freq[i], amino[i], datafreq[i], freq1[i] );
588	#endif
589	average = 0.0;
590	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
591	average += n_distmp[i][j] * freq1[i] * freq1[j];
592	}
593	#if TEST
594	fprintf( stdout, "####### average2 = %f\n", average );
595	#endif
596
597	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
598	n_distmp[i][j] -= average;
599	#if TEST
600	fprintf( stdout, "average2 = %f\n", average );
601	fprintf( stdout, "after average substruction : \n" );
602	for( i=0; i<20; i++ )
603	{
604	for( j=0; j<20; j++ )
605	{
606	fprintf( stdout, "%6.2f", n_distmp[i][j] );
607	}
608	fprintf( stdout, "\n" );
609	}
610	#endif
611
612	average = 0.0;
613	for( i=0; i<20; i++ )
614	average += n_distmp[i][i] * freq1[i];
615	#if TEST
616	fprintf( stdout, "####### average1 = %f\n", average );
617	#endif
618
619	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
620	n_distmp[i][j] *= 600.0 / average;
621	#if TEST
622	fprintf( stdout, "after average division : \n" );
623	for( i=0; i<20; i++ )
624	{
625	for( j=0; j<=i; j++ )
626	{
627	fprintf( stdout, "%7.1f", n_distmp[i][j] );
628	}
629	fprintf( stdout, "\n" );
630	}
631	#endif
632
633	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
634	n_distmp[i][j] -= offset;
635	#if TEST
636	fprintf( stdout, "after offset substruction (offset = %d): \n", offset );
637	for( i=0; i<20; i++ )
638	{
639	for( j=0; j<=i; j++ )
640	{
641	fprintf( stdout, "%7.1f", n_distmp[i][j] );
642	}
643	fprintf( stdout, "\n" );
644	}
645	#endif
646	#if 0
647	/* Ãí°Õ ¡ª¡ª¡ª¡ª¡ª¡ª¡ª¡ª¡ª¡ª */
648	penalty -= offset;
649	#endif
650
651
652	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
653	n_distmp[i][j] = shishagonyuu( n_distmp[i][j] );
654
655	if( disp )
656	{
657	fprintf( stdout, " scoring matrix \n" );
658	for( i=0; i<20; i++ )
659	{
660	fprintf( stdout, "%c ", amino[i] );
661	for( j=0; j<20; j++ )
662	fprintf( stdout, "%5.0f", n_distmp[i][j] );
663	fprintf( stdout, "\n" );
664	}
665	fprintf( stdout, " " );
666	for( i=0; i<20; i++ )
667	fprintf( stdout, " %c", amino[i] );
668
669	average = 0.0;
670	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
671	average += n_distmp[i][j] * freq1[i] * freq1[j];
672	fprintf( stdout, "average = %f\n", average );
673
674	average = 0.0;
675	for( i=0; i<20; i++ )
676	average += n_distmp[i][i] * freq1[i];
677	fprintf( stdout, "itch average = %f\n", average );
678	fprintf( stderr, "parameters: %d, %d, %d\n", penalty, penalty_ex, offset );
679
680
681	exit( 1 );
682	}
683
684	for( i=0; i<26; i++ ) for( j=0; j<26; j++ ) n_dis[i][j] = 0;
685	for( i=0; i<20; i++ ) for( j=0; j<20; j++ ) n_dis[i][j] = (int)n_distmp[i][j];
686
687	FreeDoubleMtx( n_distmp );
688	FreeDoubleVec( datafreq );
689	FreeDoubleVec( freq );
690
691	fprintf( stderr, "done.\n" );
692
693	}
694	else if( dorp == 'p' && scoremtx == 2 ) /* Miyata-Yasunaga */
695	{
696	fprintf( stderr, "Not supported\n" );
697	exit( 1 );
698	for( i=0; i<26; i++ ) for( j=0; j<26; j++ ) n_dis[i][j] = locn_dism[i][j];
699	for( i=0; i<26; i++ ) if( i != 24 ) n_dis[i][24] = n_dis[24][i] = exgpm;
700	n_dis[24][24] = 0;
701	if( ppenalty == NOTSPECIFIED ) ppenalty = locpenaltym;
702	if( poffset == NOTSPECIFIED ) poffset = -20;
703	if( pamN == NOTSPECIFIED ) pamN = 0;
704	if( kimuraR == NOTSPECIFIED ) kimuraR = 1;
705
706	penalty = ppenalty;
707	offset = poffset;
708
709	sprintf( modelname, "Miyata-Yasunaga, %6.3f, %6.3f", -(double)ppenalty/1000, -(double)poffset/1000 );
710	for( i=0; i<26; i++ ) amino[i] = locaminom[i];
711	for( i=0; i<26; i++ ) amino_grp[(int)amino[i]] = locgrpm[i];
712	#if DEBUG
713	fprintf( stdout, "scoreing matrix : \n" );
714	for( i=0; i<26; i++ )
715	{
716	for( j=0; j<26; j++ )
717	{
718	fprintf( stdout, "%#5d", n_dis[i][j] );
719	}
720	fprintf( stdout, "\n" );
721	}
722	#endif
723	}
724	else /* JTT */
725	{
726	double **rsr;
727	double **pam1;
728	double **pamx;
729	double *freq;
730	double *freq1;
731	double *mutab;
732	double *datafreq;
733	double average;
734	double tmp;
735	double delta;
736
737	rsr = AllocateDoubleMtx( 20, 20 );
738	pam1 = AllocateDoubleMtx( 20, 20 );
739	pamx = AllocateDoubleMtx( 20, 20 );
740	freq = AllocateDoubleVec( 20 );
741	mutab = AllocateDoubleVec( 20 );
742	datafreq = AllocateDoubleVec( 20 );
743
744	if( ppenalty == NOTSPECIFIED ) ppenalty = DEFAULTGOP_J;
745	if( ppenalty_OP == NOTSPECIFIED ) ppenalty_OP = DEFAULTGOP_J;
746	if( ppenalty_ex == NOTSPECIFIED ) ppenalty_ex = DEFAULTGEP_J;
747	if( ppenalty_EX == NOTSPECIFIED ) ppenalty_EX = DEFAULTGEP_J;
748	if( poffset == NOTSPECIFIED ) poffset = DEFAULTOFS_J;
749	if( pamN == NOTSPECIFIED ) pamN = DEFAULTPAMN;
750	if( kimuraR == NOTSPECIFIED ) kimuraR = 1;
751
752	penalty = (int)( 600.0 / 1000.0 * ppenalty + 0.5 );
753	penalty_OP = (int)( 600.0 / 1000.0 * ppenalty_OP + 0.5 );
754	penalty_ex = (int)( 600.0 / 1000.0 * ppenalty_ex + 0.5 );
755	penalty_EX = (int)( 600.0 / 1000.0 * ppenalty_EX + 0.5 );
756	offset = (int)( 600.0 / 1000.0 * poffset + 0.5 );
757	offsetFFT = (int)( 600.0 / 1000.0 * (-0) + 0.5 );
758	offsetLN = (int)( 600.0 / 1000.0 * 100 + 0.5);
759	penaltyLN = (int)( 600.0 / 1000.0 * -2000 + 0.5);
760	penalty_exLN = (int)( 600.0 / 1000.0 * -100 + 0.5);
761
762	sprintf( modelname, "%s %dPAM, %6.3f, %6.3f", (TMorJTT==TM)?"Transmembrane":"JTT", pamN, -(double)ppenalty/1000, -(double)poffset/1000 );
763
764	JTTmtx( rsr, freq, amino, amino_grp, (int)(TMorJTT==TM) );
765
766	for( i=0; i<0x80; i++ ) amino_n[i] = -1;
767	for( i=0; i<26; i++ ) amino_n[(int)amino[i]] = i;
768	if( fmodel == 1 )
769	{
770	calcfreq( nseq, seq, datafreq );
771	freq1 = datafreq;
772	}
773	else
774	freq1 = freq;
775
776
777	#if TEST
778	fprintf( stdout, "rsr = \n" );
779	for( i=0; i<20; i++ )
780	{
781	for( j=0; j<20; j++ )
782	{
783	fprintf( stdout, "%9.2f ", rsr[i][j] );
784	}
785	fprintf( stdout, "\n" );
786	}
787	#endif
788
789
790	fprintf( stderr, "generating %dPAM %s scoring matrix for amino acids ... ", pamN, (TMorJTT==TM)?"Transmembrane":"JTT" );
791
792	tmp = 0.0;
793	for( i=0; i<20; i++ )
794	{
795	mutab[i] = 0.0;
796	for( j=0; j<20; j++ )
797	mutab[i] += rsr[i][j] * freq1[j];
798	tmp += mutab[i] * freq1[i];
799	}
800	#if TEST
801	fprintf( stdout, "mutability = \n" );
802	for( i=0; i<20; i++ )
803	fprintf( stdout, "%5.3f\n", mutab[i] );
804
805	fprintf( stdout, "tmp = %f\n", tmp );
806	#endif
807	delta = 0.01 / tmp;
808	for( i=0; i<20; i++ )
809	{
810	for( j=0; j<20; j++ )
811	{
812	if( i != j )
813	pam1[i][j] = delta * rsr[i][j] * freq1[j];
814	else
815	pam1[i][j] = 1.0 - delta * mutab[i];
816	}
817	}
818
819	if( disp )
820	{
821	fprintf( stdout, "pam1 = \n" );
822	for( i=0; i<20; i++ )
823	{
824	for( j=0; j<20; j++ )
825	{
826	fprintf( stdout, "%9.6f ", pam1[i][j] );
827	}
828	fprintf( stdout, "\n" );
829	}
830	}
831
832	MtxuntDouble( pamx, 20 );
833	for( x=0; x < pamN; x++ ) MtxmltDouble( pamx, pam1, 20 );
834
835	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
836	pamx[i][j] /= freq1[j];
837
838	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
839	{
840	if( pamx[i][j] == 0.0 )
841	{
842	fprintf( stderr, "WARNING: pamx[%d][%d] = 0.0?\n", i, j );
843	pamx[i][j] = 0.00001; /* by J. Thompson */
844	}
845	pamx[i][j] = log10( pamx[i][j] ) * 1000.0;
846	}
847
848	#if TEST
849	fprintf( stdout, "raw scoring matrix : \n" );
850	for( i=0; i<20; i++ )
851	{
852	for( j=0; j<20; j++ )
853	{
854	fprintf( stdout, "%5.0f", pamx[i][j] );
855	}
856	fprintf( stdout, "\n" );
857	}
858	average = tmp = 0.0;
859	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
860	{
861	average += pamx[i][j] * freq1[i] * freq1[j];
862	tmp += freq1[i] * freq1[j];
863	}
864	average /= tmp;
865	fprintf( stdout, "Zenbu average = %f, tmp = %f \n", average, tmp );
866	average = tmp = 0.0;
867	for( i=0; i<20; i++ ) for( j=i; j<20; j++ )
868	{
869	average += pamx[i][j] * freq1[i] * freq1[j];
870	tmp += freq1[i] * freq1[j];
871	}
872	average /= tmp;
873	fprintf( stdout, "Zenbu average2 = %f, tmp = %f \n", average, tmp );
874	average = tmp = 0.0;
875	for( i=0; i<20; i++ )
876	{
877	average += pamx[i][i] * freq1[i];
878	tmp += freq1[i];
879	}
880	average /= tmp;
881	fprintf( stdout, "Itch average = %f, tmp = %f \n", average, tmp );
882	#endif
883
884	#if NORMALIZE1
885	if( fmodel == -1 )
886	average = 0.0;
887	else
888	{
889	#if TEST
890	for( i=0; i<20; i++ )
891	fprintf( stdout, "freq[%c] = %f, datafreq[%c] = %f, freq1[] = %f\n", amino[i], freq[i], amino[i], datafreq[i], freq1[i] );
892	#endif
893	average = 0.0;
894	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
895	average += pamx[i][j] * freq1[i] * freq1[j];
896	}
897	#if TEST
898	fprintf( stdout, "####### average2 = %f\n", average );
899	#endif
900
901	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
902	pamx[i][j] -= average;
903	#if TEST
904	fprintf( stdout, "average2 = %f\n", average );
905	fprintf( stdout, "after average substruction : \n" );
906	for( i=0; i<20; i++ )
907	{
908	for( j=0; j<20; j++ )
909	{
910	fprintf( stdout, "%5.0f", pamx[i][j] );
911	}
912	fprintf( stdout, "\n" );
913	}
914	#endif
915
916	average = 0.0;
917	for( i=0; i<20; i++ )
918	average += pamx[i][i] * freq1[i];
919	#if TEST
920	fprintf( stdout, "####### average1 = %f\n", average );
921	#endif
922
923	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
924	pamx[i][j] *= 600.0 / average;
925	#if TEST
926	fprintf( stdout, "after average division : \n" );
927	for( i=0; i<20; i++ )
928	{
929	for( j=0; j<=i; j++ )
930	{
931	fprintf( stdout, "%5.0f", pamx[i][j] );
932	}
933	fprintf( stdout, "\n" );
934	}
935	#endif
936
937	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
938	pamx[i][j] -= offset;
939	#if TEST
940	fprintf( stdout, "after offset substruction (offset = %d): \n", offset );
941	for( i=0; i<20; i++ )
942	{
943	for( j=0; j<=i; j++ )
944	{
945	fprintf( stdout, "%5.0f", pamx[i][j] );
946	}
947	fprintf( stdout, "\n" );
948	}
949	#endif
950	#if 0
951	/* Ãí°Õ ¡ª¡ª¡ª¡ª¡ª¡ª¡ª¡ª¡ª¡ª */
952	penalty -= offset;
953	#endif
954
955
956	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
957	pamx[i][j] = shishagonyuu( pamx[i][j] );
958
959	#else
960
961	average = 0.0;
962	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
963	average += pamx[i][j];
964	average /= 400.0;
965
966	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
967	{
968	pamx[i][j] -= average;
969	pamx[i][j] = shishagonyuu( pamx[i][j] );
970	}
971	#endif
972	if( disp )
973	{
974	fprintf( stdout, " scoring matrix \n" );
975	for( i=0; i<20; i++ )
976	{
977	fprintf( stdout, "%c ", amino[i] );
978	for( j=0; j<20; j++ )
979	fprintf( stdout, "%5.0f", pamx[i][j] );
980	fprintf( stdout, "\n" );
981	}
982	fprintf( stdout, " " );
983	for( i=0; i<20; i++ )
984	fprintf( stdout, " %c", amino[i] );
985
986	average = 0.0;
987	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
988	average += pamx[i][j] * freq1[i] * freq1[j];
989	fprintf( stdout, "average = %f\n", average );
990
991	average = 0.0;
992	for( i=0; i<20; i++ )
993	average += pamx[i][i] * freq1[i];
994	fprintf( stdout, "itch average = %f\n", average );
995	fprintf( stderr, "parameters: %d, %d, %d\n", penalty, penalty_ex, offset );
996
997
998	exit( 1 );
999	}
1000
1001	for( i=0; i<26; i++ ) for( j=0; j<26; j++ ) n_dis[i][j] = 0;
1002	for( i=0; i<20; i++ ) for( j=0; j<20; j++ ) n_dis[i][j] = (int)pamx[i][j];
1003
1004	fprintf( stderr, "done.\n" );
1005	FreeDoubleMtx( rsr );
1006	FreeDoubleMtx( pam1 );
1007	FreeDoubleMtx( pamx );
1008	FreeDoubleVec( freq );
1009	FreeDoubleVec( mutab );
1010	FreeDoubleVec( datafreq );
1011	}
1012	fprintf( stderr, "scoremtx = %d\n", scoremtx );
1013
1014	#if DEBUG
1015	fprintf( stderr, "scoremtx = %d\n", scoremtx );
1016	fprintf( stderr, "amino[] = %s\n", amino );
1017	#endif
1018
1019	for( i=0; i<0x80; i++ )amino_n[i] = -1;
1020	for( i=0; i<26; i++) amino_n[(int)amino[i]] = i;
1021	for( i=0; i<0x80; i++ ) for( j=0; j<0x80; j++ ) amino_dis[i][j] = 0;
1022	for( i=0; i<0x80; i++ ) for( j=0; j<0x80; j++ ) amino_disLN[i][j] = 0;
1023	for( i=0; i<0x80; i++ ) for( j=0; j<0x80; j++ ) amino_dis_consweight_multi[i][j] = 0.0;
1024	for( i=0; i<26; i++) for( j=0; j<26; j++ )
1025	{
1026	amino_dis[(int)amino[i]][(int)amino[j]] = n_dis[i][j];
1027	n_dis_consweight_multi[i][j] = (float)n_dis[i][j] * consweight_multi;
1028	amino_dis_consweight_multi[(int)amino[i]][(int)amino[j]] = (double)n_dis[i][j] * consweight_multi;
1029	}
1030
1031	if( dorp == 'd' ) /* DNA */
1032	{
1033	for( i=0; i<5; i++) for( j=0; j<5; j++ )
1034	amino_disLN[(int)amino[i]][(int)amino[j]] = n_dis[i][j] + offset - offsetLN;
1035	for( i=5; i<10; i++) for( j=5; j<10; j++ )
1036	amino_disLN[(int)amino[i]][(int)amino[j]] = n_dis[i][j] + offset - offsetLN;
1037	for( i=0; i<5; i++) for( j=0; j<5; j++ )
1038	n_disFFT[i][j] = n_dis[i][j] + offset - offsetFFT;
1039	for( i=5; i<10; i++) for( j=5; j<10; j++ )
1040	n_disFFT[i][j] = n_dis[i][j] + offset - offsetFFT;
1041	}
1042	else // protein
1043	{
1044	for( i=0; i<20; i++) for( j=0; j<20; j++ )
1045	amino_disLN[(int)amino[i]][(int)amino[j]] = n_dis[i][j] + offset - offsetLN;
1046	for( i=0; i<20; i++) for( j=0; j<20; j++ )
1047	n_disFFT[i][j] = n_dis[i][j] + offset - offsetFFT;
1048	}
1049
1050	#if 0
1051	fprintf( stderr, "amino_dis (offset = %d): \n", offset );
1052	for( i=0; i<20; i++ )
1053	{
1054	for( j=0; j<20; j++ )
1055	{
1056	fprintf( stderr, "%5d", amino_dis[(int)amino[i]][(int)amino[j]] );
1057	}
1058	fprintf( stderr, "\n" );
1059	}
1060
1061	fprintf( stderr, "amino_disLN (offsetLN = %d): \n", offsetLN );
1062	for( i=0; i<20; i++ )
1063	{
1064	for( j=0; j<20; j++ )
1065	{
1066	fprintf( stderr, "%5d", amino_disLN[(int)amino[i]][(int)amino[j]] );
1067	}
1068	fprintf( stderr, "\n" );
1069	}
1070
1071	fprintf( stderr, "n_dis (offset = %d): \n", offset );
1072	for( i=0; i<26; i++ )
1073	{
1074	for( j=0; j<26; j++ )
1075	{
1076	fprintf( stderr, "%5d", n_dis[i][j] );
1077	}
1078	fprintf( stderr, "\n" );
1079	}
1080
1081	fprintf( stderr, "n_disFFT (offsetFFT = %d): \n", offsetFFT );
1082	for( i=0; i<26; i++ )
1083	{
1084	for( j=0; j<26; j++ )
1085	{
1086	fprintf( stderr, "%5d", n_disFFT[i][j] );
1087	}
1088	fprintf( stderr, "\n" );
1089	}
1090	exit( 1 );
1091	#endif
1092
1093
1094	ppid = 0;
1095
1096
1097	if( fftThreshold == NOTSPECIFIED )
1098	{
1099	fftThreshold = FFT_THRESHOLD;
1100	}
1101	if( fftWinSize == NOTSPECIFIED )
1102	{
1103	if( dorp == 'd' )
1104	fftWinSize = FFT_WINSIZE_D;
1105	else
1106	fftWinSize = FFT_WINSIZE_P;
1107	}
1108
1109
1110	if( fftscore )
1111	{
1112	double av, sd;
1113
1114	for( i=0; i<20; i++ ) polarity[i] = polarity_[i];
1115	for( av=0.0, i=0; i<20; i++ ) av += polarity[i];
1116	av /= 20.0;
1117	for( sd=0.0, i=0; i<20; i++ ) sd += ( polarity[i]-av ) * ( polarity[i]-av );
1118	sd /= 20.0; sd = sqrt( sd );
1119	for( i=0; i<20; i++ ) polarity[i] -= av;
1120	for( i=0; i<20; i++ ) polarity[i] /= sd;
1121
1122	for( i=0; i<20; i++ ) volume[i] = volume_[i];
1123	for( av=0.0, i=0; i<20; i++ ) av += volume[i];
1124	av /= 20.0;
1125	for( sd=0.0, i=0; i<20; i++ ) sd += ( volume[i]-av ) * ( volume[i]-av );
1126	sd /= 20.0; sd = sqrt( sd );
1127	for( i=0; i<20; i++ ) volume[i] -= av;
1128	for( i=0; i<20; i++ ) volume[i] /= sd;
1129
1130	#if 0
1131	for( i=0; i<20; i++ ) fprintf( stdout, "amino=%c, pol = %f<-%f, vol = %f<-%f\n", amino[i], polarity[i], polarity_[i], volume[i], volume_[i] );
1132	for( i=0; i<20; i++ ) fprintf( stdout, "%c %+5.3f %+5.3f\n", amino[i], volume[i], polarity[i] );
1133	#endif
1134	}
1135	}

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

Context Navigation

source: branches/items/GDE/MAFFT/mafft-7.055-with-extensions/core/constants.c

Download in other formats: