1 #include "mltaln.h"
2 
3 static int upperCase = 0;
4 
5 #define DEBUG   0
6 #define IODEBUG 0
7 
creverse(char f)8 char creverse( char f )
9 {
10 	static char *table = NULL;
11 	if( table == NULL )
12 	{
13 		int i;
14 		table = AllocateCharVec(0x80);
15 		for( i=0; i<0x80; i++ ) table[i] = i;
16 		table['A'] = 'T';
17 		table['C'] = 'G';
18 		table['G'] = 'C';
19 		table['T'] = 'A';
20 		table['U'] = 'A';
21 		table['M'] = 'K';
22 		table['R'] = 'Y';
23 		table['W'] = 'W';
24 		table['S'] = 'S';
25 		table['Y'] = 'R';
26 		table['K'] = 'M';
27 		table['V'] = 'B';
28 		table['H'] = 'D';
29 		table['D'] = 'H';
30 		table['B'] = 'V';
31 		table['N'] = 'N';
32 		table['a'] = 't';
33 		table['c'] = 'g';
34 		table['g'] = 'c';
35 		table['t'] = 'a';
36 		table['u'] = 'a';
37 		table['m'] = 'k';
38 		table['r'] = 'y';
39 		table['w'] = 'w';
40 		table['s'] = 's';
41 		table['y'] = 'r';
42 		table['k'] = 'm';
43 		table['v'] = 'b';
44 		table['h'] = 'd';
45 		table['d'] = 'h';
46 		table['b'] = 'v';
47 		table['n'] = 'n';
48 //		table['-'] = '-';
49 //		table['.'] = '.';
50 //		table['*'] = '*';
51 	}
52 	return( table[(int)f] );
53 }
54 
sreverse(char * r,char * s)55 void sreverse( char *r, char *s )
56 {
57 	r += strlen( s );
58 	*r-- = 0;
59 	while( *s )
60 		*r-- = creverse( *s++ );
61 //		*r-- = ( *s++ );
62 }
63 
gappick_samestring(char * seq)64 void gappick_samestring( char *seq )
65 {
66 	char *aseq = seq;
67 
68 	for( ; *seq != 0; seq++ )
69 	{
70 		if( *seq != '-' )
71 			*aseq++ = *seq;
72 	}
73 	*aseq = 0;
74 }
75 
76 #if 0
77 
78 static int addlocalhom2( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa, int skip )
79 {
80 	int pos1, pos2, start1, start2, end1, end2;
81 	char *pt1, *pt2;
82 	int iscore;
83 	int isumscore;
84 	int sumoverlap;
85 	LocalHom *tmppt;
86 	int st;
87 	int nlocalhom = 0;
88 	pt1 = al1; pt2 = al2;
89 	pos1 = off1; pos2 = off2;
90 
91 	isumscore = 0;
92 	sumoverlap = 0;
93 
94 #if 0
95 	fprintf( stderr, "nlocalhom = %d in addlocalhom\n", nlocalhom );
96 	fprintf( stderr, "al1 = %s, al2 = %s\n", al1, al2 );
97 	fprintf( stderr, "off1 = %d, off2 = %d\n", off1, off2 );
98 	fprintf( stderr, "localhopt = %p, skip = %d\n", localhompt, skip );
99 	fprintf( stderr, "pt1 = \n%s\n, pt2 = \n%s\n", pt1, pt2 );
100 #endif
101 
102 	if( skip )
103 	{
104 		while( --skip > 0 ) localhompt = localhompt->next;
105 		localhompt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
106 		localhompt = localhompt->next;
107 //		fprintf( stderr, "tmppt = %p, localhompt = %p\n", tmppt, localhompt );
108 	}
109 	tmppt = localhompt;
110 
111 	st = 0;
112 	iscore = 0;
113 	while( *pt1 != 0 )
114 	{
115 //		fprintf( stderr, "In in while loop\n" );
116 //		fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
117 		if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
118 		{
119 			end1 = pos1 - 1;
120 			end2 = pos2 - 1;
121 
122 			if( nlocalhom++ > 0 )
123 			{
124 //				fprintf( stderr, "reallocating ...\n" );
125 				tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
126 //				fprintf( stderr, "done\n" );
127 				tmppt = tmppt->next;
128 				tmppt->next = NULL;
129 			}
130 			tmppt->start1 = start1;
131 			tmppt->start2 = start2;
132 			tmppt->end1   = end1  ;
133 			tmppt->end2   = end2  ;
134 
135 #if 1
136 			isumscore += iscore;
137 			sumoverlap += end2-start2+1;
138 #else
139 			tmppt->overlapaa   = end2-start2+1;
140 			tmppt->opt = iscore * 5.8 / 600;
141 			tmppt->overlapaa   = overlapaa;
142 			tmppt->opt = (double)opt;
143 #endif
144 
145 #if 0
146 			fprintf( stderr, "iscore (1)= %d\n", iscore );
147 			fprintf( stderr, "al1: %d - %d\n", start1, end1 );
148 			fprintf( stderr, "al2: %d - %d\n", start2, end2 );
149 #endif
150 			iscore = 0;
151 			st = 0;
152 		}
153 		else if( *pt1 != '-' && *pt2 != '-' )
154 		{
155 			if( st == 0 )
156 			{
157 				start1 = pos1; start2 = pos2;
158 				st = 1;
159 			}
160 			iscore += n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]];
161 //			fprintf( stderr, "%c-%c, score(0) = %d\n", *pt1, *pt2, iscore );
162 		}
163 		if( *pt1++ != '-' ) pos1++;
164 		if( *pt2++ != '-' ) pos2++;
165 	}
166 
167 	if( st )
168 	{
169 		if( nlocalhom++ > 0 )
170 		{
171 //			fprintf( stderr, "reallocating ...\n" );
172 			tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
173 //			fprintf( stderr, "done\n" );
174 			tmppt = tmppt->next;
175 			tmppt->next = NULL;
176 		}
177 		end1 = pos1 - 1;
178 		end2 = pos2 - 1;
179 		tmppt->start1 = start1;
180 		tmppt->start2 = start2;
181 		tmppt->end1   = end1  ;
182 		tmppt->end2   = end2  ;
183 
184 #if 1
185 		isumscore += iscore;
186 		sumoverlap += end2-start2+1;
187 #else
188 		tmppt->overlapaa   = end2-start2+1;
189 		tmppt->opt = (double)iscore * 5.8 / 600;
190 		tmppt->overlapaa   = overlapaa;
191 		tmppt->opt = (double)opt;
192 #endif
193 #if 0
194 		fprintf( stderr, "score (2)= %d\n", iscore );
195 		fprintf( stderr, "al1: %d - %d\n", start1, end1 );
196 		fprintf( stderr, "al2: %d - %d\n", start2, end2 );
197 #endif
198 	}
199 
200 	for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
201 	{
202 		tmppt->overlapaa = sumoverlap;
203 		tmppt->opt = (double)sumscore * 5.8 / 600 / sumoverlap;
204 	}
205 	return( nlocalhom );
206 }
207 
208 #endif
209 
210 
211 
addlocalhom_r(char * al1,char * al2,LocalHom * localhompt,int off1,int off2,int opt,int overlapaa,int skip)212 static int addlocalhom_r( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa, int skip )
213 {
214 	int pos1, pos2, start1, start2, end1, end2;
215 	char *pt1, *pt2;
216 	double score;
217 	double sumscore;
218 	int sumoverlap;
219 	LocalHom *tmppt = NULL; // by D.Mathog, a guess
220 	int st;
221 	int nlocalhom = 0;
222 	pt1 = al1; pt2 = al2;
223 	pos1 = off1; pos2 = off2;
224 
225 	sumscore = 0.0;
226 	sumoverlap = 0;
227 	start1 = 0; // by D.Mathog, a guess
228 	start2 = 0; // by D.Mathog, a guess
229 
230 #if 0
231 	fprintf( stderr, "nlocalhom = %d in addlocalhom\n", nlocalhom );
232 	fprintf( stderr, "al1 = %s, al2 = %s\n", al1, al2 );
233 	fprintf( stderr, "off1 = %d, off2 = %d\n", off1, off2 );
234 	fprintf( stderr, "localhopt = %p, skip = %d\n", localhompt, skip );
235 #endif
236 	fprintf( stderr, "pt1 = \n%s\n, pt2 = \n%s\n", pt1, pt2 );
237 
238 	if( skip )
239 	{
240 		while( --skip > 0 ) localhompt = localhompt->next;
241 		localhompt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
242 		localhompt = localhompt->next;
243 		fprintf( stderr, "tmppt = %p, localhompt = %p\n", (void *)tmppt, (void *)localhompt );
244 	}
245 	tmppt = localhompt;
246 
247 	st = 0;
248 	score = 0.0;
249 	while( *pt1 != 0 )
250 	{
251 		fprintf( stderr, "In in while loop\n" );
252 		fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
253 		if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
254 		{
255 			end1 = pos1 - 1;
256 			end2 = pos2 - 1;
257 
258 			if( nlocalhom++ > 0 )
259 			{
260 //				fprintf( stderr, "reallocating ...\n" );
261 				tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
262 //				fprintf( stderr, "done\n" );
263 				tmppt = tmppt->next;
264 				tmppt->next = NULL;
265 			}
266 			tmppt->start1 = start1;
267 			tmppt->start2 = start2;
268 			tmppt->end1   = end1  ;
269 			tmppt->end2   = end2  ;
270 
271 #if 1
272 			sumscore += score;
273 			sumoverlap += end2-start2+1;
274 #else
275 			tmppt->overlapaa   = end2-start2+1;
276 			tmppt->opt = score * 5.8 / 600;
277 			tmppt->overlapaa   = overlapaa;
278 			tmppt->opt = (double)opt;
279 #endif
280 
281 			fprintf( stderr, "score (1)= %f\n", score );
282 			fprintf( stderr, "al1: %d - %d\n", start1, end1 );
283 			fprintf( stderr, "al2: %d - %d\n", start2, end2 );
284 			score = 0.0;
285 			st = 0;
286 		}
287 		else if( *pt1 != '-' && *pt2 != '-' )
288 		{
289 			if( st == 0 )
290 			{
291 				start1 = pos1; start2 = pos2;
292 				st = 1;
293 			}
294 			score += (double)n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]];
295 //			fprintf( stderr, "%c-%c, score(0) = %f\n", *pt1, *pt2, score );
296 		}
297 		if( *pt1++ != '-' ) pos1++;
298 		if( *pt2++ != '-' ) pos2++;
299 	}
300 	if( nlocalhom++ > 0 )
301 	{
302 //		fprintf( stderr, "reallocating ...\n" );
303 		tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
304 //		fprintf( stderr, "done\n" );
305 		tmppt = tmppt->next;
306 		tmppt->next = NULL;
307 	}
308 	end1 = pos1 - 1;
309 	end2 = pos2 - 1;
310 	tmppt->start1 = start1;
311 	tmppt->start2 = start2;
312 	tmppt->end1   = end1  ;
313 	tmppt->end2   = end2  ;
314 
315 #if 1
316 	sumscore += score;
317 	sumoverlap += end2-start2+1;
318 #else
319 	tmppt->overlapaa   = end2-start2+1;
320 	tmppt->opt = score * 5.8 / 600;
321 	tmppt->overlapaa   = overlapaa;
322 	tmppt->opt = (double)opt;
323 #endif
324 
325 	fprintf( stderr, "score (2)= %f\n", score );
326 	fprintf( stderr, "al1: %d - %d\n", start1, end1 );
327 	fprintf( stderr, "al2: %d - %d\n", start2, end2 );
328 
329 	for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
330 	{
331 		tmppt->overlapaa = sumoverlap;
332 		tmppt->opt = sumscore * 5.8 / 600 / sumoverlap;
333 	}
334 	return( nlocalhom );
335 }
putlocalhom3(char * al1,char * al2,LocalHom * localhompt,int off1,int off2,int opt,int overlapaa)336 void putlocalhom3( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
337 {
338 	int pos1, pos2, start1, start2, end1, end2;
339 	char *pt1, *pt2;
340 	double score;
341 	double sumscore;
342 	int sumoverlap;
343 	LocalHom *tmppt;
344 	LocalHom *subnosento;
345 	int st;
346 	int saisho;
347 
348 	pt1 = al1; pt2 = al2;
349 	pos1 = off1; pos2 = off2;
350 
351 	sumscore = 0.0;
352 	sumoverlap = 0;
353 	start1 = 0; // by Mathog, a guess
354 	start2 = 0; // by Mathog, a guess
355 
356 	subnosento = localhompt;
357 	while( subnosento->next ) subnosento = subnosento->next;
358 	tmppt = subnosento;
359 
360 	saisho = ( localhompt->nokori == 0 );
361 
362 	fprintf( stderr, "localhompt = %p\n", (void *)localhompt );
363 	fprintf( stderr, "tmppt = %p\n", (void *)tmppt );
364 	fprintf( stderr, "subnosento = %p\n", (void *)subnosento );
365 
366 	st = 0;
367 	score = 0.0;
368 	while( *pt1 != 0 )
369 	{
370 //		fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
371 		if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
372 		{
373 			end1 = pos1 - 1;
374 			end2 = pos2 - 1;
375 
376 			if( localhompt->nokori++ > 0 )
377 			{
378 //				fprintf( stderr, "reallocating ...\n" );
379 				tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
380 //				fprintf( stderr, "done\n" );
381 				tmppt = tmppt->next;
382 				tmppt->next = NULL;
383 			}
384 			tmppt->start1 = start1;
385 			tmppt->start2 = start2;
386 			tmppt->end1   = end1  ;
387 			tmppt->end2   = end2  ;
388 
389 #if 1
390 			if( divpairscore )
391 			{
392 				tmppt->overlapaa   = end2-start2+1;
393 				tmppt->opt = score / tmppt->overlapaa * 5.8 / 600;
394 			}
395 			else
396 			{
397 				sumscore += score;
398 				sumoverlap += end2-start2+1;
399 			}
400 #else
401 			tmppt->overlapaa   = overlapaa;
402 			tmppt->opt = (double)opt;
403 #endif
404 
405 #if 0
406 			fprintf( stderr, "score (1)= %f\n", score );
407 			fprintf( stderr, "al1: %d - %d\n", start1, end1 );
408 			fprintf( stderr, "al2: %d - %d\n", start2, end2 );
409 #endif
410 			score = 0.0;
411 			st = 0;
412 		}
413 		else if( *pt1 != '-' && *pt2 != '-' )
414 		{
415 			if( st == 0 )
416 			{
417 				start1 = pos1; start2 = pos2;
418 				st = 1;
419 			}
420 			score += (double)n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]]; // - offset �Ϥ���ʤ�����
421 //			fprintf( stderr, "%c-%c, score(0) = %f\n", *pt1, *pt2, score );
422 		}
423 		if( *pt1++ != '-' ) pos1++;
424 		if( *pt2++ != '-' ) pos2++;
425 	}
426 	if( *(pt1-1) != '-' && *(pt2-1) != '-'  )
427 	{
428 		if( localhompt->nokori++ > 0 )
429 		{
430 //			fprintf( stderr, "reallocating ...\n" );
431 			tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
432 //			fprintf( stderr, "done\n" );
433 			tmppt = tmppt->next;
434 			tmppt->next = NULL;
435 		}
436 
437 		end1 = pos1 - 1;
438 		end2 = pos2 - 1;
439 		tmppt->start1 = start1;
440 		tmppt->start2 = start2;
441 		tmppt->end1   = end1  ;
442 		tmppt->end2   = end2  ;
443 
444 
445 #if 1
446 		if( divpairscore )
447 		{
448 			tmppt->overlapaa   = end2-start2+1;
449 			tmppt->opt = score / tmppt->overlapaa * 5.8 / 600;
450 		}
451 		else
452 		{
453 			sumscore += score;
454 			sumoverlap += end2-start2+1;
455 		}
456 #else
457 		tmppt->overlapaa   = overlapaa;
458 		tmppt->opt = (double)opt;
459 #endif
460 
461 #if 0
462 		fprintf( stderr, "score (2)= %f\n", score );
463 		fprintf( stderr, "al1: %d - %d\n", start1, end1 );
464 		fprintf( stderr, "al2: %d - %d\n", start2, end2 );
465 #endif
466 	}
467 
468 	fprintf( stderr, "sumscore = %f\n", sumscore );
469 	if( !divpairscore )
470 	{
471 
472 		if( !saisho ) subnosento = subnosento->next;
473 		for( tmppt=subnosento; tmppt; tmppt=tmppt->next )
474 		{
475 			tmppt->overlapaa = sumoverlap;
476 			tmppt->opt = sumscore * 5.8 / 600 / sumoverlap;
477 			fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
478 		}
479 	}
480 }
putlocalhom_ext(char * al1,char * al2,LocalHom * localhompt,int off1,int off2,int opt,int overlapaa)481 void putlocalhom_ext( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
482 {
483 	int pos1, pos2, start1, start2, end1, end2;
484 	char *pt1, *pt2;
485 	int iscore;
486 	int isumscore;
487 	int sumoverlap;
488 	LocalHom *tmppt = localhompt;
489 	int nlocalhom = 0;
490 	int st;
491 	pt1 = al1; pt2 = al2;
492 	pos1 = off1; pos2 = off2;
493 
494 
495 	isumscore = 0;
496 	sumoverlap = 0;
497 	start1 = 0; // by D.Mathog, a guess
498 	start2 = 0; // by D.Mathog, a guess
499 
500 	st = 0;
501 	iscore = 0;
502 	while( *pt1 != 0 )
503 	{
504 //		fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
505 		if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
506 		{
507 			end1 = pos1 - 1;
508 			end2 = pos2 - 1;
509 
510 			if( nlocalhom++ > 0 )
511 			{
512 //				fprintf( stderr, "reallocating ...\n" );
513 				tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
514 //				fprintf( stderr, "done\n" );
515 				tmppt = tmppt->next;
516 				tmppt->next = NULL;
517 			}
518 			tmppt->start1 = start1;
519 			tmppt->start2 = start2;
520 			tmppt->end1   = end1  ;
521 			tmppt->end2   = end2  ;
522 
523 #if 1
524 			if( divpairscore )
525 			{
526 				tmppt->overlapaa   = end2-start2+1;
527 				tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
528 			}
529 			else
530 			{
531 				isumscore += iscore;
532 				sumoverlap += end2-start2+1;
533 			}
534 #else
535 			tmppt->overlapaa   = overlapaa;
536 			tmppt->opt = (double)opt;
537 #endif
538 
539 #if 0
540 			fprintf( stderr, "iscore (1)= %d\n", iscore );
541 			fprintf( stderr, "al1: %d - %d\n", start1, end1 );
542 			fprintf( stderr, "al2: %d - %d\n", start2, end2 );
543 #endif
544 			iscore = 0;
545 			st = 0;
546 		}
547 		else if( *pt1 != '-' && *pt2 != '-' )
548 		{
549 			if( st == 0 )
550 			{
551 				start1 = pos1; start2 = pos2;
552 				st = 1;
553 			}
554 			iscore += n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]]; // - offset �Ϥ���ʤ�����
555 //			fprintf( stderr, "%c-%c, iscore(0) = %d\n", *pt1, *pt2, iscore );
556 		}
557 		if( *pt1++ != '-' ) pos1++;
558 		if( *pt2++ != '-' ) pos2++;
559 	}
560 	if( *(pt1-1) != '-' && *(pt2-1) != '-'  )
561 	{
562 		if( nlocalhom++ > 0 )
563 		{
564 //			fprintf( stderr, "reallocating ...\n" );
565 			tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
566 //			fprintf( stderr, "done\n" );
567 			tmppt = tmppt->next;
568 			tmppt->next = NULL;
569 		}
570 		end1 = pos1 - 1;
571 		end2 = pos2 - 1;
572 		tmppt->start1 = start1;
573 		tmppt->start2 = start2;
574 		tmppt->end1   = end1  ;
575 		tmppt->end2   = end2  ;
576 
577 #if 1
578 		if( divpairscore )
579 		{
580 			tmppt->overlapaa   = end2-start2+1;
581 			tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
582 		}
583 		else
584 		{
585 			isumscore += iscore;
586 			sumoverlap += end2-start2+1;
587 		}
588 #else
589 		tmppt->overlapaa   = overlapaa;
590 		tmppt->opt = (double)opt;
591 #endif
592 
593 #if 0
594 		fprintf( stderr, "iscore (2)= %d\n", iscore );
595 		fprintf( stderr, "al1: %d - %d\n", start1, end1 );
596 		fprintf( stderr, "al2: %d - %d\n", start2, end2 );
597 #endif
598 	}
599 
600 	if( !divpairscore )
601 	{
602 		for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
603 		{
604 			tmppt->overlapaa = sumoverlap;
605 //			tmppt->opt = (double)isumscore * 5.8 / ( 600 * sumoverlap );
606 			tmppt->opt = (double)600 * 5.8 / 600;
607 //			fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
608 		}
609 	}
610 }
611 
putlocalhom_str(char * al1,char * al2,double * equiv,double scale,LocalHom * localhompt,int off1,int off2,int opt,int overlapaa)612 void putlocalhom_str( char *al1, char *al2, double *equiv, double scale, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
613 {
614 	int posinaln, pos1, pos2, start1, start2, end1, end2;
615 	char *pt1, *pt2;
616 	int isumscore;
617 	int sumoverlap;
618 	LocalHom *tmppt = localhompt;
619 	int nlocalhom = 0;
620 //	int st;
621 	pt1 = al1; pt2 = al2;
622 	pos1 = off1; pos2 = off2;
623 
624 	isumscore = 0;
625 	sumoverlap = 0;
626 	start1 = 0; // by D.Mathog, a guess
627 	start2 = 0; // by D.Mathog, a guess
628 
629 	posinaln = 0;
630 	while( *pt1 != 0 )
631 	{
632 		if( *pt1 != '-' && *pt2 != '-' && equiv[posinaln] > 0.0 )
633 		{
634 			start1 = end1 = pos1; start2 = end2 = pos2;
635 			if( nlocalhom++ > 0 )
636 			{
637 //				fprintf( stderr, "reallocating ... (posinaln=%d)\n", posinaln );
638 				tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
639 //				fprintf( stderr, "done\n" );
640 				tmppt = tmppt->next;
641 				tmppt->next = NULL;
642 			}
643 			tmppt->start1 = start1;
644 			tmppt->start2 = start2;
645 			tmppt->end1   = end1  ;
646 			tmppt->end2   = end2  ;
647 
648 			tmppt->overlapaa   = 1;
649 //			tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
650 			tmppt->opt = equiv[posinaln] * scale;
651 //			fprintf( stdout, "*pt1=%c, *pt2=%c, equiv=%f\n", *pt1, *pt2, equiv[posinaln] );
652 
653 		}
654 		if( *pt1++ != '-' ) pos1++;
655 		if( *pt2++ != '-' ) pos2++;
656 		posinaln++;
657 	}
658 }
659 
660 
putlocalhom2(char * al1,char * al2,LocalHom * localhompt,int off1,int off2,int opt,int overlapaa)661 void putlocalhom2( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
662 {
663 	int pos1, pos2, start1, start2, end1, end2;
664 	char *pt1, *pt2;
665 	int iscore;
666 	int isumscore;
667 	int sumoverlap;
668 	LocalHom *tmppt = localhompt;
669 	int nlocalhom = 0;
670 	int st;
671 	pt1 = al1; pt2 = al2;
672 	pos1 = off1; pos2 = off2;
673 
674 
675 	isumscore = 0;
676 	sumoverlap = 0;
677 	start1 = 0; // by D.Mathog, a guess
678 	start2 = 0; // by D.Mathog, a guess
679 
680 	st = 0;
681 	iscore = 0;
682 	while( *pt1 != 0 )
683 	{
684 //		fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
685 		if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
686 		{
687 			end1 = pos1 - 1;
688 			end2 = pos2 - 1;
689 
690 			if( nlocalhom++ > 0 )
691 			{
692 //				fprintf( stderr, "reallocating ...\n" );
693 				tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
694 //				fprintf( stderr, "done\n" );
695 				tmppt = tmppt->next;
696 				tmppt->next = NULL;
697 			}
698 			tmppt->start1 = start1;
699 			tmppt->start2 = start2;
700 			tmppt->end1   = end1  ;
701 			tmppt->end2   = end2  ;
702 
703 #if 1
704 			if( divpairscore )
705 			{
706 				tmppt->overlapaa   = end2-start2+1;
707 				tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
708 			}
709 			else
710 			{
711 				isumscore += iscore;
712 				sumoverlap += end2-start2+1;
713 			}
714 #else
715 			tmppt->overlapaa   = overlapaa;
716 			tmppt->opt = (double)opt;
717 #endif
718 
719 #if 0
720 			fprintf( stderr, "iscore (1)= %d\n", iscore );
721 			fprintf( stderr, "al1: %d - %d\n", start1, end1 );
722 			fprintf( stderr, "al2: %d - %d\n", start2, end2 );
723 #endif
724 			iscore = 0;
725 			st = 0;
726 		}
727 		else if( *pt1 != '-' && *pt2 != '-' )
728 		{
729 			if( st == 0 )
730 			{
731 				start1 = pos1; start2 = pos2;
732 				st = 1;
733 			}
734 			iscore += n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]]; // - offset �Ϥ���ʤ�����
735 //			fprintf( stderr, "%c-%c, iscore(0) = %d\n", *pt1, *pt2, iscore );
736 		}
737 		if( *pt1++ != '-' ) pos1++;
738 		if( *pt2++ != '-' ) pos2++;
739 	}
740 	if( *(pt1-1) != '-' && *(pt2-1) != '-'  )
741 	{
742 		if( nlocalhom++ > 0 )
743 		{
744 //			fprintf( stderr, "reallocating ...\n" );
745 			tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
746 //			fprintf( stderr, "done\n" );
747 			tmppt = tmppt->next;
748 			tmppt->next = NULL;
749 		}
750 		end1 = pos1 - 1;
751 		end2 = pos2 - 1;
752 		tmppt->start1 = start1;
753 		tmppt->start2 = start2;
754 		tmppt->end1   = end1  ;
755 		tmppt->end2   = end2  ;
756 
757 #if 1
758 		if( divpairscore )
759 		{
760 			tmppt->overlapaa   = end2-start2+1;
761 			tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
762 		}
763 		else
764 		{
765 			isumscore += iscore;
766 			sumoverlap += end2-start2+1;
767 		}
768 #else
769 		tmppt->overlapaa   = overlapaa;
770 		tmppt->opt = (double)opt;
771 #endif
772 
773 #if 0
774 		fprintf( stderr, "iscore (2)= %d\n", iscore );
775 		fprintf( stderr, "al1: %d - %d\n", start1, end1 );
776 		fprintf( stderr, "al2: %d - %d\n", start2, end2 );
777 #endif
778 	}
779 
780 	if( !divpairscore )
781 	{
782 		for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
783 		{
784 			tmppt->overlapaa = sumoverlap;
785 			tmppt->opt = (double)isumscore * 5.8 / ( 600 * sumoverlap );
786 //			fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
787 		}
788 	}
789 }
putlocalhom(char * al1,char * al2,LocalHom * localhompt,int off1,int off2,int opt,int overlapaa)790 void putlocalhom( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
791 {
792 	int pos1, pos2, start1, start2, end1, end2;
793 	char *pt1, *pt2;
794 	double score;
795 	double sumscore;
796 	int sumoverlap;
797 	LocalHom *tmppt = localhompt;
798 	int nlocalhom = 0;
799 	int st;
800 	pt1 = al1; pt2 = al2;
801 	pos1 = off1; pos2 = off2;
802 
803 
804 	sumscore = 0.0;
805 	sumoverlap = 0;
806 	start1 = 0; // by D.Mathog, a guess
807 	start2 = 0; // by D.Mathog, a guess
808 
809 	st = 0;
810 	score = 0.0;
811 	while( *pt1 != 0 )
812 	{
813 //		fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
814 		if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
815 		{
816 			end1 = pos1 - 1;
817 			end2 = pos2 - 1;
818 
819 			if( nlocalhom++ > 0 )
820 			{
821 //				fprintf( stderr, "reallocating ...\n" );
822 				tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
823 //				fprintf( stderr, "done\n" );
824 				tmppt = tmppt->next;
825 				tmppt->next = NULL;
826 			}
827 			tmppt->start1 = start1;
828 			tmppt->start2 = start2;
829 			tmppt->end1   = end1  ;
830 			tmppt->end2   = end2  ;
831 
832 #if 1
833 			if( divpairscore )
834 			{
835 				tmppt->overlapaa   = end2-start2+1;
836 				tmppt->opt = score / tmppt->overlapaa * 5.8 / 600;
837 			}
838 			else
839 			{
840 				sumscore += score;
841 				sumoverlap += end2-start2+1;
842 			}
843 #else
844 			tmppt->overlapaa   = overlapaa;
845 			tmppt->opt = (double)opt;
846 #endif
847 
848 #if 0
849 			fprintf( stderr, "score (1)= %f\n", score );
850 			fprintf( stderr, "al1: %d - %d\n", start1, end1 );
851 			fprintf( stderr, "al2: %d - %d\n", start2, end2 );
852 #endif
853 			score = 0.0;
854 			st = 0;
855 		}
856 		else if( *pt1 != '-' && *pt2 != '-' )
857 		{
858 			if( st == 0 )
859 			{
860 				start1 = pos1; start2 = pos2;
861 				st = 1;
862 			}
863 			score += (double)n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]]; // - offset �Ϥ���ʤ�����
864 //			fprintf( stderr, "%c-%c, score(0) = %f\n", *pt1, *pt2, score );
865 		}
866 		if( *pt1++ != '-' ) pos1++;
867 		if( *pt2++ != '-' ) pos2++;
868 	}
869 	if( nlocalhom++ > 0 )
870 	{
871 //		fprintf( stderr, "reallocating ...\n" );
872 		tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
873 //		fprintf( stderr, "done\n" );
874 		tmppt = tmppt->next;
875 		tmppt->next = NULL;
876 	}
877 	end1 = pos1 - 1;
878 	end2 = pos2 - 1;
879 	tmppt->start1 = start1;
880 	tmppt->start2 = start2;
881 	tmppt->end1   = end1  ;
882 	tmppt->end2   = end2  ;
883 
884 #if 1
885 	if( divpairscore )
886 	{
887 		tmppt->overlapaa   = end2-start2+1;
888 		tmppt->opt = score / tmppt->overlapaa * 5.8 / 600;
889 	}
890 	else
891 	{
892 		sumscore += score;
893 		sumoverlap += end2-start2+1;
894 	}
895 #else
896 	tmppt->overlapaa   = overlapaa;
897 	tmppt->opt = (double)opt;
898 #endif
899 
900 #if 0
901 	fprintf( stderr, "score (2)= %f\n", score );
902 	fprintf( stderr, "al1: %d - %d\n", start1, end1 );
903 	fprintf( stderr, "al2: %d - %d\n", start2, end2 );
904 #endif
905 
906 	if( !divpairscore )
907 	{
908 		for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
909 		{
910 			tmppt->overlapaa = sumoverlap;
911 			tmppt->opt = sumscore * 5.8 / 600 / sumoverlap;
912 //			fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
913 		}
914 	}
915 }
916 
cutal(char * al,int al_display_start,int start,int end)917 char *cutal( char *al, int al_display_start, int start, int end )
918 {
919 	int pos;
920 	char *pt = al;
921 	char *val = NULL;
922 
923 	pos = al_display_start;
924 	do
925 	{
926 		if( start == pos ) val = pt;
927 		if( end == pos ) break;
928 //		fprintf( stderr, "pos=%d, *pt=%c, val=%p\n", pos, *pt, val );
929 		if( *pt != '-' ) pos++;
930 	} while( *pt++ != 0 );
931 	*(pt+1) = 0;
932 	return( val );
933 }
934 
ErrorExit(char * message)935 void ErrorExit( char *message )
936 {
937 	fprintf( stderr, "%s\n", message );
938 	exit( 1 );
939 }
940 
strncpy_caseC(char * str1,char * str2,int len)941 void strncpy_caseC( char *str1, char *str2, int len )
942 {
943 	if( dorp == 'd' && upperCase > 0 )
944 	{
945 		while( len-- )
946 			*str1++ = toupper( *str2++ );
947 	}
948 	else strncpy( str1, str2, len );
949 }
950 
seqUpper(int nseq,char ** seq)951 void seqUpper( int nseq, char **seq )
952 {
953 	int i, j, len;
954 	for( i=0; i<nseq; i++ )
955 	{
956 		len = strlen( seq[i] );
957 		for( j=0; j<len; j++ )
958 			seq[i][j] = toupper( seq[i][j] );
959 	}
960 }
961 
seqLower(int nseq,char ** seq)962 void seqLower( int nseq, char **seq )
963 {
964 	int i, j, len;
965 	for( i=0; i<nseq; i++ )
966 	{
967 		len = strlen( seq[i] );
968 		for( j=0; j<len; j++ )
969 			seq[i][j] = tolower( seq[i][j] );
970 	}
971 }
972 
getaline_fp_eof(char * s,int l,FILE * fp)973 int getaline_fp_eof( char *s, int l, FILE *fp )  /* end of file -> return 1 */
974 {
975     int c, i = 0 ;
976     int noteofflag = 0;
977     for( i=0; i<l && ( noteofflag = ( (c=getc(fp)) != EOF ) ) && c != '\n'; i++ )
978     	*s++ = c;
979     *s = '\0' ;
980      return( !noteofflag );
981 }
982 
getaline_fp_eof_new(s,l,fp)983 int getaline_fp_eof_new(s, l, fp)  /* end of file -> return 1 */
984 char    s[] ; int l ; FILE *fp ;
985 {
986         int c = 0, i = 0 ;
987 		int noteofflag = 0;
988 
989 		if( feof( fp ) ) return( 1 );
990 
991 		for( i=0; i<l && ( noteofflag = ( (c=getc(fp)) != EOF ) ) && c != '\n'; i++ )
992         { *s++ = c ; }
993         *s = '\0' ;
994 		if( c != '\n' && c != EOF ) while( getc(fp) != '\n' )
995 			;
996 		return( !noteofflag );
997 }
998 
myfgets(s,l,fp)999 int myfgets(s, l, fp)  /* l�ʾ�ϡ������ޤ��ɤ����Ф� */
1000 char    s[] ; int l ; FILE *fp ;
1001 {
1002         int     c = 0, i = 0 ;
1003 
1004 		if( feof( fp ) ) return( 1 );
1005 
1006 		for( i=0; i<l && ( c=getc( fp ) ) != '\n'; i++ )
1007         	*s++ = c;
1008         *s = '\0' ;
1009 		if( c != '\n' )
1010 			while( getc(fp) != '\n' )
1011 				;
1012 		return( 0 );
1013 }
1014 
input_new(FILE * fp,int d)1015 double input_new( FILE *fp, int d )
1016 {
1017 	char mojiretsu[10];
1018 	int i, c;
1019 
1020 	c = getc( fp );
1021 	if( c != '\n' )
1022 		ungetc( c, fp );
1023 
1024 	for( i=0; i<d; i++ )
1025 		mojiretsu[i] = getc( fp );
1026 	mojiretsu[i] = 0;
1027 
1028 	return( atof( mojiretsu ) );
1029 }
1030 
1031 
PreRead(FILE * fp,int * locnjob,int * locnlenmax)1032 void PreRead( FILE *fp, int *locnjob, int *locnlenmax )
1033 {
1034 	int i, nleni;
1035 	char b[B];
1036 
1037 	fgets( b, B-1, fp ); *locnjob = atoi( b );
1038 	*locnlenmax = 0;
1039 	i=0;
1040 	while( i<*locnjob )
1041 	{
1042 		fgets( b, B-1, fp );
1043 		if( !strncmp( b, "=", 1 ) )
1044 		{
1045 			i++;
1046 			fgets( b, B-1, fp ); nleni = atoi( b );
1047 			if( nleni > *locnlenmax ) *locnlenmax = nleni;
1048 		}
1049 	}
1050 	if( *locnlenmax > N )
1051 	{
1052 		fprintf( stderr, "TOO LONG SEQUENCE!\n" );
1053 		exit( 1 );
1054 	}
1055 	if( njob > M  )
1056 	{
1057 		fprintf( stderr, "TOO MANY SEQUENCE!\n" );
1058 		fprintf( stderr, "%d > %d\n", njob, M );
1059 		exit( 1 );
1060 	}
1061 }
1062 
allSpace(char * str)1063 int allSpace( char *str )
1064 {
1065 	int value = 1;
1066 	while( *str ) value *= ( !isdigit( *str++ ) );
1067 	return( value );
1068 }
1069 
Read(char name[M][B],int nlen[M],char ** seq)1070 void Read( char name[M][B], int nlen[M], char **seq )
1071 {
1072 	extern void FRead( FILE *x, char y[M][B], int z[M], char **w );
1073 	FRead( stdin, name, nlen, seq );
1074 }
1075 
1076 
FRead(FILE * fp,char name[][B],int nlen[],char ** seq)1077 void FRead( FILE *fp, char name[][B], int nlen[], char **seq )
1078 {
1079 	int i, j;
1080 	char b[B];
1081 
1082 	fgets( b, B-1, fp );
1083 #if DEBUG
1084 	fprintf( stderr, "b = %s\n", b );
1085 #endif
1086 
1087     if( strstr( b, "onnet" ) ) scoremtx = 1;
1088     else if( strstr( b, "DnA" ) )
1089 	{
1090 		scoremtx = -1;
1091 		upperCase = -1;
1092 	}
1093     else if( strstr( b, "dna" ) )
1094 	{
1095 		scoremtx = -1;
1096 		upperCase = 0;
1097 	}
1098 	else if( strstr( b, "DNA" ) )
1099 	{
1100 		scoremtx = -1;
1101 		upperCase = 1;
1102 	}
1103     else if( strstr( b, "M-Y" ) || strstr( b, "iyata" ) ) scoremtx = 2;
1104     else scoremtx = 0;
1105 #if DEBUG
1106 	fprintf( stderr, " %s->scoremtx = %d\n", b, scoremtx );
1107 #endif
1108 
1109 	geta2 = GETA2;
1110 
1111 #if 0
1112 	if( strlen( b ) >=25 )
1113 	{
1114 		b[25] = 0;
1115 	#if DEBUG
1116 		fprintf( stderr, "kimuraR = %s\n", b+20 );
1117 	#endif
1118 		kimuraR = atoi( b+20 );
1119 
1120 		if( kimuraR < 0 || 20 < kimuraR ) ErrorExit( "Illeagal kimuraR value.\n" );
1121 		if( allSpace( b+20 ) ) kimuraR = NOTSPECIFIED;
1122 	}
1123 	else kimuraR = NOTSPECIFIED;
1124 	#if DEBUG
1125 	fprintf( stderr, "kimuraR = %d\n", kimuraR );
1126 	#endif
1127 
1128 	if( strlen( b ) >=20 )
1129 	{
1130 		b[20] = 0;
1131 	#if DEBUG
1132 		fprintf( stderr, "pamN = %s\n", b+15 );
1133 	#endif
1134 		pamN = atoi( b+15 );
1135 		if( pamN < 0 || 400 < pamN ) ErrorExit( "Illeagal pam value.\n" );
1136 		if( allSpace( b+15 ) ) pamN = NOTSPECIFIED;
1137 	}
1138 	else pamN = NOTSPECIFIED;
1139 
1140 	if( strlen( b ) >= 15 )
1141 	{
1142 		b[15] = 0;
1143 	#if DEBUG
1144 		fprintf( stderr, "poffset = %s\n", b+10 );
1145 	#endif
1146 		poffset = atoi( b+10 );
1147 		if( poffset > 500 ) ErrorExit( "Illegal extending gap ppenalty\n" );
1148 		if( allSpace( b+10 ) ) poffset = NOTSPECIFIED;
1149 	}
1150 	else poffset = NOTSPECIFIED;
1151 
1152 	if( strlen( b ) >= 10 )
1153 	{
1154 		b[10] = 0;
1155 	#if DEBUG
1156 		fprintf( stderr, "ppenalty = %s\n", b+5 );
1157 	#endif
1158 		ppenalty = atoi( b+5 );
1159 		if( ppenalty > 0 ) ErrorExit( "Illegal opening gap ppenalty\n" );
1160 		if( allSpace( b+5 ) ) ppenalty = NOTSPECIFIED;
1161 	}
1162 	else ppenalty = NOTSPECIFIED;
1163 #endif
1164 
1165 	for( i=0; i<njob; i++ )
1166 	{
1167 		getaline_fp_eof_new( b, B-1, fp );
1168 		strcpy( name[i], b );
1169 #if DEBUG
1170 		fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1171 #endif
1172 		fgets( b, B-1, fp ); nlen[i] = atoi( b );      /* seq i no nagasa */
1173 		seq[i][0] = 0;
1174 		if( nlen[i] ) for( j=0; j <= (nlen[i]-1)/C; j++ )
1175 		{
1176 			getaline_fp_eof_new( b, B-1, fp );
1177 			/*	b[C] = 0;  */
1178 			strcat( seq[i], b );
1179 		}
1180 		seq[i][nlen[i]] = 0;
1181 	}
1182 	if( scoremtx == -1 && upperCase != -1 ) seqLower( njob, seq );
1183 }
1184 
1185 
countKUorWA(FILE * fp)1186 static int countKUorWA( FILE *fp )
1187 {
1188 	int value;
1189 	int c, b;
1190 
1191 	value= 0;
1192 	b = '\n';
1193 	while( ( c = getc( fp ) ) != EOF )
1194 	{
1195 		if( b == '\n' && ( c == '>' ) )
1196 			value++;
1197 		b = c;
1198 	}
1199 	rewind( fp );
1200 	return( value );
1201 }
1202 
searchKUorWA(FILE * fp)1203 void searchKUorWA( FILE *fp )
1204 {
1205 	int c, b;
1206 	b = '\n';
1207 	while( !( ( ( c = getc( fp ) ) == '>' || c == EOF ) && b == '\n' ) )
1208 		b = c;
1209 	ungetc( c, fp );
1210 }
1211 
onlyGraph(char * str)1212 static int onlyGraph( char *str )
1213 {
1214 	char tmp;
1215 	char *res = str;
1216 	char *bk = str;
1217 
1218 //	while( (tmp=*str++) ) if( isgraph( tmp ) ) *res++ = tmp;
1219 	while( (tmp=*str++) )
1220 	{
1221 		if( 0x20 < tmp && tmp < 0x7f ) *res++ = tmp;
1222 		if( tmp == '>' || tmp == '(' )
1223 		{
1224 			fprintf( stderr, "========================================================\n" );
1225 			fprintf( stderr, "========================================================\n" );
1226 			fprintf( stderr, "=== \n" );
1227 			fprintf( stderr, "=== ERROR!! \n" );
1228 //			fprintf( stderr, "=== In the '--anysymbol' and '--preservecase' modes, \n" );
1229 			fprintf( stderr, "=== '>' and '(' are acceptable only in title lines.\n" );
1230 			fprintf( stderr, "=== \n" );
1231 			fprintf( stderr, "========================================================\n" );
1232 			fprintf( stderr, "========================================================\n" );
1233 			exit( 1 );
1234 		}
1235 	}
1236 	*res = 0;
1237 	return( res - bk );
1238 }
1239 
charfilter(char * str)1240 static int charfilter( char *str )
1241 {
1242 	char tmp;
1243 	char *res = str;
1244 	char *bk = str;
1245 
1246 	while( (tmp=*str++) )
1247 	{
1248 		if( tmp == '=' || tmp == '*' || tmp == '<' || tmp == '>' || tmp == '(' || tmp == ')' )
1249 		{
1250 			fprintf( stderr, "\n" );
1251 			fprintf( stderr, "Characters '= * < > ( )' are not accepted in the --text mode, \nalthough most printable characters are ok.\n" );
1252 			fprintf( stderr, "\n" );
1253 			exit( 1 );
1254 		}
1255 		if( 0x20 < tmp && tmp < 0x7f )
1256 //		if( tmp != '\n' && tmp != ' ' && tmp != '\t' ) // unprintable characters mo ok.
1257 			*res++ = tmp;
1258 	}
1259 	*res = 0;
1260 	return( res - bk );
1261 }
onlyAlpha_lower(char * str)1262 static int onlyAlpha_lower( char *str )
1263 {
1264 	char tmp;
1265 	char *res = str;
1266 	char *bk = str;
1267 
1268 	while( (tmp=*str++) )
1269 		if( isalpha( tmp ) || tmp == '-' || tmp == '*' || tmp == '.' )
1270 			*res++ = tolower( tmp );
1271 	*res = 0;
1272 	return( res - bk );
1273 }
onlyAlpha_upper(char * str)1274 static int onlyAlpha_upper( char *str )
1275 {
1276 	char tmp;
1277 	char *res = str;
1278 	char *bk = str;
1279 
1280 	while( (tmp=*str++) )
1281 		if( isalpha( tmp ) || tmp == '-' || tmp == '*' || tmp == '.' )
1282 			*res++ = toupper( tmp );
1283 	*res = 0;
1284 	return( res - bk );
1285 }
1286 
kake2hiku(char * str)1287 void kake2hiku( char *str )
1288 {
1289 	do
1290 		if( *str == '*' ) *str = '-';
1291 	while( *str++ );
1292 }
1293 
load1SeqWithoutName_realloc_casepreserve(FILE * fpp)1294 char *load1SeqWithoutName_realloc_casepreserve( FILE *fpp )
1295 {
1296 	int c, b;
1297 	char *cbuf;
1298 	int size = N;
1299 	char *val;
1300 
1301 	val = malloc( (size+1) * sizeof( char ) );
1302 	cbuf = val;
1303 
1304 	b = '\n';
1305 	while( ( c = getc( fpp ) ) != EOF &&
1306           !( ( c == '>' || c == '(' || c == EOF ) && b == '\n' ) )
1307 	{
1308 		*cbuf++ = (char)c;  /* Ĺ�����Ƥ⤷��ʤ� */
1309 		if( cbuf - val == size )
1310 		{
1311 			size += N;
1312 			fprintf( stderr, "reallocating...\n" );
1313 			val = (char *)realloc( val, (size+1) * sizeof( char ) );
1314 			if( !val )
1315 			{
1316 				fprintf( stderr, "Allocation error in load1SeqWithoutName_realloc \n" );
1317 				exit( 1 );
1318 			}
1319 			fprintf( stderr, "done.\n" );
1320 			cbuf = val + size-N;
1321 		}
1322 		b = c;
1323 	}
1324 	ungetc( c, fpp );
1325 	*cbuf = 0;
1326 	onlyGraph( val );
1327 //	kake2hiku( val );
1328 	return( val );
1329 }
1330 
load1SeqWithoutName_realloc(FILE * fpp)1331 char *load1SeqWithoutName_realloc( FILE *fpp )
1332 {
1333 	int c, b;
1334 	char *cbuf;
1335 	int size = N;
1336 	char *val;
1337 
1338 	val = malloc( (size+1) * sizeof( char ) );
1339 	cbuf = val;
1340 
1341 	b = '\n';
1342 	while( ( c = getc( fpp ) ) != EOF &&
1343           !( ( c == '>' || c == '(' || c == EOF ) && b == '\n' ) )
1344 	{
1345 		*cbuf++ = (char)c;  /* Ĺ�����Ƥ⤷��ʤ� */
1346 		if( cbuf - val == size )
1347 		{
1348 			size += N;
1349 			fprintf( stderr, "reallocating...\n" );
1350 			val = (char *)realloc( val, (size+1) * sizeof( char ) );
1351 			if( !val )
1352 			{
1353 				fprintf( stderr, "Allocation error in load1SeqWithoutName_realloc \n" );
1354 				exit( 1 );
1355 			}
1356 			fprintf( stderr, "done.\n" );
1357 			cbuf = val + size-N;
1358 		}
1359 		b = c;
1360 	}
1361 	ungetc( c, fpp );
1362 	*cbuf = 0;
1363 
1364 	if( nblosum == -2 )
1365 	{
1366 		charfilter( val );
1367 	}
1368 	else
1369 	{
1370 		if( dorp == 'd' )
1371 			onlyAlpha_lower( val );
1372 		else
1373 			onlyAlpha_upper( val );
1374 		kake2hiku( val );
1375 	}
1376 	return( val );
1377 }
1378 
load1SeqWithoutName_new(FILE * fpp,char * cbuf)1379 int load1SeqWithoutName_new( FILE *fpp, char *cbuf )
1380 {
1381 	int c, b;
1382 	char *bk = cbuf;
1383 
1384 	b = '\n';
1385 	while( ( c = getc( fpp ) ) != EOF &&                    /* by T. Nishiyama */
1386           !( ( c == '>' || c == '(' || c == EOF ) && b == '\n' ) )
1387 	{
1388 		*cbuf++ = (char)c;  /* Ĺ�����Ƥ⤷��ʤ� */
1389 		b = c;
1390 	}
1391 	ungetc( c, fpp );
1392 	*cbuf = 0;
1393 	if( dorp == 'd' )
1394 		onlyAlpha_lower( bk );
1395 	else
1396 		onlyAlpha_upper( bk );
1397 	kake2hiku( bk );
1398 	return( 0 );
1399 }
1400 
1401 
readDataforgaln(FILE * fp,char ** name,int * nlen,char ** seq)1402 void readDataforgaln( FILE *fp, char **name, int *nlen, char **seq )
1403 {
1404 	int i;
1405 	static char *tmpseq = NULL;
1406 
1407 #if 0
1408 	if( !tmpseq )
1409 	{
1410 		tmpseq = AllocateCharVec( N );
1411 	}
1412 #endif
1413 
1414 	rewind( fp );
1415 	searchKUorWA( fp );
1416 
1417 	for( i=0; i<njob; i++ )
1418 	{
1419 		name[i][0] = '='; getc( fp );
1420 #if 0
1421 		fgets( name[i]+1, B-2, fp );
1422 		j = strlen( name[i] );
1423 		if( name[i][j-1] != '\n' )
1424 			ErrorExit( "Too long name\n" );
1425 		name[i][j-1] = 0;
1426 #else
1427 		myfgets( name[i]+1, B-2, fp );
1428 #endif
1429 #if 0
1430 		fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1431 #endif
1432 		tmpseq = load1SeqWithoutName_realloc( fp );
1433 		strcpy( seq[i], tmpseq );
1434 		nlen[i] = strlen( seq[i] );
1435 		free( tmpseq );
1436 	}
1437 	if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1438 #if 0
1439 	free( tmpseq );
1440 #endif
1441 }
1442 
readData_varlen(FILE * fp,char ** name,int * nlen,char ** seq)1443 void readData_varlen( FILE *fp, char **name, int *nlen, char **seq )
1444 {
1445 	int i;
1446 	static char *tmpseq = NULL;
1447 
1448 	rewind( fp );
1449 	searchKUorWA( fp );
1450 
1451 	for( i=0; i<njob; i++ )
1452 	{
1453 		name[i][0] = '='; getc( fp );
1454 #if 0
1455 		fgets( name[i]+1, B-2, fp );
1456 		j = strlen( name[i] );
1457 		if( name[i][j-1] != '\n' )
1458 			ErrorExit( "Too long name\n" );
1459 		name[i][j-1] = 0;
1460 #else
1461 		myfgets( name[i]+1, B-2, fp );
1462 #endif
1463 #if 0
1464 		fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1465 #endif
1466 		tmpseq = load1SeqWithoutName_realloc( fp );
1467 		nlen[i] = strlen( tmpseq );
1468 //		fprintf( stderr, "nlen[%d] = %d\n", i+1, nlen[i] );
1469 		seq[i] = calloc( nlen[i]+1, sizeof( char ) );
1470 		strcpy( seq[i], tmpseq );
1471 		free( tmpseq );
1472 	}
1473 	if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1474 #if 0
1475 	free( tmpseq );
1476 #endif
1477 }
1478 
readData_pointer2(FILE * fp,int nseq,char ** name,int * nlen,char ** seq)1479 void readData_pointer2( FILE *fp, int nseq, char **name, int *nlen, char **seq )
1480 {
1481 	int i;
1482 	static char *tmpseq = NULL;
1483 
1484 #if 0
1485 	if( !tmpseq )
1486 	{
1487 		tmpseq = AllocateCharVec( N );
1488 	}
1489 #endif
1490 
1491 	rewind( fp );
1492 	searchKUorWA( fp );
1493 
1494 	for( i=0; i<nseq; i++ )
1495 	{
1496 		name[i][0] = '='; getc( fp );
1497 #if 0
1498 		fgets( name[i]+1, B-2, fp );
1499 		j = strlen( name[i] );
1500 		if( name[i][j-1] != '\n' )
1501 			ErrorExit( "Too long name\n" );
1502 		name[i][j-1] = 0;
1503 #else
1504 		myfgets( name[i]+1, B-2, fp );
1505 #endif
1506 #if 0
1507 		fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1508 #endif
1509 		tmpseq = load1SeqWithoutName_realloc( fp );
1510 		strcpy( seq[i], tmpseq );
1511 		free( tmpseq );
1512 		nlen[i] = strlen( seq[i] );
1513 	}
1514 	if( dorp == 'd' && upperCase != -1 ) seqLower( nseq, seq );
1515 #if 0
1516 	free( tmpseq );
1517 #endif
1518 	if( outnumber )
1519 	{
1520 		char *namebuf;
1521 		char *cptr;
1522 		namebuf = calloc( B+100, sizeof( char ) );
1523 		for( i=0; i<nseq; i++ )
1524 		{
1525 			namebuf[0] = '=';
1526 			cptr = strstr( name[i], "_numo_e_" );
1527 			if( cptr )
1528 				sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, cptr+8 );
1529 			else
1530 				sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, name[i]+1 );
1531 			strncpy( name[i], namebuf, B );
1532 			name[i][B-1] = 0;
1533 		}
1534 		free( namebuf );
1535 //		exit( 1 );
1536 	}
1537 }
1538 
1539 
readData_pointer_casepreserve(FILE * fp,char ** name,int * nlen,char ** seq)1540 void readData_pointer_casepreserve( FILE *fp, char **name, int *nlen, char **seq )
1541 {
1542 	int i;
1543 	static char *tmpseq = NULL;
1544 
1545 #if 0
1546 	if( !tmpseq )
1547 	{
1548 		tmpseq = AllocateCharVec( N );
1549 	}
1550 #endif
1551 
1552 	rewind( fp );
1553 	searchKUorWA( fp );
1554 
1555 	for( i=0; i<njob; i++ )
1556 	{
1557 		name[i][0] = '='; getc( fp );
1558 #if 0
1559 		fgets( name[i]+1, B-2, fp );
1560 		j = strlen( name[i] );
1561 		if( name[i][j-1] != '\n' )
1562 			ErrorExit( "Too long name\n" );
1563 		name[i][j-1] = 0;
1564 #else
1565 		myfgets( name[i]+1, B-2, fp );
1566 #endif
1567 #if 0
1568 		fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1569 #endif
1570 		tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1571 		strcpy( seq[i], tmpseq );
1572 		free( tmpseq );
1573 		nlen[i] = strlen( seq[i] );
1574 	}
1575 }
1576 
1577 
copydatafromgui(char ** namegui,char ** seqgui,char ** name,int * nlen,char ** seq)1578 int copydatafromgui( char **namegui, char **seqgui, char **name, int *nlen, char **seq )
1579 {
1580 	int i;
1581 
1582 
1583 	for( i=0; i<njob; i++ )
1584 	{
1585 		name[i][0] = '=';
1586 		strncpy( name[i]+1, namegui[i], B-2 );
1587 		name[i][B-1] = 0;
1588 
1589 		strcpy( seq[i], seqgui[i] );
1590 		nlen[i] = strlen( seq[i] );
1591 	}
1592 	if( dorp == 'd' )
1593 		seqLower( njob, seq );
1594 	else if( dorp == 'p' )
1595 		seqUpper( njob, seq );
1596 	else
1597 	{
1598 		reporterr( "DNA or Protein?\n" );
1599 		return( 1 );
1600 	}
1601 #if 0
1602 	free( tmpseq );
1603 #endif
1604 	if( outnumber )
1605 	{
1606 		char *namebuf;
1607 		char *cptr;
1608 		namebuf = calloc( B+100, sizeof( char ) );
1609 		for( i=0; i<njob; i++ )
1610 		{
1611 			namebuf[0] = '=';
1612 			cptr = strstr( name[i], "_numo_e_" );
1613 			if( cptr )
1614 				sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, cptr+8 );
1615 			else
1616 				sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, name[i]+1 );
1617 			strncpy( name[i], namebuf, B );
1618 			name[i][B-1] = 0;
1619 		}
1620 		free( namebuf );
1621 	}
1622 	return( 0 );
1623 }
1624 
readData_pointer(FILE * fp,char ** name,int * nlen,char ** seq)1625 void readData_pointer( FILE *fp, char **name, int *nlen, char **seq )
1626 {
1627 	int i;
1628 	static char *tmpseq = NULL;
1629 
1630 #if 0
1631 	if( !tmpseq )
1632 	{
1633 		tmpseq = AllocateCharVec( N );
1634 	}
1635 #endif
1636 
1637 	rewind( fp );
1638 	searchKUorWA( fp );
1639 
1640 	for( i=0; i<njob; i++ )
1641 	{
1642 		name[i][0] = '='; getc( fp );
1643 #if 0
1644 		fgets( name[i]+1, B-2, fp );
1645 		j = strlen( name[i] );
1646 		if( name[i][j-1] != '\n' )
1647 			ErrorExit( "Too long name\n" );
1648 		name[i][j-1] = 0;
1649 #else
1650 		myfgets( name[i]+1, B-2, fp );
1651 #endif
1652 #if 0
1653 		fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1654 #endif
1655 		tmpseq = load1SeqWithoutName_realloc( fp );
1656 		strcpy( seq[i], tmpseq );
1657 		free( tmpseq );
1658 		nlen[i] = strlen( seq[i] );
1659 	}
1660 	if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1661 #if 0
1662 	free( tmpseq );
1663 #endif
1664 	if( outnumber )
1665 	{
1666 		char *namebuf;
1667 		char *cptr;
1668 		namebuf = calloc( B+100, sizeof( char ) );
1669 		for( i=0; i<njob; i++ )
1670 		{
1671 			namebuf[0] = '=';
1672 			cptr = strstr( name[i], "_numo_e_" );
1673 			if( cptr )
1674 				sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, cptr+8 );
1675 			else
1676 				sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, name[i]+1 );
1677 			strncpy( name[i], namebuf, B );
1678 			name[i][B-1] = 0;
1679 		}
1680 		free( namebuf );
1681 //		exit( 1 );
1682 	}
1683 }
1684 
readData(FILE * fp,char name[][B],int nlen[],char ** seq)1685 void readData( FILE *fp, char name[][B], int nlen[], char **seq )
1686 {
1687 	int i;
1688 	static char *tmpseq = NULL;
1689 
1690 #if 0
1691 	if( !tmpseq )
1692 	{
1693 		tmpseq = AllocateCharVec( N );
1694 	}
1695 #endif
1696 
1697 	rewind( fp );
1698 	searchKUorWA( fp );
1699 
1700 	for( i=0; i<njob; i++ )
1701 	{
1702 		name[i][0] = '='; getc( fp );
1703 #if 0
1704 		fgets( name[i]+1, B-2, fp );
1705 		j = strlen( name[i] );
1706 		if( name[i][j-1] != '\n' )
1707 			ErrorExit( "Too long name\n" );
1708 		name[i][j-1] = 0;
1709 #else
1710 		myfgets( name[i]+1, B-2, fp );
1711 #endif
1712 #if 0
1713 		fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1714 #endif
1715 		tmpseq = load1SeqWithoutName_realloc( fp );
1716 		strcpy( seq[i], tmpseq );
1717 		nlen[i] = strlen( seq[i] );
1718 		free( tmpseq );
1719 	}
1720 	if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1721 #if 0
1722 	free( tmpseq );
1723 #endif
1724 }
1725 
cutAlignment(FILE * fp,int ** regtable,char ** revtable,int * outtable,char ** name,char ** outseq)1726 void cutAlignment( FILE *fp, int **regtable, char **revtable, int *outtable, char **name, char **outseq )
1727 {
1728 	int i, j;
1729 	int outlen;
1730 	static char *tmpseq = NULL;
1731 	static char *dumname = NULL;
1732 	char *fs, *rs;
1733 	int npos, lpos;
1734 	int startpos, endpos, seqlen;
1735 
1736 	if( dumname == NULL )
1737 	{
1738 		dumname = AllocateCharVec( N );
1739 	}
1740 
1741 	rewind( fp );
1742 	searchKUorWA( fp );
1743 
1744 
1745 	npos = 0;
1746 	for( i=0; i<njob; i++ )
1747 	{
1748 		dumname[0] = '>'; getc( fp );
1749 		myfgets( dumname+1, B-1, fp );
1750 		tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1751 
1752 		if( outtable[i] )
1753 		{
1754 //			putc( '>', stdout );
1755 //			puts( dumname+1 );
1756 
1757 
1758 			strncat( name[npos], dumname, B-1 );
1759 			name[npos][B-1] = 0;
1760 
1761 			if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1762 			seqlen = strlen( tmpseq );
1763 			lpos = 0;
1764 			for( j=0; j<5; j++ )
1765 			{
1766 				if( regtable[0][j*2] == -1 && regtable[0][j*2+1] == -1 ) continue;
1767 
1768 				startpos = regtable[0][j*2];
1769 				endpos   = regtable[0][j*2+1];
1770 				if( startpos > endpos )
1771 				{
1772 					endpos   = regtable[0][j*2];
1773 					startpos = regtable[0][j*2+1];
1774 				}
1775 
1776 				if( startpos < 0 ) startpos = 0;
1777 				if( endpos   < 0 ) endpos   = 0;
1778 				if( endpos   >= seqlen ) endpos   = seqlen-1;
1779 				if( startpos >= seqlen ) startpos = seqlen-1;
1780 
1781 //				fprintf( stderr, "startpos = %d, endpos = %d\n", startpos, endpos );
1782 
1783 				outlen = endpos - startpos+1;
1784 				if( revtable[0][j] == 'f' )
1785 				{
1786 //					fprintf( stderr, "regtable[%d][st] = %d\n", i, regtable[0][j*2+0] );
1787 //					fprintf( stderr, "regtable[%d][en] = %d\n", i, regtable[0][j*2+1] );
1788 //					fprintf( stderr, "outlen = %d\n", outlen );
1789 //					fprintf( stdout, "%.*s\n", outlen, tmpseq+regtable[0][j*2] );
1790 					strncpy( outseq[npos] + lpos, tmpseq+startpos, outlen );
1791 					lpos += outlen;
1792 				}
1793 				else
1794 				{
1795 					fs = AllocateCharVec( outlen+1 );
1796 					rs = AllocateCharVec( outlen+1 );
1797 
1798 					fs[outlen] = 0;
1799 					strncpy( fs, tmpseq+startpos, outlen );
1800 					sreverse( rs, fs );
1801 //					fprintf( stdout, "%s\n", rs );
1802 					strncpy( outseq[npos] + lpos, rs, outlen );
1803 					lpos += outlen;
1804 					free( fs );
1805 					free( rs );
1806 				}
1807 				outseq[npos][lpos] = 0;
1808 			}
1809 			npos++;
1810 		}
1811 		free( tmpseq );
1812 	}
1813 }
1814 
cutData(FILE * fp,int ** regtable,char ** revtable,int * outtable)1815 void cutData( FILE *fp, int **regtable, char **revtable, int *outtable )
1816 {
1817 	int i, j;
1818 	int outlen, seqlen, startpos, endpos;
1819 	static char *tmpseq = NULL;
1820 	static char *dumname = NULL;
1821 	char *fs, *rs;
1822 
1823 	if( dumname == NULL )
1824 	{
1825 		dumname = AllocateCharVec( N );
1826 	}
1827 
1828 	rewind( fp );
1829 	searchKUorWA( fp );
1830 
1831 	for( i=0; i<njob; i++ )
1832 	{
1833 		dumname[0] = '='; getc( fp );
1834 		myfgets( dumname+1, B-2, fp );
1835 		tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1836 
1837 		if( outtable[i] )
1838 		{
1839 			gappick_samestring( tmpseq );
1840 			putc( '>', stdout );
1841 			puts( dumname+1 );
1842 
1843 			seqlen = strlen( tmpseq );
1844 
1845 			if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1846 			if( outtable[i] == 2 )
1847 			{
1848 				startpos = 0;
1849 				endpos   = seqlen-1;
1850 				outlen = endpos - startpos + 1;
1851 				fprintf( stdout, "%.*s\n", outlen, tmpseq+startpos );
1852 			}
1853 			else
1854 			{
1855 				for( j=0; j<5; j++ )
1856 				{
1857 					if( regtable[i][j*2] == -1 && regtable[i][j*2+1] == -1 ) continue;
1858 
1859 					startpos = regtable[i][j*2];
1860 					endpos   = regtable[i][j*2+1];
1861 
1862 					if( startpos > endpos )
1863 					{
1864 						endpos   = regtable[i][j*2];
1865 						startpos = regtable[i][j*2+1];
1866 					}
1867 
1868 					if( startpos < 0 ) startpos = 0;
1869 					if( endpos   < 0 ) endpos   = 0;
1870 					if( endpos   >= seqlen ) endpos   = seqlen-1;
1871 					if( startpos >= seqlen ) startpos = seqlen-1;
1872 
1873 					outlen = endpos - startpos + 1;
1874 					if( revtable[i][j] == 'f' )
1875 					{
1876 						fprintf( stderr, "startpos = %d\n", startpos );
1877 						fprintf( stderr, "endpos   = %d\n", endpos );
1878 						fprintf( stderr, "outlen = %d\n", outlen );
1879 						fprintf( stdout, "%.*s\n", outlen, tmpseq+startpos );
1880 					}
1881 					else
1882 					{
1883 						fs = AllocateCharVec( outlen+1 );
1884 						rs = AllocateCharVec( outlen+1 );
1885 
1886 						fs[outlen] = 0;
1887 						strncpy( fs, tmpseq+startpos, outlen );
1888 						sreverse( rs, fs );
1889 						fprintf( stdout, "%s\n", rs );
1890 						free( fs );
1891 						free( rs );
1892 					}
1893 				}
1894 			}
1895 		}
1896 		free( tmpseq );
1897 	}
1898 }
1899 
catData(FILE * fp)1900 void catData( FILE *fp )
1901 {
1902 	int i;
1903 	static char *tmpseq = NULL;
1904 	static char *dumname = NULL;
1905 //	char *cptr;
1906 
1907 	if( dumname == NULL )
1908 	{
1909 		dumname = AllocateCharVec( N );
1910 	}
1911 
1912 	rewind( fp );
1913 	searchKUorWA( fp );
1914 
1915 	for( i=0; i<njob; i++ )
1916 	{
1917 		dumname[0] = '='; getc( fp );
1918 		myfgets( dumname+1, B-2, fp );
1919 		if( outnumber )
1920 		{
1921 			fprintf( stdout, ">_numo_s_%08d_numo_e_", i+1 );
1922 		}
1923 		else
1924 		{
1925 			putc( '>', stdout );
1926 		}
1927 		puts( dumname+1 );
1928 		tmpseq = load1SeqWithoutName_realloc( fp );
1929 		if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1930 		puts( tmpseq );
1931 		free( tmpseq );
1932 	}
1933 }
1934 
countATGC(char * s,int * total)1935 int countATGC( char *s, int *total )
1936 {
1937 	int nATGC;
1938 	int nChar;
1939 	char c;
1940 	nATGC = nChar = 0;
1941 
1942 	if( *s == 0 )
1943 	{
1944 		*total = 0;
1945 		return( 0 );
1946 	}
1947 
1948 	do
1949 	{
1950 		c = tolower( *s );
1951 		if( isalpha( c ) )
1952 		{
1953 			nChar++;
1954 			if( c == 'a' || c == 't' || c == 'g' || c == 'c' || c == 'u' || c == 'n' )
1955 				nATGC++;
1956 		}
1957 	}
1958 	while( *++s );
1959 
1960 	*total = nChar;
1961 	return( nATGC );
1962 }
1963 
countATGCbk(char * s)1964 double countATGCbk( char *s )
1965 {
1966 	int nATGC;
1967 	int nChar;
1968 	char c;
1969 	nATGC = nChar = 0;
1970 
1971 	do
1972 	{
1973 		c = tolower( *s );
1974 		if( isalpha( c ) )
1975 		{
1976 			nChar++;
1977 			if( c == 'a' || c == 't' || c == 'g' || c == 'c' || c == 'u' || c == 'n' )
1978 				nATGC++;
1979 		}
1980 	}
1981 	while( *++s );
1982 	return( (double)nATGC / nChar );
1983 }
1984 
1985 
countnogaplen(char * seq)1986 int countnogaplen( char *seq )
1987 {
1988 	int val = 0;
1989 	while( *seq )
1990 		if( *seq++ != '-' ) val++;
1991 	return( val );
1992 }
1993 
countnormalletters(char * seq,char * ref)1994 int countnormalletters( char *seq, char *ref )
1995 {
1996 	int val = 0;
1997 	while( *seq )
1998 		if( strchr( ref, *seq++ ) ) val++;
1999 	return( val );
2000 }
2001 
getnumlen_casepreserve(FILE * fp,int * nlenminpt)2002 void getnumlen_casepreserve( FILE *fp, int *nlenminpt )
2003 {
2004 	int total;
2005 	int nsite = 0;
2006 	int atgcnum;
2007 	int i, tmp;
2008 	char *tmpseq, *tmpname;
2009 	double atgcfreq;
2010 	tmpname = AllocateCharVec( N );
2011 	njob = countKUorWA( fp );
2012 	searchKUorWA( fp );
2013 	nlenmax = 0;
2014 	*nlenminpt = 99999999;
2015 	atgcnum = 0;
2016 	total = 0;
2017 	for( i=0; i<njob; i++ )
2018 	{
2019 		myfgets( tmpname, N-1, fp );
2020 		tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
2021 		tmp = strlen( tmpseq );
2022 		if( tmp > nlenmax ) nlenmax  = tmp;
2023 		if( tmp < *nlenminpt ) *nlenminpt  = tmp;
2024 		atgcnum += countATGC( tmpseq, &nsite );
2025 		total += nsite;
2026 		free( tmpseq );
2027 	}
2028 	free( tmpname );
2029 	atgcfreq = (double)atgcnum / total;
2030 //	fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2031 	if( dorp == NOTSPECIFIED )
2032 	{
2033 		if( atgcfreq > 0.75 )
2034 		{
2035 			dorp = 'd';
2036 			upperCase = -1;
2037 		}
2038 		else
2039 		{
2040 			dorp = 'p';
2041 			upperCase = 0;
2042 		}
2043 	}
2044 }
2045 
getnumlen_nogap(FILE * fp,int * nlenminpt)2046 void getnumlen_nogap( FILE *fp, int *nlenminpt )
2047 {
2048 	int total;
2049 	int nsite = 0;
2050 	int atgcnum;
2051 	int i, tmp;
2052 	char *tmpseq, *tmpname;
2053 	double atgcfreq;
2054 	tmpname = AllocateCharVec( N );
2055 	njob = countKUorWA( fp );
2056 	searchKUorWA( fp );
2057 	nlenmax = 0;
2058 	*nlenminpt = 99999999;
2059 	atgcnum = 0;
2060 	total = 0;
2061 	for( i=0; i<njob; i++ )
2062 	{
2063 		myfgets( tmpname, N-1, fp );
2064 		tmpseq = load1SeqWithoutName_realloc( fp );
2065 		tmp = countnogaplen( tmpseq );
2066 		if( tmp > nlenmax ) nlenmax  = tmp;
2067 		if( tmp < *nlenminpt ) *nlenminpt  = tmp;
2068 		atgcnum += countATGC( tmpseq, &nsite );
2069 		total += nsite;
2070 		free( tmpseq );
2071 	}
2072 	free( tmpname );
2073 	atgcfreq = (double)atgcnum / total;
2074 	fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2075 	if( dorp == NOTSPECIFIED )
2076 	{
2077 		if( atgcfreq > 0.75 )
2078 		{
2079 			dorp = 'd';
2080 			upperCase = -1;
2081 		}
2082 		else
2083 		{
2084 			dorp = 'p';
2085 			upperCase = 0;
2086 		}
2087 	}
2088 }
2089 
2090 
getnumlen_nogap_outallreg(FILE * fp,int * nlenminpt)2091 void getnumlen_nogap_outallreg( FILE *fp, int *nlenminpt )
2092 {
2093 	int total;
2094 	int nsite = 0;
2095 	int atgcnum;
2096 	int i, tmp;
2097 	char *tmpseq, *tmpname;
2098 	double atgcfreq;
2099 	tmpname = AllocateCharVec( N );
2100 	njob = countKUorWA( fp );
2101 	searchKUorWA( fp );
2102 	nlenmax = 0;
2103 	*nlenminpt = 99999999;
2104 	atgcnum = 0;
2105 	total = 0;
2106 	for( i=0; i<njob; i++ )
2107 	{
2108 		myfgets( tmpname, N-1, fp );
2109 		fprintf( stdout, "%s\n", tmpname );
2110 		tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
2111 		tmp = countnogaplen( tmpseq );
2112 		fprintf( stdout, "%d\n", tmp );
2113 		if( tmp > nlenmax ) nlenmax  = tmp;
2114 		if( tmp < *nlenminpt ) *nlenminpt  = tmp;
2115 		atgcnum += countATGC( tmpseq, &nsite );
2116 		total += nsite;
2117 		free( tmpseq );
2118 	}
2119 	free( tmpname );
2120 	atgcfreq = (double)atgcnum / total;
2121 	fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2122 	if( dorp == NOTSPECIFIED )
2123 	{
2124 		if( atgcfreq > 0.75 )
2125 		{
2126 			dorp = 'd';
2127 			upperCase = -1;
2128 		}
2129 		else
2130 		{
2131 			dorp = 'p';
2132 			upperCase = 0;
2133 		}
2134 	}
2135 }
2136 
escapehtml(char * res,char * ori,int maxlen)2137 static void escapehtml( char *res, char *ori, int maxlen )
2138 {
2139 	char *res0 = res;
2140 	while( *ori )
2141 	{
2142 		if( *ori == '<' )
2143 		{
2144 			strcpy( res, "&lt;" );
2145 			res += 3;
2146 		}
2147 		else if( *ori == '>' )
2148 		{
2149 			strcpy( res, "&gt;" );
2150 			res += 3;
2151 		}
2152 		else if( *ori == '&' )
2153 		{
2154 			strcpy( res, "&amp;" );
2155 			res += 4;
2156 		}
2157 		else if( *ori == '"' )
2158 		{
2159 			strcpy( res, "&quot;" );
2160 			res += 5;
2161 		}
2162 		else if( *ori == ' ' )
2163 		{
2164 			strcpy( res, "&nbsp;" );
2165 			res += 5;
2166 		}
2167 		else
2168 		{
2169 			*res = *ori;
2170 		}
2171 		res++;
2172 		ori++;
2173 
2174 		if( res - res0 -10 > N ) break;
2175 	}
2176 	*res = 0;
2177 }
2178 
getnumlen_nogap_outallreg_web(FILE * fp,FILE * ofp,int * nlenminpt,int * isalignedpt)2179 void getnumlen_nogap_outallreg_web( FILE *fp, FILE *ofp, int *nlenminpt, int *isalignedpt )
2180 {
2181 	int total;
2182 	int nsite = 0;
2183 	int atgcnum;
2184 	int alnlen = 0, alnlen_prev;
2185 	int i, tmp, lennormalchar;
2186 	char *tmpseq, *tmpname, *tmpname2;
2187 	double atgcfreq;
2188 	tmpname = AllocateCharVec( N );
2189 	tmpname2 = AllocateCharVec( N );
2190 	njob = countKUorWA( fp );
2191 	searchKUorWA( fp );
2192 	nlenmax = 0;
2193 	*nlenminpt = 99999999;
2194 	atgcnum = 0;
2195 	total = 0;
2196 	alnlen_prev = -1;
2197 	*isalignedpt = 1;
2198 	for( i=0; i<njob; i++ )
2199 	{
2200 		myfgets( tmpname, N-1, fp );
2201 		tmpname2[0] = tmpname[0];
2202 		escapehtml( tmpname2+1, tmpname+1, N );
2203 //		fprintf( stdout, "%s\n", tmpname );
2204 //		fprintf( stdout, "%s\n", tmpname2 );
2205 //		exit(1);
2206 		tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
2207 		tmp = countnogaplen( tmpseq );
2208 //		fprintf( stdout, "%d\n", tmp );
2209 		if( tmp > nlenmax ) nlenmax  = tmp;
2210 		if( tmp < *nlenminpt ) *nlenminpt  = tmp;
2211 		atgcnum += countATGC( tmpseq, &nsite );
2212 		total += nsite;
2213 
2214 		alnlen = strlen( tmpseq );
2215 //		fprintf( stdout, "##### alnlen, alnlen_prev = %d, %d\n", alnlen, alnlen_prev );
2216 		if( i>0 && alnlen_prev != alnlen ) *isalignedpt = 0;
2217 		alnlen_prev = alnlen;
2218 
2219 		atgcfreq = (double)atgcnum / total;
2220 //		fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2221 //		if( dorp == NOTSPECIFIED ) // you kentou
2222 		{
2223 			if( atgcfreq > 0.75 )
2224 			{
2225 				dorp = 'd';
2226 				upperCase = -1;
2227 			}
2228 			else
2229 			{
2230 				dorp = 'p';
2231 				upperCase = 0;
2232 			}
2233 		}
2234 
2235 		if( dorp == 'd' ) lennormalchar = countnormalletters( tmpseq, "atgcuATGCU" );
2236 		else              lennormalchar = countnormalletters( tmpseq, "ARNDCQEGHILKMFPSTWYVarndcqeghilkmfpstwyv" );
2237 		free( tmpseq );
2238 
2239 		fprintf( ofp, " <label for='s%d'><span id='ss%d'><input type='checkbox' id='s%d' name='s%d' checked></span> <input type='text' class='ll' id='ll%d' style='display:none' size='6' value='%d' readonly='readonly'>%s</label>\n", i, i, i, i, i, lennormalchar, tmpname2 );
2240 		fprintf( ofp, "<span id='t%d-0' style='display:none'>", i );
2241 		fprintf( ofp, " <a href='javascript:void(0)' onclick='ddcycle(this.form,\"t%d\")'>+reg</a>", i );
2242 		fprintf( ofp, " Begin:<input type='text' name='b%d-0' size='8' value='1' class='ie'> End:<input type='text' name='e%d-0' size='8' value='%d' class='ie'>", i, i, tmp );
2243 		if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-0'><input type='checkbox' name='r%d-0' id='r%d-0'>Reverse</label>", i, i, i );
2244 //		fprintf( ofp, "  Sequence Length:<input type='text' name='l%d' size='8' value='%d' readonly='readonly'>", i, tmp );
2245 		fprintf( ofp, "\n</span>" );
2246 		fprintf( ofp, "<span id='t%d-1' style='display:none'>", i );
2247 		fprintf( ofp, "      Begin:<input type='text' name='b%d-1' size='8' value='' class='ie'> End:<input type='text' name='e%d-1' size='8' value='' class='ie'>", i, i );
2248 		if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-1'><input type='checkbox' name='r%d-1' id='r%d-1'>Reverse</label>", i, i, i );
2249 		fprintf( ofp, "\n</span>" );
2250 		fprintf( ofp, "<span id='t%d-2' style='display:none'>", i );
2251 		fprintf( ofp, "      Begin:<input type='text' name='b%d-2' size='8' value='' class='ie'> End:<input type='text' name='e%d-2' size='8' value='' class='ie'>", i, i );
2252 		if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-2'><input type='checkbox' name='r%d-2' id='r%d-2'>Reverse</label>", i, i, i );
2253 		fprintf( ofp, "\n</span>" );
2254 		fprintf( ofp, "<span id='t%d-3' style='display:none'>", i );
2255 		fprintf( ofp, "      Begin:<input type='text' name='b%d-3' size='8' value='' class='ie'> End:<input type='text' name='e%d-3' size='8' value='' class='ie'>", i, i );
2256 		if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-3'><input type='checkbox' name='r%d-3' id='r%d-3'>Reverse</label>", i, i, i );
2257 		fprintf( ofp, "\n</span>" );
2258 		fprintf( ofp, "<span id='t%d-4' style='display:none'>", i );
2259 		fprintf( ofp, "      Begin:<input type='text' name='b%d-4' size='8' value='' class='ie'> End:<input type='text' name='e%d-4' size='8' value='' class='ie'>", i, i );
2260 		if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-4'><input type='checkbox' name='r%d-4' id='r%d-4'>Reverse</label>", i, i, i );
2261 		fprintf( ofp, "\n</span>" );
2262 	}
2263 	free( tmpname );
2264 	free( tmpname2 );
2265 	atgcfreq = (double)atgcnum / total;
2266 	fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2267 //	if( dorp == NOTSPECIFIED ) // you kentou
2268 	{
2269 		if( atgcfreq > 0.75 )
2270 		{
2271 			dorp = 'd';
2272 			upperCase = -1;
2273 		}
2274 		else
2275 		{
2276 			dorp = 'p';
2277 			upperCase = 0;
2278 		}
2279 	}
2280 	fprintf( ofp, "\n" );
2281 	if( *isalignedpt )
2282 	{
2283 		fprintf( ofp, "<span id='tall-0' style='display:none'>" );
2284 		fprintf( ofp, "Cut the alignment\n" );
2285 		fprintf( ofp, " <a href='javascript:void(0)' onclick='ddcycle(this.form,\"tall\")'>+reg</a>" );
2286 		fprintf( ofp, " Begin:<input type='text' name='ball-0' size='8' value='1'> End:<input type='text' name='eall-0' size='8' value='%d'>", alnlen );
2287 		if( dorp == 'd' ) fprintf( ofp, " <label for='rall-0'><input type='checkbox' name='rall-0' id='rall-0'>Reverse</label>" );
2288 		fprintf( ofp, "  Alignment length:<input type='text' name='lall' size='8' value='%d' readonly='readonly'>", alnlen );
2289 		fprintf( ofp, "\n</span>" );
2290 		fprintf( ofp, "<span id='tall-1' style='display:none'>" );
2291 		fprintf( ofp, "      Begin:<input type='text' name='ball-1' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2292 		if( dorp == 'd' ) fprintf( ofp, " <label for='rall-1'><input type='checkbox' name='rall-1' id='rall-1'>Reverse</label>" );
2293 		fprintf( ofp, "\n</span>" );
2294 		fprintf( ofp, "<span id='tall-2' style='display:none'>" );
2295 		fprintf( ofp, "      Begin:<input type='text' name='ball-2' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2296 		if( dorp == 'd' ) fprintf( ofp, " <label for='rall-2'><input type='checkbox' name='rall-2' id='rall-2'>Reverse</label>" );
2297 		fprintf( ofp, "\n</span>" );
2298 		fprintf( ofp, "<span id='tall-3' style='display:none'>" );
2299 		fprintf( ofp, "      Begin:<input type='text' name='ball-3' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2300 		if( dorp == 'd' ) fprintf( ofp, " <label for='rall-3'><input type='checkbox' name='rall-3' id='rall-3'>Reverse</label>" );
2301 		fprintf( ofp, "\n</span>" );
2302 		fprintf( ofp, "<span id='tall-4' style='display:none'>" );
2303 		fprintf( ofp, "      Begin:<input type='text' name='ball-4' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2304 		if( dorp == 'd' ) fprintf( ofp, " <label for='rall-4'><input type='checkbox' name='rall-4' id='rall-4'>Reverse</label>" );
2305 		fprintf( ofp, "\n</span>" );
2306 	}
2307 
2308 }
2309 
getnumlen(FILE * fp)2310 void getnumlen( FILE *fp )
2311 {
2312 	int total;
2313 	int nsite = 0;
2314 	int atgcnum;
2315 	int i, tmp;
2316 	char *tmpseq;
2317 	char *tmpname;
2318 	double atgcfreq;
2319 	tmpname = AllocateCharVec( N );
2320 	njob = countKUorWA( fp );
2321 	searchKUorWA( fp );
2322 	nlenmax = 0;
2323 	atgcnum = 0;
2324 	total = 0;
2325 	for( i=0; i<njob; i++ )
2326 	{
2327 		myfgets( tmpname, N-1, fp );
2328 		tmpseq = load1SeqWithoutName_realloc( fp );
2329 		tmp = strlen( tmpseq );
2330 		if( tmp > nlenmax ) nlenmax  = tmp;
2331 		atgcnum += countATGC( tmpseq, &nsite );
2332 		total += nsite;
2333 //		fprintf( stderr, "##### total = %d\n", total );
2334 		free( tmpseq );
2335 	}
2336 
2337 
2338 	atgcfreq = (double)atgcnum / total;
2339 //	fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2340 	if( dorp == NOTSPECIFIED )
2341 	{
2342 		if( atgcfreq > 0.75 )
2343 		{
2344 			dorp = 'd';
2345 			upperCase = -1;
2346 		}
2347 		else
2348 		{
2349 			dorp = 'p';
2350 			upperCase = 0;
2351 		}
2352 	}
2353 	free( tmpname );
2354 }
2355 
2356 
2357 
WriteGapFill(FILE * fp,int locnjob,char name[][B],int nlen[M],char ** aseq)2358 void WriteGapFill( FILE *fp, int locnjob, char name[][B], int nlen[M], char **aseq )
2359 {
2360 	static char b[N];
2361 	int i, j;
2362 	int nalen[M];
2363 	static char gap[N];
2364 	static char buff[N];
2365 
2366 #if IODEBUG
2367 	fprintf( stderr, "IMAKARA KAKU\n" );
2368 #endif
2369 	nlenmax = 0;
2370 	for( i=0; i<locnjob; i++ )
2371 	{
2372 		int len = strlen( aseq[i] );
2373 		if( nlenmax < len ) nlenmax = len;
2374 	}
2375 
2376 	for( i=0; i<nlenmax; i++ ) gap[i] = '-';
2377 	gap[nlenmax] = 0;
2378 
2379 	fprintf( fp, "%5d", locnjob );
2380 	fprintf( fp, "\n" );
2381 
2382 	for( i=0; i<locnjob; i++ )
2383 	{
2384 		strcpy( buff, aseq[i] );
2385 		strncat( buff, gap, nlenmax-strlen( aseq[i] ) );
2386 		buff[nlenmax] = 0;
2387 		nalen[i] = strlen( buff );
2388 		fprintf( fp, "%s\n", name[i] );
2389 		fprintf( fp, "%5d\n", nalen[i] );
2390 		for( j=0; j<nalen[i]; j=j+C )
2391 		{
2392 			strncpy_caseC( b, buff+j, C ); b[C] = 0;
2393 			fprintf( fp, "%s\n",b );
2394 		}
2395 	}
2396 #if DEBUG
2397 	fprintf( stderr, "nalen[0] = %d\n", nalen[0] );
2398 #endif
2399 #if IODEBUG
2400 	fprintf( stderr, "KAKIOWATTA\n" );
2401 #endif
2402 }
2403 
writeDataforgaln(FILE * fp,int locnjob,char ** name,int * nlen,char ** aseq)2404 void writeDataforgaln( FILE *fp, int locnjob, char **name, int *nlen, char **aseq )
2405 {
2406 	int i, j;
2407 	int nalen;
2408 
2409 	for( i=0; i<locnjob; i++ )
2410 	{
2411 		nalen = strlen( aseq[i] );
2412 		fprintf( fp, ">%s\n", name[i]+1 );
2413 		for( j=0; j<nalen; j=j+C )
2414 		{
2415 #if 0
2416 			strncpy( b, aseq[i]+j, C ); b[C] = 0;
2417 			fprintf( fp, "%s\n",b );
2418 #else
2419 			fprintf( fp, "%.*s\n", C, aseq[i]+j );
2420 #endif
2421 		}
2422 	}
2423 }
2424 
writeData_pointer(FILE * fp,int locnjob,char ** name,int * nlen,char ** aseq)2425 void writeData_pointer( FILE *fp, int locnjob, char **name, int *nlen, char **aseq )
2426 {
2427 	int i, j;
2428 	int nalen;
2429 
2430 	for( i=0; i<locnjob; i++ )
2431 	{
2432 #if DEBUG
2433 		fprintf( stderr, "i = %d in writeData\n", i );
2434 #endif
2435 		nalen = strlen( aseq[i] );
2436 		fprintf( fp, ">%s\n", name[i]+1 );
2437 		for( j=0; j<nalen; j=j+C )
2438 		{
2439 #if 0
2440 			strncpy( b, aseq[i]+j, C ); b[C] = 0;
2441 			fprintf( fp, "%s\n",b );
2442 #else
2443 			fprintf( fp, "%.*s\n", C, aseq[i]+j );
2444 #endif
2445 		}
2446 	}
2447 }
2448 
writeData(FILE * fp,int locnjob,char name[][B],int nlen[],char ** aseq)2449 void writeData( FILE *fp, int locnjob, char name[][B], int nlen[], char **aseq )
2450 {
2451 	int i, j;
2452 	int nalen;
2453 
2454 	for( i=0; i<locnjob; i++ )
2455 	{
2456 #if DEBUG
2457 		fprintf( stderr, "i = %d in writeData\n", i );
2458 #endif
2459 		nalen = strlen( aseq[i] );
2460 		fprintf( fp, ">%s\n", name[i]+1 );
2461 		for( j=0; j<nalen; j=j+C )
2462 		{
2463 #if 0
2464 			strncpy( b, aseq[i]+j, C ); b[C] = 0;
2465 			fprintf( fp, "%s\n",b );
2466 #else
2467 			fprintf( fp, "%.*s\n", C, aseq[i]+j );
2468 #endif
2469 		}
2470 	}
2471 }
2472 
2473 
write1seq(FILE * fp,char * aseq)2474 void write1seq( FILE *fp, char *aseq )
2475 {
2476 	int j;
2477 	int nalen;
2478 
2479 	nalen = strlen( aseq );
2480 	for( j=0; j<nalen; j=j+C )
2481 		fprintf( fp, "%.*s\n", C, aseq+j );
2482 }
2483 
readhat2_doublehalf_part_pointer(FILE * fp,int nseq,int nadd,char ** name,double ** mtx)2484 void readhat2_doublehalf_part_pointer( FILE *fp, int nseq, int nadd, char **name, double **mtx )
2485 {
2486     int i, j, nseq0, norg;
2487     char b[B];
2488 
2489     fgets( b, B, fp );
2490     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 )
2491 	{
2492 		fprintf( stderr, "%d != %d\n", nseq, nseq0 );
2493 		ErrorExit( "hat2 is wrong." );
2494 	}
2495     fgets( b, B, fp );
2496     for( i=0; i<nseq; i++ )
2497     {
2498 #if 0
2499         getaline_fp_eof( b, B, fp );
2500 #else
2501 		myfgets( b, B-2, fp );
2502 #endif
2503 #if 0
2504 		j = MIN( strlen( b+6 ), 10 );
2505         if( strncmp( name[i], b+6 , j ) )
2506 		{
2507 			fprintf( stderr, "Error in hat2\n" );
2508 			fprintf( stderr, "%s != %s\n", b, name[i] );
2509 			exit( 1 );
2510 		}
2511 #endif
2512     }
2513 	norg = nseq-nadd;
2514     for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ )
2515     {
2516         mtx[i][j] = ( input_new( fp, D ) );
2517     }
2518 }
2519 
readhat2_doublehalf_pointer(FILE * fp,int nseq,char ** name,double ** mtx)2520 void readhat2_doublehalf_pointer( FILE *fp, int nseq, char **name, double **mtx )
2521 {
2522     int i, j, nseq0;
2523     char b[B];
2524 
2525     fgets( b, B, fp );
2526     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 )
2527 	{
2528 		fprintf( stderr, "%d != %d\n", nseq, nseq0 );
2529 		ErrorExit( "hat2 is wrong." );
2530 	}
2531     fgets( b, B, fp );
2532     for( i=0; i<nseq; i++ )
2533     {
2534 #if 0
2535         getaline_fp_eof( b, B, fp );
2536 #else
2537 		myfgets( b, B-2, fp );
2538 #endif
2539 #if 0
2540 		j = MIN( strlen( b+6 ), 10 );
2541         if( strncmp( name[i], b+6 , j ) )
2542 		{
2543 			fprintf( stderr, "Error in hat2\n" );
2544 			fprintf( stderr, "%s != %s\n", b, name[i] );
2545 			exit( 1 );
2546 		}
2547 #endif
2548     }
2549     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2550     {
2551         mtx[i][j-i] = ( input_new( fp, D ) );
2552     }
2553 }
readhat2_doublehalf(FILE * fp,int nseq,char name[M][B],double ** mtx)2554 void readhat2_doublehalf( FILE *fp, int nseq, char name[M][B], double **mtx )
2555 {
2556     int i, j, nseq0;
2557     char b[B];
2558 
2559     fgets( b, B, fp );
2560     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2561     fgets( b, B, fp );
2562     for( i=0; i<nseq; i++ )
2563     {
2564 #if 0
2565         getaline_fp_eof( b, B, fp );
2566 #else
2567 		myfgets( b, B-2, fp );
2568 #endif
2569 #if 0
2570 		j = MIN( strlen( b+6 ), 10 );
2571         if( strncmp( name[i], b+6 , j ) )
2572 		{
2573 			fprintf( stderr, "Error in hat2\n" );
2574 			fprintf( stderr, "%s != %s\n", b, name[i] );
2575 			exit( 1 );
2576 		}
2577 #endif
2578     }
2579     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2580     {
2581         mtx[i][j-i] = ( input_new( fp, D ) );
2582     }
2583 }
readhat2_double(FILE * fp,int nseq,char name[M][B],double ** mtx)2584 void readhat2_double( FILE *fp, int nseq, char name[M][B], double **mtx )
2585 {
2586     int i, j, nseq0;
2587     char b[B];
2588 
2589     fgets( b, B, fp );
2590     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2591     fgets( b, B, fp );
2592     for( i=0; i<nseq; i++ )
2593     {
2594 #if 0
2595         getaline_fp_eof( b, B, fp );
2596 #else
2597 		myfgets( b, B-2, fp );
2598 #endif
2599 #if 0
2600 		j = MIN( strlen( b+6 ), 10 );
2601         if( strncmp( name[i], b+6 , j ) )
2602 		{
2603 			fprintf( stderr, "Error in hat2\n" );
2604 			fprintf( stderr, "%s != %s\n", b, name[i] );
2605 			exit( 1 );
2606 		}
2607 #endif
2608     }
2609     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2610     {
2611         mtx[i][j] = ( input_new( fp, D ) );
2612     }
2613 }
readhat2_int(FILE * fp,int nseq,char name[M][B],int ** mtx)2614 void readhat2_int( FILE *fp, int nseq, char name[M][B], int **mtx )
2615 {
2616     int i, j, nseq0;
2617     char b[B];
2618 
2619     fgets( b, B, fp );
2620     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2621     fgets( b, B, fp );
2622     for( i=0; i<nseq; i++ )
2623     {
2624 #if 0
2625         getaline_fp_eof( b, B, fp );
2626 #else
2627 		myfgets( b, B-2, fp );
2628 #endif
2629 #if 0
2630 		j = MIN( strlen( b+6 ), 10 );
2631         if( strncmp( name[i], b+6 , j ) )
2632 		{
2633 			fprintf( stderr, "Error in hat2\n" );
2634 			fprintf( stderr, "%s != %s\n", b, name[i] );
2635 			exit( 1 );
2636 		}
2637 #endif
2638     }
2639     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2640     {
2641         mtx[i][j] = (int)( input_new( fp, D ) * INTMTXSCALE + 0.5 );
2642     }
2643 }
2644 
readhat2_pointer(FILE * fp,int nseq,char ** name,double ** mtx)2645 void readhat2_pointer( FILE *fp, int nseq, char **name, double **mtx )
2646 {
2647     int i, j, nseq0;
2648     char b[B];
2649 
2650     fgets( b, B, fp );
2651     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2652     fgets( b, B, fp );
2653     for( i=0; i<nseq; i++ )
2654     {
2655 #if 0
2656         getaline_fp_eof( b, B, fp );
2657 #else
2658 		myfgets( b, B-2, fp );
2659 #endif
2660 #if 0
2661 		j = MIN( strlen( b+6 ), 10 );
2662         if( strncmp( name[i], b+6 , j ) )
2663 		{
2664 			fprintf( stderr, "Error in hat2\n" );
2665 			fprintf( stderr, "%s != %s\n", b, name[i] );
2666 			exit( 1 );
2667 		}
2668 #endif
2669     }
2670     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2671     {
2672         mtx[i][j] = (double)input_new( fp, D);
2673     }
2674 }
readhat2(FILE * fp,int nseq,char name[M][B],double ** mtx)2675 void readhat2( FILE *fp, int nseq, char name[M][B], double **mtx )
2676 {
2677     int i, j, nseq0;
2678     char b[B];
2679 
2680     fgets( b, B, fp );
2681     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2682     fgets( b, B, fp );
2683     for( i=0; i<nseq; i++ )
2684     {
2685 #if 0
2686         getaline_fp_eof( b, B, fp );
2687 #else
2688 		myfgets( b, B-2, fp );
2689 #endif
2690 #if 0
2691 		j = MIN( strlen( b+6 ), 10 );
2692         if( strncmp( name[i], b+6 , j ) )
2693 		{
2694 			fprintf( stderr, "Error in hat2\n" );
2695 			fprintf( stderr, "%s != %s\n", b, name[i] );
2696 			exit( 1 );
2697 		}
2698 #endif
2699     }
2700     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2701     {
2702         mtx[i][j] = (double)input_new( fp, D);
2703     }
2704 }
2705 
WriteFloatHat2_pointer_halfmtx(FILE * hat2p,int locnjob,char ** name,double ** mtx)2706 void WriteFloatHat2_pointer_halfmtx( FILE *hat2p, int locnjob, char **name, double **mtx )
2707 {
2708 	int i, j, ijsa;
2709 	double max = 0.0;
2710 	for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2711 
2712 	fprintf( hat2p, "%5d\n", 1 );
2713 	fprintf( hat2p, "%5d\n", locnjob );
2714 	fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2715 
2716 	for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2717 	for( i=0; i<locnjob; i++ )
2718 	{
2719 		for( j=i+1; j<njob; j++ )
2720 		{
2721 			fprintf( hat2p, DFORMAT, mtx[i][j-i] );
2722 			ijsa = j-i;
2723 			if( ijsa % 12 == 0 || ijsa == locnjob-i-1 ) fprintf( hat2p, "\n" );
2724 		}
2725 	}
2726 }
2727 
WriteFloatHat2_pointer(FILE * hat2p,int locnjob,char ** name,double ** mtx)2728 void WriteFloatHat2_pointer( FILE *hat2p, int locnjob, char **name, double **mtx )
2729 {
2730 	int i, j;
2731 	double max = 0.0;
2732 	for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2733 
2734 	fprintf( hat2p, "%5d\n", 1 );
2735 	fprintf( hat2p, "%5d\n", locnjob );
2736 	fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2737 
2738 	for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2739 	for( i=0; i<locnjob; i++ )
2740 	{
2741 		for( j=1; j<locnjob-i; j++ )
2742 		{
2743 			fprintf( hat2p, DFORMAT, mtx[i][j] );
2744 			if( j % 12 == 0 || j == locnjob-i-1 ) fprintf( hat2p, "\n" );
2745 		}
2746 	}
2747 }
2748 
WriteFloatHat2(FILE * hat2p,int locnjob,char name[M][B],double ** mtx)2749 void WriteFloatHat2( FILE *hat2p, int locnjob, char name[M][B], double **mtx )
2750 {
2751 	int i, j;
2752 	double max = 0.0;
2753 	for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2754 
2755 	fprintf( hat2p, "%5d\n", 1 );
2756 	fprintf( hat2p, "%5d\n", locnjob );
2757 	fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2758 
2759 	for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2760 	for( i=0; i<locnjob; i++ )
2761 	{
2762 		for( j=1; j<locnjob-i; j++ )
2763 		{
2764 			fprintf( hat2p, DFORMAT, mtx[i][j] );
2765 			if( j % 12 == 0 || j == locnjob-i-1 ) fprintf( hat2p, "\n" );
2766 		}
2767 	}
2768 }
2769 
WriteHat2_int(FILE * hat2p,int locnjob,char name[M][B],int ** mtx)2770 void WriteHat2_int( FILE *hat2p, int locnjob, char name[M][B], int **mtx )
2771 {
2772 	int i, j;
2773 	double max = 0.0;
2774 	for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2775 	max /= INTMTXSCALE;
2776 
2777 	fprintf( hat2p, "%5d\n", 1 );
2778 	fprintf( hat2p, "%5d\n", locnjob );
2779 	fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2780 
2781 	for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2782 	for( i=0; i<locnjob-1; i++ )
2783 	{
2784 		for( j=i+1; j<locnjob; j++ )
2785 		{
2786 			fprintf( hat2p, DFORMAT, (double)mtx[i][j] / INTMTXSCALE );
2787 			if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2788 		}
2789 	}
2790 }
2791 
WriteHat2_part_pointer(FILE * hat2p,int locnjob,int nadd,char ** name,double ** mtx)2792 void WriteHat2_part_pointer( FILE *hat2p, int locnjob, int nadd, char **name, double **mtx )
2793 {
2794 	int i, j;
2795 	int norg = locnjob-nadd;
2796 	double max = 0.0;
2797 //	for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2798 
2799 	fprintf( hat2p, "%5d\n", 1 );
2800 	fprintf( hat2p, "%5d\n", locnjob );
2801 	fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2802 
2803 	for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2804 	for( i=0; i<norg; i++ )
2805 	{
2806 		for( j=0; j<nadd; j++ )
2807 		{
2808 			fprintf( hat2p, DFORMAT, mtx[i][j] );
2809 			if( (j+1) % 12 == 0 || j == nadd-1 ) fprintf( hat2p, "\n" );
2810 		}
2811 	}
2812 }
2813 
WriteHat2_pointer(FILE * hat2p,int locnjob,char ** name,double ** mtx)2814 void WriteHat2_pointer( FILE *hat2p, int locnjob, char **name, double **mtx )
2815 {
2816 	int i, j;
2817 	double max = 0.0;
2818 	for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2819 
2820 	fprintf( hat2p, "%5d\n", 1 );
2821 	fprintf( hat2p, "%5d\n", locnjob );
2822 	fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2823 
2824 	for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2825 	for( i=0; i<locnjob-1; i++ )
2826 	{
2827 		for( j=i+1; j<locnjob; j++ )
2828 		{
2829 			fprintf( hat2p, DFORMAT, mtx[i][j] );
2830 			if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2831 		}
2832 	}
2833 }
2834 
WriteHat2(FILE * hat2p,int locnjob,char name[M][B],double ** mtx)2835 void WriteHat2( FILE *hat2p, int locnjob, char name[M][B], double **mtx )
2836 {
2837 	int i, j;
2838 	double max = 0.0;
2839 	for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2840 
2841 	fprintf( hat2p, "%5d\n", 1 );
2842 	fprintf( hat2p, "%5d\n", locnjob );
2843 	fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2844 
2845 	for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2846 	for( i=0; i<locnjob-1; i++ )
2847 	{
2848 		for( j=i+1; j<locnjob; j++ )
2849 		{
2850 			fprintf( hat2p, DFORMAT, mtx[i][j] );
2851 			if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2852 		}
2853 	}
2854 }
2855 
2856 #if 0
2857 void WriteHat2plain( FILE *hat2p, int locnjob, double **mtx )
2858 {
2859 	int i, j, ilim;
2860 
2861 	ilim = locnjob-1;
2862 	for( i=0; i<ilim; i++ )
2863 	{
2864 		fprintf( hat2p, "%d-%d d=%.3f\n", i+1, i+1, 0.0 );
2865 		for( j=i+1; j<locnjob; j++ )
2866 		{
2867 			fprintf( hat2p, "%d-%d d=%.3f\n", i+1, j+1, mtx[i][j] );
2868 		}
2869 	}
2870 }
2871 #endif
2872 
ReadFasta_sub(FILE * fp,double * dis,int nseq,char name[M][B])2873 int ReadFasta_sub( FILE *fp, double *dis, int nseq, char name[M][B] )
2874 {
2875     int i, count=0;
2876     char b[B];
2877     int junban[M];
2878 
2879     count = 0;
2880     for( i=0; i<10000000 && count<nseq; i++ )
2881     {
2882         fgets( b, B-1, fp );
2883         if( !strncmp( "+==========+", b, 12 ) )
2884         {
2885             junban[count] = atoi( b+12 );
2886             count++;
2887         }
2888     }
2889 
2890 	for( i=0; i<nseq; i++ ) dis[i] = 0.0;
2891     count = 0;
2892     for( i=0; i<100000 && count<nseq; i++ )
2893     {
2894 		if( fgets( b, B-1, fp ) ) break;
2895         if( !strncmp( name[junban[count]], b, 20  ) )
2896         {
2897             fgets( b, B-1, fp );
2898             dis[junban[count]] = atof( b );
2899             count++;
2900         }
2901     }
2902     return 0;
2903 }
2904 
2905 
ReadSsearch(FILE * fp,double * dis,int nseq,char name[M][B])2906 int ReadSsearch( FILE *fp, double *dis, int nseq, char name[M][B] )
2907 {
2908     int i, count=0;
2909     char b[B];
2910     int junban[M];
2911 	int opt;
2912 
2913     count = 0;
2914     for( i=0; i<10000000 && count<nseq; i++ )
2915     {
2916         fgets( b, B-1, fp );
2917         if( !strncmp( "+==========+", b, 12 ) )
2918         {
2919             junban[count] = atoi( b+12 );
2920 			sscanf( b+75, "%d", &opt );
2921             dis[junban[count]] = (double)opt;
2922             count++;
2923         }
2924     }
2925 
2926 /*
2927     count = 0;
2928     for( i=0; i<100000 && count<nseq; i++ )
2929     {
2930         fgets( b, B-1, fp );
2931         if( !strncmp( name[junban[count]], b, 20  ) )
2932         {
2933             dis[junban[count]] = atof( b+65 );
2934             count++;
2935         }
2936     }
2937 */
2938     return 0;
2939 }
2940 
ReadBlastm7_avscore(FILE * fp,double * dis,int nin)2941 int ReadBlastm7_avscore( FILE *fp, double *dis, int nin )
2942 {
2943     int count=0;
2944     char b[B];
2945 	char *pt;
2946     int *junban;
2947 	double score, sumscore;
2948 	double len, sumlen;
2949 	int qstart, qend, tstart, tend;
2950 	double scorepersite;
2951 	static char qal[N], tal[N], al[N];
2952 	int nlocalhom;
2953 
2954 	junban = calloc( nin, sizeof( int ) );
2955 
2956 	count = 0;
2957 	sumscore = 0.0;
2958 	sumlen = 0.0;
2959 	score = 0.0;
2960 	len = 0.0;
2961 	scorepersite = 0.0; // by D.Mathog, a guess
2962     while( 1 )
2963 	{
2964 
2965 		if( feof( fp ) ) break;
2966 
2967 		while( fgets( b, B-1, fp ) )
2968 		{
2969 			if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
2970 		}
2971 
2972 		if( !strncmp( "          <Hit_def>", b, 19 ) )
2973 		{
2974 			junban[count] = atoi( b+31 );
2975 			nlocalhom = 0;
2976 		}
2977 
2978 
2979 		while( fgets( b, B-1, fp ) )
2980 			if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
2981 		pt = b + 25;
2982 		score = atof( pt );
2983 		sumscore += score;
2984 
2985 
2986 		while( fgets( b, B-1, fp ) )
2987 			if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
2988 		pt = b + 30;
2989 		qstart = atoi( pt ) - 1;
2990 
2991 
2992 		while( fgets( b, B-1, fp ) )
2993 			if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
2994 		pt = b + 28;
2995 		qend = atoi( pt ) - 1;
2996 
2997 
2998 		while( fgets( b, B-1, fp ) )
2999 			if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
3000 		pt = b + 28;
3001 		tstart = atoi( pt ) - 1;
3002 
3003 
3004 		while( fgets( b, B-1, fp ) )
3005 			if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
3006 		pt = b + 26;
3007 		tend = atoi( pt ) - 1;
3008 
3009 
3010 		while( fgets( b, B-1, fp ) )
3011 			if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
3012 		pt = b + 29;
3013 		len = atoi( pt );
3014 		sumlen += len;
3015 
3016 
3017 		while( fgets( al, N-100, fp ) )
3018 			if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
3019 
3020 		strcpy( qal, al+24 );
3021 		pt = qal;
3022 		while( *++pt != '<' )
3023 			;
3024 		*pt = 0;
3025 
3026 
3027 		while( fgets( al, N-100, fp ) )
3028 			if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
3029 
3030 		strcpy( tal, al+24 );
3031 		pt = tal;
3032 		while( *++pt != '<' )
3033 			;
3034 		*pt = 0;
3035 
3036 
3037 //		fprintf( stderr, "t=%d, score = %f, qstart=%d, qend=%d, tstart=%d, tend=%d, overlapaa=%d\n", junban[count], score, qstart, qend, tstart, tend, overlapaa );
3038 
3039 
3040 		while( fgets( b, B-1, fp ) )
3041 			if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
3042 
3043 
3044 		fgets( b, B-1, fp );
3045 
3046 
3047 		if( !strncmp( "          </Hit_hsps>", b, 21 ) )
3048 		{
3049 			dis[junban[count++]] = sumscore;
3050 			sumscore = 0.0;
3051 			fgets( b, B-1, fp );
3052 			fgets( b, B-1, fp );
3053 			scorepersite = sumscore / sumlen;
3054 			if( scorepersite != (int)scorepersite )
3055 			{
3056 				fprintf( stderr, "ERROR! sumscore=%f, sumlen=%f, and scorepersite=%f\n", sumscore, sumlen, scorepersite );
3057 				exit( 1 );
3058 			}
3059 
3060 			if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
3061 		}
3062 	}
3063 
3064 	free( junban );
3065 
3066     return (int)scorepersite;
3067 }
ReadBlastm7_scoreonly(FILE * fp,double * dis,int nin)3068 int ReadBlastm7_scoreonly( FILE *fp, double *dis, int nin )
3069 {
3070     int count=0;
3071     char b[B];
3072 	char *pt;
3073     int *junban;
3074 	int overlapaa;
3075 	double score, sumscore;
3076 	int qstart, qend, tstart, tend;
3077 	static char qal[N], tal[N], al[N];
3078 	int nlocalhom;
3079 
3080 	junban = calloc( nin, sizeof( int ) );
3081 
3082 	count = 0;
3083 	sumscore = 0.0;
3084 	score = 0.0;
3085     while( 1 )
3086 	{
3087 
3088 		if( feof( fp ) ) break;
3089 
3090 		while( fgets( b, B-1, fp ) )
3091 		{
3092 			if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
3093 		}
3094 
3095 		if( !strncmp( "          <Hit_def>", b, 19 ) )
3096 		{
3097 			junban[count] = atoi( b+31 );
3098 			nlocalhom = 0;
3099 		}
3100 
3101 
3102 		while( fgets( b, B-1, fp ) )
3103 			if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
3104 		pt = b + 25;
3105 		score = atof( pt );
3106 		sumscore += score;
3107 
3108 
3109 		while( fgets( b, B-1, fp ) )
3110 			if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
3111 		pt = b + 30;
3112 		qstart = atoi( pt ) - 1;
3113 
3114 
3115 		while( fgets( b, B-1, fp ) )
3116 			if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
3117 		pt = b + 28;
3118 		qend = atoi( pt ) - 1;
3119 
3120 
3121 		while( fgets( b, B-1, fp ) )
3122 			if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
3123 		pt = b + 28;
3124 		tstart = atoi( pt ) - 1;
3125 
3126 
3127 		while( fgets( b, B-1, fp ) )
3128 			if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
3129 		pt = b + 26;
3130 		tend = atoi( pt ) - 1;
3131 
3132 
3133 		while( fgets( b, B-1, fp ) )
3134 			if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
3135 		pt = b + 29;
3136 		overlapaa = atoi( pt );
3137 
3138 
3139 		while( fgets( al, N-100, fp ) )
3140 			if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
3141 
3142 		strcpy( qal, al+24 );
3143 		pt = qal;
3144 		while( *++pt != '<' )
3145 			;
3146 		*pt = 0;
3147 
3148 
3149 		while( fgets( al, N-100, fp ) )
3150 			if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
3151 
3152 		strcpy( tal, al+24 );
3153 		pt = tal;
3154 		while( *++pt != '<' )
3155 			;
3156 		*pt = 0;
3157 
3158 
3159 //		fprintf( stderr, "t=%d, score = %f, qstart=%d, qend=%d, tstart=%d, tend=%d, overlapaa=%d\n", junban[count], score, qstart, qend, tstart, tend, overlapaa );
3160 
3161 //		nlocalhom += addlocalhom_r( qal, tal, localhomlist+junban[count], qstart, tstart, score, overlapaa, nlocalhom );
3162 
3163 		while( fgets( b, B-1, fp ) )
3164 			if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
3165 
3166 
3167 		fgets( b, B-1, fp );
3168 
3169 
3170 		if( !strncmp( "          </Hit_hsps>", b, 21 ) )
3171 		{
3172 			dis[junban[count++]] = sumscore;
3173 			sumscore = 0.0;
3174 			fgets( b, B-1, fp );
3175 			fgets( b, B-1, fp );
3176 			if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
3177 		}
3178 	}
3179 
3180 	free( junban );
3181 
3182     return count;
3183 }
3184 
ReadBlastm7(FILE * fp,double * dis,int qmem,char ** name,LocalHom * localhomlist)3185 int ReadBlastm7( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3186 {
3187     int count=0;
3188     char b[B];
3189 	char *pt;
3190     static int junban[M];
3191 	int overlapaa;
3192 	double score, sumscore;
3193 	int qstart, qend, tstart, tend;
3194 	static char qal[N], tal[N], al[N];
3195 	int nlocalhom;
3196 
3197 
3198 
3199 	count = 0;
3200 	sumscore = 0.0;
3201 	score = 0.0;
3202 	nlocalhom = 0;
3203     while( 1 )
3204 	{
3205 
3206 		if( feof( fp ) ) break;
3207 
3208 		while( fgets( b, B-1, fp ) )
3209 		{
3210 			if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
3211 		}
3212 
3213 		if( !strncmp( "          <Hit_def>", b, 19 ) )
3214 		{
3215 			junban[count] = atoi( b+31 );
3216 			nlocalhom = 0;
3217 		}
3218 
3219 
3220 		while( fgets( b, B-1, fp ) )
3221 			if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
3222 		pt = b + 25;
3223 		score = atof( pt );
3224 		sumscore += score;
3225 
3226 
3227 		while( fgets( b, B-1, fp ) )
3228 			if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
3229 		pt = b + 30;
3230 		qstart = atoi( pt ) - 1;
3231 
3232 
3233 		while( fgets( b, B-1, fp ) )
3234 			if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
3235 		pt = b + 28;
3236 		qend = atoi( pt ) - 1;
3237 
3238 
3239 		while( fgets( b, B-1, fp ) )
3240 			if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
3241 		pt = b + 28;
3242 		tstart = atoi( pt ) - 1;
3243 
3244 
3245 		while( fgets( b, B-1, fp ) )
3246 			if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
3247 		pt = b + 26;
3248 		tend = atoi( pt ) - 1;
3249 
3250 
3251 		while( fgets( b, B-1, fp ) )
3252 			if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
3253 		pt = b + 29;
3254 		overlapaa = atoi( pt );
3255 
3256 
3257 		while( fgets( al, N-100, fp ) )
3258 			if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
3259 
3260 		strcpy( qal, al+24 );
3261 		pt = qal;
3262 		while( *++pt != '<' )
3263 			;
3264 		*pt = 0;
3265 
3266 
3267 		while( fgets( al, N-100, fp ) )
3268 			if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
3269 
3270 		strcpy( tal, al+24 );
3271 		pt = tal;
3272 		while( *++pt != '<' )
3273 			;
3274 		*pt = 0;
3275 
3276 
3277 //		fprintf( stderr, "t=%d, score = %f, qstart=%d, qend=%d, tstart=%d, tend=%d, overlapaa=%d\n", junban[count], score, qstart, qend, tstart, tend, overlapaa );
3278 
3279 		nlocalhom += addlocalhom_r( qal, tal, localhomlist+junban[count], qstart, tstart, score, overlapaa, nlocalhom );
3280 
3281 		while( fgets( b, B-1, fp ) )
3282 			if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
3283 
3284 
3285 		fgets( b, B-1, fp );
3286 
3287 
3288 		if( !strncmp( "          </Hit_hsps>", b, 21 ) )
3289 		{
3290 			dis[junban[count++]] = sumscore;
3291 			sumscore = 0.0;
3292 			fgets( b, B-1, fp );
3293 			fgets( b, B-1, fp );
3294 			if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
3295 		}
3296 	}
3297     return count;
3298 }
3299 
ReadFasta34noalign(FILE * fp,double * dis,int qmem,char ** name,LocalHom * localhomlist)3300 int ReadFasta34noalign( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3301 {
3302     int count=0;
3303     char b[B];
3304 	char *pt;
3305     static int junban[M];
3306 	int opt;
3307 	double z, bits;
3308 
3309 
3310     count = 0;
3311 #if 0
3312     for( i=0; i<10000000 && count<nseq; i++ )
3313 #else
3314     while( !feof( fp ) )
3315 #endif
3316     {
3317         fgets( b, B-1, fp );
3318         if( !strncmp( "+==========+", b, 12 ) )
3319         {
3320             junban[count] = atoi( b+12 );
3321 
3322 			pt = strchr( b, ')' ) + 1;
3323 			sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3324             dis[junban[count]] = (double)opt;
3325             count++;
3326 
3327         }
3328     }
3329 
3330     return count;
3331 }
ReadFasta34m10_nuc(FILE * fp,double * dis,int qmem,char ** name,LocalHom * localhomlist)3332 int ReadFasta34m10_nuc( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3333 {
3334     int count=0;
3335     char b[B];
3336 	char *pt;
3337     static int junban[M];
3338 	int overlapaa;
3339 	int opt, qstart, qend, tstart, tend;
3340 	double z, bits;
3341 	int qal_display_start, tal_display_start;
3342 	static char qal[N], tal[N];
3343 	char *qal2, *tal2;
3344 	int c;
3345 
3346 
3347     count = 0;
3348 #if 0
3349     for( i=0; i<10000000 && count<nseq; i++ )
3350 #else
3351     while( !feof( fp ) )
3352 #endif
3353     {
3354         fgets( b, B-1, fp );
3355         if( !strncmp( "+==========+", b, 12 ) )
3356         {
3357             junban[count] = atoi( b+12 );
3358 
3359 			if( strchr( b, 'r' ) ) continue;
3360 
3361 			pt = strchr( b, ']' ) + 1;
3362 			sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3363             dis[junban[count]] = (double)opt;
3364             count++;
3365 
3366         }
3367 		else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3368 		{
3369 			break;
3370 		}
3371 
3372     }
3373 	if( !count ) return -1;
3374 
3375 	count = 0;
3376     while( 1 )
3377 	{
3378 		if( strncmp( ">>+==========+", b, 14 ) )
3379 		{
3380 			fgets( b, B-1, fp );
3381 			if( feof( fp ) ) break;
3382 			continue;
3383 		}
3384 		junban[count++] = atoi( b+14 );
3385 //		fprintf( stderr, "t = %d\n", atoi( b+14 ) );
3386 		while( fgets( b, B-1, fp ) )
3387 			if( !strncmp( "; fa_opt:", b, 9 ) || !strncmp( "; sw_s-w opt:", b, 13 ) ) break;
3388 		pt = strstr( b, ":" ) +1;
3389 		opt = atoi( pt );
3390 
3391 
3392 		while( fgets( b, B-1, fp ) )
3393 			if( !strncmp( "_overlap:", b+4, 9 ) ) break;
3394 		pt = strstr( b, ":" ) +1;
3395 		overlapaa = atoi( pt );
3396 
3397 		while( fgets( b, B-1, fp ) )
3398 			if( !strncmp( "_start:", b+4, 7 ) ) break;
3399 		pt = strstr( b, ":" ) +1;
3400 		qstart = atoi( pt ) - 1;
3401 
3402 		while( fgets( b, B-1, fp ) )
3403 			if( !strncmp( "_stop:", b+4, 6 ) ) break;
3404 		pt = strstr( b, ":" ) +1;
3405 		qend = atoi( pt ) - 1;
3406 
3407 		while( fgets( b, B-1, fp ) )
3408 			if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3409 		pt = strstr( b, ":" ) +1;
3410 		qal_display_start = atoi( pt ) - 1;
3411 
3412 		pt = qal;
3413 		while( (c = fgetc( fp )) )
3414 		{
3415 			if( c == '>' )
3416 			{
3417 				ungetc( c, fp );
3418 				break;
3419 			}
3420 			if( isalpha( c ) || c == '-' )
3421 			*pt++ = c;
3422 		}
3423 		*pt = 0;
3424 
3425 		while( fgets( b, B-1, fp ) )
3426 			if( !strncmp( "_start:", b+4, 7 ) ) break;
3427 		pt = strstr( b, ":" ) + 1;
3428 		tstart = atoi( pt ) - 1;
3429 
3430 		while( fgets( b, B-1, fp ) )
3431 			if( !strncmp( "_stop:", b+4, 6 ) ) break;
3432 		pt = strstr( b, ":" ) + 1;
3433 		tend = atoi( pt ) - 1;
3434 
3435 		while( fgets( b, B-1, fp ) )
3436 			if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3437 		pt = strstr( b, ":" ) + 1;
3438 		tal_display_start = atoi( pt ) - 1;
3439 
3440 		pt = tal;
3441 		while( ( c = fgetc( fp ) ) )
3442 		{
3443 			if( c == '>' )
3444 			{
3445 				ungetc( c, fp );
3446 				break;
3447 			}
3448 			if( isalpha( c ) || c == '-' )
3449 			*pt++ = c;
3450 		}
3451 		*pt = 0;
3452 
3453 //		fprintf( stderr, "(%d-%d:%d-%d)\n", qstart, qend, tstart, tend );
3454 //		fprintf( stderr, "qal_display_start = %d, tal_display_start = %d\n", qal_display_start, tal_display_start );
3455 
3456 //		fprintf( stderr, "qal = %s\n", qal );
3457 //		fprintf( stderr, "tal = %s\n", tal );
3458 
3459 		qal2 = cutal( qal, qal_display_start, qstart, qend );
3460 		tal2 = cutal( tal, tal_display_start, tstart, tend );
3461 
3462 //		fprintf( stderr, "qal2 = %s\n", qal2 );
3463 //		fprintf( stderr, "tal2 = %s\n", tal2 );
3464 
3465 //		fprintf( stderr, "putting   %d - %d, opt = %d\n", qmem, junban[count-1], opt );
3466 		putlocalhom( qal2, tal2, localhomlist+junban[count-1], qstart, tstart, opt, overlapaa );
3467 	}
3468 //	fprintf( stderr, "count = %d\n", count );
3469     return count;
3470 }
ReadFasta34m10(FILE * fp,double * dis,int qmem,char ** name,LocalHom * localhomlist)3471 int ReadFasta34m10( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3472 {
3473     int count=0;
3474     char b[B];
3475 	char *pt;
3476     static int junban[M];
3477 	int overlapaa;
3478 	int opt, qstart, qend, tstart, tend;
3479 	double z, bits;
3480 	int qal_display_start, tal_display_start;
3481 	static char qal[N], tal[N];
3482 	char *qal2, *tal2;
3483 	int c;
3484 
3485 
3486     count = 0;
3487 #if 0
3488     for( i=0; i<10000000 && count<nseq; i++ )
3489 #else
3490     while( !feof( fp ) )
3491 #endif
3492     {
3493         fgets( b, B-1, fp );
3494         if( !strncmp( "+==========+", b, 12 ) )
3495         {
3496             junban[count] = atoi( b+12 );
3497 
3498 			pt = strchr( b, ')' ) + 1;
3499 			sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3500             dis[junban[count]] = (double)opt;
3501             count++;
3502 
3503         }
3504 		else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3505 		{
3506 			break;
3507 		}
3508 
3509     }
3510 	if( !count ) return -1;
3511 
3512 	count = 0;
3513     while( 1 )
3514 	{
3515 		if( strncmp( ">>+==========+", b, 14 ) )
3516 		{
3517 			fgets( b, B-1, fp );
3518 			if( feof( fp ) ) break;
3519 			continue;
3520 		}
3521 		junban[count++] = atoi( b+14 );
3522 //		fprintf( stderr, "t = %d\n", atoi( b+14 ) );
3523 		while( fgets( b, B-1, fp ) )
3524 			if( !strncmp( "; fa_opt:", b, 9 ) || !strncmp( "; sw_s-w opt:", b, 13 ) ) break;
3525 		pt = strstr( b, ":" ) +1;
3526 		opt = atoi( pt );
3527 
3528 
3529 		while( fgets( b, B-1, fp ) )
3530 			if( !strncmp( "_overlap:", b+4, 9 ) ) break;
3531 		pt = strstr( b, ":" ) +1;
3532 		overlapaa = atoi( pt );
3533 
3534 		while( fgets( b, B-1, fp ) )
3535 			if( !strncmp( "_start:", b+4, 7 ) ) break;
3536 		pt = strstr( b, ":" ) +1;
3537 		qstart = atoi( pt ) - 1;
3538 
3539 		while( fgets( b, B-1, fp ) )
3540 			if( !strncmp( "_stop:", b+4, 6 ) ) break;
3541 		pt = strstr( b, ":" ) +1;
3542 		qend = atoi( pt ) - 1;
3543 
3544 		while( fgets( b, B-1, fp ) )
3545 			if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3546 		pt = strstr( b, ":" ) +1;
3547 		qal_display_start = atoi( pt ) - 1;
3548 
3549 		pt = qal;
3550 		while( (c = fgetc( fp )) )
3551 		{
3552 			if( c == '>' )
3553 			{
3554 				ungetc( c, fp );
3555 				break;
3556 			}
3557 			if( isalpha( c ) || c == '-' )
3558 			*pt++ = c;
3559 		}
3560 		*pt = 0;
3561 
3562 		while( fgets( b, B-1, fp ) )
3563 			if( !strncmp( "_start:", b+4, 7 ) ) break;
3564 		pt = strstr( b, ":" ) + 1;
3565 		tstart = atoi( pt ) - 1;
3566 
3567 		while( fgets( b, B-1, fp ) )
3568 			if( !strncmp( "_stop:", b+4, 6 ) ) break;
3569 		pt = strstr( b, ":" ) + 1;
3570 		tend = atoi( pt ) - 1;
3571 
3572 		while( fgets( b, B-1, fp ) )
3573 			if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3574 		pt = strstr( b, ":" ) + 1;
3575 		tal_display_start = atoi( pt ) - 1;
3576 
3577 		pt = tal;
3578 		while( ( c = fgetc( fp ) ) )
3579 		{
3580 			if( c == '>' )
3581 			{
3582 				ungetc( c, fp );
3583 				break;
3584 			}
3585 			if( isalpha( c ) || c == '-' )
3586 			*pt++ = c;
3587 		}
3588 		*pt = 0;
3589 
3590 //		fprintf( stderr, "(%d-%d:%d-%d)\n", qstart, qend, tstart, tend );
3591 //		fprintf( stderr, "qal_display_start = %d, tal_display_start = %d\n", qal_display_start, tal_display_start );
3592 
3593 //		fprintf( stderr, "qal = %s\n", qal );
3594 //		fprintf( stderr, "tal = %s\n", tal );
3595 
3596 		qal2 = cutal( qal, qal_display_start, qstart, qend );
3597 		tal2 = cutal( tal, tal_display_start, tstart, tend );
3598 
3599 //		fprintf( stderr, "qal2 = %s\n", qal2 );
3600 //		fprintf( stderr, "tal2 = %s\n", tal2 );
3601 
3602 //		fprintf( stderr, "putting   %d - %d, opt = %d\n", qmem, junban[count-1], opt );
3603 		putlocalhom( qal2, tal2, localhomlist+junban[count-1], qstart, tstart, opt, overlapaa );
3604 	}
3605 //	fprintf( stderr, "count = %d\n", count );
3606     return count;
3607 }
ReadFasta34m10_scoreonly_nucbk(FILE * fp,double * dis,int nin)3608 int ReadFasta34m10_scoreonly_nucbk( FILE *fp, double *dis, int nin )
3609 {
3610     int count=0;
3611     char b[B];
3612 	char *pt;
3613     int pos;
3614 	int opt;
3615 	double z, bits;
3616 
3617     count = 0;
3618     while( !feof( fp ) )
3619     {
3620         fgets( b, B-1, fp );
3621         if( !strncmp( "+===========+", b, 13 ) )
3622         {
3623             pos = atoi( b+13 );
3624 
3625 			if( strchr( b, 'r' ) ) continue;
3626 
3627 //			pt = strchr( b, ')' ) + 1;
3628 			pt = strchr( b, ']' ) + 1;
3629 			sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3630             dis[pos] += (double)opt;
3631             count++;
3632 #if 0
3633 			fprintf( stderr, "b=%s\n", b );
3634 			fprintf( stderr, "opt=%d\n", opt );
3635 			fprintf( stderr, "pos=%d\n", pos );
3636 			fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3637 #endif
3638 
3639         }
3640 		else if( 0 == strncmp( ">>><<<", b, 6 ) )
3641 		{
3642 			break;
3643 		}
3644 
3645     }
3646 	if( !count ) return -1;
3647 
3648     return count;
3649 }
3650 
ReadFasta34m10_scoreonly_nuc(FILE * fp,double * dis,int nin)3651 int ReadFasta34m10_scoreonly_nuc( FILE *fp, double *dis, int nin )
3652 {
3653     int count=0;
3654     char b[B];
3655 	char *pt;
3656     int pos;
3657 	int opt;
3658 	double z, bits;
3659 	int c;
3660 	int *yonda;
3661 
3662 
3663 	yonda = AllocateIntVec( nin );
3664 	for( c=0; c<nin; c++ ) yonda[c] = 0;
3665 	for( c=0; c<nin; c++ ) dis[c] = 0.0;
3666 
3667     count = 0;
3668     while( !feof( fp ) )
3669     {
3670         fgets( b, B-1, fp );
3671         if( !strncmp( "+===========+", b, 13 ) )
3672         {
3673             pos = atoi( b+13 );
3674 
3675 			if( strchr( b, 'r' ) ) continue;
3676 
3677 //			pt = strchr( b, ')' ) + 1;
3678 			pt = strchr( b, ']' ) + 1;
3679 			sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3680 			if( yonda[pos] == 0 )
3681 			{
3682 	            dis[pos] += (double)opt;
3683 				yonda[pos] = 1;
3684 			}
3685             count++;
3686 #if 0
3687 			fprintf( stderr, "b=%s\n", b );
3688 			fprintf( stderr, "opt=%d\n", opt );
3689 			fprintf( stderr, "pos=%d\n", pos );
3690 			fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3691 #endif
3692 
3693         }
3694         else if( !strncmp( ">>>", b, 3 ) )
3695 		{
3696 			for( c=0; c<nin; c++ ) yonda[c] = 0;
3697 		}
3698 		else if( 0 == strncmp( ">>><<<", b, 6 ) )
3699 		{
3700 			break;
3701 		}
3702 
3703     }
3704 
3705 	free( yonda );
3706 
3707 	if( !count ) return -1;
3708 
3709     return count;
3710 }
3711 
ReadFasta34m10_scoreonly(FILE * fp,double * dis,int nin)3712 int ReadFasta34m10_scoreonly( FILE *fp, double *dis, int nin )
3713 {
3714     int count=0;
3715     char b[B];
3716 	char *pt;
3717     int pos;
3718 	int opt;
3719 	double z, bits;
3720 	int c;
3721 	int *yonda;
3722 
3723 
3724 	yonda = AllocateIntVec( nin );
3725 	for( c=0; c<nin; c++ ) yonda[c] = 0;
3726 	for( c=0; c<nin; c++ ) dis[c] = 0.0;
3727 
3728     count = 0;
3729     while( !feof( fp ) )
3730     {
3731         fgets( b, B-1, fp );
3732         if( !strncmp( "+===========+", b, 13 ) )
3733         {
3734             pos = atoi( b+13 );
3735 
3736 			pt = strchr( b, ')' ) + 1;
3737 			sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3738 			if( yonda[pos] == 0 )
3739 			{
3740 	            dis[pos] += (double)opt;
3741 				yonda[pos] = 1;
3742 			}
3743             count++;
3744 #if 0
3745 			fprintf( stderr, "b=%s\n", b );
3746 			fprintf( stderr, "opt=%d\n", opt );
3747 			fprintf( stderr, "pos=%d\n", pos );
3748 			fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3749 #endif
3750 
3751         }
3752         else if( !strncmp( ">>>", b, 3 ) )
3753 		{
3754 			for( c=0; c<nin; c++ ) yonda[c] = 0;
3755 		}
3756 		else if( 0 == strncmp( ">>><<<", b, 6 ) )
3757 		{
3758 			break;
3759 		}
3760 
3761     }
3762 
3763 	free( yonda );
3764 
3765 	if( !count ) return -1;
3766 
3767     return count;
3768 }
ReadFasta34(FILE * fp,double * dis,int nseq,char name[M][B],LocalHom * localhomlist)3769 int ReadFasta34( FILE *fp, double *dis, int nseq, char name[M][B], LocalHom *localhomlist )
3770 {
3771     int count=0;
3772     char b[B];
3773 	char *pt;
3774     static int junban[M];
3775 	int overlapaa;
3776 	int opt, qstart, qend, tstart, tend;
3777 	double z, bits;
3778 
3779 
3780     count = 0;
3781 #if 0
3782     for( i=0; i<10000000 && count<nseq; i++ )
3783 #else
3784     while( !feof( fp ) )
3785 #endif
3786     {
3787         fgets( b, B-1, fp );
3788         if( !strncmp( "+==========+", b, 12 ) )
3789         {
3790             junban[count] = atoi( b+12 );
3791 
3792 			pt = strchr( b, ')' ) + 1;
3793 			sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3794             dis[junban[count]] = (double)opt;
3795             count++;
3796 
3797         }
3798 		else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3799 		{
3800 			break;
3801 		}
3802 
3803     }
3804 	if( !count ) return -1;
3805 
3806 	count = 0;
3807     while( !feof( fp ) )
3808 	{
3809 		if( !strncmp(">>+==========+", b, 14 ) )
3810 		{
3811             junban[count] = atoi( b+14 );
3812             count++;
3813         	fgets( b, B-1, fp ); // initn:
3814 			pt = strstr( b, "opt: " ) + 5;
3815 			localhomlist[junban[count-1]].opt = atof( pt );
3816         	fgets( b, B-1, fp ); // Smith-Waterman score
3817 			pt = strstr( b, "ungapped) in " ) + 13;
3818 			sscanf( pt, "%d", &overlapaa );
3819 			fprintf( stderr, "pt = %s, overlapaa = %d\n", pt, overlapaa );
3820 			pt = strstr( b, "overlap (" ) + 8;
3821 			sscanf( pt, "(%d-%d:%d-%d)", &qstart, &qend, &tstart, &tend );
3822 			localhomlist[junban[count-1]].overlapaa = overlapaa;
3823 			localhomlist[junban[count-1]].start1 = qstart-1;
3824 			localhomlist[junban[count-1]].end1   = qend-1;
3825 			localhomlist[junban[count-1]].start2 = tstart-1;
3826 			localhomlist[junban[count-1]].end2   = tend-1;
3827 		}
3828         fgets( b, B-1, fp );
3829 	}
3830 	fprintf( stderr, "count = %d\n", count );
3831     return count;
3832 }
3833 
ReadFasta3(FILE * fp,double * dis,int nseq,char name[M][B])3834 int ReadFasta3( FILE *fp, double *dis, int nseq, char name[M][B] )
3835 {
3836     int count=0;
3837     char b[B];
3838 	char *pt;
3839     int junban[M];
3840 	int initn, init1, opt;
3841 	double z;
3842 
3843     count = 0;
3844 #if 0
3845     for( i=0; i<10000000 && count<nseq; i++ )
3846 #else
3847     while( !feof( fp ) )
3848 #endif
3849     {
3850         fgets( b, B-1, fp );
3851         if( !strncmp( "+==========+", b, 12 ) )
3852         {
3853             junban[count] = atoi( b+12 );
3854 
3855 			pt = strchr( b, ')' ) + 1;
3856 			sscanf( pt, "%d %d %d %lf", &initn, &init1, &opt, &z );
3857             dis[junban[count]] = (double)opt;
3858             count++;
3859         }
3860     }
3861     return 0;
3862 }
3863 
ReadFasta(FILE * fp,double * dis,int nseq,char name[M][B])3864 int ReadFasta( FILE *fp, double *dis, int nseq, char name[M][B] )
3865 {
3866     int i, count=0;
3867     char b[B];
3868     int junban[M];
3869 	int initn, init1, opt;
3870 
3871     count = 0;
3872 	for( i=0; i<nseq; i++ ) dis[i] = 0.0;
3873     for( i=0; !feof( fp ) && count<nseq; i++ )
3874     {
3875         fgets( b, B-1, fp );
3876         if( !strncmp( "+==========+", b, 12 ) )
3877         {
3878             junban[count] = atoi( b+12 );
3879 
3880 			sscanf( b+50, "%d %d %d", &initn, &init1, &opt );
3881             dis[junban[count]] = (double)opt;
3882             count++;
3883         }
3884     }
3885 
3886 /*
3887     count = 0;
3888     for( i=0; i<100000 && count<nseq; i++ )
3889     {
3890         fgets( b, B-1, fp );
3891         if( !strncmp( name[junban[count]], b, 20  ) )
3892         {
3893             dis[junban[count]] = atof( b+65 );
3894             count++;
3895         }
3896     }
3897 */
3898     return 0;
3899 }
3900 
3901 
ReadOpt(FILE * fp,int opt[M],int nseq,char name[M][B])3902 int ReadOpt( FILE *fp, int opt[M], int nseq, char name[M][B] )
3903 {
3904     int i, count=0;
3905     char b[B];
3906     int junban[M];
3907 	int optt, initn, init1;
3908 
3909     count = 0;
3910     for( i=0; i<10000000 && count<nseq; i++ )
3911     {
3912         fgets( b, B-1, fp );
3913         if( !strncmp( "+==========+", b, 12 ) )
3914         {
3915             junban[count] = atoi( b+12 );
3916 			sscanf( b+50, "%d %d %d", &initn, &init1, &optt );
3917             opt[junban[count]] = (double)optt;
3918             count++;
3919         }
3920     }
3921     return 0;
3922 }
3923 
ReadOpt2(FILE * fp,int opt[M],int nseq,char name[M][B])3924 int ReadOpt2( FILE *fp, int opt[M], int nseq, char name[M][B] )
3925 {
3926     int i, count=0;
3927     char b[B];
3928     int junban[M];
3929 
3930     count = 0;
3931     for( i=0; i<10000000 && count<nseq; i++ )
3932     {
3933         fgets( b, B-1, fp );
3934         if( !strncmp( "+==========+", b, 12 ) )
3935         {
3936             junban[count] = atoi( b+12 );
3937             opt[junban[count]] = atoi( b+65 );
3938             count++;
3939         }
3940     }
3941     return 0;
3942 }
3943 
3944 
3945 
writePre(int nseq,char ** name,int nlen[M],char ** aseq,int force)3946 int writePre( int nseq, char **name, int nlen[M], char **aseq, int force )
3947 {
3948 #if USE_XCED
3949 	int i, value;
3950 	if( !signalSM )
3951 	{
3952 		if( force )
3953 		{
3954 			rewind( prep_g );
3955 			if( devide ) dvWrite( prep_g, nseq, name, nlen, aseq );
3956 #if 0
3957 			else    WriteGapFill( prep_g, nseq, name, nlen, aseq );
3958 #else
3959 			else    writeData( prep_g, nseq, name, nlen, aseq );
3960 #endif
3961 		}
3962 		return( 0 );
3963 	}
3964 	for( i=0; i<10; i++ )
3965 	{
3966 #if IODEBUG
3967 		fprintf( stderr, "SEMAPHORE = %d\n", signalSM[SEMAPHORE] );
3968 #endif
3969 		if( signalSM[SEMAPHORE]-- > 0 )
3970 		{
3971 #if 0 /* /tmp/pre �δط��ǤϤ����� */
3972 			if( ferror( prep_g ) ) prep_g = fopen( "pre", "w" );
3973 			if( !prep_g ) ErrorExit( "Cannot re-open pre." );
3974 #endif
3975 			rewind( prep_g );
3976 			signalSM[STATUS] = IMA_KAITERU;
3977 #if IODEBUG
3978 			if( force ) fprintf( stderr, "FINAL " );
3979 #endif
3980 			if( devide ) dvWrite( prep_g, nseq, name, nlen, aseq );
3981 			else    WriteGapFill( prep_g, nseq, name, nlen, aseq );
3982 			/*
3983 			fprintf( prep_g, '\EOF' );
3984 			*/
3985 			fflush( prep_g );
3986 			if( force ) signalSM[STATUS] = OSHIMAI;
3987 			else        signalSM[STATUS] = KAKIOWATTA;
3988 			value = 1;
3989 			signalSM[SEMAPHORE]++;
3990 #if IODEBUG
3991 			fprintf( stderr, "signalSM[STATUS] = %c\n", signalSM[STATUS] );
3992 #endif
3993 			break;
3994 		}
3995 		else
3996 		{
3997 #if IODEBUG
3998 			fprintf( stderr, "YONDERUKARA_AKIRAMERU\n" );
3999 #endif
4000 			value = 0;
4001 			signalSM[SEMAPHORE]++;
4002 			if( !force ) break;
4003 #if IODEBUG
4004 			fprintf( stderr, "MATSU\n" );
4005 #endif
4006 			sleep( 1 );
4007 		}
4008 	}
4009 	if( force && !value ) ErrorExit( "xced ga pre wo hanasanai \n" );
4010 	return( value );
4011 #else
4012 	if( force )
4013 	{
4014 		rewind( prep_g );
4015 		writeData_pointer( prep_g, nseq, name, nlen, aseq );
4016 	}
4017 #endif
4018 	return( 0 );
4019 }
4020 
4021 
readOtherOptions(int * ppidptr,int * fftThresholdptr,int * fftWinSizeptr)4022 void readOtherOptions( int *ppidptr, int *fftThresholdptr, int *fftWinSizeptr )
4023 {
4024 	if( calledByXced )
4025 	{
4026 		FILE *fp = fopen( "pre", "r" );
4027 		char b[B];
4028 		if( !fp ) ErrorExit( "Cannot open pre.\n" );
4029 		fgets( b, B-1, fp );
4030 		sscanf( b, "%d %d %d", ppidptr, fftThresholdptr, fftWinSizeptr );
4031 		fclose( fp );
4032 #if IODEBUG
4033 	fprintf( stderr, "b = %s\n", b );
4034 	fprintf( stderr, "ppid = %d\n", ppid );
4035 	fprintf( stderr, "fftThreshold = %d\n", fftThreshold );
4036 	fprintf( stderr, "fftWinSize = %d\n", fftWinSize );
4037 #endif
4038 	}
4039 	else
4040 	{
4041 		*ppidptr = 0;
4042 		*fftThresholdptr = FFT_THRESHOLD;
4043 		if( dorp == 'd' )
4044 			*fftWinSizeptr = FFT_WINSIZE_D;
4045 		else
4046 			*fftWinSizeptr = FFT_WINSIZE_P;
4047 	}
4048 #if 0
4049 	fprintf( stderr, "fftThresholdptr=%d\n", *fftThresholdptr );
4050 	fprintf( stderr, "fftWinSizeptr=%d\n", *fftWinSizeptr );
4051 #endif
4052 }
4053 
initSignalSM(void)4054 void initSignalSM( void )
4055 {
4056 //	int signalsmid;
4057 
4058 #if IODEBUG
4059 	if( ppid ) fprintf( stderr, "PID of xced = %d\n", ppid );
4060 #endif
4061 	if( !ppid )
4062 	{
4063 		signalSM = NULL;
4064 		return;
4065 	}
4066 
4067 #if 0
4068 	signalsmid = shmget( (key_t)ppid, 3, IPC_ALLOC | 0666 );
4069 	if( signalsmid == -1 ) ErrorExit( "Cannot get Shared memory for signal.\n" );
4070 	signalSM = shmat( signalsmid, 0, 0 );
4071 	if( (int)signalSM == -1 ) ErrorExit( "Cannot attatch Shared Memory for signal!\n" );
4072 	signalSM[STATUS] = IMA_KAITERU;
4073 	signalSM[SEMAPHORE] = 1;
4074 #endif
4075 }
4076 
initFiles(void)4077 void initFiles( void )
4078 {
4079 	char pname[100];
4080 	if( ppid )
4081 		sprintf( pname, "/tmp/pre.%d", ppid );
4082 	else
4083 		sprintf( pname, "pre" );
4084 	prep_g = fopen( pname, "w" );
4085 	if( !prep_g ) ErrorExit( "Cannot open pre" );
4086 
4087 	trap_g = fopen( "trace", "w" );
4088 	if( !trap_g ) ErrorExit( "cannot open trace" );
4089 	fprintf( trap_g, "PID = %d\n", getpid() );
4090 	fflush( trap_g );
4091 }
4092 
closeFiles(void)4093 void closeFiles( void )
4094 {
4095 	fclose( prep_g );
4096 	fclose( trap_g );
4097 }
4098 
4099 
WriteForFasta(FILE * fp,int locnjob,char ** name,int nlen[M],char ** aseq)4100 void WriteForFasta( FILE *fp, int locnjob, char **name, int nlen[M], char **aseq )
4101 {
4102     static char b[N];
4103     int i, j;
4104     int nalen[M];
4105 
4106     for( i=0; i<locnjob; i++ )
4107     {
4108         nalen[i] = strlen( aseq[i] );
4109         fprintf( fp, ">%s\n", name[i] );
4110         for( j=0; j<nalen[i]; j=j+C )
4111         {
4112             strncpy( b, aseq[i]+j, C ); b[C] = 0;
4113             fprintf( fp, "%s\n",b );
4114         }
4115     }
4116 }
4117 
readlocalhomtable2(FILE * fp,int njob,LocalHom ** localhomtable,char * kozoarivec)4118 void readlocalhomtable2( FILE*fp, int njob, LocalHom **localhomtable, char *kozoarivec )
4119 {
4120 	double opt;
4121 	static char buff[B];
4122 	char infor[100];
4123 	int i, j, overlapaa, start1, end1, start2, end2;
4124 	LocalHom *tmpptr1, *tmpptr2;
4125 
4126 //	for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) nlocalhom[i][j] = 0;
4127 
4128 	while ( NULL != fgets( buff, B-1, fp ) )
4129 	{
4130 //		fprintf( stderr, "\n" );
4131 		sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4132 		if( *infor == 'k' ) kozoarivec[i] = kozoarivec[j] = 1;
4133 
4134 #if 0
4135 		if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4136 #endif
4137 
4138 //		if( i < j )
4139 		{
4140 			if( localhomtable[i][j].nokori++ > 0 )
4141 			{
4142 				tmpptr1 = localhomtable[i][j].last;
4143 //				fprintf( stderr, "reallocating, localhomtable[%d][%d].nokori = %d\n", i, j, localhomtable[i][j].nokori );
4144 				tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4145 				tmpptr1 = tmpptr1->next;
4146 				tmpptr1->extended = -1;
4147 				tmpptr1->next = NULL;
4148 				localhomtable[i][j].last = tmpptr1;
4149 //				fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", i, j, localhomtable[i][j].nokori );
4150 			}
4151 			else
4152 			{
4153 				tmpptr1 = localhomtable[i]+j;
4154 //				fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", i, j, localhomtable[i][j].nokori );
4155 			}
4156 
4157 			tmpptr1->start1 = start1;
4158 			tmpptr1->start2 = start2;
4159 			tmpptr1->end1 = end1;
4160 			tmpptr1->end2 = end2;
4161 //			tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4162 //			tmpptr1->opt = opt;
4163 			tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4164 			tmpptr1->overlapaa = overlapaa;
4165 			tmpptr1->korh = *infor;
4166 		}
4167 //		else
4168 		{
4169 			if( localhomtable[j][i].nokori++ > 0 )
4170 			{
4171 				tmpptr2 = localhomtable[j][i].last;
4172 				tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4173 				tmpptr2 = tmpptr2->next;
4174 				tmpptr2->extended = -1;
4175 				tmpptr2->next = NULL;
4176 				localhomtable[j][i].last = tmpptr2;
4177 //				fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", j, i, localhomtable[j][i].nokori );
4178 			}
4179 			else
4180 			{
4181 				tmpptr2 = localhomtable[j]+i;
4182 //				fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", j, i, localhomtable[j][i].nokori );
4183 			}
4184 
4185 			tmpptr2->start2 = start1;
4186 			tmpptr2->start1 = start2;
4187 			tmpptr2->end2 = end1;
4188 			tmpptr2->end1 = end2;
4189 //			tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4190 //			tmpptr2->opt = opt;
4191 			tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
4192 			tmpptr2->overlapaa = overlapaa;
4193 			tmpptr2->korh = *infor;
4194 
4195 //			fprintf( stderr, "i=%d, j=%d, st1=%d, en1=%d, opt = %f\n", i, j, tmpptr1->start1, tmpptr1->end1, opt );
4196 		}
4197 
4198 	}
4199 }
readlocalhomtable(FILE * fp,int njob,LocalHom ** localhomtable,char * kozoarivec)4200 void readlocalhomtable( FILE*fp, int njob, LocalHom **localhomtable, char *kozoarivec )
4201 {
4202 	double opt;
4203 	static char buff[B];
4204 	char infor[100];
4205 	int i, j, overlapaa, start1, end1, start2, end2;
4206 	int **nlocalhom = NULL;
4207 	LocalHom *tmpptr1=NULL, *tmpptr2=NULL; // by D.Mathog, a guess
4208 
4209 	nlocalhom = AllocateIntMtx( njob, njob );
4210 	for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) nlocalhom[i][j] = 0;
4211 
4212 	while ( NULL != fgets( buff, B-1, fp ) )
4213 	{
4214 //		fprintf( stderr, "\n" );
4215 		sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4216 		if( *infor == 'k' ) kozoarivec[i] = kozoarivec[j] = 1;
4217 
4218 #if 0
4219 		if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4220 #endif
4221 
4222 
4223 //		if( i < j )
4224 		{
4225 			if( nlocalhom[i][j]++ > 0 )
4226 			{
4227 //				fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4228 				tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4229 				tmpptr1 = tmpptr1->next;
4230 				tmpptr1->next = NULL;
4231 			}
4232 			else
4233 			{
4234 				tmpptr1 = localhomtable[i]+j;
4235 //				fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4236 			}
4237 
4238 			tmpptr1->start1 = start1; // CHUUI!!!!
4239 			tmpptr1->start2 = start2;
4240 			tmpptr1->end1 = end1; // CHUUI!!!!
4241 			tmpptr1->end2 = end2;
4242 //			tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4243 //			tmpptr1->opt = opt;
4244 			tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4245 			tmpptr1->overlapaa = overlapaa;
4246 			tmpptr1->korh = *infor;
4247 
4248 //			fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
4249 		}
4250 //		else
4251 		{
4252 			if( nlocalhom[j][i]++ > 0 )
4253 			{
4254 				tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4255 				tmpptr2 = tmpptr2->next;
4256 				tmpptr2->next = NULL;
4257 			}
4258 			else
4259 				tmpptr2 = localhomtable[j]+i;
4260 
4261 			tmpptr2->start2 = start1; // CHUUI!!!!
4262 			tmpptr2->start1 = start2;
4263 			tmpptr2->end2 = end1; // CHUUI!!!!
4264 			tmpptr2->end1 = end2;
4265 //			tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4266 //			tmpptr2->opt = opt;
4267 			tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
4268 			tmpptr2->overlapaa = overlapaa;
4269 			tmpptr2->korh = *infor;
4270 
4271 //			fprintf( stderr, "j=%d, i=%d, opt = %f\n", j, i, opt );
4272 		}
4273 
4274 	}
4275 	FreeIntMtx( nlocalhom );
4276 }
4277 
4278 
readlocalhomtable_two(FILE * fp,int norg,int nadd,LocalHom ** localhomtable,LocalHom ** localhomtablex,char * kozoarivec)4279 void readlocalhomtable_two( FILE*fp, int norg, int nadd, LocalHom **localhomtable, LocalHom **localhomtablex, char *kozoarivec ) // for test only
4280 {
4281 	double opt;
4282 	static char buff[B];
4283 	char infor[100];
4284 	int i, j, overlapaa, start1, end1, start2, end2;
4285 	int **nlocalhom = NULL;
4286 	int **nlocalhomx = NULL;
4287 	LocalHom *tmpptr1=NULL, *tmpptr2=NULL; // by D.Mathog, a guess
4288 
4289 	nlocalhom = AllocateIntMtx( norg, nadd );
4290 	for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ ) nlocalhom[i][j] = 0;
4291 	nlocalhomx = AllocateIntMtx( nadd, norg );
4292 	for( i=0; i<nadd; i++ ) for( j=0; j<norg; j++ ) nlocalhomx[i][j] = 0;
4293 
4294 	while ( NULL != fgets( buff, B-1, fp ) )
4295 	{
4296 //		fprintf( stderr, "\n" );
4297 		sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4298 		if( *infor == 'k' )
4299 		{
4300 			fprintf( stderr, "Not supported!\n" );
4301 			exit( 1 );
4302 		}
4303 		j -= norg;
4304 
4305 #if 0
4306 		if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4307 #endif
4308 
4309 
4310 //		if( i < j )
4311 		{
4312 			if( nlocalhom[i][j]++ > 0 )
4313 			{
4314 //				fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4315 				tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4316 				tmpptr1 = tmpptr1->next;
4317 				tmpptr1->next = NULL;
4318 			}
4319 			else
4320 			{
4321 				tmpptr1 = localhomtable[i]+j;
4322 //				fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4323 			}
4324 
4325 			tmpptr1->start1 = start1; // CHUUI!!!!
4326 			tmpptr1->start2 = start2;
4327 			tmpptr1->end1 = end1; // CHUUI!!!!
4328 			tmpptr1->end2 = end2;
4329 //			tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4330 //			tmpptr1->opt = opt;
4331 			tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4332 			tmpptr1->overlapaa = overlapaa;
4333 			tmpptr1->korh = *infor;
4334 
4335 //			fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
4336 		}
4337 
4338 		{
4339 			if( nlocalhomx[j][i]++ > 0 )
4340 			{
4341 				tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4342 				tmpptr2 = tmpptr2->next;
4343 				tmpptr2->next = NULL;
4344 			}
4345 			else
4346 				tmpptr2 = localhomtablex[j]+i;
4347 
4348 			tmpptr2->start2 = start1+1; // CHUUI!!!!
4349 			tmpptr2->start1 = start2;
4350 			tmpptr2->end2 = end1+1; // CHUUI!!!!
4351 			tmpptr2->end1 = end2;
4352 //			tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4353 //			tmpptr2->opt = opt;
4354 			tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
4355 			tmpptr2->overlapaa = overlapaa;
4356 			tmpptr2->korh = *infor;
4357 
4358 //			fprintf( stderr, "j=%d, i=%d, opt = %f\n", j, i, opt );
4359 		}
4360 
4361 	}
4362 	FreeIntMtx( nlocalhom );
4363 	FreeIntMtx( nlocalhomx );
4364 }
4365 
readlocalhomtable_one(FILE * fp,int norg,int nadd,LocalHom ** localhomtable,char * kozoarivec)4366 void readlocalhomtable_one( FILE*fp, int norg, int nadd, LocalHom **localhomtable, char *kozoarivec ) // for test only
4367 {
4368 	double opt;
4369 	static char buff[B];
4370 	char infor[100];
4371 	int i, j, overlapaa, start1, end1, start2, end2;
4372 	int **nlocalhom = NULL;
4373 	LocalHom *tmpptr1=NULL; // by D.Mathog, a guess
4374 
4375 	nlocalhom = AllocateIntMtx( norg, nadd );
4376 	for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ ) nlocalhom[i][j] = 0;
4377 
4378 	while ( NULL != fgets( buff, B-1, fp ) )
4379 	{
4380 //		fprintf( stderr, "\n" );
4381 		sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4382 		if( *infor == 'k' )
4383 		{
4384 			fprintf( stderr, "Not supported!\n" );
4385 			exit( 1 );
4386 		}
4387 		j -= norg;
4388 
4389 #if 0
4390 		if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4391 #endif
4392 
4393 
4394 //		if( i < j )
4395 		{
4396 			if( nlocalhom[i][j]++ > 0 )
4397 			{
4398 //				fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4399 				tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4400 				tmpptr1 = tmpptr1->next;
4401 				tmpptr1->next = NULL;
4402 			}
4403 			else
4404 			{
4405 				tmpptr1 = localhomtable[i]+j;
4406 //				fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4407 			}
4408 
4409 			tmpptr1->start1 = start1; // CHUUI!!!!
4410 			tmpptr1->start2 = start2;
4411 			tmpptr1->end1 = end1; // CHUUI!!!!
4412 			tmpptr1->end2 = end2;
4413 //			tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4414 //			tmpptr1->opt = opt;
4415 			tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4416 			tmpptr1->overlapaa = overlapaa;
4417 			tmpptr1->korh = *infor;
4418 
4419 //			fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
4420 		}
4421 
4422 	}
4423 	FreeIntMtx( nlocalhom );
4424 }
4425 
outlocalhom_part(LocalHom ** localhom,int norg,int nadd)4426 void outlocalhom_part( LocalHom **localhom, int norg, int nadd )
4427 {
4428 	int i, j;
4429 	LocalHom *tmpptr;
4430 	for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ )
4431 	{
4432 		tmpptr = localhom[i]+j;
4433 		fprintf( stderr, "%d-%d\n", i, j+norg );
4434 		do
4435 		{
4436 			fprintf( stderr, "reg1=%d-%d, reg2=%d-%d, imp=%f, opt=%f\n", tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, tmpptr->importance, tmpptr->opt );
4437 		}
4438 		while( (tmpptr=tmpptr->next) );
4439 	}
4440 }
4441 
outlocalhom(LocalHom ** localhom,int nseq)4442 void outlocalhom( LocalHom **localhom, int nseq )
4443 {
4444 	int i, j;
4445 	LocalHom *tmpptr;
4446 	for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ )
4447 	{
4448 		tmpptr = localhom[i]+j;
4449 		fprintf( stderr, "%d-%d\n", i, j );
4450 		do
4451 		{
4452 			fprintf( stderr, "reg1=%d-%d, reg2=%d-%d, imp=%f, opt=%f\n", tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, tmpptr->importance, tmpptr->opt );
4453 		}
4454 		while( (tmpptr=tmpptr->next) );
4455 	}
4456 }
4457 
outlocalhompt(LocalHom *** localhom,int n1,int n2)4458 void outlocalhompt( LocalHom ***localhom, int n1, int n2 )
4459 {
4460 	int i, j;
4461 	LocalHom *tmpptr;
4462 	for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
4463 	{
4464 		tmpptr = localhom[i][j];
4465 //		fprintf( stdout, "%d-%d\n", i, j );
4466 		do
4467 		{
4468 			fprintf( stdout, "%d-%d, reg1=%d-%d, reg2=%d-%d, imp=%f, opt=%f, wimp=%f\n", i, j, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, tmpptr->importance, tmpptr->opt, tmpptr->wimportance );
4469 		}
4470 		while( (tmpptr=tmpptr->next) );
4471 	}
4472 }
4473 
FreeLocalHomTable_part(LocalHom ** localhomtable,int n,int m)4474 void FreeLocalHomTable_part( LocalHom **localhomtable, int n, int m )
4475 {
4476 	int i, j;
4477 	LocalHom *ppp, *tmpptr;
4478 	for( i=0; i<n; i++ )
4479 	{
4480 		for( j=0; j<m; j++ )
4481 		{
4482 			tmpptr=localhomtable[i]+j;
4483 			ppp = tmpptr->next;
4484 			for( ; tmpptr; tmpptr=ppp )
4485 			{
4486 #if DEBUG
4487 				fprintf( stderr, "i=%d, j=%d\n", i, j );
4488 #endif
4489 				ppp = tmpptr->next;
4490 				if( tmpptr!=localhomtable[i]+j )
4491 				{
4492 #if DEBUG
4493 					fprintf( stderr, "freeing %p\n", tmpptr );
4494 #endif
4495 					free( tmpptr );
4496 				}
4497 			}
4498 		}
4499 #if DEBUG
4500 		fprintf( stderr, "freeing localhomtable[%d]\n", i );
4501 #endif
4502 		free( localhomtable[i] );
4503 	}
4504 #if DEBUG
4505 	fprintf( stderr, "freeing localhomtable\n" );
4506 #endif
4507 	free( localhomtable );
4508 #if DEBUG
4509 	fprintf( stderr, "freed\n" );
4510 #endif
4511 }
4512 
FreeLocalHomTable_two(LocalHom ** localhomtable,int n,int m)4513 void FreeLocalHomTable_two( LocalHom **localhomtable, int n, int m )
4514 {
4515 	int i, j;
4516 	LocalHom *ppp, *tmpptr;
4517 	for( i=0; i<n; i++ )
4518 	{
4519 		for( j=0; j<m; j++ )
4520 		{
4521 			tmpptr=localhomtable[i]+j;
4522 			ppp = tmpptr->next;
4523 			for( ; tmpptr; tmpptr=ppp )
4524 			{
4525 #if DEBUG
4526 				fprintf( stderr, "i=%d, j=%d\n", i, j );
4527 #endif
4528 				ppp = tmpptr->next;
4529 				if( tmpptr!=localhomtable[i]+j )
4530 				{
4531 #if DEBUG
4532 					fprintf( stderr, "freeing %p\n", tmpptr );
4533 #endif
4534 					free( tmpptr );
4535 				}
4536 			}
4537 		}
4538 #if DEBUG
4539 		fprintf( stderr, "freeing localhomtable[%d]\n", i );
4540 #endif
4541 		free( localhomtable[i] );
4542 	}
4543 
4544 	for( i=n; i<n+m; i++ )
4545 	{
4546 		for( j=0; j<n; j++ )
4547 		{
4548 			tmpptr=localhomtable[i]+j;
4549 			ppp = tmpptr->next;
4550 			for( ; tmpptr; tmpptr=ppp )
4551 			{
4552 #if DEBUG
4553 				fprintf( stderr, "i=%d, j=%d\n", i, j );
4554 #endif
4555 				ppp = tmpptr->next;
4556 				if( tmpptr!=localhomtable[i]+j )
4557 				{
4558 #if DEBUG
4559 					fprintf( stderr, "freeing %p\n", tmpptr );
4560 #endif
4561 					free( tmpptr );
4562 				}
4563 			}
4564 		}
4565 #if DEBUG
4566 		fprintf( stderr, "freeing localhomtable[%d]\n", i );
4567 #endif
4568 		free( localhomtable[i] );
4569 	}
4570 #if DEBUG
4571 	fprintf( stderr, "freeing localhomtable\n" );
4572 #endif
4573 	free( localhomtable );
4574 #if DEBUG
4575 	fprintf( stderr, "freed\n" );
4576 #endif
4577 }
4578 
FreeLocalHomTable_one(LocalHom ** localhomtable,int n,int m)4579 void FreeLocalHomTable_one( LocalHom **localhomtable, int n, int m )
4580 {
4581 	int i, j;
4582 	LocalHom *ppp, *tmpptr;
4583 	for( i=0; i<n; i++ )
4584 	{
4585 		for( j=0; j<m; j++ )
4586 		{
4587 			tmpptr=localhomtable[i]+j;
4588 			ppp = tmpptr->next;
4589 			for( ; tmpptr; tmpptr=ppp )
4590 			{
4591 #if DEBUG
4592 				fprintf( stderr, "i=%d, j=%d\n", i, j );
4593 #endif
4594 				ppp = tmpptr->next;
4595 				if( tmpptr!=localhomtable[i]+j )
4596 				{
4597 #if DEBUG
4598 					fprintf( stderr, "freeing %p\n", tmpptr );
4599 #endif
4600 					free( tmpptr );
4601 				}
4602 			}
4603 		}
4604 #if DEBUG
4605 		fprintf( stderr, "freeing localhomtable[%d]\n", i );
4606 #endif
4607 		free( localhomtable[i] );
4608 	}
4609 
4610 #if DEBUG
4611 	fprintf( stderr, "freeing localhomtable\n" );
4612 #endif
4613 	free( localhomtable );
4614 #if DEBUG
4615 	fprintf( stderr, "freed\n" );
4616 #endif
4617 }
4618 
FreeLocalHomTable(LocalHom ** localhomtable,int n)4619 void FreeLocalHomTable( LocalHom **localhomtable, int n )
4620 {
4621 	int i, j;
4622 	LocalHom *ppp, *tmpptr;
4623 	for( i=0; i<n; i++ )
4624 	{
4625 		for( j=0; j<n; j++ )
4626 		{
4627 			tmpptr=localhomtable[i]+j;
4628 			ppp = tmpptr->next;
4629 			for( ; tmpptr; tmpptr=ppp )
4630 			{
4631 #if DEBUG
4632 				fprintf( stderr, "i=%d, j=%d\n", i, j );
4633 #endif
4634 				ppp = tmpptr->next;
4635 				if( tmpptr!=localhomtable[i]+j )
4636 				{
4637 #if DEBUG
4638 					fprintf( stderr, "freeing %p\n", tmpptr );
4639 #endif
4640 					free( tmpptr );
4641 				}
4642 			}
4643 		}
4644 #if DEBUG
4645 		fprintf( stderr, "freeing localhomtable[%d]\n", i );
4646 #endif
4647 		free( localhomtable[i] );
4648 	}
4649 #if DEBUG
4650 	fprintf( stderr, "freeing localhomtable\n" );
4651 #endif
4652 	free( localhomtable );
4653 #if DEBUG
4654 	fprintf( stderr, "freed\n" );
4655 #endif
4656 }
4657 
progName(char * str)4658 char *progName( char *str )
4659 {
4660     char *value;
4661     if( ( value = strrchr( str, '/' ) ) != NULL )
4662         return( value+1 );
4663     else
4664         return( str );
4665 }
4666 
tabtospace(char * str)4667 static void tabtospace( char *str )
4668 {
4669 	char *p;
4670 //	fprintf( stderr, "before = %s\n", str );
4671 	while( NULL != ( p = strchr( str , '\t' ) ) )
4672 	{
4673 		*p = ' ';
4674 	}
4675 //	fprintf( stderr, "after = %s\n", str );
4676 }
4677 
extractfirstword(char * str)4678 static char *extractfirstword( char *str )
4679 {
4680 	char *val = str;
4681 
4682 	tabtospace( str );
4683 	while( *str )
4684 	{
4685 		if( val == str && *str == ' ' )
4686 		{
4687 			val++; str++;
4688 		}
4689 		else if( *str != ' ' )
4690 		{
4691 			str++;
4692 		}
4693 		else if( *str == ' ' )
4694 		{
4695 			*str = 0;
4696 		}
4697 	}
4698 	return( val );
4699 }
4700 
phylipout_pointer(FILE * fp,int nseq,int maxlen,char ** seq,char ** name,int * order,int namelen)4701 void phylipout_pointer( FILE *fp, int nseq, int maxlen, char **seq, char **name, int *order, int namelen )
4702 {
4703 	int pos, pos2, j;
4704 	if( namelen == -1 ) namelen = 10;
4705 	pos = 0;
4706 
4707 	fprintf( fp, " %d %d\n", nseq, maxlen );
4708 
4709 	while( pos < maxlen )
4710 	{
4711 		for( j=0; j<nseq; j++ )
4712 		{
4713 			if( pos == 0 )
4714 				fprintf( fp, "%-*.*s", namelen, namelen, extractfirstword( name[order[j]]+1 ) );
4715 			else
4716 				fprintf( fp, "%-*.*s", namelen, namelen, "" );
4717 
4718 			pos2 = pos;
4719 			while( pos2 < pos+41 && pos2 < maxlen )
4720 			{
4721 				fprintf( fp, " %.10s", seq[order[j]]+pos2 );
4722 				pos2 += 10;
4723 			}
4724 			fprintf( fp, "\n" );
4725 		}
4726 		fprintf( fp, "\n" );
4727 		pos += 50;
4728 	}
4729 }
4730 
clustalout_pointer(FILE * fp,int nseq,int maxlen,char ** seq,char ** name,char * mark,char * comment,int * order,int namelen)4731 void clustalout_pointer( FILE *fp, int nseq, int maxlen, char **seq, char **name, char *mark, char *comment, int *order, int namelen )
4732 {
4733 	int pos, j;
4734 	if( namelen == -1 ) namelen = 15;
4735 	pos = 0;
4736 	if( comment == NULL )
4737 		fprintf( fp, "CLUSTAL format alignment by MAFFT (v%s)\n\n", VERSION );
4738 	else
4739 		fprintf( fp, "CLUSTAL format alignment by MAFFT %s (v%s)\n\n", comment, VERSION );
4740 
4741 	while( pos < maxlen )
4742 	{
4743 		fprintf( fp, "\n" );
4744 		for( j=0; j<nseq; j++ )
4745 		{
4746 			fprintf( fp, "%-*.*s ", namelen, namelen, extractfirstword( name[order[j]]+1 ) );
4747 			fprintf( fp, "%.60s\n", seq[order[j]]+pos ); // Ĺ�����㤦�Ȥ��ᡣ
4748 		}
4749 		if( mark )
4750 		{
4751 			fprintf( fp, "%-*.*s ", namelen, namelen, "" );
4752 			fprintf( fp, "%.60s\n", mark + pos ); // Ĺ�����㤦�Ȥ��ᡣ
4753 		}
4754 		pos += 60;
4755 	}
4756 }
4757 
4758 
writeData_reorder_pointer(FILE * fp,int locnjob,char ** name,int * nlen,char ** aseq,int * order)4759 void writeData_reorder_pointer( FILE *fp, int locnjob, char **name, int *nlen, char **aseq, int *order )
4760 {
4761 	int i, j, k;
4762 	int nalen;
4763 
4764 	for( i=0; i<locnjob; i++ )
4765 	{
4766 		k = order[i];
4767 #if DEBUG
4768 		fprintf( stderr, "i = %d in writeData\n", i );
4769 #endif
4770 		nalen = strlen( aseq[k] );
4771 		fprintf( fp, ">%s\n", name[k]+1 );
4772 		for( j=0; j<nalen; j=j+C )
4773 		{
4774 #if 0
4775 			strncpy( b, aseq[k]+j, C ); b[C] = 0;
4776 			fprintf( fp, "%s\n",b );
4777 #else
4778 			fprintf( fp, "%.*s\n", C, aseq[k]+j );
4779 #endif
4780 		}
4781 	}
4782 }
writeData_reorder(FILE * fp,int locnjob,char name[][B],int nlen[],char ** aseq,int * order)4783 void writeData_reorder( FILE *fp, int locnjob, char name[][B], int nlen[], char **aseq, int *order )
4784 {
4785 	int i, j, k;
4786 	int nalen;
4787 
4788 	for( i=0; i<locnjob; i++ )
4789 	{
4790 		k = order[i];
4791 #if DEBUG
4792 		fprintf( stderr, "i = %d in writeData\n", i );
4793 #endif
4794 		nalen = strlen( aseq[k] );
4795 		fprintf( fp, ">%s\n", name[k]+1 );
4796 		for( j=0; j<nalen; j=j+C )
4797 		{
4798 #if 0
4799 			strncpy( b, aseq[k]+j, C ); b[C] = 0;
4800 			fprintf( fp, "%s\n",b );
4801 #else
4802 			fprintf( fp, "%.*s\n", C, aseq[k]+j );
4803 #endif
4804 		}
4805 	}
4806 }
showaamtxexample()4807 static void showaamtxexample()
4808 {
4809 	fprintf( stderr, "Format error in aa matrix\n" );
4810 	fprintf( stderr, "# Example:\n" );
4811 	fprintf( stderr, "# comment\n" );
4812 	fprintf( stderr, "   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V\n" );
4813 	fprintf( stderr, "A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0\n" );
4814 	fprintf( stderr, "R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3\n" );
4815 	fprintf( stderr, "...\n" );
4816 	fprintf( stderr, "V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4\n" );
4817 	fprintf( stderr, "frequency 0.07 0.05 0.04 0.05 0.02 .. \n" );
4818 	fprintf( stderr, "# Example end\n" );
4819 	fprintf( stderr, "Only the lower half is loaded\n" );
4820 	fprintf( stderr, "The last line (frequency) is optional.\n" );
4821 	exit( 1 );
4822 }
4823 
loadaamtx(void)4824 double *loadaamtx( void )
4825 {
4826 	int i, j, k, ii, jj;
4827 	double *val;
4828 	double **raw;
4829 	int *map;
4830 	char *aaorder = "ARNDCQEGHILKMFPSTWYV";
4831 	char *inorder;
4832 	char *line;
4833 	char *ptr1;
4834 	char *ptr2;
4835 	char *mtxfname = "_aamtx";
4836 	FILE *mf;
4837 
4838 	raw = AllocateDoubleMtx( 21, 20 );
4839 	val = AllocateDoubleVec( 420 );
4840 	map = AllocateIntVec( 20 );
4841 
4842 	if( dorp != 'p' )
4843 	{
4844 		fprintf( stderr, "User-defined matrix is not supported for DNA\n" );
4845 		exit( 1 );
4846 	}
4847 
4848 	mf = fopen( mtxfname, "r" );
4849 	if( mf == NULL )
4850 	{
4851 		fprintf( stderr, "Cannot open the _aamtx file\n" );
4852 		exit( 1 );
4853 	}
4854 
4855 	inorder = calloc( 1000, sizeof( char ) );
4856 	line = calloc( 1000, sizeof( char ) );
4857 
4858 
4859 	while( !feof( mf ) )
4860 	{
4861 		fgets( inorder, 999, mf );
4862 		if( inorder[0] != '#' ) break;
4863 	}
4864 	ptr1 = ptr2 = inorder;
4865 	while( *ptr2 )
4866 	{
4867 		if( isalpha( *ptr2 ) )
4868 		{
4869 			*ptr1 = toupper( *ptr2 );
4870 			ptr1++;
4871 		}
4872 		ptr2++;
4873 	}
4874 	inorder[20] = 0;
4875 
4876 	for( i=0; i<20; i++ )
4877 	{
4878 		ptr2 = strchr( inorder, aaorder[i] );
4879 		if( ptr2 == NULL )
4880 		{
4881 			fprintf( stderr, "%c: not found in the first 20 letters.\n", aaorder[i] );
4882 			showaamtxexample();
4883 		}
4884 		else
4885 		{
4886 			map[i] = ptr2 - inorder;
4887 		}
4888 	}
4889 
4890 	i = 0;
4891 	while( !feof( mf ) )
4892 	{
4893 		fgets( line, 999, mf );
4894 //		fprintf( stderr, "line = %s\n", line );
4895 		if( line[0] == '#' ) continue;
4896 		ptr1 = line;
4897 //		fprintf( stderr, "line = %s\n", line );
4898 		for( j=0; j<=i; j++ )
4899 		{
4900 			while( !isdigit( *ptr1 ) && *ptr1 != '-' && *ptr1 != '.' )
4901 				ptr1++;
4902 
4903 			raw[i][j] = atof( ptr1 );
4904 //			fprintf( stderr, "raw[][]=%f, %c-%c %d-%d\n", raw[i][j], inorder[i], inorder[j], i, j );
4905 			ptr1 = strchr( ptr1, ' ' );
4906 			if( ptr1 == NULL && j<i) showaamtxexample();
4907 		}
4908 		i++;
4909 		if( i > 19 ) break;
4910 	}
4911 
4912 	for( i=0; i<20; i++ ) raw[20][i] = -1.0;
4913 	while( !feof( mf ) )
4914 	{
4915 		fgets( line, 999, mf );
4916 		if( line[0] == 'f' )
4917 		{
4918 //			fprintf( stderr, "line = %s\n", line );
4919 			ptr1 = line;
4920 			for( j=0; j<20; j++ )
4921 			{
4922 				while( !isdigit( *ptr1 ) && *ptr1 != '-' && *ptr1 != '.' )
4923 					ptr1++;
4924 
4925 				raw[20][j] = atof( ptr1 );
4926 //				fprintf( stderr, "raw[20][]=%f, %c %d\n", raw[20][j], inorder[i], j );
4927 				ptr1 = strchr( ptr1, ' ' );
4928 				if( ptr1 == NULL && j<19) showaamtxexample();
4929 			}
4930 			break;
4931 		}
4932 	}
4933 
4934 	k = 0;
4935 	for( i=0; i<20; i++ )
4936 	{
4937 		for( j=0; j<=i; j++ )
4938 		{
4939 			if( i != j )
4940 			{
4941 				ii = MAX( map[i], map[j] );
4942 				jj = MIN( map[i], map[j] );
4943 			}
4944 			else ii = jj = map[i];
4945 			val[k++] = raw[ii][jj];
4946 //			fprintf( stderr, "%c-%c, %f\n", aaorder[i], aaorder[j], val[k-1] );
4947 		}
4948 	}
4949 	for( i=0; i<20; i++ ) val[400+i] = raw[20][map[i]];
4950 
4951 	fprintf( stderr, "inorder = %s\n", inorder );
4952 	fclose( mf );
4953 	free( inorder );
4954 	free( line );
4955 	FreeDoubleMtx( raw );
4956 	free( map );
4957 	return( val );
4958 }
4959 
tab2space(char * s)4960 static void tab2space( char *s ) // nen no tame
4961 {
4962 	while( *s )
4963 	{
4964 		if( *s == '\t' ) *s = ' ';
4965 		s++;
4966 	}
4967 }
4968 
readasubalignment(char * s,int * t,int * preservegaps)4969 static int readasubalignment( char *s, int *t, int *preservegaps )
4970 {
4971 	int v = 0;
4972 	char status = 's';
4973 	char *pt = s;
4974 	*preservegaps = 0;
4975 	tab2space( s );
4976 	while( *pt )
4977 	{
4978 		if( *pt == ' ' )
4979 		{
4980 			status = 's';
4981 		}
4982 		else
4983 		{
4984 			if( status == 's' )
4985 			{
4986 				if( *pt == '\n' || *pt == '#' ) break;
4987 				status = 'n';
4988 				t[v] = atoi( pt );
4989 				if( t[v] == 0 )
4990 				{
4991 					fprintf( stderr, "Format error? Sequences must be specified as 1, 2, 3...\n" );
4992 					exit( 1 );
4993 				}
4994 				if( t[v] < 0 ) *preservegaps = 1;
4995 				t[v] = abs( t[v] );
4996 				t[v] -= 1;
4997 				v++;
4998 			}
4999 		}
5000 		pt++;
5001 	}
5002 	t[v] = -1;
5003 	return( v );
5004 }
5005 
countspace(char * s)5006 static int countspace( char *s )
5007 {
5008 	int v = 0;
5009 	char status = 's';
5010 	char *pt = s;
5011 	tab2space( s );
5012 	while( *pt )
5013 	{
5014 		if( *pt == ' ' )
5015 		{
5016 			status = 's';
5017 		}
5018 		else
5019 		{
5020 			if( status == 's' )
5021 			{
5022 				if( *pt == '\n' || *pt == '#' ) break;
5023 				v++;
5024 				status = 'n';
5025 				if( atoi( pt ) == 0 )
5026 				{
5027 					fprintf( stderr, "Format error? Sequences should be specified as 1, 2, 3...\n" );
5028 					exit( 1 );
5029 				}
5030 			}
5031 		}
5032 		pt++;
5033 	}
5034 	return( v );
5035 }
5036 
5037 
readsubalignmentstable(int nseq,int ** table,int * preservegaps,int * nsubpt,int * maxmempt)5038 void readsubalignmentstable( int nseq, int **table, int *preservegaps, int *nsubpt, int *maxmempt )
5039 {
5040 	FILE *fp;
5041 	char *line;
5042 	int linelen = 1000000;
5043 	int nmem;
5044 	int lpos;
5045 	int i, p;
5046 	int *tab01;
5047 
5048 	line = calloc( linelen, sizeof( char ) );
5049 	fp = fopen( "_subalignmentstable", "r" );
5050 	if( !fp )
5051 	{
5052 		fprintf( stderr, "Cannot open _subalignmentstable\n" );
5053 		exit( 1 );
5054 	}
5055 	if( table == NULL )
5056 	{
5057 		*nsubpt = 0;
5058 		*maxmempt = 0;
5059 		while( 1 )
5060 		{
5061 			fgets( line, linelen-1, fp );
5062 			if( feof( fp ) ) break;
5063 			if( line[strlen(line)-1] != '\n' )
5064 			{
5065 				fprintf( stderr, "too long line? \n" );
5066 				exit( 1 );
5067 			}
5068 			if( line[0] == '#' ) continue;
5069 			if( atoi( line ) == 0 ) continue;
5070 			nmem = countspace( line );
5071 			if( nmem > *maxmempt ) *maxmempt = nmem;
5072 			(*nsubpt)++;
5073 		}
5074 	}
5075 	else
5076 	{
5077 		tab01 = calloc( nseq, sizeof( int ) );
5078 		for( i=0; i<nseq; i++ ) tab01[i] = 0;
5079 		lpos = 0;
5080 		while( 1 )
5081 		{
5082 			fgets( line, linelen-1, fp );
5083 			if( feof( fp ) ) break;
5084 			if( line[strlen(line)-1] != '\n' )
5085 			{
5086 				fprintf( stderr, "too long line? \n" );
5087 				exit( 1 );
5088 			}
5089 			if( line[0] == '#' ) continue;
5090 			if( atoi( line ) == 0 ) continue;
5091 			readasubalignment( line, table[lpos], preservegaps+lpos );
5092 			for( i=0; (p=table[lpos][i])!=-1; i++ )
5093 			{
5094 				if( tab01[p] )
5095 				{
5096 					fprintf( stderr, "\nSequence %d appears in different groups.\n", p+1 );
5097 					fprintf( stderr, "Hierarchical grouping is not supported.\n\n" );
5098 					exit( 1 );
5099 				}
5100 				tab01[p] = 1;
5101 				if( p > nseq-1 )
5102 				{
5103 					fprintf( stderr, "Sequence %d does not exist in the input sequence file.\n", p+1 );
5104 					exit( 1 );
5105 				}
5106 			}
5107 			lpos++;
5108 		}
5109 		free( tab01 );
5110 	}
5111 	fclose( fp );
5112 	free( line );
5113 }
5114 
5115 
readmccaskill(FILE * fp,RNApair ** pairprob,int length)5116 void readmccaskill( FILE *fp, RNApair **pairprob, int length )
5117 {
5118 	char gett[1000];
5119 	int *pairnum;
5120 	int i;
5121 	int left, right;
5122 	double prob;
5123 	int c;
5124 
5125 	pairnum = (int *)calloc( length, sizeof( int ) );
5126 	for( i=0; i<length; i++ ) pairnum[i] = 0;
5127 
5128 	c = getc( fp );
5129 	{
5130 		if( c != '>' )
5131 		{
5132 			fprintf( stderr, "format error in hat4 - 1\n" );
5133 			exit( 1 );
5134 		}
5135 	}
5136 	fgets( gett, 999, fp );
5137 	while( 1 )
5138 	{
5139 		if( feof( fp ) ) break;
5140 		c = getc( fp );
5141 		ungetc( c, fp );
5142 		if( c == '>' || c == EOF )
5143 		{
5144 			break;
5145 		}
5146 		fgets( gett, 999, fp );
5147 //		fprintf( stderr, "gett = %s\n", gett );
5148 		sscanf( gett, "%d %d %lf", &left, &right, &prob );
5149 
5150 		if( left >= length || right >= length )
5151 		{
5152 			fprintf( stderr, "format error in hat4 - 2\n" );
5153 			exit( 1 );
5154 		}
5155 
5156 		if( prob < 0.01 ) continue; // 080607, mafft ni dake eikyou
5157 
5158 		if( left != right && prob > 0.0 )
5159 		{
5160 			pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
5161 			pairprob[left][pairnum[left]].bestscore = prob;
5162 			pairprob[left][pairnum[left]].bestpos = right;
5163 			pairnum[left]++;
5164 			pairprob[left][pairnum[left]].bestscore = -1.0;
5165 			pairprob[left][pairnum[left]].bestpos = -1;
5166 //			fprintf( stderr, "%d-%d, %f\n", left, right, prob );
5167 
5168 			pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
5169 			pairprob[right][pairnum[right]].bestscore = prob;
5170 			pairprob[right][pairnum[right]].bestpos = left;
5171 			pairnum[right]++;
5172 			pairprob[right][pairnum[right]].bestscore = -1.0;
5173 			pairprob[right][pairnum[right]].bestpos = -1;
5174 //			fprintf( stderr, "%d-%d, %f\n", right, left, prob );
5175 		}
5176 	}
5177 	free( pairnum );
5178 }
5179 
readpairfoldalign(FILE * fp,char * s1,char * s2,char * aln1,char * aln2,int q1,int q2,int * of1,int * of2,int sumlen)5180 void readpairfoldalign( FILE *fp, char *s1, char *s2, char *aln1, char *aln2, int q1, int q2, int *of1, int *of2, int sumlen )
5181 {
5182 	char gett[1000];
5183 	int *maptoseq1;
5184 	int *maptoseq2;
5185 	char dumc;
5186 	int dumi;
5187 	char sinseq[100], sinaln[100];
5188 	int posinseq, posinaln;
5189 	int alnlen;
5190 	int i;
5191 	int pos1, pos2;
5192 	char *pa1, *pa2;
5193 	char qstr[1000];
5194 
5195 	*of1 = -1;
5196 	*of2 = -1;
5197 
5198 	maptoseq1 = AllocateIntVec( sumlen+1 );
5199 	maptoseq2 = AllocateIntVec( sumlen+1 );
5200 
5201 	posinaln = 0; // foldalign ga alingment wo kaesanaitok no tame.
5202 
5203 	while( !feof( fp ) )
5204 	{
5205 		fgets( gett, 999, fp );
5206 		if( !strncmp( gett, "; ALIGNING", 10 ) ) break;
5207 	}
5208 	sprintf( qstr, "; ALIGNING            %d against %d\n", q1+1, q2+1 );
5209 	if( strcmp( gett, qstr ) )
5210 	{
5211 		fprintf( stderr, "Error in FOLDALIGN\n" );
5212 		fprintf( stderr, "qstr = %s, but gett = %s\n", qstr, gett );
5213 		exit( 1 );
5214 	}
5215 
5216 	while( !feof( fp ) )
5217 	{
5218 		fgets( gett, 999, fp );
5219 		if( !strncmp( gett, "; --------", 10 ) ) break;
5220 	}
5221 
5222 
5223 	while( !feof( fp ) )
5224 	{
5225 		fgets( gett, 999, fp );
5226 		if( !strncmp( gett, "; ********", 10 ) ) break;
5227 //		fprintf( stderr, "gett = %s\n", gett );
5228 		sscanf( gett, "%c %c %s %s %d %d", &dumc, &dumc, sinseq, sinaln, &dumi, &dumi );
5229 		posinaln = atoi( sinaln );
5230 		posinseq = atoi( sinseq );
5231 //		fprintf( stderr, "posinseq = %d\n", posinseq );
5232 //		fprintf( stderr, "posinaln = %d\n", posinaln );
5233 		maptoseq1[posinaln-1] = posinseq-1;
5234 	}
5235 	alnlen = posinaln;
5236 
5237 	while( !feof( fp ) )
5238 	{
5239 		fgets( gett, 999, fp );
5240 		if( !strncmp( gett, "; --------", 10 ) ) break;
5241 	}
5242 
5243 	while( !feof( fp ) )
5244 	{
5245 		fgets( gett, 999, fp );
5246 		if( !strncmp( gett, "; ********", 10 ) ) break;
5247 //		fprintf( stderr, "gett = %s\n", gett );
5248 		sscanf( gett, "%c %c %s %s %d %d", &dumc, &dumc, sinseq, sinaln, &dumi, &dumi );
5249 		posinaln = atof( sinaln );
5250 		posinseq = atof( sinseq );
5251 //		fprintf( stderr, "posinseq = %d\n", posinseq );
5252 //		fprintf( stderr, "posinaln = %d\n", posinaln );
5253 		maptoseq2[posinaln-1] = posinseq-1;
5254 	}
5255 	if( alnlen != posinaln )
5256 	{
5257 		fprintf( stderr, "Error in foldalign?\n" );
5258 		exit( 1 );
5259 	}
5260 
5261 	pa1 = aln1;
5262 	pa2 = aln2;
5263 	for( i=0; i<alnlen; i++ )
5264 	{
5265 		pos1 = maptoseq1[i];
5266 		pos2 = maptoseq2[i];
5267 
5268 		if( pos1 > -1 )
5269 			*pa1++ = s1[pos1];
5270 		else
5271 			*pa1++ = '-';
5272 
5273 		if( pos2 > -1 )
5274 			*pa2++ = s2[pos2];
5275 		else
5276 			*pa2++ = '-';
5277 	}
5278 	*pa1 = 0;
5279 	*pa2 = 0;
5280 
5281 	*of1 = 0;
5282 	for( i=0; i<alnlen; i++ )
5283 	{
5284 		*of1 = maptoseq1[i];
5285 		if( *of1 > -1 ) break;
5286 	}
5287 	*of2 = 0;
5288 	for( i=0; i<alnlen; i++ )
5289 	{
5290 		*of2 = maptoseq2[i];
5291 		if( *of2 > -1 ) break;
5292 	}
5293 
5294 //	fprintf( stderr, "*of1=%d, aln1 = :%s:\n", *of1, aln1 );
5295 //	fprintf( stderr, "*of2=%d, aln2 = :%s:\n", *of2, aln2 );
5296 
5297 	free( maptoseq1 );
5298 	free( maptoseq2 );
5299 }
5300 
myatoi(char * in)5301 int myatoi( char *in )
5302 {
5303 	if( in == NULL )
5304 	{
5305 		fprintf( stderr, "Error in myatoi()\n" );
5306 		exit( 1 );
5307 	}
5308 	return( atoi( in ) );
5309 }
5310 
myatof(char * in)5311 double myatof( char *in )
5312 {
5313 	if( in == NULL )
5314 	{
5315 		fprintf( stderr, "Error in myatof()\n" );
5316 		exit( 1 );
5317 	}
5318 	return( atof( in ) );
5319 }
5320 
reporterr(const char * str,...)5321 void reporterr( const char *str, ... )
5322 {
5323 //	static int loglen = 0;
5324 	va_list args;
5325 
5326 	if( gmsg )
5327 	{
5328 # if 1  // ato de sakujo
5329 		static FILE *errtmpfp = NULL;
5330 		if( errtmpfp == NULL )
5331 			errtmpfp = fopen( "maffterr", "w" );
5332 		else
5333 			errtmpfp = fopen( "maffterr", "a" );
5334 		va_start( args, str );
5335 		vfprintf( errtmpfp, str, args );
5336 		va_end( args );
5337 		fclose( errtmpfp );
5338 #endif
5339 
5340 #if 0
5341 		char *tmpptr;
5342 		tmpptr = (char *)realloc( *gmsg, (loglen+10000) * sizeof( char ) );
5343 		if( tmpptr == NULL )
5344 		{
5345 			fprintf( stderr, "Cannot relloc *gmsg\n" );
5346 			exit( 1 );
5347 		}
5348 		*gmsg = tmpptr;
5349 		va_start( args, str );
5350 		loglen += vsprintf( *gmsg + loglen, str, args );
5351 		va_end( args );
5352 
5353 
5354 		va_start( args, str );
5355 		loglen += vsprintf( *gmsg + loglen, str, args );
5356 		va_end( args );
5357 		*(*gmsg + loglen) = 0;
5358 		if( loglen > gmsglen - 100 ) loglen = 0; // tekitou
5359 #endif
5360 
5361 	}
5362 	else
5363 	{
5364 		va_start( args, str );
5365 		vfprintf( stderr, str, args );
5366 		va_end( args );
5367 //		fflush( stderr ); // iru?
5368 	}
5369 	return;
5370 }
5371